经典功率谱估计

86 阅读 0 评论 57 点赞

我是靠谱客的博主从容薯片，最近开发中收集的这篇文章主要介绍经典功率谱估计，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

功率谱估计在信号处理中很重要，最近在学习《现代信号处理》，老师让用经典功率谱估计来估计信号的功率谱；记录下来，这也是一种学习。功率谱密度(Power Spectral Density,PSD)

1、简要介绍

首先，利用给定的一组样本数据估计一个平稳随机信号的功率谱密度被陈伟功率谱估计。在许多工程应用中，功率谱的估计是十分重要的，因为它能给出功率谱密度、谱峰的宽度、高度和位置；可以确定运动目标的位置、辐射强度和运动速度。在生物医学工程中，PSD的峰形和波形显示类癫痫病发作的周期，在目标识别中，PSD可以作为目标的特征之一。

其次，在数字信号处理中，一个连续时间的随机过程必须先进行采样，变成离散序列后在进行有关处理，这一处理包括：从连续函数变化为离散序列、从模拟系统变化为离散系统、从Fourier积分变化为Fourier级数；

最后，经典谱估计也是一种非参数化的谱估计；假定随机过程有N个样本x(0),x(1),x(2),...x(N-1)。不失一般性，假定这些数据已零均值化。估计离散信号x( n)的谱估计分别为非参数方法和参数化方法。非参数化方法也称为经典谱估计，他是以Fourier变化为基础的，有两种主要方法，即直接法和间接法。关系如图所示：

总结：

加窗函数虽然能够减小周期图的偏差，改善功率谱曲线的光滑性，但作为非参数化谱估计，周期图具有分辨率低的固有缺陷，不能适应高分辨功率谱估计的需要。与之相比，参数化谱估计可以提过功能比周期图高很多的频率分辨率，故常称参数化谱估计为高分辨率谱估计。

2、实际例子

为什么要把功率谱进行对数化？？

. 因为信号的傅里叶变换所得的数值的动态范围很大，但绝大部分数值都集中在低频区域，如果直接显示的话，显示效果不太好，取对数是为了压缩数值的动态范围，改善显示效果。

.对数坐标能够放大细节，可以提取出一般坐标刻度下提取不出的信号。

. DB 是一个纯计数单位：dB = 10logX。dB的意义其实再简单不过了，就是把一个很大（后面跟一长串0的）或者很小（前面有一长串0的）的数比较简短地表示出来。如：

X = 1000000000000000 = 10logX = 150 dB

X = 0.000000000000001 = 10logX = -150 dB 在不同的使环境下有不同的物理含义。

1）直接法（给出了没有取对数的情况！）通过分别取128/256/512/1024个点进行测试

不同N值的功率谱的方差值
N	128	256	512	1024
方差值	75.0426	83.9851	370.5139	870.3023

结论:功率谱随N的增大，分辨率变大，方差也变大。

matlab代码：

clc;
clear all;

%% 
%------------（直接法)----------
% 对N个点分别进行傅里叶变换，之后求和；用求得的和的模的平方除以N
%代码说明：
%构造随机信号:频率分别为;f1,f2,f3;
%             幅值分别为：A1,A2,A3；
%             噪声采用：均值为0，方差为1的高斯白噪声
%             采样频率为：Fs

Fs=1000;
f1=50;  f2=125;  f3=135;
A1=1;   A2=1.5;  A3=1;
N=128;  %采样点数
% N=256;
% N=512;
% N=1024;
Nfft=N;
n=0:N-1;
t=n/Fs; %采样的时间序列
xn=A1*cos(2*pi*f1*t)+A2*sin(2*pi*f2*t)+A3*cos(2*pi*f3*t)+1.5*randn(size(n));
P=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB
P1=abs(fft(xn,Nfft).^2)/N;%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(P)-1)*Fs/length(P);%给出频率序列
figure;
subplot(221);plot(n,xn);grid on;title('时域信号');ylabel('幅值');
xlabel('1024个点');
subplot(222);plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');
ylabel('功率谱/dB'); xlabel('频率/Hz');
subplot(223);plot(f(1:N/2),P1(1:N/2));grid on;title('功率谱（没有取对数）');%绘制功率谱曲线
ylabel('功率谱');    xlabel('频率/Hz');
std = std(P1)^2;

2） s1 = sin(2*pi*ff*t) s1 = exp(1j*2*pi*ff*t)

直接法、间接法、直接法改进、间接法改进

最后的代码：

clc;
clear all;
%% ------------------------------作业要求--------------
% 对信号x(n)=sin?(ωt)+n(t)和x(n)=exp?(jωt)+n(t)使用上述4种方法进行功率谱估计。
% 叠加高斯白噪声n(t)后的信噪比为-10dB、-3dB、3dB、10dB四种情况。
% 编程不使用fft、xcorr等MATLAB提供的函数。
%% 
%% ----------------------------（直接法)--------------------------------
%
% 对N个点分别进行傅里叶变换，之后求和；用求得的和的模的平方除以N
% fs = 1;   %将采样频率定为信号频率的10倍 fs/f=10;
% N = 1024;
% Nfft=N;
% ff = 0.1; %信号频率
% t = (0:N-1)/fs;
% % s1 = sin(2*pi*ff*t); %信号1
% s1 = exp(1j*2*pi*ff*t); %信号2
% snr1 = -10;%不同的信噪比
% snr2 = -3;
% snr3 = 10;
% snr4 = 3;
% %----------------------四种信噪比的信号------------
% S1 = awgn(s1, snr1, 'measured');
% S2 = awgn(s1, snr2, 'measured');
% S3 = awgn(s1, snr3, 'measured');
% S4 = awgn(s1, snr4, 'measured');
% %-----------四种信噪比的信号功率谱估计----------------
% P1=10*log10(abs(fft(S1,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB
% P2=10*log10(abs(fft(S2,Nfft).^2)/N);
% P3=10*log10(abs(fft(S3,Nfft).^2)/N);
% P4=10*log10(abs(fft(S4,Nfft).^2)/N);
% %-----------------------------------------------
% % P=abs(fft(S,Nfft).^2)/N;%Fourier振幅谱平方的平均值，没有转化为dB
% f=(0:length(P1)-1)*fs/length(P1);%给出频率序列
% %----------展示一个信号的功率谱-------------
% % figure;
% % subplot(221);plot((0:N-1),S);grid on;title('时域信号');ylabel('幅值');
% % xlabel('1024个点');
% % subplot(222);plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');
% % ylabel('功率谱/dB'); xlabel('频率/Hz');
% % subplot(223);plot(f(1:N/2),P1(1:N/2));grid on;title('功率谱（没有取对数）');%绘制功率谱曲线
% % ylabel('功率谱');    xlabel('频率/Hz');
% %----------展示一个信号的不同信噪比snr的功率谱-------------
% figure;
% % subplot(221);plot(f(1:N/2),P1(1:N/2),'k');grid on;title('直接法功率谱（snr=-10）信号为三角函数');ylabel('功率谱/dB'); xlabel('频率/Hz');%不同的信号要选择不同的标题
% subplot(221);plot(f(1:N/2),P1(1:N/2),'k');grid on;title('直接法功率谱（snr=-10）信号为指数函数');ylabel('功率谱/dB'); xlabel('频率/Hz');
% 
% subplot(222);plot(f(1:N/2),P2(1:N/2),'k');grid on;title('功率谱（snr=-3）');
% ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(223);plot(f(1:N/2),P3(1:N/2),'k');grid on;title('功率谱（snr=10）');
% ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(224);plot(f(1:N/2),P4(1:N/2),'k');grid on;title('功率谱（snr=3）');
% ylabel('功率谱/dB'); xlabel('频率/Hz');
%%
%% ----------------------------（间接法）--------------------------------------
%  
% 自相关函数法实现功率谱估计,通过计算自相关函数的傅里叶变换得到的；
% fs = 1;   %将采样频率定为信号频率的10倍 fs/f=10;
% N = 1024;
% NFFT=N;
% ff = 0.1; %信号频率
% t = (0:N-1)/fs;
% s1 = sin(2*pi*ff*t); %信号1
% % s1 = exp(1j*2*pi*ff*t); %信号2
% snr1 = -10;%不同的信噪比
% snr2 = -3;
% snr3 = 10;
% snr4 = 3;
% % %----------------------四种信噪比的信号-----------------------
% S1 = awgn(s1, snr1, 'measured');
% S2 = awgn(s1, snr2, 'measured');
% S3 = awgn(s1, snr3, 'measured');
% S4 = awgn(s1, snr4, 'measured');
% % ------------------------四种不同信号的功率谱------------------
% C1 = xcorr(S1,'unbiased');
% P1 = abs(fft(C1,NFFT));
% P11 = 10*log10(P1((0:round(NFFT/2 -1))+1));
% 
% C2 = xcorr(S2,'unbiased');
% P2 = abs(fft(C1,NFFT));
% P22 = 10*log10(P2((0:round(NFFT/2 -1))+1));
% 
% C3 = xcorr(S3,'unbiased');
% P3 = abs(fft(C3,NFFT));
% P33 = 10*log10(P3((0:round(NFFT/2 -1))+1));
% 
% C4 = xcorr(S4,'unbiased');
% P4 = abs(fft(C4,NFFT));
% P44 = 10*log10(P4((0:round(NFFT/2 -1))+1));
% % -------------------------------------------------------------
% tt = (0:round(NFFT/2 -1))*fs/NFFT;
% % ---------------------显示四种信号-----------------------------
% figure();
% subplot(221);plot(tt,P11,'k');grid on;title('间接法（snr=-10）信号为三角函数');ylabel('功率谱/dB'); xlabel('频率/Hz');%不同的信号要选择不同的标题
% % subplot(221);plot(tt,P11,'k');grid on;title('间接法（snr=-10）信号为指数函数');ylabel('功率谱/dB'); xlabel('频率/Hz');
% 
% subplot(222);plot(tt,P22,'k');grid on;title('（snr=-3）');ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(223);plot(tt,P33,'k');grid on;title('（snr=10）');ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(224);plot(tt,P44,'k');grid on;title('（snr=3）');ylabel('功率谱/dB'); xlabel('频率/Hz');
%% 
%% ----------------------------（直接法改进）--------------------
% %---------------给数据直接加窗------------
% fs = 1;   %将采样频率定为信号频率的10倍 fs/f=10;
% N = 1024;
% Nfft=N;
% ff = 0.1; %信号频率
% t = (0:N-1)/fs;
% % s1 = sin(2*pi*ff*t); %信号1
% s1 = exp(1j*2*pi*ff*t); %信号2
% window=hamming(length(t));%汉明窗
% 
% snr1 = -10;%不同的信噪比
% snr2 = -3;
% snr3 = 10;
% snr4 = 3;
% %----------------------四种信噪比的信号------------
% S1 = awgn(s1, snr1, 'measured');
% S2 = awgn(s1, snr2, 'measured');
% S3 = awgn(s1, snr3, 'measured');
% S4 = awgn(s1, snr4, 'measured');
% %-----------四种信噪比的信号功率谱估计----------------
% [P1,f1]=periodogram(S1,window,Nfft,fs);  %结果显示的时候要转换为dB
% [P2,f2]=periodogram(S2,window,Nfft,fs); 
% [P3,f3]=periodogram(S3,window,Nfft,fs); 
% [P4,f4]=periodogram(S4,window,Nfft,fs); 
% %-----------------------------------------------
% % P=abs(fft(S,Nfft).^2)/N;%Fourier振幅谱平方的平均值，没有转化为dB
% f=(0:length(P1)-1)*fs/length(P1);%给出频率序列
% %----------展示一个信号的功率谱-------------
% % figure;
% % subplot(221);plot((0:N-1),S);grid on;title('时域信号');ylabel('幅值');
% % xlabel('1024个点');
% % subplot(222);plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');
% % ylabel('功率谱/dB'); xlabel('频率/Hz');
% % subplot(223);plot(f(1:N/2),P1(1:N/2));grid on;title('功率谱（没有取对数）');%绘制功率谱曲线
% % ylabel('功率谱');    xlabel('频率/Hz');
% %----------展示一个信号的不同信噪比snr的功率谱-------------
% figure;
% % subplot(221);plot(f1,10*log10(P1),'k');grid on;title('加汉明窗直接法功率谱（snr=-10）信号为三角函数');ylabel('功率谱/dB'); xlabel('频率/Hz');%不同的信号要选择不同的标题
% subplot(221);plot(f1,10*log10(P1),'k');grid on;title('加汉明窗直接法功率谱（snr=-10）信号为指数函数');ylabel('功率谱/dB'); xlabel('频率/Hz');%不同的信号要选择不同的标题
% 
% subplot(222);plot(f2,10*log10(P2),'k');grid on;title('（snr=-3）');ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(223);plot(f3,10*log10(P3),'k');grid on;title('（snr=10）');ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(224);plot(f4,10*log10(P4),'k');grid on;title('（snr=3）');ylabel('功率谱/dB'); xlabel('频率/Hz');
%%
%% -----------------------------间接法改进-----------------------------
%------给自相关函数加窗-----------------
fs = 1;   %将采样频率定为信号频率的10倍 fs/f=10;
N = 1024;
ff = 0.1; %信号频率
t = (0:N-1)/fs;
% s1 = sin(2*pi*ff*t); %信号1
s1 = exp(1j*2*pi*ff*t); %信号2

snr1 = -10;%不同的信噪比
snr2 = -3;
snr3 = 10;
snr4 = 3;
%----------------------四种信噪比的信号------------
S1 = awgn(s1, snr1, 'measured');
S2 = awgn(s1, snr2, 'measured');
S3 = awgn(s1, snr3, 'measured');
S4 = awgn(s1, snr4, 'measured');


%-----------四种信噪比的信号自相关函数----------------
C1 = xcorr(S1,'unbiased');
C2 = xcorr(S2,'unbiased');
C3 = xcorr(S3,'unbiased');
C4 = xcorr(S4,'unbiased');
window=hamming(length(C1));%汉明窗
Nfft = length(C1);
 % ------------四种自相关函数加窗----------------
[P1,f1]=periodogram(C1,window,Nfft,fs);  %结果显示的时候要转换为dB
[P2,f2]=periodogram(C2,window,Nfft,fs); 
[P3,f3]=periodogram(C3,window,Nfft,fs); 
[P4,f4]=periodogram(C4,window,Nfft,fs); 
%-----------------------------------------------
% P=abs(fft(S,Nfft).^2)/N;%Fourier振幅谱平方的平均值，没有转化为dB
f=(0:length(P1)-1)*fs/length(P1);%给出频率序列
%----------展示一个信号的功率谱-------------
% figure;
% subplot(221);plot((0:N-1),S);grid on;title('时域信号');ylabel('幅值');
% xlabel('1024个点');
% subplot(222);plot(f(1:N/2),P(1:N/2));grid on;title('功率谱（dB图）');
% ylabel('功率谱/dB'); xlabel('频率/Hz');
% subplot(223);plot(f(1:N/2),P1(1:N/2));grid on;title('功率谱（没有取对数）');%绘制功率谱曲线
% ylabel('功率谱');    xlabel('频率/Hz');
%----------展示一个信号的不同信噪比snr的功率谱-------------
figure;
% subplot(221);plot(f1,10*log10(P1),'k');grid on;title('加汉明窗直接法功率谱（snr=-10）信号为三角函数');ylabel('功率谱/dB'); xlabel('频率/Hz');%不同的信号要选择不同的标题
subplot(221);plot(f1,10*log10(P1),'k');grid on;title('加汉明窗间接法功率谱（snr=-10）信号为指数函数');ylabel('功率谱/dB'); xlabel('频率/Hz');%不同的信号要选择不同的标题

subplot(222);plot(f2,10*log10(P2),'k');grid on;title('（snr=-3）');ylabel('功率谱/dB'); xlabel('频率/Hz');
subplot(223);plot(f3,10*log10(P3),'k');grid on;title('（snr=10）');ylabel('功率谱/dB'); xlabel('频率/Hz');
subplot(224);plot(f4,10*log10(P4),'k');grid on;title('（snr=3）');ylabel('功率谱/dB'); xlabel('频率/Hz');

以上内容有不对的地方请指正，刚开始学习，可能有不对的地方！！！