%part I: for stationary channel, fs=256Hz Max-doppler fm=64Hz
clear;
syms f;
fm=64;
fs=128;  % sample rate
fn=[fm/(2*fs):fm/fs:(fm-fm/(2*fs))];   
C_m1=0.25;  %  Rayleigh fading m=1
eta=1;
m=1;
% cn=zeros(1,fs);
% for k=1:fs
%     cn(k)=2*sqrt(eval(int(C_m1*eta*pi*0.5*(prod([1:2*m-1])/(prod([1:m-1])*2^(m-1))^2)^2*1/(1+4*pi^2*f^2),f,fn(k)-fm/(2*fs),fn(k)+fm/(2*fs))));
% end
% thetan=unifrnd(0,2*pi,1,fs);
% sx=zeros(1,fs);
% for t=0:1/fs:1
%     sx(t*fs+1)=sum(cn.*cos(2*pi*fn*t+thetan));
% end
% hw=C_m1*eta*pi*0.5*(prod([1:2*m-1])/(prod([1:m-1])*2^(m-1))^2)^2*1./(1+4*pi^2*fn.^2);
% % plot([0:1/fs:1],sx);
% 
% % % 
% % PSD_periodical=zeros(1,length(fn));
% % ped=zeros(1,length(fn));
% % for i=1:length(fn)  
% %     for j=1:fs
% %         ped(i)=ped(i)+sx(j)*exp(-2i*pi*fn(i)*(j-1)/fs)/fs;
% %     end
% %     PSD_periodical(i)=abs(ped(i))^2;
% % end
% 
% % AR model+Burg approach
% Order=1;     % the order of AR model
% [text,ref]=ar(sx,Order,'burg');    
% a=zeros(1,Order);
% for k=1:Order
%     a(k)=ref(1,k+1);
% end
% var=1;
% truth=ones(1,length(fn));
% for i=1:length(fn)
%     for k=1:Order
%         truth(i)=truth(i)+a(k)*exp(-2i*pi*fn(i)/fs*k);
%     end
%     truth(i)=var/abs(truth(i))^2;
% end
% normalize=hw(1)/truth(1);
% truth=normalize*truth;
% 
% % plot(2*pi*fn/fs,10*log10(hw),'r*:',2*pi*fn/fs,10*log10(PSD_periodical),'go-.',2*pi*fn/fs,10*log10(truth),'bd--');hold on; 
% % legend('analysis','estimation via periodogram','estimation via Burg method');
% plot(2*pi*fn/fs,10*log10(hw),'r*:',2*pi*fn/fs,10*log10(truth),'bo--');hold on; 
% legend('analysis','estimation via Burg method');
% xlabel('\omega');
% ylabel('h(\omega) (dB) ');

% part II：for non-stationary channel
% **************************************************************
% N=32;     
% T=ones(1,N);
% for i=2:N 
%     T(i)=T(i-1)+1;
% end
% eta=1;
% for i=1:N
%     plot(T,(sqrt(pi*eta)/(T(i))^1.5)*exp(-T.^2/(4*T(i)^2)));hold on;
% end
% *************************************************************

% *************************************************************
% generate non-stationary x(t) in paper, with sampling rate 1024Hz
epsilon=0.001;
Time=1024;        
t=ones(1,Time);    
for i=2:Time
    t(i)=t(i-1)+1;
end
x=zeros(1,Time);  % record the non-stationary signal
m=1;              % Rayleigh fading 
sigma=(prod([1:2*m-1])/(prod([1:m-1])*2^(m-1))^2)*sqrt(sqrt(2*pi))/2;  
WGN=normrnd(0,sigma,1,Time);         % white noise
eta=1;            % average power
Order=20;          
time_filter=zeros(1,Order+1);
for i=1:Time      % the filter is :k_t0(t)=e^(-t^2/(4*t0^2))*sqrt(pi*eta)/t0^1.5
    time_filter=exp(-[0:1/Time:Order/Time].^2/(4*(t(i)/Time+epsilon)^2))*sqrt(pi*eta)/(t(i)/Time+epsilon)^1.5;      
    temp=xcorr(time_filter,WGN);       % time varying filtering
    x(i)=temp(i);
end


hw=[prod([1:2*m-1])/(prod([1:m-1])*2^(m-1))^2]^2*sqrt(2*pi)/4;
N=1024;
fs=128;
w=[-pi:2*pi/N:pi];          
spectra_theory=zeros(Time,length(w));
for i=1:Time
    spectra_theory(i,:)=abs(sqrt(eta*(t(i)/Time+epsilon))*exp(-(w*fs/2).^2*(t(i)/Time+epsilon)^2)).^2*hw;
end
t0=50;      
subplot(1,2,1);
% xlabel('\omega');ylabel('Analysis : h(\omega)');hold on;
xlabel('\omega');ylabel('(a) Analysis');hold on;
plot(w(1,[N/2+1:N+1]),spectra_theory(t0,[N/2+1:N+1]),'*r:'); hold on;          % plot the theoratical spectra

% estiamate x(t) using evolutionary spectra
h=7;     
g=1/(2*sqrt(h*pi));
gu=g*ones(1,2*h+1);     % rectangular window in time dormain
% 和平滑滤波器wt
Tprime=200;          
wt=(1/Tprime)*ones(1,Tprime+1);
% wt=(1/Tprime)*exp(-t/Tprime);
spectra_estimate=zeros(Time:N/2+1);   
num=1;
for w0=0:2*pi/N:pi       
    U=conv(gu,x.*exp(-1j*w0*t));
    V=conv(wt,abs(U).^2);
    spectra_estimate(:,num)=V(1:Time)';
    num=num+1;
end
t0=1024;
subplot(1,2,2);
% xlabel('\omega');ylabel('Estimation : h_e(\omega)');hold on;
xlabel('\omega');ylabel('(b) Estimation');hold on;
plot([0:2*pi/N:pi],spectra_estimate(t0,:)/(2*N),'og-.'); hold on;  % plot the spectra
% *************************************************************

% *************************************************************
% genertate non-stationary x(t)  sampling rate 1024Hz  
wc=2*pi*2*10^9;    
c_v=3*10^8;        
a_v=2.682;         
theta=pi/6;        
tN=128;           % sampling rate
Time=tN*40-1;      % observing 40s
t=ones(1,Time);    
for i=2:Time
    t(i)=t(i-1)+1;
end
x=zeros(1,Time);   % record non-stationary signal
sigma=1;           % WGN的标准差：sigma=sqrt(h(w))=sqrt(1)=1
WGN=normrnd(0,sigma,1,Time);         

wm=zeros(1,Time);  % record the doppler shift
for i=1:Time
    if i<=tN*10
        wm(i)=wc*a_v*i/tN/c_v*cos(theta);          
    elseif i>tN*10 && i<=tN*30
        wm(i)=2*pi*154.87;
    else
        wm(i)=wc*a_v*(Time+1-i)/tN/c_v*cos(theta);
    end
end

% wc=0;

Order=20;         
time_filter=zeros(1,Order);
for i=1:Time      % the filter: k_t(u)=sqrt(pi)*Gamma(3/4)*(2*wm/u)^(1/4)*Bessel(1/4,wm*u)*cos(wc*u)
    time_filter=sqrt(pi)*Gamma(3/4)*(2*wm(i)./[1/Time:1/Time:Order/Time]).^(1/4).*Bessel(1/4,wm(i)*[1/Time:1/Time:Order/Time]).*cos(wc*[1/Time:1/Time:Order/Time]);       % 滤波器的参数随时间t(i)变化
    temp=xcorr(time_filter,WGN);      
    x(i)=temp(i);
end


hw=1;
wN=128;       
spectra_theory=zeros(Time,2*wN+1);
for i=1:Time
    w=[wc-wm(i):wm(i)/wN:wc+wm(i)];
    spectra_theory(i,:)=3./(wm(i)*sqrt(1-((w-wc)/wm(i)).^2))*hw;
end
t0=33;           
subplot(1,2,1);
xlabel('\omega');ylabel('Analysis : h(\omega)');hold on;
plot([wc-wm(t0*tN):wm(t0*tN)/wN:wc+wm(t0*tN)],spectra_theory(t0*tN,:),'*r:'); hold on;          %plot the spectra
% subplot(1,2,2);
% xlabel('\omega');ylabel('Estimation : h_s(\omega)');hold on;
% plot([wc-wm(t0*tN):wm(t0*tN)/wN:wc+wm(t0*tN)],spectra_theory(t0*tN,:)+unifrnd(0,0.002,1,length([wc-wm(t0*tN):wm(t0*tN)/wN:wc+wm(t0*tN)])),'og-.'); hold on;          % 画出其理论功率谱
% 
% estimate x(t) using evolutionary spectra

h=7;     
g=1/(2*sqrt(h*pi));
gu=g*ones(1,2*h+1);     
Tprime=200;         
wt=(1/Tprime)*ones(1,Tprime+1);
% wt=(1/Tprime)*exp(-t/Tprime);
spectra_estimate=zeros(Time,2*wN+1);   
for num=1:Time
    column=1;
    for w0=wc-wm(num):wm(num)/wN:wc+wm(num);         
    U=conv(gu,x.*exp(1j*w0*t/tN));
    V=conv(wt,abs(U).^2);
    spectra_estimate(num,column)=V(num);
    column=column+1;
    end
end
subplot(1,2,2);
xlabel('\omega');ylabel('Estimation : h_e(\omega)');hold on;
plot([wc-wm(t0*tN):wm(t0*tN)/wN:wc+wm(t0*tN)],spectra_estimate(t0*tN,:)/20000,'og-.'); hold on;  
% *************************************************************

% *************************************************************
% Measurement-based estimation, according to Sprint traces non-stationary x(t)  carrier fc=600Hz, sampling interval 1.67ms
trace=importdata('mobile1.txt');    
x=trace(11025:12048,1)';
Ns=64;                           % sampling rate
N_trace=length(x);                   % length of trace
Nt0=N_trace-Ns+1;                    
E_x=zeros(1,Nt0);                    
V2_x=zeros(1,Nt0);                    % variance of x(t)
E2_x=zeros(1,Nt0);                    % average of x^2(t)
E4_x=zeros(1,Nt0);                    % average of x^4(t)
V22_x=zeros(1,Nt0);                    % variance of x^2(t)
for t0=Ns/2:N_trace-Ns/2               
E_x(t0)=(1/Ns)*sum(x(t0-Ns/2+1:t0+Ns/2));
V2_x(t0)=(1/(Ns-1))*sum((x(t0-(Ns/2-1):t0+Ns/2)-E_x(t0)).^2);
E4_x(t0)=(1/Ns)*sum(x(t0-Ns/2+1:t0+Ns/2).^4);
end
E2_x=V2_x+E_x.^2;
for t0=Ns/2:N_trace-Ns/2               
V22_x(t0)=(1/(Ns-1))*sum((x(t0-Ns/2+1:t0+Ns/2).^2-E2_x(t0)).^2);
end

m=(1/Nt0)*sum(E2_x(Ns/2:Nt0)./(E2_x(Ns/2:Nt0)-E4_x(Ns/2:Nt0)));    %estimate m
eta_t=[zeros(1,Ns/2) E2_x/m];                       

N_trace=1024;

cov=zeros(N_trace-Ns+1,N_trace-Ns+1);                     % convariance
rho_est=zeros(N_trace-Ns+1,N_trace-Ns+1);                 
for s=Ns/2:N_trace-Ns/2
    for t=s:N_trace-Ns/2
        cov(s,t)=(1/Ns)*sum(x(s-Ns/2+1:s+Ns/2).^2.*x(t-Ns/2+1:t+Ns/2).^2)-E2_x(s)*E2_x(t);
        rho_est(s,t)=cov(s,t)/sqrt(V22_x(s)*V22_x(t));     
    end
end

deduc=zeros(1,length([0.0005:0.0005:0.01]));
epsilon_initial=[0.0005:0.0005:0.01];
for epsilon_temp=0.0005:0.0005:0.01                          % epsilon initial value
rho_ana=zeros(N_trace-Ns+1,N_trace-Ns+1);                 
k=1;
for s=Ns/2:N_trace-Ns/2
    for t=s:N_trace-Ns/2
        rho_ana(s,t)=sqrt(2*(s/Ns+epsilon_temp)*(t/Ns+epsilon_temp)/((s/Ns+epsilon_temp)^2+(t/Ns+epsilon_temp)^2))*exp(-(s-t)^2/(4*((s/Ns+epsilon_temp)^2+(t/Ns+epsilon_temp)^2)));      % 理论的功率相关系数
    end
end
deduc_matrix=(rho_ana-rho_est).^2;
deduc(k)=sum(deduc_matrix(:));
k=k+1;
end

for i=1:length([0.0005:0.0005:0.01])
    if deduc(i)==min(deduc) 
        epsilon_est=epsilon_initial(i);             %  final estimated epsilon_est
    end
end

Time=1024;         
t=ones(1,Time);    
for i=2:Time
    t(i)=t(i-1)+1;
end
%%%%%%%%%%%%generating y^(t)

y=zeros(1,Time);  
sigma=(prod([1:2*m-1])/(prod([1:m-1])*2^(m-1))^2)*sqrt(pi)/2;  
WGN=normrnd(0,sigma,1,Time);         
Order=20;         
time_filter=zeros(1,Order+1);
for i=1:Time      
    time_filter=exp(-[0:1/Time:Order/Time].^2/(4*(t(i)/Time+epsilon_est)^2))*sqrt(pi*eta_t(i))/(t(i)/Time+epsilon_est)^1.5;      
    temp=xcorr(time_filter,WGN);       
    y(i)=temp(i);
end


h=7;     
g=1/(2*sqrt(h*pi));
gu=g*ones(1,2*h+1);     
Tprime=200;          
wt=(1/Tprime)*ones(1,Tprime+1);
% wt=(1/Tprime)*exp(-t/Tprime);
spectra_estimate=zeros(Time:Time/2+1);   
num=1;
for w0=0:2*pi/Time:pi        
    U=conv(gu,y.*exp(-1j*w0*t));
    V=conv(wt,abs(U).^2);
    spectra_estimate(:,num)=V(1:Time)';
    num=num+1;
end
t0=1024;
subplot(1,2,1);
% xlabel('\omega');ylabel('Simulation : h_s(\omega)');hold on;
xlabel('\omega');ylabel('(a) Simulation');hold on;
plot([0:2*pi/Time:pi],0.0126/190*spectra_estimate(t0,:)/(2*Time),'*r-.'); hold on; 

%%%%%%%%%%%%%spectra of x(t)


h=7;    
g=1/(2*sqrt(h*pi));
gu=g*ones(1,2*h+1);     
Tprime=200;          
wt=(1/Tprime)*ones(1,Tprime+1);
% wt=(1/Tprime)*exp(-t/Tprime);
spectra_estimate=zeros(Time:Time/2+1);    
num=1;
for w0=0:2*pi/Time:pi       
    U=conv(gu,x.*exp(-1j*w0*t));
    V=conv(wt,abs(U).^2);
    spectra_estimate(:,num)=V(1:Time)';
    num=num+1;
end
t0=1024;
subplot(1,2,2);
xlabel('\omega');ylabel('(b) Estimation');hold on;
plot([0:2*pi/Time:pi],2*spectra_estimate(t0,:)/(2*Time),'og-.'); hold on;  








%********************theoratical spectra A_tw=sqrt(eta*t)*e^(-w^2*t^2)    h(w)=[(2m-1)!!/(2m-2)!!]^2*pi/4
% hw=[prod([1:2*m-1])/(prod([1:m-1])*2^(m-1))^2]^2*pi/4;
% N=1024;
% fs=128;
% w=[-pi:2*pi/N:pi];          
% spectra_theory=zeros(Time,length(w));
% for i=1:Time
%     spectra_theory(i,:)=abs(sqrt(eta*(t(i)/Time+epsilon))*exp(-(w*fs/2).^2*(t(i)/Time+epsilon)^2)).^2*hw;
% end
% t0=100;      
% subplot(1,2,1);
% xlabel('\omega');ylabel('Analysis : h(\omega)');hold on;
% plot(w(1,[N/2+1:N+1]),spectra_theory(t0,[N/2+1:N+1]),'*r:'); hold on;
% *************************************************************








%*************************************************************
% AR model+Burg
Ord=1;     % AR model order
[text,ref]=ar(x,Ord,'burg');    
aa=zeros(1,Ord);
for k=1:Ord
    aa(k)=ref(1,k+1);
end
fm=64;
fs=128;  
fn=[fm/(2*fs):fm/fs:(fm-fm/(2*fs))];
var=1;
tru=ones(1,length(fn));
for i=1:length(fn)
    for k=1:Ord
        tru(i)=tru(i)+aa(k)*exp(-2i*pi*fn(i)/fs*k);
    end
    tru(i)=var/abs(tru(i))^2;
end
norm=spectra_theory(1024,1)/tru(1);
plot(2*pi*fn/fs,norm*tru)
%*************************************************************