MATLAB处理信号得到频谱相谱功率谱.docx
《MATLAB处理信号得到频谱相谱功率谱.docx》由会员分享，可在线阅读，更多相关《MATLAB处理信号得到频谱相谱功率谱.docx（11页珍藏版）》请在冰豆网上搜索。
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=[43267890];
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'）;gridon;
subplot（2,2,2）,plot（f（1:
N/2）,mag（1:
N/2））;%绘出Nyquist频率之前随频率变化的振幅
xlabel（'频率/Hz'）;
ylabel（'振幅'）;title（'N=128'）;gridon;
%对信号采样数据为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'）;gridon;
subplot（2,2,4）
plot（f（1:
N/2）,mag（1:
N/2））;%绘出Nyquist频率之前随频率变化的振幅
xlabel（'频率/Hz'）;
ylabel（'振幅'）;title（'N=1024'）;gridon;
运行结果：
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=32Nfft=32'）;gridon;
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=32Nfft=128'）;gridon;
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=136Nfft=128'）;gridon;
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=136Nfft=512'）;gridon;
结论：
（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（'相位谱'）;
gridon;
期中的 ph=2*angle（Y（1:
N/2））;ph=ph*180/pi;是利用angle函数求出每个点的角度，并由弧度转化成角度！
angle函数解释：
Phaseangle
Syntax
P=angle（Z）
Description
P=angle（Z）returnsthephaseangles,inradians,foreachelementofcomplexarrayZ.Theanglesliebetween±π.
ForcomplexZ,themagnitudeRandphaseanglethetaaregivenby
R=abs（Z）
theta=angle（Z）
andthestatement
Z=R.*exp（i*theta）
convertsbacktotheoriginalcomplexZ.
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
Theanglefunctioncanbeexpressedasangle（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_Pxxplot_Pxx-plot_Pxxcplot_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';%频率间隔为[0Fs/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））;%oruse:
yy=fliplr（y）;
z=conv（x,yy）;
pause;
subplot（3,1,3）;
plot（b*dt,z,'r'）;
即在xcorr中不使用scaling。
升级会员 