下面的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 matrixnoverlap=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 cspdoverlap_rate=0.997;% define the window overlap rate.overlap_number=round(overlap_rate*window_size);% define the 85% overlap of the windowsig1_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 signalwindow_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 defaultpxy=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']);