%用matlab的例子來驗證我的STFT的正確性
%越界賦0值
clear ;

sample_rate=500 ;
dt=1.0/sample_rate ;
sample_number=512*2 ;
df=1/(2*dt)/(sample_number/2) ;

%-4 7 40
Ap=-30.0 ;
Bp=85 ;
Cp=30 ;

signal_t=zeros(sample_number,1) ;
signal_y=zeros(sample_number,1) ;

for m=1:sample_number
    signal_t(m)=(m-1)*dt ;
end
for n=1:sample_number
    signal_w(n)=(n-1)/dt/sample_number ;%為頻率，單位為HZ
end


for m=1:sample_number
    %signal_y(m)=12*exp(i*2*pi*(Ap*(m-1)*dt*(m-1)*dt*(m-1)*dt+Bp*(m-1)*dt*(m-1)*dt+Cp*(m-1)*dt)) ;
    signal_y(m)=2*12*(m-1)*dt*exp(-3*(m-1)*dt)*exp(i*2*pi*(Ap*(m-1)*dt*(m-1)*dt*(m-1)*dt+Bp*(m-1)*dt*(m-1)*dt+Cp*(m-1)*dt)) ;
end


figure ;
plot(signal_t,real(signal_y)) ;
title('信號的實部時程圖','FontSize',8) ;
xlabel('時間(s)','FontSize',8) ;
ylabel('幅值','FontSize',8) ;
hold off ;
figure ;
plot(signal_t,imag(signal_y)) ;
title('信號的虛部時程圖','FontSize',8) ;
xlabel('時間(s)','FontSize',8) ;
ylabel('幅值','FontSize',8) ;
hold off ;


%求信號的功率譜
fourier_spec=fft(signal_y,sample_number) ;
power_spec=abs(fourier_spec) ;
figure ;
plot(signal_w,power_spec) ;
title('解析信號的功率圖','FontSize',8) ;
xlabel('時間(s)','FontSize',8) ;
ylabel('功率','FontSize',8) ;
hold off ;

signal_phase=unwrap(angle(signal_y)) ;
figure ;
plot(signal_t,signal_phase) ;
title('解析信號的相位圖') ;
hold off ;

for n=1:sample_number-1
    diff_phase(n)=(signal_phase(n+1)-signal_phase(n)) ;
end
figure ;
plot(signal_t(1:sample_number-1),diff_phase) ;
title('解析信號的相位差圖') ;
hold off ;

instantaneous_frequency=diff_phase./dt./(2*pi) ;
figure ;
plot(signal_t(1:sample_number-1),instantaneous_frequency) ;
title('解析信號的瞬時頻率圖') ;
hold off ;

%用連續小波分析CWT
wavelet_name='morl' ;%使用小波的名稱
center_freq=centfrq(wavelet_name) ;%小波的中心頻率
freq_vector=(1:sample_number/4)*df ;%為頻率向量，單位為HZ
scale_vector=center_freq./freq_vector/dt ; %由頻率向量算出尺度向量
disp([freq_vector',scale_vector']) ;
wavelet_coef = cwt(signal_y,scale_vector,wavelet_name,'plot') ;
%colormap(pink(64)) ;

figure ;
contour(signal_t,freq_vector,(abs(wavelet_coef)).^2,150) ;
colorbar ;
hold off 



iiiiiiiiiiiiiii=0 ;%****暫停處***






%用STFT公式(1)計算信號的時頻分布(用外文論文上的公式)
%窗用高斯窗exp(-at*t),a為窗的時間分辯率參數，a越大，窗就越窄，窗的時間分辨率就越好，而頻率分辨率就越差
win_parameter=1 ;%win_parameter為高斯窗的時間分辯率參數
for n=1:sample_number
    stft_time(n)=(n-1)*dt ;
    for m=1:sample_number
        if (n+m-1)>sample_number
            temp_stft_signal=0 ;
        else
            temp_stft_signal=signal_y(n+m-1) ;
        end
        temp_stft(m)=temp_stft_signal*exp(-win_parameter*(m-1)*dt*(m-1)*dt) ;
    end
    temp_stft_matrix=fft(temp_stft,sample_number) ;
    stft_matrix(:,n)=temp_stft_matrix' ;
end

power_stft=abs(stft_matrix) ;

for n=1:sample_number
    stft_w(n)=(n-1)/dt/sample_number ;
end

figure ;
contour(stft_time,stft_w(1:sample_number/4),power_stft(1:sample_number/4,:),100) ;
title('用STFT方法計算信號的時頻功率譜－等高線圖') ;
xlabel('時間(s)','FontSize',30) ;
ylabel('頻率(Hz)','FontSize',30) ;
hold off ;

%用維格納分布（WD）計算信號的時頻分布
for n=1:sample_number   %令t為一個常數(n-1)*dt
    wd_time(n)=(n-1)*dt ;
    for m=1:sample_number  %形成WD序列temp(m)=x(t+m/2)*x'(t-m/2)
        if mod(m,2)==0
            if (n+m/2-1)>sample_number
                temp_x1=0 ;
            else
                temp_x1=signal_y(n+m/2-1) ;
            end
            if (n-m/2-1)<1
                temp_x2=0 ;
            else
                temp_x2=signal_y(n-m/2-1) ;
            end         
        else
            if (n+(m+1)/2)>sample_number
                temp_x11=0 ;
            else
                temp_x11=signal_y(n+(m+1)/2) ;
            end
            if (n+(m-1)/2)>sample_number
                temp_x12=0 ;
            else
                temp_x12=signal_y(n+(m-1)/2) ;
            end
            temp_x1=(temp_x11+temp_x12)/2 ;
            
            
            if (n-(m+1)/2)<1
                temp_x21=0 ;
            else
                temp_x21=signal_y(n-(m+1)/2) ;
            end
            if (n-(m-1)/2)<1
                temp_x22=0 ;
            else
                temp_x22=signal_y(n-(m-1)/2) ;
            end
            temp_x2=(temp_x21+temp_x22)/2 ;
            
        end    
        x_new(m)=temp_x1*conj(temp_x2)  ;   
    end
    temp_matrix=fft(x_new,sample_number) ;
    wd_matrix(:,n)=temp_matrix' ;
end
power_wd=abs(wd_matrix) ;

for n=1:sample_number
    wd_w(n)=(n-1)/dt/sample_number ;
end

figure ;
contour(wd_time,wd_w(1:sample_number/4),power_wd(1:sample_number/4,:),100) ;
title('用WD方法計算信號的時頻功率譜－等高線圖') ;
xlabel('時間(s)','FontSize',30) ;
ylabel('頻率(Hz)','FontSize',30) ;
hold off ;