离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
. v& M/ R3 S- o+ }2 W%这里的向量都用行向量，假设被测变量是速度，单位为m/s
4 s# k: m4 C) Qclear;
; q  Q* \2 a1 t) L4 kclose all;

load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
7 r* c# |, E3 kA=data;                                        %将测量数据赋给A，此时A为N×2的数组
x=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
y=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
* i/ q. ^9 H) p& N- P, L; O( N5 X( v8 T% X9 ]y=y';                                           %将列向量变成行向量

%显示数据基本信息
' p& T8 w7 e  |' gfprintf('\n数据基本信息：\n')
2 i$ H' G5 A/ _6 `4 dfprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
fprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
3 ]) d7 A7 R2 j# T  Sfprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
7 `3 p/ e( h4 k% y" q" s4 i! G2 wfprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
fprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
+ [, c4 }8 q0 m0 h* `fprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
fprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
, A" q* i3 M4 K2 u0 j0 ^: |: Dfprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
$ A, o. O* A) n" o- ifprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
  
%显示原始数据曲线图（时域）
subplot(2,1,1);
/ R7 D3 Q2 Z6 z" f; i! Eplot(x,y)                                                                                %显示原始数据曲线图
axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
$ J3 L* x$ o1 u4 G& K" A% k0 uxlabel('时间 (s)');
& [0 v* b* M6 k# dylabel('被测变量y');
# }( E& w, z; B! ?title('原始信号(时域)');
grid on;

7 J1 ], K2 d! K  j# d%傅立叶变换
y=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
Fs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
+ h9 k/ n. g5 i7 t4 F5 ]N=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
z=fft(y);

%频谱分析
9 y/ k6 Q. x% k, P, n3 y& d: If=(0:N-1)*Fs/N;
- {4 n1 x( U2 z. o$ yMag=2*abs(z)/N;                                        %幅值，单位同被测变量y
7 k9 T0 Z# A; B' DPyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式

5 p1 P. {( `, E3 |! Z%显示频谱图(频域)
3 H- m+ R3 u2 a+ e8 V/ r& Isubplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
" {" @! u5 @3 |7 e9 S! V. J%                 |
4 {( J) x' g( U& U8 B1 s%             将这里的Pyy改成Mag就是 幅值-频率图了
, a* ~# Z) ]- _6 @axis([min(f(1:N/2)) max(f(1:N/2)) 1.1*floor(min(Pyy(1:N/2))) 1.1*ceil(max(Pyy(1:N/2)))])
1 B. `* d0 M. l9 g* ~2 zxlabel('频率 (Hz)')
& x+ Z, ?) w9 V7 |5 ?ylabel('能量')
title('频谱图(频域)')
) q/ V0 n& ~6 ?& M5 E& w& @grid on;

( \  \4 @3 x( o3 W- f( z%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
fprintf('\n傅立叶变换结果：\n')
fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
7 L, C: |# R" ^3 U0 T4 cfprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期

yin123

离散信号MATLAB频谱分析程序