离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
- j+ |- o% p/ ]) L8 Rclear;
. l0 I: L/ H- E3 w0 p4 {close all;
! z- m* X$ [# _9 B" M9 s9 t0 {7 ^8 V* J
. _7 C3 m" V; C' Wload data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
A=data;                                        %将测量数据赋给A，此时A为N×2的数组
x=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
y=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
9 T6 e! ?$ G/ E; my=y';                                           %将列向量变成行向量
! Z/ p) l9 X2 C' {/ o# o
%显示数据基本信息
fprintf('\n数据基本信息：\n')
fprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
' Q. Y& T0 z( B1 ~1 J& a0 gfprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
fprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
& j5 F- u, u$ u! gfprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
* L! P( O/ `8 Bfprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
  R% X6 g- _! ofprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
9 P6 _* H) m! w1 t, r! lfprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
fprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
fprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
+ G& X/ ^0 A4 h6 w3 mfprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
  
9 ?! {/ \2 v7 x2 D$ k- ^: I%显示原始数据曲线图（时域）
$ @5 }" t9 v3 Y5 g8 p: n4 @subplot(2,1,1);
plot(x,y)                                                                                %显示原始数据曲线图
/ `7 P1 s+ W3 w# Q, F5 ^axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
xlabel('时间 (s)');
4 v; x# \# Y! S! t6 \  \( u4 Jylabel('被测变量y');
4 v8 F  i  x+ s: U  b1 Q0 E% ktitle('原始信号(时域)');
! |/ D( U* e# V, Vgrid on;

%傅立叶变换
' X0 C* ]+ h0 q" By=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
7 v# S) V' U# U: z9 t9 rFs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
N=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
8 p4 Y6 u/ L8 L8 o  Uz=fft(y);

%频谱分析
7 _2 ~* c0 M4 Q- w: f2 X8 v' Vf=(0:N-1)*Fs/N;
  K  D  `$ {) W) UMag=2*abs(z)/N;                                        %幅值，单位同被测变量y
1 i7 {7 U1 G$ W2 QPyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式

1 a6 D2 F$ D: U. \%显示频谱图(频域)
8 y4 F+ e" G( wsubplot(2,1,2)
& b- l& w# m# l1 a! b7 q% oplot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
%                 |
8 J3 _" ^5 _8 S9 \- m%             将这里的Pyy改成Mag就是 幅值-频率图了
) a7 @! D& x4 uaxis([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)))])
: ~/ u! c6 g; c/ fxlabel('频率 (Hz)')
ylabel('能量')
4 c% ]4 D) ?/ u# Ltitle('频谱图(频域)')
, d- q$ S+ w$ q4 w  o0 |grid on;
- W/ f0 n8 ]' g3 t- I
%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
2 Q; T  f0 f4 O( H" ^! a! Nfprintf('\n傅立叶变换结果：\n')
5 q  C* m; [- X" Xfprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期
* Q# }2 B5 j7 H. j

yin123

离散信号MATLAB频谱分析程序