能量谱与功率谱理解与编程

也可以这么写 

能量E、功率P的公式中，只和T周期也就是时间和f(t)信号本身有关

一、周期信号：

无限时间的正弦波，能求出他的面积吗，不能的。那再求出它平方的面积也是不能的，能量是无穷的。那么什么情况能量有限啊，肯定是能求出f(t)面积啊，只有它是无限趋近于0，才能求出来，所以它必须是非周期的，则能量有限，称为能量信号。

二、非周期信号

 

 

 

 

能量谱=能量谱密度：定义指单位频率的信号能量

帕塞瓦尔定理：

时域能量=频域能量

 

称为能量信号的能量谱密度，它表示在频率f处宽度为df的频带内的信号能量，单位：J或者可以看做是单位频带内的信号能量。

 

功率谱：

功率谱=功率谱密度，它表示单位频率的信号功率。

怎么理解能量谱与功率谱

能量信号和功率信号是随时间变换的，画图就是横轴是时间，纵轴是能量和功率幅度

由帕塞瓦尔定理，时域的能量等于频域的能量，把能量和功率取傅里叶变换，把随时间变换换成随频率变换的图

横轴是频率，纵轴是能量和功率幅度。

 

信号的功率谱显示的是信号在不同频段能量的强弱，从频域上反映信号的规律。功率谱估计的作用是把时域中幅度随时间变化的趋势变换

为信号能量随频率变化的图像，这样可以直接看到脑电信号不同频率成分在信号中所占据的比例。 

能量谱则是从能量方面考虑，每个频段的能量的规律。

 

编程：

直接计算功率谱密度，也叫周期图法

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); %直接法，直接计算功率谱

subplot(1,2,1);
plot(xn);
subplot(1,2,2);
plot(f,10*log10(Pxx));

 

维纳辛钦定理：

 

相关函数法：

fs=500;
n=0:1/fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));
cxn=xcorr(xn,'unbiased');

len=length(xn);
bx=fft(cxn,len);
px=abs(bx);%自相关函数的傅里叶

index=0:round(len/2-1);
k=index*fs/len; 
plot_Pxx=10*log10(px(index+1));
subplot(1,2,1);
plot(n,xn);
subplot(1,2,2);
plot(k,plot_Pxx);

 

经典谱估计：直接法=周期图法

现代谱估计

改进方法----分段平均周期图法（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

Welch法
Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n) w(n)w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小
 

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,f1]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f2]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);

subplot(1,3,1);
plot(f,plot_Pxx);

subplot(1,3,2);
plot(f1,plot_Pxx1);

subplot(1,3,3);
plot(f2,plot_Pxx2);