利用经典谱估计法估计信号的功率谱(随机信号)

plot_Pxx=10*log(Pxx(index+1))/(log(10));
axes(handles.axes2);
plot(k,plot_Pxx);
xlabel('频率（hz）');
ylabel('功率谱密度（Db/Hz)');
title('BT法');
BT法所得功率谱波形如下：
(3)周期图法
xlabel('频率（hz）');
ylabel('功率谱密度（Db/Hz)');
title('平均法');
functionpushbutton5_Callback(hObject, eventdata, handles)
globaly Fs nfft n;
cxn=xcorr(y,'unbiased');
functionpushbutton1_Callback(hObject, eventdata, handles)
globaly Fs nfft n;
Fs=600;
nfft=1024;
n=0:1/Fs:1;
y=cos(2*pi*30*n)+3*cos(2*pi*100*n)+randn(size(n));
用有限长样本序列的DFT来表示随机序列的功率谱只是一种估计或近似，不可避免存在误差。为了减少误差，使功率谱估计更加平滑，可采用平均法。
平均法采用了分段的功率谱估计，较之于周期图法，平均法的估计曲线较为平滑。在程序中加入窗函数，使得谱分辨率不会由于分段而下降。
Welch法就是用改进的平均法来求取随机信号的功率谱密度估计。Welch法采用信号重叠分段，加窗函数和FFT算法等计算一个信号序列的功率谱。
由于上式中的离散傅里叶变换可以用快速傅里叶变换计算，因此就可以估计出功率谱。
3.平均法
理论基础：
平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少，达到一致估计的目的，有一个定理：如果 是不相关的随机变量，且都有个均值 及其方差 ，则可以证明它们的算术平均的均值为 。
由定理可见：具有 ，当 时，方差 ，可以达到一致估计的目的。因此，将 个独立的估计量经过算术平均后得到的估计量的方差也是原估计量方差的 。
六
郭 敏 （学号：3222008054）完成了平均法以及Welch法的程序编写以及实验报告的统筹,编写.
七
[1]常建平,李林海.随机信号分析.科学出版社.2010年，第5章第3节：184－189；
[2]黄文梅,熊桂林,杨勇.信号分析与处理.国防科技大学出版社.2000年,第6章第3节:222—229
'gui_OutputFcn', @ming_OutputFcn,...
'gui_LayoutFcn', [] ,...
'gui_Callback', []);
ifnargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log(Pxx(index+1))/(log(10));
4.Welch
理论基础：
Welch法又称修正周期法，其步骤为：
先将N个数据分成L段，每数据段M个数据，N=ML；
选择适当的窗函数w（n），并用该w（n）依次对每段数据做相应的加权，然后确定每段的周期图
U为归一化因子
对每段周期图进行平均得到功率谱估计：
当数据量一定时，若分段数L增加，M点数减少，则分辨率下降；若L减少，虽M增加，但方差增大。解决这一矛盾的方法是，让数据段间适当“重叠”。
axes(handles.axes5);
plot(f,Qxx1);
xlabel('频率（hz）');
ylabel('功率谱密度（Db/Hz)');
title('Welch法');
Welch法所得功率谱波形如下：
四
综述
从给出一段信号y=cos(2*pi*30*n)+3*cos(2*pi*100*n),利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数，并利用GUI界面呈现出不同谱估计方法所得的结果。
三．算法设计与实现
1
（1）BT法的流程图
（2）周期图法的流程图
（3）平均法的流程图
（4）WHale Waihona Puke Baidulch法的流程图
2
(1)产生原始信号
Fs=600;
nfft=1024;
n=0:1/Fs:1;
y=cos(2*pi*30*n)+3*cos(2*pi*100*n)+randn(size(n));
axes(handles.axes1);
平均法所得功率谱波形如下：
(5)Welch法
global y Fs nfft n;
window1=hamming(100);
noverlap=20;
range='half';
[Pxx1,f]=pwelch(y,window1,noverlap,nfft,Fs,range);
Qxx1=10*log10(Pxx1);
guidata(hObject, handles);
functionvarargout = ming_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
functionfigure1_CreateFcn(hObject, eventdata, handles)
plot(n,y);
title(‘原始信号波形‘);
原始信号波形如下：
(2)BT法
global y Fs nfft n;
cxn=xcorr(y,'unbiased');
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
global y Fs nfft n;
XF=fft(y,nfft);
Pxx=abs(XF).^2/length(n);
index=0:round(nfft/2-1);
f=index*Fs/nfft;
axes(handles.axes4);
plot(f,10*log(Pxx(index+1)));
xlabel('频率（hz）');
ylabel('功率谱密度（Db/Hz)');
title('周期图法');
周期图法所得功率谱波形如下：
(4)平均法
global y Fs nfft n;
window=hamming(nfft);
noverlap=0;
p=0.9;
[Pxx,Pxxc]=psd(y,nfft,Fs,window,noverlap,p);
八．附录
程序源代码：
functionvarargout = ming(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename,...
'gui_Singleton', gui_Singleton,...
'gui_OpeningFcn', @ming_OpeningFcn,...
平均法
平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少，达到一致估计的目的。
获取x(n)后，将x(n)分为10段，对每段用 计算出周期图，对以上10个周期图加以平均得出功率谱。
Welch法
Welch法又称修正周期法。
获取x(n)后，先将N个数据分成L段，选择汉明窗w(n)，并用该w(n)依次对每段数据做相应的加权，确定每段的周期图
随机信号
利用经典谱估计法估计信号的功率谱
作业综述
给出一段信号“asd.wav”,利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数。采用MATLAB语言，利用MATLAB语言强大的数据处理和数据可视化能力，通过GUI的对话框模板，使操作更为简便！在一个GUI界面中，同时呈现出不同方法产生出的功率谱。
axes(handles.axes1);
plot(n,y);
title(‘原始信号波形‘);
functionpushbutton3_Callback(hObject, eventdata, handles)
globaly Fs nfft n;
window=hamming(nfft);
noverlap=0;
BT法
BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值 。
由式 = 估算出 ，再对 作FFT变换，得到 。
周期图法
周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。
获取x(n)后，对x(n)作FFT得到X(W)，再由 得出功率谱。
由 得出每段谱估计。
GUI界面如下：
GUI界面
五
由图可以看出，在频率30hz和100hz处功率谱有两个峰值，说明信号中有30hz和100hz的周期成分。
通过BT法能观察到两个峰值，但是所呈现的波形不能准确表达出信号的功率谱变化情况。
通过周期图法求出的功率谱密度在很大范围内波动，而且容易证明，即使增加信号取样点数N，实验效果依然没有明显改进。
1.
理论基础：
(1)随机序列的维纳-辛钦定理
由于随机序列{X(n)}的自相关函数Rx(m)=E[X(n)X(n+m)]定义在离散点m上，设取样间隔为 ，则可将随机序列的自相关函数用连续时间函数表示为
等式两边取傅里叶变换，则随机序列的功率谱密度
(2)谱估计
BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值 。即
这里给出了几种不同的方法：BT法，周期图法，平均法以及Welch法。把几种不同方法所得到的功率谱都呈现在一个界面中，便于对几种不同方法得到的功率谱作对比。
一．
给出一段信号及采样率，利用经典谱估计法估计出信号的功率谱。
二．
经典谱估计的方法，实质上依赖于传统的傅里叶变换法。它是将数据工作区外的未知数据假设为零，相当于数据加窗，主要方法有BT法，周期图法，平均法以及Welch法。
其中 可有式得到。
2.
理论基础：
周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。在前面我们已知，各态历经的连续随机过程的功率谱密度满足
式中 是连续随机过程第i个样本的截取函数 的频谱。对应在随机序列中则有
由于随机序列中观测数据 仅在 的点上存在，则 的N点离散傅里叶变换为：
因此有随机信号的观测数据 的功率谱估计值（称“周期图”）如下：
index=0:round(nfft/2-1);
k=index*Fs/nfft;
Qxx=10*log10(Pxx(index+1));
axes(handles.axes6);
plot(k,Qxx );
xlabel('频率（hz）');
ylabel('功率谱密度（Db/Hz)');
title('平均法');
end
ifnargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
functionming_OpeningFcn(hObject, eventdata, handles, varargin)
D[ ]= D[ ]
即：平均法将 的N个数据分成L段（N=ML），若各数据段相互独立，则平方后估计量的方差是原来不分段估计量方差的 。所以当 时，估计量的方差趋于0，达到一致估计的目的。但是，随着分段数L的增加，M点数减少，分辨率减少，使估计变成有偏估计。相反，若L减少，M增加，虽偏差减少，但方差增大。所以，在应用中，必须兼顾分辨率和方差的要求来适当选择M和L的值。
p=0.9;
[Pxx,Pxxc]=psd(y,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
Qxx=10*log10(Pxx(index+1));
axes(handles.axes6);
plot(k,Qxx );
平均图法即是将数据 分段求周期图法后再平均。例如，给定N=1000个数据样本（平均法适用于数据量大的场合），则可以将它分成10个长度为100的小段，分别计算每一段的周期图
然后将这10个周期图加以平均得谱估计值：
由于这10小段的周期图取决于同一个过程，因而其均值相同。若这10个小段的周期图是统计独立的，则这10个小段平均之后的方差却是单段方差的 。