功率谱一点介绍

2024-03-23 18:24:16

本篇为《信号处理》系列博客的第六篇，该系列博客主要记录信号处理相关知识的学习过程和自己的理解，方便以后查阅。

功率谱一点介绍

经典功率谱估计
现代谱估计法
三种方法对比结果

信号的功率谱密度描述随机信号的功率在频域随频率的分布。

利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计，功率谱密度一般简称功率谱。

谱估计方法分为参数化方法和非参数化方法。非参数化方法又叫经典谱估计，如周期图法、自相关法等，其主要缺点是描述功率谱波动的数字特征方差性能较差，频率分辨率低；而参数化谱估计又叫做现代谱估计，如AR模型法、MA模型法、自回归移动平均模型法（ARMA模型法）等。

功率谱的单位是W/Hz,如果做了对数处理10log，就是分贝（dB）。

经典功率谱估计

经典功率谱估计是截取较长的数据链中的一段作为工作区，相当于将数据加一个矩形窗函数。

根据截取的N个样本数据用估计出其功率谱。其中可以利用相关函数法估计功率谱、也可以利用周期图法估计出功率谱。

根据自相关函数计算功率谱

先计算出自相关函数
R x ( m ) = 1 N ∑ n = 0 N − m − 1 x ( n ) x ( n + m ) \large R_{x}(m)=\frac{1}{N}\sum_{n=0}^{N-m-1}x(n)x(n+m) Rx(m)=N1n=0∑N−m−1x(n)x(n+m)
对自相关函数做傅里叶变换，得到功率谱

Matlab代码

clear all;
close all;
clc;

Fs=1024;       % 采样频率
nfft = 1024;   % fft计算点数大于采样数据点时，补零

%产生含有噪声的序列
n=(0:Fs-1)/Fs;
xn=cos(2*pi*10*n)+3*cos(2*pi*20*n)+(2*randn(size(n)));

subplot(2,1,1);
plot(xn);
title('加噪信号');
grid on

% 自相关法
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
CXk=fft(cxn,nfft);
psd2=abs(CXk);

index=0:round(nfft/2-1);
psd2 = psd2 / max(psd2);
psd2=10*log10(psd2(index+1)+0.000001);% 换成分贝作单位

k=index*Fs/nfft;

subplot(2,1,2);
plot(k,psd2(index+1));
title('自相关法');
grid on

据周期图法计算功率谱

周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。

周期图法是把随机序列 x ( n ) \large x(n) x(n)的 N N N个观测数据视为一能量有限的序列，直接计算 x ( n ) \large x(n) x(n)的离散傅里叶变换，得 x ( k ) \large x(k) x(k),然后再取其幅值的平方，并除以 N N N。

X N ( e − j w ) = ∑ n = 0 N − 1 x ( n ) e − j w n ， S ( w ) = 1 N ∣ X N ( e − j w ) ∣ 2 \large X_{N}(e^{-jw})=\sum_{n=0}^{N-1}x(n)e^{-jwn}，\large S(w)=\frac{1}{N}\left | X_{N}(e^{-jw}) \right |^{2} XN(e−jw)=n=0∑N−1x(n)e−jwn，S(w)=N1∣∣∣XN(e−jw)∣∣∣2

Matlab代码

close all;
clear all;
clc;

Fs=1024;       % 采样频率
nfft = 1024;   % fft计算点数大于采样数据点时，补零

%产生含有噪声的序列
n=(0:Fs-1)/Fs;
xn=cos(2*pi*10*n)+3*cos(2*pi*20*n)+(2*randn(size(n)));

subplot(2,1,1);
plot(xn);
title('加噪信号');
grid on

% 周期图法
window=boxcar(length(xn)); %矩形窗
[psd1,f]=periodogram(xn,window,nfft,Fs); %直接法,周期图功率谱密度估计
psd1 = psd1 / max(psd1);

subplot(2,1,2);
plot(f,10*log10(psd1+0.000001));
title('周期图法');
grid on

现代谱估计法

自相关法是AR模型参数估计中较为简单的一种功率谱估计方法，按照模型阶数从小到大的顺序进行计算。

close all;
clear all;
clc;

Fs=1024;       % 采样频率
nfft = 1024;   % fft计算点数大于采样数据点时，补零

%产生含有噪声的序列
n=(0:Fs-1)/Fs;
xn=cos(2*pi*10*n)+3*cos(2*pi*20*n)+(2*randn(size(n)));

subplot(2,1,1);
plot(xn);
title('加噪信号');
grid on
 
index=0:round(nfft/2-1);
k=index*Fs/nfft;

% AR谱
psd3 = pyulear(xn, Fs, nfft);% 自回归功率谱密度估计-Yule-Walker方法

psd3=psd3/max(psd3);
index=0:round(nfft/2-1);
psd3=10*log10(psd3(index+1)+0.000001);

subplot(2,1,2);
plot(k, psd3);
title('AR谱估计');
grid on;

码农公寓