MATLAB处理信号得到频谱、相谱、功率谱全解
第一：频谱
一.调用方法
X=FFT(x)；
X=FFT(x，N)；
x=IFFT(X);
x=IFFT(X,N)
用MATLAB进行谱分析时注意：
（1）函数FFT返回值的数据结构具有对称性。
例：
N=8;
n=0:N-1;
xn=[4 3 2 6 7 8 9 0];
Xk=fft(xn)
→
Xk =
39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 -
7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i
Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。
（2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。
二.FFT应用举例
例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
clf;
fs=100;N=128; %采样频率和数据点数
n=0:N-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y=fft(x,N); %对信号进行快速Fourier变换
mag=abs(y); %求得Fourier变换后的振幅
f=n*fs/N; %频率序列
subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz');
ylabel('振幅');title('N=128');grid on;
%对信号采样数据为1024点的处理
fs=100;N=1024;n=0:N-1;t=n/fs;
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y=fft(x,N); %对信号进行快速Fourier变换
mag=abs(y); %求取Fourier变换的振幅
f=n*fs/N;
subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;
subplot(2,2,4)
plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅
xlabel('频率/Hz');
ylabel('振幅');title('N=1024');grid on;
运行结果：
fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。
例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制：
（1）数据个数N=32，FFT所用的采样点数NFFT=32；
（2）N=32，NFFT=128；
（3）N=136，NFFT=128；
（4）N=136，NFFT=512。
clf;fs=100; %采样频率
Ndata=32; %数据长度
N=32; %FFT的数据长度
n=0:Ndata-1;t=n/fs; %数据对应的时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %时间域信号
y=fft(x,N); %信号的Fourier变换
mag=abs(y); %求取振幅
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,1),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=32');grid on;
Ndata=32; %数据个数
N=128; %FFT采用的数据长度
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅');
title('Ndata=32 Nfft=128');grid on;
Ndata=136; %数据个数
N=128; %FFT采用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,3),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=128');grid on;
Ndata=136; %数据个数
N=512; %FFT所用的数据个数
n=0:Ndata-1;t=n/fs; %时间序列
x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);
y=fft(x,N);
mag=abs(y);
f=(0:N-1)*fs/N; %真实频率
subplot(2,2,4),plot(f(1:N/2),mag(1:N/2)*2/N); %绘出Nyquist频率之前的振幅xlabel('频率/Hz');ylabel('振幅');
title('Ndata=136 Nfft=512');grid on;
结论：
（1）当数据个数和FFT采用的数据个数均为32时，频率分辨率较低，但没有由于添零而导致的其他频率成分。
（2）由于在时间域内信号加零，致使振幅谱中出现很多其他成分，这是加零造成的。其振幅由于加了多个零而明显减小。
（3）FFT程序将数据截断，这时分辨率较高。
（4）也是在数据的末尾补零，但由于含有信号的数据个数足够多，FFT振幅谱也基本不受影响。
对信号进行频谱分析时，数据样本应有足够的长度，一般FFT程序中所用数据点数与原含有信号数据点数相同，这样的频谱图具有较高的质量，可减小因补零或截断而产生的影响。
例3：x=cos(2*pi*0.24*n)+cos(2*pi*0.26*n)
（1）数据点过少，几乎无法看出有关信号频谱的详细信息；
（2）中间的图是将x(n)补90个零，幅度频谱的数据相当密，称为高密度频谱图。但从图中很难看出信号的频谱成分。
（3）信号的有效数据很长，可以清楚地看出信号的频率成分，一个是0.24Hz，一个是0.26Hz，称为高分辨率频谱。
可见，采样数据过少，运用FFT变换不能分辨出其中的频率成分。添加零后可增加频谱中的数据个数，谱的密度增高了，但仍不能分辨其中的频率成分，即谱的分辨率没有提高。只有数据点数足够多时才能分辨其中的频率成分。
第二：相谱
（相位谱和频率普是回事儿，想着把频谱中的幅值部分换成相角就可以了）
由于没有找到具体的理论，就举几个例子说明一下。
比如要求y=sin(2*pi*60*t)的相位谱，
程序如下：
fs=200;N=1024;n=0:N-1;t=n/fs;y=sin(2*pi*60*t);
Y=fft(y,N);
A=abs(Y);f=n*fs/N;
ph=2*angle(Y(1:N/2));
ph=ph*180/pi;
plot(f(1:N/2),ph(1:N/2));
xlabel('频率/hz'),ylabel('相角'),title('相位谱');
grid on;
期中的 ph=2*angle(Y(1:N/2));ph=ph*180/pi;是利用angle函数求出每个点的角度，并由弧度转化成角度！
angle函数解释：
Phase angle
Syntax
P = angle(Z)
Description
P = angle(Z) returns the phase angles, in radians, for each element of complex array Z. The angles lie between ±π.
For complex Z, the magnitude R and phase angle theta are given by
R = abs(Z)
theta = angle(Z)
and the statement
Z = R.*exp(i*theta)
converts back to the original complex Z.
Examples
Z = [ 1 - 1i 2 + 1i 3 - 1i 4 + 1i
1 + 2i
2 - 2i
3 + 2i
4 - 2i
1 - 3i
2 + 3i
3 - 3i
4 + 3i
P = angle(Z)
P =
-0.7854 0.4636 -0.3218 0.2450
1.1071 -0.7854 0.5880 -0.4636
-1.2490 0.9828 -0.7854 0.6435
1.3258 -1.1071 0.9273 -0.7854
Algorithms
The angle function can be expressed as angle(z) = imag(log(z)) = atan2(imag(z),real(z)).
第三：功率谱
matlab实现经典功率谱估计
fft做出来是频谱，psd做出来是功率谱；功率谱丢失了频谱的相位信息；频谱不同的信号其功率谱是可能相同的；功率谱是幅度取模后平方，结果是个实数
matlab中自功率谱密度直接用psd函数就可以求，按照matlab的说法，psd能实现Welch 法估计，即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。
1、直接法：
直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
Matlab代码示例：
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); %直接法
plot(f,10*log10(Pxx));
2、间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
Matlab代码示例：
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;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);
3、改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
3.1、Bartlett法
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。
Matlab代码示例：
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(length(n)); %矩形窗
noverlap=0; %数据无重叠
p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot_Pxxc=10*log10(Pxxc(index+1));
figure(1)
plot(k,plot_Pxx);
pause;
figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);
3.2、Welch法
Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。
Matlab代码示例：
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,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);
[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);
plot_Pxx=10*log10(Pxx);
plot_Pxx1=10*log10(Pxx1);
plot_Pxx2=10*log10(Pxx2);
figure(1)
plot(f,plot_Pxx);
pause;
figure(2)
plot(f,plot_Pxx1);
pause;
figure(3)
plot(f,plot_Pxx2);
第四：相关性分析
1. 首先说说自相关和互相关的概念。
这个是信号分析里的概念，他们分别表示的是两个时间序列之间和同一个时间序列在任意两个不同时刻的取值之间的相关程度，即互相关函数是描述随机信号x(t),y(t)在任意两个不同时刻t1，t2的取值之间的相关程度，自相关函数是描述随机信号x(t)在任意两个不同时刻t1，t2的取值之间的相关程度。
自相关函数是描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度；互相关函数给出了在频域内两个信号是否相关的一个
判断指标，把两测点之间信号的互谱与各自的自谱联系了起来。它能用来确定输出信号有多大程度来自输入信号，对修正测量中接入噪声源而产生
的误差非常有效.
事实上，在图象处理中，自相关和互相关函数的定义如下：设原函数是f(t)，则自相关函数定义为R(u)=f(t)*f(-t)，其中*表示卷积；设两个函数分别是f(t)和g(t)，则互相关函数定义为R(u)=f(t)*g(-t)，它反映的是两个函数在不同的相对位置上互相匹配的程度。
那么，如何在matlab中实现这两个相关并用图像显示出来呢？
dt=.1;
t=[0:dt:100];
x=cos(t);
[a,b]=xcorr(x,'unbiased');
plot(b*dt,a)
上面代码是求自相关函数并作图，对于互相关函数，稍微修改一下就可以了，即把
[a,b]=xcorr(x,'unbiased');改为[a,b]=xcorr(x,y,'unbiased');便可。
2. 实现过程：
在Matalb中，求解xcorr的过程事实上是利用Fourier变换中的卷积定理进行的，即R(u)=ifft(fft(f)×fft(g))，其中×表示乘法，注：此公式仅表示形式计算，并非实际计算所用的公式。当然也可以直接采用卷积进行计算，但是结果会与xcorr的不同。事实上，两者既然有定理保证，那么结果一定是相同的，只是没有用对公式而已。下面是检验两者结果相同的代码：
dt=.1;
t=[0:dt:100];
x=3*sin(t);
y=cos(3*t);
subplot(3,1,1);
plot(t,x);
subplot(3,1,2);
plot(t,y);
[a,b]=xcorr(x,y);
subplot(3,1,3);
plot(b*dt,a);
yy=cos(3*fliplr(t)); % or use: yy=fliplr(y);
z=conv(x,yy);
pause;
subplot(3,1,3);
plot(b*dt,z,'r');
即在xcorr中不使用scaling。
3. 其他相关问题：
（1）相关程度与相关函数的取值有什么联系？
相关系数只是一个比率，不是等单位量度，无什么单位名称，也不是相关的百分数，一般取小数点后两位来表示。相关系数的正负号只表
示相关的方向，绝对值表示相关的程度。因为不是等单位的度量，因而不能说相关系数0.7是0.35两倍，只能说相关系数为0.7的二列变量相关程度
比相关系数为0.35的二列变量相关程度更为密切和更高。也不能说相关系数从0.70到0.80与相关系数从0.30到0.40增加的程度一样大。
对于相关系数的大小所表示的意义目前在统计学界尚不一致，但通常按下是这样认为的：相关系数相关程度
0.00-±0.30 微相关
±0.30-±0.50 实相关
±0.50-±0.80 显著相关
±0.80-±1.00 高度相关
（2）matlab计算自相关函数autocorr和xcorr有什么不一样的？
分别用这两个函数对同一个序列计算，为什么结果不太一样？因为xcorr是没有将均值减掉做的相关，autocorr则是减掉了均值的。而且，用离散信号做自相关时，信号截取长度（采样点N）不一样，自相关函数就不一样。
（3）xcorr是计算互相关函数，带有一个option的参数:
a=xcorr(x,y,'option')
option=baised时，是计算互相关函数的有偏估计；
option=unbaised时，是计算互相关函数的无偏估计；
option=coeff时，是计算归一化的互相关函数，即为互相关系数，在-1至1之间；option=none，是缺省的情况。
所以想要计算互相关系数，可用'coeff'参数。
************************************************************************ *
用这个xcorr函数作离散互相关运算时要注意，当x, y是不等长向量时，短的向量会自动填0与长的对齐，运算结果是行向量还是列向量就与x一样。
互相关运算计算的是x,y两组随机数据的相关程度，使用参数coeff时，结果就是互相关系数，在-1至1之间，否则结果不一定在这范围，有可能很大也有可能很小，这视乎x, y数据的大小，所以一般要计算两组数据的相关程度，一般选择coeff参数，对结果进行归一化。所谓归一化简单理解就是将数据系列缩放到-1到1范围，正式的就是一种简化计算的方式，即将有量纲的表达式，经过变换，化为无量纲的表达式，成为纯量。变换式为X=(X实测--Xmin)/(Xmax-Xmin)。
一般来说选择归一化进行互相关运算后，得到结果绝对值越大，两组数据相关程度就越高。
matlab功率谱估计
功率谱估计及其MATLAB仿真 1经典功率谱估计 经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。 1.1相关函数法(BT法) 该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab代码示例1(Btfangfa.M)： Fs=500;%采样频率 n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));%产生含有噪声的序列 nfft=512; cxn=xcorr(xn,'unbiased');%计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); %Round towards nearest integer. k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); figure(1); plot(k,plot_Pxx); 结果如下： 1.2周期图法(periodogram) 周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换得X(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。Matlab代码示例2(PEriod.M)： Fs=600; n=0:1/Fs:1;
matlab实现功率谱密度分析psd
matlab实现功率谱密度分析psd及详细解说 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！ 直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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); %直接法 plot(f,10*log10(Pxx)); 间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk);
(完整版)功率谱估计性能分析及Matlab仿真
功率谱估计性能分析及Matlab 仿真 1 引言 随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。 信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。 2 经典功率谱估计 经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。 周期图法( Periodogram ) Schuster 首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。 取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换 1 ()()N j j n N n X e x n e ω ω---==∑ 然后进行谱估计
功率谱估计介绍(介绍了matlab函数)
功率谱估计介绍 谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以用函数periodogram实现。 周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。 种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT 等技术来计算功率谱。与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。matlab中，welch法用函数psd实现。调用格式如下： [Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP) X：输入样本数据 NFFT：FFT点数 Fs:采样率 WINDOW：窗类型 NOVERLAP，重叠长度 现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。由于涉及的问题太多，这里不再详述，可以参考有关资料。matlab中，现代谱估计的很多方法都可以实现。music方法用pmusic命令实现；pburg函数利用burg法实现功率谱估计；pyulear函数利用yule-walker算法实现功率谱估计等等。 另外，sptool工具箱也具有功率谱估计的功能。窗口化的操作界面很方便，而且有多种方法可以选择 在海杂波抑制的研究中，对海杂波谱分析一定要用到谱估计理论，一定得花时间学好！
功率谱密度估计方法的MATLAB实现
功率谱密度估计方法的MATLAB实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。 以下程序运行平台：Matlab R2015a（8.5.0.197613） 一、周期图法谱估计程序 1、源程序 Fs=100000; %采样频率100kHz N=1024; %数据长度N=1024 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声 subplot(2,1,1); plot(n,Y) title('信号') xlabel('时间');ylabel('幅度');
AR功率谱估计MatlAB
AR模型的谱估计是现代谱估计的主要内容 AR模型的谱估计是现代谱估计的主要内容。 1．AR 模型的Yule—Walker方程和Levinson-Durbin递推算法：在MATLAB中，函数levinson和aryule都采用 Levinson-Durbin递推算法来求解AR模型的参数a1,a2,……,ap及白噪声序列的方差，只是两者的输入参数不同，它们的格式为： A=LEVINSON(R,ORDER) A=ARYULE(x,ORDER) 两函数均为定阶ORDER的求解，但是函数levinson的输入参数要求是序列的自相关函数，而函数aryule的输入参数为采样序列。 下面语句说明函数levinson和函数aryule的功能是相同的： 例子： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; r=xcorr(x,'biased'); r(1:length(x)-1)=[]; A=levinson(r,5) B=aryule(x,5) 2．Burg算法： 格式为：A=ARBURG(x,ORDER); 其中x为有限长序列，参数ORDER用于指定AR 模型的阶数。以上面的例子为例： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; A=arburg(x,5)
3.改进的协方差法： 格式为：A=ARMCOV(x,ORDER); 该函数用来计算有限长序列x(n)的ORDER阶AR 模型的参数。例如：输入下面语句： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; A=armcov(x,5) AR模型阶数P的选择： AR 模型阶数P一般事先是不知道的，需要事先选定一个较大的值，在递推的过程中确定。在使用Levinson—Durbin递推方法时，可以给出由低阶到高阶的每一组参数，且模型的最小预测误差功率Pmin（相当于白噪声序列的方差）是递减的。直观上讲，当预测误差功率P达到指定的希望值时，或是不再发生变化时，这时的阶数即是应选的正确阶数。 因为预测误差功率P是单调下降的，因此，该值降到多少才合适，往往不好选择。比较常见的准则是： 最终预测误差准则：FPE(r)=Pr{[N+(r+1)]/ [N-(r+1)]} 信息论准则：AIC(r)=N*log(Pr)+2*r 上面的N为有限长序列x(n)的长度，当阶数r由1增加时，FPE(r) 和AIC(r)都将在某一r处取得极小值。将此时的r定为最合适的阶数p。 MATLAB中AR模型的谱估计的函数说明： 1． Pyulear函数： 功能：利用Yule--Walker方法进行功率谱估计. 格式： Pxx=Pyulear(x,ORDER,NFFT) [Pxx,W]=Pyulear(x,ORDER,NFFT) [Pxx,W]=Pyulear(x,ORDER,NFFT,Fs) Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
功率谱估计的MATLAB实现
实验功率谱估计 实验目的: 1、掌握最大熵谱估计的基本原理。 2、了解最终预测误差（FPE）准则。 3、掌握周期图谱估计的基本原理。 4、掌握传统谱估计中直接法与间接法之间的关系。 5、复习快速傅里叶变换与离散傅里叶变换之间关系。 实验内容: 1、设两正弦信号的归一化频率分别为0.175和0.20，用最大熵法编程计算信噪比S/N=30dB、N=32点时该信号的最大熵谱估计结果。 2、用周期图法编程计算上述信号的谱估计结果。 程序示例: 1、最大熵谱估计 clc; N=32; SNR=30; fs=1; t=1:N; t=t/fs; y=sin(2*pi*0.175*t)+sin(2*pi*0.20*t); x = awgn(y,SNR); M=1; P(M)=0; Rx(M)=0; for n=1:N P(M)=P(M)+(abs(x(n)))^2; ef(1,n)=x(n); eb(1,n)=x(n); end P(M)=P(M)/N; Rx(M)=P(M); M=2;
A=0; D=0; for n=M:N A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; end xishu=-2*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(xishu))^2); FPE(M-1)=P(M)*(N+M)/(N-M); TH=FPE(M-1); for n=M:N ef(M,n)=ef(M-1,n)+xishu*eb(M-1,n-1); eb(M,n)=eb(M-1,n-1)+xishu*ef(M-1,n); end M=M+1; A=0; D=0; for n=M:N A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; end xishu=-2*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(xishu))^2); FPE(M-1)=P(M)*(N+M)/(N-M); for m=1:M-2 a(M-1,m)=a(M-2,m)+xishu*a(M-2,M-1-m); end while FPE(M-1)
matlab求功率谱
matlab实现经典功率谱估计 fft做出来是频谱，psd做出来是功率谱；功率谱丢失了频谱的相位信息；频谱不同的信号其功率谱是可能相同的；功率谱是幅度取模后平方，结果是个实数 matlab中自功率谱密度直接用psd函数就可以求，按照matlab的说法，psd能实现Welch法估计，即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。 1、直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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); %直接法 plot(f,10*log10(Pxx)); 2、间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft;
(完整word版)自己编写算法的功率谱密度的三种matlab实现方法
功率谱密度的三种matlab 实现方法 一：实验目的： (1)掌握三种算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用； (3) 能够利用burg 法对信号作谱估计，对信号的特点加以分析。 二;实验内容： （1）简单说明三种方法的原理。 （2）用三种方法编写程序，在matlab 中实现。 （3）将计算结果表示成图形的形式，给出三种情况的功率谱图。 （4）比较三种方法的特性。 （5）写出自己的心得体会。 三：实验原理： 1.周期图法： 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。这当然必然带来误差。由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
2.相关法（间接法): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列 )(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧ 序列。 )()(1)(1 m n x n x N m R N n N N x += ∑-=∧ (2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即 jwm M M m X jw x e m R e S ----=∧∧ ∑= )()(1) 1( 以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M< [matlab实现经典功率谱估计]matlab功率 谱估计 1、直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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); %直接法 plot(f,10*log10(Pxx)); 2、间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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; cxn=xcorr(xn,”unbiased”); %计算序列的自相关函数CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot(k,plot_Pxx); 3、改进的直接法： 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。3.1、Bartlett法 Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。 Matlab代码示例： 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(length(n)); %矩形窗 noverlap=0; %数据无重叠 第一：频谱 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[4 3 2 6 7 8 9 0]; Xk=fft(xn) → Xk = 39.0000 -10.7782 + 6.2929i 0 - 5.0000i 4.7782 - 7.7071i 5.0000 4.7782 + 7.7071i 0 + 5.0000i -10.7782 - 6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。 clf; fs=100;N=128; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求得Fourier变换后的振幅 f=n*fs/N; %频率序列 subplot(2,2,1),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号 y=fft(x,N); %对信号进行快速Fourier变换 mag=abs(y); %求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag); %绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2)); %绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果： fs=100Hz，Nyquist频率为fs/2=50Hz。整个频谱图是以Nyquist频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz和40Hz。由此可以知道FFT变换数据的对称性。因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：1，与真实振幅0.5：2是一致的。为了与真实振幅对应，需要将变换后结果乘以2除以N。 例2：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t),fs=100Hz,绘制： （1）数据个数N=32，FFT所用的采样点数NFFT=32； （2）N=32，NFFT=128； （3）N=136，NFFT=128； （4）N=136，NFFT=512。 功率谱估计性能分析及其MATLAB实现 一、经典功率谱估计分类简介 1．间接法 根据维纳-辛钦定理，1958年Blackman和Turkey给出了这一方法的具体实现，即先由N个观察值，估计出自相关函数，求自相关函数傅里叶变换，以此变换结果作为对功率谱的估计。 2．直接法 直接法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机信号的N 个观察值直接进行傅里叶变换，得到，然后取其幅值的平方，再除以N，作为对功率谱的估计。 3．改进的周期图法 将N点的观察值分成L个数据段，每段的数据为M，然后计算L个数据段的周期图的平均 ，作为功率谱的估计，以此来改善用N点观察数据直接计算的周期图的方差特性。根据分段方法的不同，又可以分为Welch法和Bartlett法。 Welch法 所分的数据段可以互相重叠，选用的数据窗可以是任意窗。 Bartlett法 所分的数据段互不重叠，选用的数据窗是矩形窗。 二、经典功率谱估计的性能比较 1．仿真结果 为了比较经典功率谱估计的性能，本文采用的信号是高斯白噪声加两个正弦信号，采样率Fs=1000Hz，两个正弦信号的频率分别为f1=200Hz，f2=210Hz。所用数据长度N=400. 仿真结果如下： Figure1(a)示出了待估计信号的时域波形； Figure2(b)示出了用该数据段直接求出的周期图，所用的数据窗为矩形窗； Figure2(c)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为矩形窗，长度M=128，数据没有加窗； Figure2(d)是用BT法（间接法）求出的功率谱曲线，对自相关函数用的平滑窗为Hamming 窗，长度M=64，数据没有加窗； Figure2(e)是用Welch平均法求出的功率谱曲线，每段数据的长度为64点，重叠32点，使用的Hamming窗； Figure2(f)是用Welch平均法求出的功率谱曲线，每段数据的长度为100点，重叠48点，使用的Hamming窗； 2．性能比较 1)直接法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现 虚假谱峰； 2)间接法由于使用了平滑窗对直接法估计的功率谱进行了平滑，因此方差性能比直接 法好，功率谱比直接法估计的要平滑，但其分辨率比直接法低。 3)Welch平均周期图法是三种经典功率谱估计方法中方差性能最好的，估计的功率谱 也最为平滑，但这是以分辨率的下降及偏差的增大为代价的。 3．关于经典功率谱估计的总结 1)功率谱估计，不论是直接法还是间接法都可以用FFT快速计算，且物理概念明确，因而 仍是目前较常用的谱估计方法。 2)谱的分辨率较低，它正比于2π/N，N是所使用的数据长度。 3)方差性能不好，不是真实功率谱的一致估计，且N增大时，功率谱起伏加剧。 4)周期图的平滑和平均是和窗函数的使用紧密关联的，平滑和平均主要是用来改善周期图 的方差性能，但往往又减小了分辨率和增加了偏差，没有一个窗函数能使估计的功率谱在方差、偏差和分辨率各个方面都得到改善，因此使用窗函数只是改进估计质量的一个技巧问题，并不能从根本上解决问题。 三、AR模型功率谱估计 1．A R模型功率谱估计简介 AR模型功率谱估计是现代谱估计中最常用的一种方法，这是因为AR模型参数的精确估计可以用解一组线性方程（Yule-Walker方程）的方法求得。其核心思想是：将信号看成是一个p 阶AR过程，通过建立Yule-Walker方程求解AR模型的参数，从而得到功率谱的估计。 由于已知的仅仅是长度有限的观测数据，因此AR模型参数的求得，通常是首先通过某种算法求得自相关函数的估计值，进而求得AR模型参数的估计值。常用的几种AR模型参数提取方法有： 1)自相关法 假定观测数据区间之外的数据为0，在均方误差意义下使得数据的前向预测误差最小。 经典功率谱估计 直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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); %直接法 plot(f,10*log10(Pxx)); 间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot(k,plot_Pxx); 改进的直接法： 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。 1. Bartlett法 Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。 Matlab代码示例： clear； Fs=1000; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); MATLAB中AR模型功率谱估计中AR阶次估计的实现 （最近看了几个关于功率谱的问题，有关AR模型的谱估计，在此分享一下，希望大家不吝指正） （声明：本文内容摘自我的毕业论文——心率变异信号的预处理及功率谱估计） （按：AR模型功率谱估计是对非平稳随机信号功率谱估计的常用方法，但是其模型阶次的估计，除了HOSA工具箱里的arorder函数外，没有现成的函数可用，arorder函数是基于矩阵SVD分解的阶次估计方法，为了比较各种阶次估计方法的区别，下面的函数使用了'FPE', 'AIC', 'MDL', 'CAT'集中准则一并估计，并采用试验方法确定那一个阶次更好。） ………………………………以上省略…………………………………………………………………… 假设原始数据序列为x，那么n阶参数使用最小二乘估计在MATLAB中实现如下： 复制内容到剪贴板 代码: Y = x; Y(1:n) = []; m = N-n; X = [];% 构造系数矩阵 for i = 1:m for j = 1:n X(i,j) = xt(n+i-j); end end beta = inv(X'*X)*X'*Y'; beta即为用最小二乘法估计出的模型参数。 此外，还有估计AR模型参数的Yule-Walker方程法、基于线性预测理论的Burg算法和修正的协方差算法等[26]。相应的参数估计方法在MATLAB中都有现成的函数，比如aryule、arburg以及arcov等。 4.3.3 AR模型阶次的选择及实验设计 文献[26]中介绍了五种不同的AR模型定阶准则，分别为矩阵奇异值分解（Singular Value Decomposition, SVD）定阶法、最小预测定误差阶准则 （Final Prediction Error Criterion, FPE）、AIC定阶准则（Akaika’s Information theoretic Criterion, AIC）、MDL定阶准则以及CAT定阶准则。文献[28]中还介绍了一种BIC定阶准则。SVD方法是对Yule-Walker方 程中的自相关矩阵进行SVD分解来实现的，在MATLAB工具箱中arorder函数就是使用的该算法。其他五种算法的基本思想都是建立目标函数，阶次估计的标准是使目标函数最小化。 以上定阶准则在MATLAB中也可以方便的实现，下面是本文实现FPE、AIC、MDL、CAT定阶准则的程序（部分）： 复制内容到剪贴板 代码: for m = 1:N-1 …… % 判断是否达到所选定阶准则的要求 if strcmp(criterion,'FPE') objectfun(m+1) = (N+(m+1))/(N-(m+1))*E(m+1); elseif strcmp(criterion,'AIC') objectfun(m+1) = N*log(E(m+1))+2*(m+1); elseif strcmp(criterion,'MDL') objectfun(m+1) = N*log(E(m+1))+(m+1)*log(N); elseif strcmp(criterion,'CAT') for index = 1:m+1 temp = temp+(N-index)/(N*E(index)); end objectfun(m+1) = 1/N*temp-(N-(m+1))/(N*E(m+1)); end if objectfun(m+1) >= objectfun(m) orderpredict = m; break; end end orderpredict变量即为使用相应准则预测的AR模型阶次。 用matlab做经典功率谱估计 经典功率谱估计 1、直接法： 直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。 Matlab代码示例： 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); %直接法 plot(f,10*log10(Pxx)); 2、间接法： 间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。 Matlab代码示例： 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; cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot(k,plot_Pxx); 3、改进的直接法： 对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N 太小，谱的分辨率又不好，因此需要改进。 3.1、Bartlett法 Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。 Matlab代码示例： 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(length(n)); %矩形窗 noverlap=0; %数据无重叠 p=0.9; %置信概率 [Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p); index=0:round(nfft/2-1); k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot_Pxxc=10*log10(Pxxc(index+1)); figure(1) plot(k,plot_Pxx); pause; figure(2) plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]); 3.2、Welch法 Welch法对Bartlett法进行了两方面的修正，一是选择适当的窗函数 w(n)，并再周期图计算前直接加进去，加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时，可使各段之间有重叠，这样会使方差减小。 Matlab代码示例： clear; Fs=1000; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); 【转】matlab的功率谱计算 功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与 现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab 中，周期图法可以用函数periodogram实现。 但是周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行 谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数 据的功率谱估计。 这2种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT等技术来计算功率谱。与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。matlab中，welch法用函数psd实现。调用格式如下：[Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP) X：输入样本数据 NFFT：FFT点数 Fs:采样率 WINDOW：窗类型 NOVERLAP，重叠长度 现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和 数字信号处理 课程实验报告 实验指导教师：黄启宏 实验名称 MATLAB 仿真实现功率谱估计（采用自相关法） 专业、班级 电子与通信工程 姓 名 张帅 实验地点 仿古楼301 实验日期 2013.11.17 一、实验内容 采用自相关法（参数模型）实现功率谱估计。 二、实验目的 （1）掌握自相关法（参数模型）估计出功率谱的步骤和方法； （2）在实验的过程中找到影响功率谱估计的因素； （3）比较自相关法（参数模型）和周期图法（直接法）估计功率谱的优点。 三、实验原理 （1）假定研究的()x n 过程是由一个输入序列()u n 激励一个线性系统()H z 的输出； （2）由已知的 ()x n ，或者其自相关函数()x r m 来估计()H z 的参数； （3）由()H z 的参数来估计()x n 的功率谱。 四、涉及实验的相关情况介绍（包含使用软件或实验设备等情况） 一台安装MATLAB 软件的电脑 五、实验记录 程序、相关的图形、相关数据记录及分析）（ %自相关法（参数模型）实现功率谱估计； fs = 1000; n = 0:1/fs:0.3; xn = cos(200*pi*n)+0.1*randn(size(n)); nfft = 1024; window=boxcar(length(xn)); %输出随机信号xn； subplot(211); plot(n,xn);xlabel('时间');ylabel('幅度');title('输入信号x（n）'); axis([0 0.3 -2 2]); grid on; %自相关法功率谱估计 p1=floor(length(xn)/3)+1;%阶数 [px1,f1] = pyulear(xn,100,nfft,fs); px1max = max(px1); px1 = px1/px1max; px1=10*log10(px1+0.000001); [Pxx,f] = periodogram(xn,window,nfft,fs);%直接法功率谱估计PxxMAX =max(Pxx); Pxx = Pxx/PxxMAX;[matlab实现经典功率谱估计]matlab功率谱估计
MATLAB处理信号得到频谱、相谱、功率谱全解
功率谱估计MATLAB实现
基于matlab的经典功率谱估计(有源程序)
MATLAB中AR模型功率谱估计中AR阶次估计的实现
用matlab做经典功率谱估计
Matlab功率谱计算
MATLAB仿真实现功率谱估计(采用自相关法)