精华内容
下载资源
问答
  • Root_Music算法估计功率谱密度matlab代码,用于学习使用。
  • 基于matlab,计算信号功率谱,调用了高阶谱分析工具箱的函数 bispecd 和 bispeci计算信号的双谱
  • MUSIC算法估计功率谱密度MATLAB程序,用于学习使用的的
  • 算法允许绘制在运动图像(MI)实验中获得的α和β脑节律的功率谱图。 该算法分为三个阶段: 预处理,α和β振荡的构象以及所述节律的功率谱密度估计。 在采集脑电信号期间,MI阶段由十次试验组成,每个阶段又由...
  • 含注释,参数模型功率谱估计AR模型自相关法仿真,原理来自《数字信号处理理论、算法与实现》第三版:P545-P547
  • AR模型功率谱估计的典型算法比较及MATLAB实现.pdf
  • burg算法估计功率谱 ,完全自编 ,没用matlab自带函数
  • ESPRIT算法估计功率谱密度 matlab程序,用于学习使用的
  • MUSIC算法估计功率谱

    2018-06-12 15:01:30
    MUSIC算法估计功率谱
  • MatLab功率谱估计

    千次阅读 2018-01-20 17:24:00
    也是一种基于矩阵特征分解的谱估计非参数方法,它主要适用于混有噪声的正弦信号的功率谱估计,此方法利用相关矩阵的特征值来对MUSIC法公式中的求和进行加权得到的。 现代模型谱估计-AR模型谱估计 Yule-...

    转载:http://blog.sina.com.cn/gjchunqiu

    随机信号处理

    • 随机变量分布特征量
    • 均值mean
    • 协方差矩阵cov
    • 相关系数矩阵corrcoef
      [R, P] = corrcoef(X),P值用于检验相关性,越小越相关,0.05以下为显著相关。

    相关函数估计

    • 相关函数估计xcorr
      [c,lags] = xcorr(x,y,maxlags,’option’)
      Maxlags可以指定计算的的延迟,为[-maxlags:maxlags];
      ‘biased’: 相关函数的无偏估计
      ‘unbiased’: 相关函数的有偏估计
      ‘coeff’: 归一化相关函数,即把0延迟处的自相关系数归一化为1
      ‘none’: 使用原始非归一化相关
    • 协方差函数估计xcov
      内部执行过程为序列减去均值,再执行xcorr
    • 相关函数mscohere

    经典功率谱估计

    • 直接法(周期图法,直接FFT)periodogram
    • 改进算法
      Bartlett法
    • 功率谱估计dspdata.psd
    • 互功率谱估计cpsd
      Welch法pwelch
    • 间接法(自相关法或BT法)
    • 基于经典谱估计的系统辨识tfestimate
      使用Welch平均周期图法计算系统的谱估计

    现代谱估计-非参数法

    • MTM法pmtm
      使用正交窗口来截取获得相互独立的改进周期图法功率谱估计,然后再把这些估计结果结合得到最终的估计。随着NW的增大,窗的个数增多,会有更多的谱估计,从而谱估计的方差得到减小,但同时带来谱泄露的增大,而且正的谱估计的结果将会有更大的偏差。
    • MUSIC法pmusic
      基于矩阵特征分解的谱估计非参数方法,它把相关数据矩阵中的信息分类,把信息分配到信号的子空间或噪声的子空间。它适合于普遍情况下的正弦信号参数估计的方法,是多信号分类法的简称。
    • 特征向量法peig
      也是一种基于矩阵特征分解的谱估计非参数方法,它主要适用于混有噪声的正弦信号的功率谱估计,此方法利用相关矩阵的特征值来对MUSIC法公式中的求和进行加权得到的。

    现代模型谱估计-AR模型谱估计

    • Yule-Walker法估计pyulear
    • Burg法
    • AR模型参数估计arburg

    • AR模型功率谱估计pburg

    • 协方差谱估计pmcov

    小结

    模型谱估计适合于估计系统模型,如滤波器的谱估计;非参数法适合于正弦波信号谱估计。
    工具箱方法
    这里写图片描述
    image
    Burg—AR模型谱估计
    Covariance—AR模型谱估计
    FFT—传统谱估计
    Mod.Covar—AR模型谱估计
    MTM—现代非参数谱估计
    MUSIC—现代非参数谱估计
    Welch—传统谱估计
    Yule AR—AR模型谱估计

    展开全文
  • AR算法估计功率谱密度matlab程序,用于学习使用.。
  • ARAM功率谱估计:AR模型功率谱估计是一个全极点的模型,要利用AR模型进行功率谱估计须通过levinson—dubin递推算法由正则方程求得AR的参数...在Matlab仿真中可调用Pburg函数直接画出基于burg算法功率谱估计的曲线图。
  • 现代数字信号处理与应用 5.24关于burg算法功率谱实现的仿真实验 自己参考相关的资料写的burg算法,可以运行,结果和课本上基本一直,有些地方处理不好,比较简单就不写注释了。学习burg算法的可以参考一下
  • 功率谱估计,直接法,welch算法求取信号功率谱,包含FFT直接法估计功率谱,求取的结果和matlab的库函数pwelch完全一致,信号默认采用hamming窗,可自行修改。
  • LOMB算法进行功率谱估计,Lomb算法可以用于非均匀采样序列,自己编写的LOMB函数,经过测试,结果和Matlab的plomb一致。
  • 包含MVDR算法,RLS算法,奇异值分解MVDR算法用于功率谱估计的MATLAB仿真,里面附有详细的注释,有不理解的地方可以评论区回复,适合初学者参考
  • 随机信号大作业,陈恩庆老师的课程。完整的报告 成绩90分 经典法功率谱估计、现代法谱估计(Burg 算法、Yule-walker法、Levison-Durbin法)含误差分析 十分详细,代码有详细备注
  • AR模型功率谱估计burg算法matlab完整,直接可运行。
  • 计算功率谱

    2017-10-23 20:09:57
    用直接法得到数据的功率谱,给出的是matlab程序,可以通过修改数据进行应用
  • matlab功率谱估计

    2013-04-26 23:27:43
    matlab功率谱估计 多种函数调用算法 例子简单易懂
  • AR模型功率谱估计的典型算法比较及MATLAB实现
  • MATLAB仿真软件平台上对AR模型参数的四种不同功率谱估计算法进行了仿真,同时对功率谱估计结果进行了分析比较并得到了预期的谱估计效果。最后从实际应用角度出发讨论了AR模型参数不同功率谱估计算法的特点,以便在...
  • MVDR算法估计功率谱密度matlab程序,用于学习使用。
  • 数字信号中功率谱估计相关方法简介及MATLAB实现

    万次阅读 多人点赞 2018-12-11 17:17:35
    数字信号功率谱估计相关方法的MATLAB实现 在参阅了其他博客关于功率谱估计Matlab程序实现方法,进行重新整理、运行调试程序以及知识点总结,最终可直接运行!主要包括pmem(最大熵法)、psd(Welch法)、pmtm(多...

                                        数字信号功率谱估计相关方法的MATLAB实现

        在参阅了其他博客关于功率谱估计Matlab程序实现方法,进行重新整理、运行调试程序以及知识点总结,最终可直接运行!主要包括pmem(最大熵法)、psd(Welch法)、pmtm(多窗口法)、pmusic(多信号分类法)、分段平均周期图法(Bartlett)、

    加窗平均周期图法等的介绍,具体使用哪种方法可以根据自己程序实现结果进行选择。


    1. 基本方法

          功率谱方法可以分为两种,直接法和间接法。直接法也称为周期图法,它是直接对有限个样本数据进行傅里叶变换来得到功率谱。样本数据越长,直接法获得的分辨率越高。间接法则是先对有限个样本数据进行自相关估计,再进行傅里叶变换,最后得到功率谱。

        直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系:

                                                                         

    式中,N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有:

                                                

    其中,FFT[x(n)]为对序列x(n)的Fourier变换,由于FFT[x(n)]的周期为N,求得的功率谱估计以N为周期。

    clear;  
    fs=1000; %采样频率  
    nfft=1024; %采样的点数
    n=0:1/fs:1; 
    xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); %产生含有噪声的序列 
    window=boxcar(length(xn)); %使用矩形窗(默认为矩形窗,窗的大小与信号长度一样)  
    [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法  
    plot(f,10*log10(Pxx));  %绘制分贝形式的功率谱
    

    用有限长样本序列的Fourier变换来表示随机序列的功率谱,只是一种估计或近似,不可避免存在误差。为了减少误差,使功率谱估计更加平滑,可采用分段平均周期图法(Bartlett法)、加窗平均周期图法(Welch法)等方法加以改进。

    2. 改进方法----分段平均周期图法(Bartlett法)

           将信号序列x(n),n=0,1,…,N-1,分成互不重叠的P个小段,每小段由m个采样值,则P*m=N。对每个小段信号序列进行功率谱估计,然后再取平均作为整个序列x(n)的功率谱估计。

    %% 分段平均周期图法(Bartlett),不重叠分段功率谱估计
    
    fs=128;sigLength=1152;                        %采样频率、数据长度
    Nsec=256;n=0:sigLength-1;t=n/fs;              %分段数据点数,分段间隔,时间序列
    pxx1=abs(fft(x_sigal(1:256),Nsec).^2)/Nsec;   %第一段功率谱
    pxx2=abs(fft(x_sigal(257:512),Nsec).^2)/Nsec; %第二段功率谱
    pxx3=abs(fft(x_sigal(515:768),Nsec).^2)/Nsec; %第三段功率谱
    pxx4=abs(fft(x_sigal(769:1024),Nsec).^2)/Nsec; %第四段功率谱
    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);         %平均得到整个序列功率谱
    f=(0:length(Pxx)-1)*fs/length(Pxx);            %给出功率谱对应的频率
    figure;
    plot(f(1:Nsec/2),Pxx(1:Nsec/2));               %绘制功率谱曲线
    xlabel('频率/Hz');ylabel('功率谱 /dB');
    title('平均周期图(无重叠) N=4*256');
    grid on
    

            平均周期图法还可以对信号x(n)进行重叠分段,如按2:1重叠分段,即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计,然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法,一般后者比前者好。

    fs=128;sigLength=1152;                    %采样频率、数据长度
    Nsec=256;n=0:sigLength-1;t=n/fs;          %分段数据点数,分段间隔,时间序列
    pxx1=abs(fft(x_sigal(1:256),Nsec).^2)/Nsec;    %第一段功率谱
    pxx2=abs(fft(x_sigal(129:384),Nsec).^2)/Nsec;  %第二段功率谱
    pxx3=abs(fft(x_sigal(257:512),Nsec).^2)/Nsec;  %第三段功率谱
    pxx4=abs(fft(x_sigal(385:640),Nsec).^2)/Nsec;  %第四段功率谱
    pxx5=abs(fft(x_sigal(513:768),Nsec).^2)/Nsec;  %第五段功率谱
    pxx6=abs(fft(x_sigal(641:896),Nsec).^2)/Nsec;  %第六段功率谱
    pxx7=abs(fft(x_sigal(769:1024),Nsec).^2)/Nsec; %第七段功率谱
    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);         %平均得到整个序列功率谱
    f=(0:length(Pxx)-1)*fs/length(Pxx);            %给出功率谱对应的频率
    figure;	
    plot(f(1:Nsec/2),Pxx(1:Nsec/2));               %绘制功率谱曲线
    xlabel('频率/Hz');ylabel('功率谱 /dB');
    title('平均周期图(重叠一半) N=4*256');
    grid on
    

    3. 改进方法----加窗平均周期图法

    加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后,用非矩形窗口对每一小段信号序列进行预处理,再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识可知,采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”,同时可增加频峰的宽度,从而提高频谱分辨率。

    %% 加窗平均周期图法,不重叠分段加窗方法功率谱估计
    
    fs=128;sigLength=1152;                         %采样频率、数据长度
    w=hanning(256);                                %采用的窗口数据
    Nsec=256;n=0:sigLength-1;t=n/fs;               %分段数据点数,分段间隔,时间序列
    pxx1=abs(fft(w.*x_sigal(1:256),Nsec).^2)/norm(w)^2;   %第一段加窗振幅谱平方
    pxx2=abs(fft(w.*x_sigal(257:512),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
    pxx3=abs(fft(w.*x_sigal(513:768),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
    pxx4=abs(fft(w.*x_sigal(769:1024),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);                 %求得平均功率谱,转换为dB
    f=(0:length(Pxx)-1)*fs/length(Pxx);                    %求得频率序列
    figure;
    plot(f(1:Nsec/2),Pxx(1:Nsec/2));                       %绘制功率谱曲线
    xlabel('频率/Hz');ylabel('功率谱/dB');
    title('加窗平均周期图(无重叠) N=4*256');
    grid on
    

         加窗平均周期图法还可以对信号x(n)进行重叠分段,如按2:1重叠分段,即前一段信号和后一段信号有一半是重叠的。

    %% 加窗平均周期图法,重叠一半分段加窗功率谱估计
    fs=128;sigLength=1152;                         %采样频率、数据长度
    w=hanning(256);                                %采用的窗口数据
    Nsec=256;n=0:sigLength-1;t=n/fs;               %分段数据点数,分段间隔,时间序列
    pxx1=abs(fft(w.*x_sigal(1:256),Nsec).^2)/norm(w)^2;   %第一段加窗振幅谱平方
    pxx2=abs(fft(w.*x_sigal(129:384),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方
    pxx3=abs(fft(w.*x_sigal(257:512),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方
    pxx4=abs(fft(w.*x_sigal(385:640),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方
    pxx5=abs(fft(w.*x_sigal(513:768),Nsec).^2)/norm(w)^2; %第五段加窗振幅谱平方
    pxx6=abs(fft(w.*x_sigal(641:896),Nsec).^2)/norm(w)^2; %第六段加窗振幅谱平方
    pxx7=abs(fft(w.*x_sigal(769:1024),Nsec).^2)/norm(w)^2; %第七段加窗振幅谱平方
    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %平均功率谱转换为dB
    f=(0:length(Pxx)-1)*fs/length(Pxx);                   %求得频率序列
    figure;
    plot(f(1:Nsec/2),Pxx(1:Nsec/2));                      %绘制功率谱曲线
    xlabel('频率/Hz');ylabel('功率谱/dB');
    title('加窗平均周期图(重叠一半) N=4*256');
    grid on
    

    4. Welch法功率谱估计及其MATLAB函数  

            Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD)和两个信号序列的互功率谱估计(CSD)。

    MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计,而不需要自己编程。

    (1) 函数psd利用Welch法估计一个信号自功率谱密度,函数调用格式为:

    [Pxx,f]=psd(x,Nfft,Fs,window,Noverlap,’dflag’)

    式中,

    x为信号序列;

    Nfft为采用的FFT长度。这一值决定了功率谱估计速度,当Nfft采用2的幂时,程序采用快速算法;

    Fs为采样频率;

    Window定义窗函数和x分段序列的长度。窗函数长度必须小于或等于Nfft,否则会给出错误信息;

    Noverlap为分段序列重叠的采样点数(长度),它应小于Nfft;

    dflag为去除信号趋势分量的选择项:’linear’,去除线性趋势分量,’mean’去除均值分量,’none’不做去除趋势处理。

    Pxx为信号x的自功率谱密度估计。

    f为返回的频率向量,它和Pxx对应,并且有相同长度。

    注意:

    ①在psd函数调用格式中,缺省值为:Nfft=min(256,length(x)),Fs=2Hz, window=hanning(Nfft),noverlap=0. 若x是实序列,函数psd仅计算频率为正的功率

    ②程序前半部分中频率向量f的创建方法。它与函数psd的输出Pxx长度的关系如下:若x为实序列,当Nfft为奇数时,f=(0:(Nfft+1)/2-1)/Nfft;当Nfft为偶数时,f=(0:Nfft/2)/Nfft。

    函数还可以计算带有置信区间的功率谱估计,调用格式为:

    [Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)

    式中,p为置信区间,0<=p<=1。

    由此可知,滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频率特性。

    %% psd(Welch法)功率谱估计
    
    fs=128;sigLength=1152;                %采样频率、数据点数
    Nfft=1024;n=0:sigLength-1;t=n/128;    %傅里叶变换点数、时间序列
    window=hanning(256);                  %选用的窗口大小	
    noverlap=128;                         %分段序列重叠的采样点数(长度)
    dflag='none';                         %不做趋势处理
    [Pxx,Pxxc,f]=psd(x_sigal(:,1),Nfft,fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间,无返回值是绘制出置信区间
    figure
    plot(f,10*log10(Pxx));  %绘制功率谱
    xlabel('频率/Hz');ylabel('功率谱/dB');
    title('PSD—Welch方法'); 
    grid on
    

    利用welch法估计两个信号的互功率谱密度CSD,函数调用格式为:

        [Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)

        [Pxy,Pxyc,f]=csd(x,y,Nfft,Fs,window,Noverlap,p)

    这里,x,y为两个信号序列;Pxy为x,y的互功率谱估计;其他参数的意义同自功率谱函数psd。

    5. 多窗口法MTM功率谱估计及其MATLAB函数  

           多窗口法利用多个正交窗口(Tapers)获得各自独立的近似功率谱估计,然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法,这种功率谱估计具有更大的自由度,并在估计精度和估计波动方面均有较好的效果。普通的功率谱估计只利用单一窗口,因此在序列始端和末端均会丢失相关信息,而且无法找回。而MTM法估计增加窗口使得丢失的信息尽量减少。

          MTM法简单地采用一个参数:时间带宽积(Time-bandwidth product)NW,这个参数用以定义计算功率谱所用窗的数目,为2*NW-1。NW越大,功率谱计算次数越多,时间域分辨率越高,而频率域分辨率降低,使得功率谱估计的波动减小。随着NW增大,每次估计中谱泄漏增多,总功率谱估计的偏差增大。对于每一个数据组,通常有一个最优的NW使得在估计偏差和估计波动两方面求得折中,这需要在程序中反复调试来获得。

    MATLAB信号处理工具箱中函数PMTM就是采用MTM法估计功率谱密度。函数调用格式为:

    [Pxx,f]=pmtm(x,nw,Nfft,Fs)

    式中,x为信号序列;

    nw为时间带宽积,缺省值为4。通常可取2,5/2,3,7/2;

    Nfft为快速傅里叶变换长度;

    Fs为采样频率。

    %% pmtm(多窗口法)功率谱估计
    fs=128;sigLength=1152;                  %采样频率、数据点数
    Nfft=1024;n=0:sigLength-1;t=n/fs;       %傅里叶变换点数,时间序列
    nw=4;                                   %为时间带宽积,缺省值为4
    [Pxx1,f]=pmtm(x_sigal(:,1),nw,Nfft,fs); %用多窗口法(NW=4)估计功率谱
    figure;
    plot(f,10*log10(Pxx1));                 %绘制功率谱
    xlabel('频率/Hz');ylabel('功率谱/dB');
    title('多窗口法(MTM) 功率谱求解');
    grid on;
    

    函数还可以通过无返回值而绘出置信区间,如pmtm(x,nw,Nfft,Fs,’option’,p)绘制带置信区间的功率谱密度估计曲线,0<=p<=1。

    6 . 最大熵法pmem功率谱估计

          周期图法功率谱估计需要对信号序列“截断”或加窗处理,其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积,导致误差的产生。最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息,使估计谱的熵最大。主要方法是以已知的自相关序列rxx(0),rxx(1),…,rxx(p)为基础,外推自相关序列rxx(p+1),rxx(p+2),…,保证信息熵最大。

    最大熵功率谱估计法假定随机过程是平稳高斯过程,可以证明,随机信号的最大熵谱与AR自回归(全极点滤波器)模型谱是等价的。

    MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem,其调用格式为:

    [Pxx,f]=pmem(x,p,Nfft,Fs)

    式中,x为输入信号序列或输入相关矩阵;

    p为全极点滤波器阶次;

    %% pmem(最大熵法)功率谱估计
    
    sigLength=1152;     %数据长度
    Nfft=1024;          %傅里叶变换点数(一般为数据点数临近的2的n次方,或者直接取数据点数)
    fs=128;             %采样频率
    p=6;                %全极点滤波器阶次
    n=0:sigLength-1;
    t=n/fs;             %时间序列
    [Pxx1,f]=pmem(x_sigal(:,1),p,Nfft,fs);   %采用最大熵法,采用滤波器阶数6,估计功率谱
    figure;	
    plot(f,10*log10(Pxx1));                  %绘制功率谱
    xlabel('频率/Hz');ylabel('功率谱/dB');
    title('最大熵法 Order=20原始信号功率谱');
    grid on
    

    7 .多信号分类法功率谱估计

    MATLAB信号处理工具箱还提供另一种功率谱估计函数pmusic。该函数执行多信号分类法(multiple signal classification, Music法)。将数据自相关矩阵看成由信号自相关矩阵和噪声自相关矩阵两部分组成,即数据自相关矩阵R包含有两个子空间信息:信号子空间和噪声子空间。这样,矩阵特征值向量(Eigen vector)也可分为两个子空间:信号子空间和噪声子空间。为了求得功率谱估计,函数pmusic计算信号子空间和噪声子空间的特征值向量函数,使得在周期信号频率处函数值最大,功率谱估计出现峰值,而在其他频率处函数值最小。其调用格式为:

    [Pxx,f,a]=pmusic(x, p ,Nfft, Fs, window, Noverlap)

    式中,x为输入信号的向量或矩阵;p为信号子空间维数;其他参数的意义与函数psd相同。

    %% pmusic(多信号分类法)功率谱估计
    fs=128;sigLength=1152;             %采样频率、数据长度
    window=hanning(256);               %选用的窗口大小
    noverlap=128;                      %分段序列重叠的采样点数(长度)
    Nfft=1024;n=0:sigLength-1;t=n/fs;  %傅里叶变换点数,时间序列
    [Pxx,f,a]=pmusic(x_sigal(:,1),[7,1.1],Nfft,fs,window,noverlap);   %采用多信号分类法估计功率谱
    figure;
    plot(f,10*log10(Pxx));                                            %绘制功率谱
    xlabel('频率/Hz'); ylabel('功率谱/dB')
    title('通过MUSIC法估计的伪谱')
    

    特别说明:x_sigal(:,1)我使用的脑电运动想象数据作为输入数据

    供大家参考网址:

     https://wenku.baidu.com/view/7046318b0b1c59eef9c7b451.html

     

     

    展开全文
  • 随机信号处理大作业,用Burg算法求解该AR模型的参数程序,并画出p取值分别为3、4、5时的功率谱曲线,效果还可以。
  • 功率谱密度相关方法的MATLAB实现

    万次阅读 2016-07-11 19:30:32
    1. 基本方法 ...周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系:   式中,N

    转自http://blog.sina.com.cn/s/blog_79ecf6980100vcqn.html


    1基本方法

    周期图法是直接将信号的采样数据x(n)进行Fourier变换求取功率谱密度估计的方法。假定有限长随机信号序列为x(n)。它的Fourier变换和功率谱密度估计存在下面的关系:

          功率谱密度相关方法的MATLAB实现                    

    式中,N为随机信号序列x(n)的长度。在离散的频率点f=kΔf,有:

             功率谱密度相关方法的MATLAB实现 

    其中,FFT[x(n)]为对序列x(n)的Fourier变换,由于FFT[x(n)]的周期为N,求得的功率谱估计以N为周期,因此这种方法称为周期图法。下面用例子说明如何采用这种方法进行功率谱

    用有限长样本序列的Fourier变换来表示随机序列的功率谱,只是一种估计或近似,不可避免存在误差。为了减少误差,使功率谱估计更加平滑,可采用分段平均周期图法(Bartlett法)、加窗平均周期图法(Welch法)等方法加以改进。

    2. 分段平均周期图法(Bartlett法)

    将信号序列x(n),n=0,1,…,N-1,分成互不重叠的P个小段,每小段由m个采样值,则P*m=N。对每个小段信号序列进行功率谱估计,然后再取平均作为整个序列x(n)的功率谱估计。

    平均周期图法还可以对信号x(n)进行重叠分段,如按2:1重叠分段,即前一段信号和后一段信号有一半是重叠的。对每一小段信号序列进行功率谱估计,然后再取平均值作为整个序列x(n)的功率谱估计。这两种方法都称为平均周期图法,一般后者比前者好。程序运行结果为图9-5,上图采用不重叠分段法的功率谱估计,下图为2:1重叠分段的功率谱估计,可见后者估计曲线较为平滑。与上例比较,平均周期图法功率谱估计具有明显效果(涨落曲线靠近0dB)。

    3.加窗平均周期图法

    加窗平均周期图法是对分段平均周期图法的改进。在信号序列x(n)分段后,用非矩形窗口对每一小段信号序列进行预处理,再采用前述分段平均周期图法进行整个信号序列x(n)的功率谱估计。由窗函数的基本知识(第7章)可知,采用合适的非矩形窗口对信号进行处理可减小“频谱泄露”,同时可增加频峰的宽度,从而提高频谱分辨率。

    其中上图采用无重叠数据分段的加窗平均周期图法进行功率谱估计,而下图采用重叠数据分段的加窗平均周期图法进行功率谱估计,显然后者是更佳的,信号谱峰加宽,而噪声谱均在0dB附近,更为平坦(注意采用无重叠数据分段噪声的最大的下降分贝数大于5dB,而重叠数据分段周期图法噪声的最大下降分贝数小于5dB)。

    4.  Welch法估计及其MATLAB函数

    Welch功率谱密度就是用改进的平均周期图法来求取随机信号的功率谱密度估计的。Welch法采用信号重叠分段、加窗函数和FFT算法等计算一个信号序列的自功率谱估计(PSD如上例中的下半部分的求法)和两个信号序列的互功率谱估计(CSD)。

    MATLAB信号处理工具箱函数提供了专门的函数PSD和CSD自动实现Welch法估计,而不需要自己编程。

    (1) 函数psd利用Welch法估计一个信号自功率谱密度,函数调用格式为:

    [Pxx[,f]]=psd(x[,Nfft,Fs,window,Noverlap,’dflag’])

    式中,x为信号序列;Nfft为采用的FFT长度。这一值决定了功率谱估计速度,当Nfft采用2的幂时,程序采用快速算法;Fs为采样频率;Window定义窗函数和x分段序列的长度。窗函数长度必须小于或等于Nfft,否则会给出错误信息;Noverlap为分段序列重叠的采样点数(长度),它应小于Nfft;dflag为去除信号趋势分量的选择项:’linear’,去除线性趋势分量,’mean’去除均值分量,’none’不做去除趋势处理。Pxx为信号x的自功率谱密度估计。f为返回的频率向量,它和Pxx对应,并且有相同长度。

    在psd函数调用格式中,缺省值为:Nfft=min(256,length(x)),Fs=2Hz, window=hanning(Nfft),noverlap=0. 若x是实序列,函数psd仅计算频率为正的功率

    注意程序前半部分中频率向量f的创建方法。它与函数psd的输出Pxx长度的关系如下:若x为实序列,当Nfft为奇数时,f=(0:(Nfft+1)/2-1)/Nfft;当Nfft为偶数时,f=(0:Nfft/2)/Nfft。

    函数还有一种缺省返回值的调用格式,用于直接绘制信号序列x的功率谱估计曲线。

    函数还可以计算带有置信区间的功率谱估计,调用格式为:

    [Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,Noverlap,p)

    式中,p为置信区间,0<=p<=1

     

    由此可知,滤波器输入白噪声序列的输出信号的功率谱或自相关可以确定滤波器的频率特性。

    (2)函数csd利用welch法估计两个信号的互功率谱密度,函数调用格式为:

        [Pxy[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,’dflag’)

        [Pxy,Pxyc[,f]]=csd(x,y,Nfft,Fs,window,Noverlap,p)

    这里,x,y为两个信号序列;Pxy为x,y的互功率谱估计;其他参数的意义同自功率谱函数psd。

    可以看到,两个白噪声信号的互功率谱(上图)杂乱无章,看不出周期成分,大部分功率谱在-5dB以下。然而白噪声与带有噪声的周期信号的功率谱在其周期(频率为1000Hz)处有一峰值,清楚地表明了周期信号的周期或频率。因此,利用未知信号与白噪声信号的互功率谱也可以检测未知信号中所含有的频率成分。

     

    5 多 窗 口 

    多窗口法(Multitaper method,简称MTM法)利用多个正交窗口(Tapers)获得各自独立的近似功率谱估计,然后综合这些估计得到一个序列的功率谱估计。相对于普通的周期图法,这种功率谱估计具有更大的自由度,并在估计精度和估计波动方面均有较好的效果。普通的功率谱估计只利用单一窗口,因此在序列始端和末端均会丢失相关信息,而且无法找回。而MTM法估计增加窗口使得丢失的信息尽量减少。

    MTM法简单地采用一个参数:时间带宽积(Time-bandwidth product)NW,这个参数用以定义计算功率谱所用窗的数目,为2*NW-1。NW越大,功率谱计算次数越多,时间域分辨率越高,而频率域分辨率降低,使得功率谱估计的波动减小。随着NW增大,每次估计中谱泄漏增多,总功率谱估计的偏差增大。对于每一个数据组,通常有一个最优的NW使得在估计偏差和估计波动两方面求得折中,这需要在程序中反复调试来获得。

    MATLAB信号处理工具箱中函数PMTM就是采用MTM法估计功率谱密度。函数调用格式为:

    [Pxx[,f]]=pmtm(x[,nw,Nfft,Fs])

    式中,x为信号序列;nw为时间带宽积,缺省值为4。通常可取2,5/2,3,7/2;Nfft为FFT长度;Fs为采样频率。

    上面的函数还可以通过无返回值而绘出置信区间,如pmtm(x,nw,Nfft,Fs,’option’,p)绘制带置信区间的功率谱密度估计曲线,0<=p<=1

    6 最 大 熵 法(Maxmum entropy method, MEM法

    如上所述,周期图法功率谱估计需要对信号序列“截断”或加窗处理,其结果是使估计的功率谱密度为信号序列真实谱和窗谱的卷积,导致误差的产生。

    最大熵功率谱估计的目的是最大限度地保留截断后丢失的“窗口”以外信号的信息,使估计谱的熵最大。主要方法是以已知的自相关序列rxx(0),rxx(1),…,rxx(p)为基础,外推自相关序列rxx(p+1),rxx(p+2),…,保证信息熵最大。

    最大熵功率谱估计法假定随机过程是平稳高斯过程,可以证明,随机信号的最大熵谱与AR自回归(全极点滤波器)模型谱是等价的。

    MATLAB信号处理工具箱提供最大熵功率谱估计函数pmem,其调用格式为:

    [Pxx,f,a]=pmem(x,p,Nfft,Fs,’xcorr’)

    式中,x为输入信号序列或输入相关矩阵;p为全极点滤波器阶次;a为全极点滤波器模型系数向量;’xcorr’是把x认为是相关矩阵。

     

    比较最大熵功率谱估计(MEM)和改进的平均周期图功率谱估计,可见,MEM法估计的功率谱曲线较光滑。在这一方法中,MEM法选定全极点滤波器的阶数取得越大,能够获得的窗口外的信息越多,但计算量也越大,需要根据情况折中考虑。

     

    7 多信号分类法

    MATLAB信号处理工具箱还提供另一种功率谱估计函数pmusic。该函数执行多信号分类法(multiple signal classification, Music法)。将数据自相关矩阵看成由信号自相关矩阵和噪声自相关矩阵两部分组成,即数据自相关矩阵R包含有两个子空间信息:信号子空间和噪声子空间。这样,矩阵特征值向量(Eigen vector)也可分为两个子空间:信号子空间和噪声子空间。为了求得功率谱估计,函数pmusic计算信号子空间和噪声子空间的特征值向量函数,使得在周期信号频率处函数值最大,功率谱估计出现峰值,而在其他频率处函数值最小。其调用格式为:

    [Pxx,f,a]=pmusic(x,p[[,thresh],Nfft,Fs,window,Noverlap])

    式中,x为输入信号的向量或矩阵;p为信号子空间维数;thresh为阈值,其他参数的意义与函数psd相同。

     

    功率谱密度相关方法的MATLAB实现

    %%

    %分段平均周期图法(Bartlett法)

    %运用信号不重叠分段估计功率谱 

    Nsec=256;n=0:sigLength-1;t=n/Fs;  %数据点数,分段间隔,时间序列

    pxx1=abs(fft(y(1:256),Nsec).^2)/Nsec;   %第一段功率谱

    pxx2=abs(fft(y(257:512),Nsec).^2)/Nsec; %第二段功率谱

    pxx3=abs(fft(y(515:768),Nsec).^2)/Nsec; %第三段功率谱

    pxx4=abs(fft(y(769:1024),Nsec).^2)/Nsec; %第四段功率谱

    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);   %平均得到整个序列功率谱

    f=(0:length(Pxx)-1)*Fs/length(Pxx);     %给出功率谱对应的频率

    %%plot(f(1:Nsec/2),Pxx(1:Nsec/2));       %绘制功率谱曲线

    figure,plot(f(1:Nsec),Pxx(1:Nsec));       %绘制功率谱曲线

    xlabel('频率/Hz');ylabel('功率谱 /dB');

    title('平均周期图(无重叠) N=4*256');

    grid on

    %%

    %运用信号重叠分段估计功率谱 

    pxx1=abs(fft(y(1:256),Nsec).^2)/Nsec;   %第一段功率谱

    pxx2=abs(fft(y(129:384),Nsec).^2)/Nsec; %第二段功率谱

    pxx3=abs(fft(y(257:512),Nsec).^2)/Nsec; %第三段功率谱

    pxx4=abs(fft(y(385:640),Nsec).^2)/Nsec; %第四段功率谱

    pxx5=abs(fft(y(513:768),Nsec).^2)/Nsec; %第五段功率谱

    pxx6=abs(fft(y(641:896),Nsec).^2)/Nsec; %第六段功率谱

    pxx7=abs(fft(y(769:1024),Nsec).^2)/Nsec; %第七段功率谱

    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7); %功率谱平均并转化为dB

    f=(0:length(Pxx)-1)*Fs/length(Pxx);  %频率序列

    figure,plot(f(1:Nsec/2),Pxx(1:Nsec/2));   %绘制功率谱曲线

    xlabel('频率/Hz');  ylabel('功率谱/dB');

    title('平均周期图(重叠一半) N=1024');

    grid on

     

    %%

    Nsec=256;n=0:sigLength-1;t=n/Fs; %数据长度,分段数据长度、时间序列

    w=hanning(256);   %采用的窗口数据

    %采用不重叠加窗方法的功率谱估计

    pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方

    pxx2=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方

    pxx3=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方

    pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方

    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4);    %求得平均功率谱,转换为dB

    f=(0:length(Pxx)-1)*Fs/length(Pxx);  %求得频率序列

    figure

    subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线

    xlabel('频率/Hz');ylabel('功率谱/dB');

    title('加窗平均周期图(无重叠) N=4*256');

    grid on

    %采用重叠加窗方法的功率谱估计

    pxx1=abs(fft(w.*y(1:256),Nsec).^2)/norm(w)^2; %第一段加窗振幅谱平方

    pxx2=abs(fft(w.*y(129:384),Nsec).^2)/norm(w)^2; %第二段加窗振幅谱平方

    pxx3=abs(fft(w.*y(257:512),Nsec).^2)/norm(w)^2; %第三段加窗振幅谱平方

    pxx4=abs(fft(w.*y(385:640),Nsec).^2)/norm(w)^2; %第四段加窗振幅谱平方

    pxx5=abs(fft(w.*y(513:768),Nsec).^2)/norm(w)^2; %第五段加窗振幅谱平方

    pxx6=abs(fft(w.*y(641:896),Nsec).^2)/norm(w)^2; %第六段加窗振幅谱平方

    pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/norm(w)^2; %第七段加窗振幅谱平方

    Pxx=10*log10((pxx1+pxx2+pxx3+pxx4+pxx5+pxx6+pxx7)/7);%平均功率谱转换为dB

    f=(0:length(Pxx)-1)*Fs/length(Pxx); %频率序列

    subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2)); %绘制功率谱曲线

    xlabel('频率/Hz');ylabel('功率谱/dB');

    title('加窗平均周期图(重叠一半)N=1024');

    grid on

    %%

    %4分段平均周期图法(hanning窗)

    Nsec=256;

    n=0:sigLength-1;

    t=n/Fs;

    w=hanning(256);

    Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;

    Pxx2=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;

    Pxx3=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;

    Pxx4=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;

    Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4)/4);

    f=(0:length(Pxx)-1)*Fs/length(Pxx);

    figure

    subplot(2,1,1)

    plot(f,Pxx);

    xlabel('频率/Hz');ylabel('功率谱/dB');

    title('Averaged Modified Periodogram (none overlap) N=4*256');

    grid

    %4分段(2:1重叠)平均周期图法(hanning窗)

    Nsec=256;

    n=0:sigLength-1;

    t=n/Fs;

    w=hanning(256);

    Pxx1=abs(fft(w.*y(1:256),Nsec).^2)/Nsec;

    Pxx2=abs(fft(w.*y(129:384),Nsec).^2)/Nsec;

    Pxx3=abs(fft(w.*y(257:512),Nsec).^2)/Nsec;

    Pxx4=abs(fft(w.*y(385:640),Nsec).^2)/Nsec;

    Pxx5=abs(fft(w.*y(513:768),Nsec).^2)/Nsec;

    Pxx6=abs(fft(w.*y(641:896),Nsec).^2)/Nsec;

    Pxx7=abs(fft(w.*y(769:1024),Nsec).^2)/Nsec;

    Pxx=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);

    f=(0:length(Pxx)-1)*Fs/length(Pxx);

    subplot(2,1,2);plot(f,Pxx);

    xlabel('频率/Hz');ylabel('Power Spectrum (dB)');

    title('Averaged Modified Periodogram (half overlap) N=1024');

    grid

    %%

    %PSD_WELCH方法

    %采样频率

    Nfft=256;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列

    window=hanning(256);   %选用的窗口

    noverlap=128;   %分段序列重叠的采样点数(长度)

    dflag='none' %不做趋势处理

    [Pxx,Pxxc,f]=psd(y,Nfft,Fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间,无返回值是绘制出置信区间

    figure

    plot(f,10*log10(Pxx));  %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');

    title('PSD—Welch方法'); grid on

        

     

    %%

    %最大熵法(MEM法)

    Nfft=256;n=0:sigLength-1;t=n/Fs;   %数据长度、分段长度和时间序列

    window=hanning(256);   %采用窗口

    [Pxx1,f]=pmem(x,20,Nfft,Fs);   %采用最大熵法,采用滤波器阶数14,估计功率谱

    figure,subplot(2,1,1),plot(f,10*log10(Pxx1));   %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');

    title('最大熵法 Order=20原始信号功率谱');grid on

    [Pxx1,f]=pmem(y0,20,Nfft,Fs);   %采用最大熵法,采用滤波器阶数14,估计功率谱

    subplot(2,1,2),plot(f,10*log10(Pxx1));   %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -20 0])

    title('最大熵法 Order=20滤波后的信号功率谱');grid on

    %%

    %%功率谱密度

    %PSD_WELCH方法

    Nfft=512;n=0:sigLength-1;t=n/Fs; %数据长度、时间序列

    window=hanning(256);   %选用的窗口

    noverlap=128;   %分段序列重叠的采样点数(长度)

    dflag='none' %不做趋势处理

    [Pxx,Pxxc,f]=psd(x,Nfft,Fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间,无返回值是绘制出置信区间

    figure;subplot(211);

    plot(f,10*log10(Pxx));  %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');

    grid on;title('PSD—Welch方法的原始信号功率谱')

    subplot(212)

    [Pxx,Pxxc,f]=psd(y0,Nfft,Fs,window,noverlap,0.95);   %功率谱估计,并以0.95的置信度给出置信区间,无返回值是绘制出置信区间

    plot(f,10*log10(Pxx));  %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -30 0])

    grid on;title('PSD—Welch方法的滤波后的信号功率谱')

    %%

    %用多窗口法(MTM)

    n=0:sigLength-1;t=n/Fs;  %数据长度、分段数据长度,时间序列

    [Pxx1,f]=pmtm(x,2,Nfft,Fs); %用多窗口法(NW=4)估计功率谱

    figure;subplot(2,1,1),plot(f,10*log10(Pxx1));  %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');

    title('多窗口法(MTM) nw=2原始信号功率谱');grid on

    [Pxx,f]=pmtm(y0,2,Nfft,Fs); %用多窗口法(NW=2)估计功率谱

    subplot(2,1,2),plot(f,10*log10(Pxx));  %绘制功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB');axis([0 4000 -80 -20])

    title('多窗口法(MTM) nw=2滤波后的信号功率谱');grid on

        

    %%

    %采用Welch方法估计功率谱

    noverlap=128;    %重叠数据

    dflag='none'Nfft=1024;

    figure;subplot(2,1,1)

    psd(x,Nfft,Fs,window,noverlap,dflag);  %采用Welch方法估计功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB')

    title('Welch方法原始信号功率谱');grid on

    subplot(2,1,2)

    psd(y0,Nfft,Fs,window,noverlap,dflag);  %采用Welch方法估计功率谱

    xlabel('频率/Hz');ylabel('功率谱/dB'),axis([0 4000 -30 0])

    title('Welch方法滤波后的信号功率谱');grid on

    pmusic(y0,[7,1.1],Nfft,Fs,32,16);   %采用多信号分类法估计功率谱

    xlabel('频率/Hz'); ylabel('功率谱/dB')

    title('通过MUSIC法估计的伪谱')


    展开全文
  • 2ask的功率谱密度matlab代码 Automatic-Modulation-Classification 使用matlab生成数据集 参考 仿真方式 MATLAB(版本2019b以上) 调制种类 28种,数字调制(25种): "BPSK", "QPSK", "8PSK","16PSK","32PSK",....
  • 功率谱估计

    2018-03-15 19:37:32
    用burg算法实现功率谱估计,附程序,可运行。现代信号处理
  • 利用周期图法进行估计,并绘制结果,窗函数采用矩形窗。利用Levinson-Durbin递推法求解Yule-walker方程,进行AR(6)的建模。与Matlab中periodogram(周期图)和pyulear(Yule-walker方程)中相应方法的结果进行比较...
  • 用SVDTLS算法实现功率谱估计,现代信号处理。可直接运行

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,371
精华内容 548
关键字:

matlab功率谱算法

matlab 订阅