MATLAB数字信号处理(1)四种经典功率谱估计方法比较

数字信号处理 专栏收录该内容
51 篇文章 87 订阅

这是我研究生课程“现代信号处理”中的作业报告,上传到blog中。


经典功率谱估计

可以采用直接法,也称周期图法,利用公式在这里插入图片描述计算功率谱密度。或者根据自相关函数和谱密度之间的傅里叶变换关系在这里插入图片描述
来计算,称为间接法或自相关函数法。

还可以先作加窗平滑处,对序列x(n)或估计的自相关函数进行加窗(如汉宁窗、汉明窗)截断,前者称作数据窗,后者称作滞后窗


MATLAB编程实现

对信号x(n)=sin⁡(ωt)+n(t)和x(n)=exp⁡(jωt)+n(t)使用上述4种方法进行功率谱估计。叠加高斯白噪声n(t)后的信噪比为-10dB、-3dB、3dB、10dB四种情况。编程不使用fft、xcorr等MATLAB提供的函数。

关键程序如下,四个函数分别使用不同方法进行功率谱估计,一个函数用于绘图。调用自己编写的函数实现功能(MATLAB R2017a版本下测试)。

N = 512; fs = 1; M=256;
t = (0:N-1)/fs;
x1n = sin(2*pi*0.1*t);           %信号1,实信号
x2n = exp(1j*2*pi*0.1*t);      %信号2,复信号

P_estimation(x1n, 10, N, M, fs);    %估计不同SNR、不同信号的功率谱
P_estimation(x1n, 3, N, M, fs);      
P_estimation(x1n, -3, N, M, fs);    
P_estimation(x1n, -10, N, M, fs);    
P_estimation(x2n, 10, N, M, fs);    
P_estimation(x2n, 3, N, M, fs);    
P_estimation(x2n, -3, N, M, fs);   
P_estimation(x2n, -10, N, M, fs);    

%   对叠加SNR的信号xn,使用方法1-4进行功率谱估计
function P_estimation(xn, snr, N, M, fs)
    sn = awgn(xn, snr, 'measured');
    f1=(0:N-1)*fs/N;                        %功率谱坐标轴
    Sx1 = direct_method(sn, N);     % 方法1:直接法  
    Sx2 = indir_method(sn,N,M);    % 方法2:间接法  
    Sx3 = add_datawin(sn, N);        % 方法3:数据加窗 
    Sx4 = add_rxwin(sn,N,M);         % 方法4:自相关函数加窗
    figure;
    subplot(221); plot(f1,10*log10(Sx1));
    title('直接法 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');
    subplot(222); plot(f1,10*log10(Sx2));
    title('间接法 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');
    subplot(223); plot(f1,10*log10(Sx3));
    title('数据加窗 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');
    subplot(224); plot(f1,10*log10(Sx4));
    title('自相关函数加窗 SNR='+string(snr)); xlabel('f/Hz'); ylabel('Sx(f)(dB/Hz)');    
end

%    方法1:直接法得到频率谱估计
function Sx = direct_method(xn, N)   %估计N点序列xn的功率谱Sx
    Sx = zeros(1,N);
    for index = 1 : N
        for n = 0 : N-1       %求Fourier变换
            Sx(index) = Sx(index) + xn(n+1)*exp(-1j*n*index*2*pi/N);  
        end
        Sx(index) = abs(Sx(index))*abs(Sx(index))/N;  %求功率谱
    end
end

%    方法2:间接法得到频率谱估计
function Sx = indir_method(xn,N,M)   %估计N点序列xn的功率谱Sx    
    Rx = zeros(1,2*M-1);   %M为估计的自相关函数的点数,1<<M<N
    xn = [xn, zeros(1,M)];
    for i = 1 : M                 %估计自相关函数
        for n = 1 : N
            Rx(i+M-1) = Rx(i+M-1) + xn(n+i)*conj(xn(n));
        end
        Rx(i+M-1) = Rx(i+M-1) / N;
    end
    for i = 1 : M-1             %根据共轭关系得到另一半
        Rx(i) = conj(Rx(2*M-i));
    end
    Sx = zeros(1,N);
    for index = 1 : N
        for k = 1 : 2*M-1        %求Fourier变换
            Sx(index) = Sx(index) + Rx(k)*exp(-1j*(k-M)*index*2*pi/N);  
        end
        Sx(index) = abs(Sx(index));
    end
end

%   方法3:数据加窗-修正周期图
	     在方法1基础上增加了加窗语句,其余重复故省略
%  方法4:自相关函数加窗-周期图平滑
在方法2基础上增加了加窗语句,其余重复故省略

结果分析

下面是不同信噪比下对x(n)=sin⁡(ωt)+n(t)信号的功率谱估计结果。
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
程序中采样频率为1,信号频率为0.1。实信号的功率谱是对称的,在图中可以看到两个尖峰在0.1和与之对称的0.9处。对于直接法而言,数据加窗后的“尖峰”更窄;对于间接法而言,自相关函数加窗后的功率谱图更平滑。这也是其称作“修正周期图”和“周期图平滑”的直观体验。

下面是不同信噪比下对x(n)=exp⁡(jωt)+n(t)信号的功率谱估计结果。
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
可以得到和上一组测试相同的结论。只不过复信号的功率谱不具有对称性。总体而言,信噪比越高,信号的功率峰也就越突出。


不同窗函数比较

下面给出直接法和间接法中,不加窗、汉明窗、汉宁窗、布莱克曼窗的差别。为了更直观地观察使用不同窗函数之间的区别,没有添加噪声。直接法中汉明窗将功率峰变得更窄;汉宁窗和布莱克曼窗效果类似,增大了功率峰值和其余功率值之间的倍数。
在这里插入图片描述
间接法中加三个窗函数都起到了周期图平滑的作用,汉宁窗和布莱克曼窗的平滑效果看起来比汉明窗要好一些。但三个窗函数都没有改变功率峰值和其余功率值之间的倍数关系。
在这里插入图片描述

3种经典功率谱估计方法MATLABA代码-功率谱代码.doc 3种MATLAB经典谱估计方法 希望对大家有用~ 附件所有代码: 直接法: 直接法又称周期图法,它是把随机序列x的N个观测数据视为一能量有限的序列,直接计算x的离散傅立叶变换,得X,然后再取其幅值的平方,并除以N,作为序列x真实功率谱估计Matlab代码: clear; Fs=1000; %采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos 3*cos randn); window=boxcar); %矩形窗 nfft=1024; [Pxx,f]=periodogram; %直接法 plot); 改进的直接法: 对于直接法的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。 1. Bartlett法 Bartlett平均周期图的方法是将N点的有限长序列x分段求周期图再平均。 Matlab代码: clear; Fs=1000; n=0:1/Fs:1; xn=cos 3*cos randn); nfft=1024; window=boxcar); %矩形窗 noverlap=0; %数据无重叠 p=0.9; %置信概率 [Pxx,Pxxc]=psd; index=0:round; k=index*Fs/nfft; plot_Pxx=10*log10); plot_Pxxc=10*log10); figure plot; pause; figure plot; 2. Welch法 Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w,并再周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。 Matlab代码: clear; Fs=1000; n=0:1/Fs:1; xn=cos 3*cos randn); nfft=1024; window=boxcar; %矩形窗 window1=hamming; %汉明窗 window2=blackman; %blackman窗 noverlap=20; %数据无重叠 range='half'; %频率间隔为[0 Fs/2],只计算一半的频率 [Pxx,f]=pwelch; [Pxx1,f]=pwelch; [Pxx2,f]=pwelch; plot_Pxx=10*log10; plot_Pxx1=10*log10; plot_Pxx2=10*log10; figure plot; pause; figure plot; pause; figure plot; 复制代码
©️2021 CSDN 皮肤主题: 酷酷鲨 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值