时频分析：窗口傅立叶变换

下面的matlab代码通过窗口傅立叶变换计算了信号的功率谱和互功率谱(power spectrum for time frequency analysis)

% this code will test the buildin cpsd function and  calculate windowed FAST Fourier Transformation
% by DocNan 
clear all; close all; clc;
tic;

t=(2.5:0.00001:3.5)';
f=(linspace(0,211000,length(t)))';
phi=cumtrapz(t,2*pi*f);

sig1=3*sin(phi)+5*randn(length(t),1);
%sig2=3*sin(phi+1*pi/8.0*0)+15*randn(length(t),1);
sig2=5*randn(length(t),1);


% prepare to buffer the original signal, this step is likely to cut original signal into many small windows and form a large matrix
noverlap=round(0.8*256);
Fs=round(1/(t(2)-t(1)));
window_size=256;%(2,4,8,16,32,64,128,256,512,1024,2048,4096,8192) define the size of window function, small section to do cspd
overlap_rate=0.997;% define the window overlap rate.
overlap_number=round(overlap_rate*window_size);% define the 85% overlap of the window
sig1_buffer=buffer(sig1,window_size,overlap_number); % get the overlapped time window matrix.
sig2_buffer=buffer(sig2,window_size,overlap_number); % get the overlapped time window matrix.


% apply the smooth window to the buffered signal
window_matrix=gausswin(window_size)*ones(1,size(sig1_buffer,2));
sig1_buffer=sig1_buffer.*window_matrix;
sig2_buffer=sig2_buffer.*window_matrix;


% calculate power spectrum, only keep the part with positive or zero frequency(half the spectrum)
pxx=fft(sig1_buffer);
pxx=2*pxx(1:size(pxx,1)/2+1,:);

% calculate cross power spectrum, cpsd will remove the negative part of spectrum by default
pxy=cpsd(sig2_buffer,sig1_buffer,window_size);


t_new=linspace(t(1),t(end),size(pxy,2));
f=linspace(0,Fs/2,size(pxy,1));


fig1=figure();
pcolor(t_new,f/1000,10*log10(abs(pxy)));
shading flat;
ylabel('f(kHz)');
xlabel('t(s)');
legend('CPSD');
colorbar;
colormap(gca,jet);


fig2=figure();
pcolor(t_new,f/1000,10*log10(abs(pxx)));
shading flat;
ylabel('f(kHz)');
xlabel('t(s)');
legend('PSD');
colorbar;
colormap(gca,jet);

toc;
disp(['cost time=',num2str(toc),'s']);