zoukankan      html  css  js  c++  java
  • 用matlab做经典功率谱估计

    经典功率谱估计
    直接法:
    直接法又称周期图法,它是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换,得X(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。
    Matlab代码示例:
    clear;
    Fs=1000; %采样频率
    n=0:1/Fs:1;
    %产生含有噪声的序列
    xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
    window=boxcar(length(xn)); %矩形窗
    nfft=1024;
    [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法
    plot(f,10*log10(Pxx));

    间接法:
    间接法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。
    Matlab代码示例:
    clear;
    Fs=1000; %采样频率
    n=0:1/Fs:1;
    %产生含有噪声的序列
    xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
    nfft=1024;
    cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
    CXk=fft(cxn,nfft);
    Pxx=abs(CXk);
    index=0:round(nfft/2-1);
    k=index*Fs/nfft;
    plot_Pxx=10*log10(Pxx(index+1));
    plot(k,plot_Pxx);

    改进的直接法:
    对于直接法的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。
    1. Bartlett法
    Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
    Matlab代码示例:
    clear;
    Fs=1000;
    n=0:1/Fs:1;
    xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
    nfft=1024;
    window=boxcar(length(n)); %矩形窗
    noverlap=0; %数据无重叠
    p=0.9; %置信概率
    [Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
    index=0:round(nfft/2-1);
    k=index*Fs/nfft;
    plot_Pxx=10*log10(Pxx(index+1));
    plot_Pxxc=10*log10(Pxxc(index+1));
    figure(1)
    plot(k,plot_Pxx);
    pause;
    figure(2)
    plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

    2. Welch法
    Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w(n),并再周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。
    Matlab代码示例:
    clear;
    Fs=1000;
    n=0:1/Fs:1;
    xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));
    nfft=1024;
    window=boxcar(100); %矩形窗
    window1=hamming(100); %海明窗
    window2=blackman(100); %blackman窗
    noverlap=20; %数据无重叠
    range='half'; %频率间隔为[0 Fs/2],只计算一半的频率
    [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);
    [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
    [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
    plot_Pxx=10*log10(Pxx);
    plot_Pxx1=10*log10(Pxx1);
    plot_Pxx2=10*log10(Pxx2);

    figure(1)
    plot(f,plot_Pxx);

    pause;

    figure(2)
    plot(f,plot_Pxx1);

    pause;

    figure(3)
    plot(f,plot_Pxx2);

  • 相关阅读:
    详解GaussDB(for MySQL)服务:复制策略与可用性分析
    华为云的研究成果又双叒叕被MICCAI收录了!
    充分释放数据价值:安全、可信6到飞起
    未来云原生世界的“领头羊”:容器批量计算项目Volcano 1.0版本发布
    一文带你掌握OBS的两种常见的鉴权方式
    数据库实践丨MySQL多表join分析
    技术贴丨教你使用华为云鲲鹏服务器部署Discuz!论坛
    python Scrapy 从零开始学习笔记(二)
    python Scrapy 从零开始学习笔记(一)
    从零开始学Electron笔记(七)
  • 原文地址:https://www.cnblogs.com/yunbo/p/747769.html
Copyright © 2011-2022 走看看