1、第一：频谱一.调用方法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.2929iXk 与 xn 的维数相同，共有 8 个元素。Xk 的第一个数对应于直流分量，即频率值为 0。（
2、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 变换
3、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(
4、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;运行
5、结果：fs=100Hz，Nyquist 频率为 fs/2=50Hz。整个频谱图是以 Nyquist 频率为对称轴的。并且可以明显识别出信号中含有两种频率成分：15Hz 和 40Hz。由此可以知道 FFT 变换数据的对称性。因此用 FFT 对信号做 谱分析，只需考察 0Nyquist 频率范围内的福频特性。若没有给出采样频率和采样间隔，则分析通常对归一化频率 01 进行。另外，振幅的大小与所用采样点数有关，采用 128 点和 1024 点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz 与 15Hz 振动幅值之比均为 4：1 ，与真实振幅 0.5：2 是一致的。为了与真实振幅对应，需要
6、将变换后结果乘以 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); %信号的 F
7、ourier 变换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)*f
8、s/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
9、/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 频
10、率之前的振幅xlabel(频率/Hz);ylabel( 振幅);title(Ndata=136 Nfft=512);grid on;结论：（1）当数据个数和 FFT 采用的数据个数均为 32 时，频率分辨率较低，但没有由于添零而导致的其他频率成分。（2）由于在 时间 域内信号加零，致使振幅谱中出现很多其他成分，这是加零造成的。其振幅由于加了多个零而明显减小。（3）FFT 程序将数据截断，这时分辨率较高。（4）也是在数据的末尾补零，但由于含有信号的数据个数足够多，FFT 振幅谱也基本不受影响。对信号进行 频谱分析时，数据样本应有足够的长度，一般 FFT 程序中所用数据点数与原含有信号数据点数相同
11、，这样的频谱图具有较高的质量，可减小因补零或截断而产生的影响。例 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 变换不能分辨出其中的频率成分。添加零后可增加频谱中的数据个数，谱的密度增高了，但仍不能分辨其中的频率成分，即谱的分辨率没有提高。只有
12、数据点数足够多时才能分辨其中的频率成分。第二： 相谱（相位谱和频率普是回事儿，想着把频谱中的幅值部分换成相角就可以了）由于没有找到具体的理论，就举几个例子说明一下。比如要求 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*a
13、ngle(Y(1:N/2);ph=ph*180/pi;是利用 angle 函数求出每个点的角度，并由弧度转化成角度！angle 函数解释：Phase angle SyntaxP = angle(Z)DescriptionP = 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 byR = abs(Z)the
14、ta = angle(Z)and the statementZ = R.*exp(i*theta)converts back to the original complex Z.ExamplesZ = 1 - 1i 2 + 1i 3 - 1i 4 + 1i1 + 2i 2 - 2i 3 + 2i 4 - 2i1 - 3i 2 + 3i 3 - 3i 4 + 3iP = angle(Z)P =-0.7854 0.4636 -0.3218 0.24501.1071 -0.7854 0.5880 -0.4636-1.2490 0.9828 -0.7854 0.64351.3258 -1.1071 0
15、.9273 -0.7854AlgorithmsThe 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 求
16、出的结果应该更光滑吧。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);
17、%直接法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*F
18、s/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);
19、 %矩形窗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 法对 Ba
20、rtlett 法进行了两方面的修正，一是 选择适当的窗函数 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; %数据无重叠r
21、ange=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,pl
22、ot_Pxx1);pause;figure(3)plot(f,plot_Pxx2); 第四： 相关性分析1. 首先说说自相关和互相关的概念。这个是信号分析里的概念，他们分别表示的是两个时间 序列之间和同一个时间序列在任意两个不同时刻的取值之间的相关程度，即互相关函数是描述随机信号 x(t),y(t)在任意两个不同时刻 t1，t2 的取值之间的相关程度，自相关函数是描述随机信号 x(t)在任意两个不同时刻 t1，t2 的取值之间的相关程度。自相关函数是描述随机信号 X(t)在任意两个不同时刻 t1，t2 的取值之间的相关程度；互相关函数给出了在频域内两个信号是否相关的一个判断指标，把两测点之间信
23、号的互谱与各自的自谱联系了起来。它能用来确定输出信号有多大程度来自输入信号，对修正测量中接入噪声源而产生的误差非常有效.事实 上，在图象处理中，自相关和互相关函数的定义 如下：设原函数是 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)上面代
24、码是求自相关函数并作图，对于互相关函数，稍微修改一下就可以了，即把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
25、;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）相关程度与相关函数的取值有什么联系？相关系数只是一个比率，不是等单位量度，无什么单 位名称，也不是相关的百分数，
26、一般取小数点后两位来表示。相关系数的正负号只表示相关的方向，绝对值表示相关的程度。因为不是等单位的度量，因而不能说相关系数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 计 算自相关函数 aut
27、ocorr 和 xcorr 有什么不一样的？分别用这两个函数对同一个序列计算，为什么结果不太一样？因为 xcorr 是没有将均值减掉做的相关，autocorr 则是减掉了均值的。而且，用离散信号做自相关时，信号截取长度（采样点 N）不一样，自相关函数就不一样。（3）xcorr 是计算互相关函数，带有一个 option 的参数:a=xcorr(x,y,option)option=baised 时，是计算互相关函数的有偏估计；option=unbaised 时，是计算互相关函数的无偏估计；option=coeff 时，是计算 归一化的互相关函数，即为互相关系数，在-1 至 1 之间；option=
28、none，是缺省的情况。所以想要计算互相关系数，可用coeff参数。*用这个 xcorr 函数作离散互相关运算时要注意，当 x, y 是不等长向量时，短的向量会自动填 0 与长 的对齐 ，运算结果是行向量还是列向量就与 x 一样。互相关运算计算的是 x,y 两组随机数据的相关程度，使用参数 coeff 时，结果就是互相关系数，在-1 至 1 之间，否则结果不一定在这范围，有可能很大也有可能很小，这视乎 x, y 数据的大小，所以一般要计算两组数据的相关程度，一般选择 coeff 参数，对结果进行归一化。所谓归一化简单理解就是将数据系列缩放到-1 到 1 范围，正式的就是一种简化计算的方式，即将有量纲的表达式，经过变换，化为无量纲的表达式，成为纯量。变换式为 X=(X 实测-Xmin)/(Xmax-Xmin)。一般来说选择归一化进行互相关运算后，得到结果绝对值越大，两组数据相关程度就越高。
