zoukankan      html  css  js  c++  java
  • 时频分析:窗口傅立叶变换

    下面的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']);
    
    
  • 相关阅读:
    [Bzoj2243][SDOI2011]染色(线段树&&树剖||LCT)
    [poj3074]Sudoku(舞蹈链)
    [Bzoj1047][HAOI2007]理想的正方形(ST表)
    [Bzoj1030][JSOI2007]文本生成器(AC自动机&dp)
    [Bzoj2431][HAOI2009]逆序对数列(前缀和优化dp)
    [Bzoj1072][SCOI2007]排列perm(状压dp)
    [Bzoj1195][HNOI2006]最短母串(AC自动机)
    Ajax解决IE浏览器兼容问题
    运行eclipse弹出“Failed to load the JNI shared”解决方法
    Java表单类双击提交
  • 原文地址:https://www.cnblogs.com/docnan/p/7050860.html
Copyright © 2011-2022 走看看