离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
7 f1 {6 W  H+ c%这里的向量都用行向量，假设被测变量是速度，单位为m/s
clear;
' c7 C; u0 D4 z9 `! H8 Pclose all;
- R6 p4 n5 }9 V. L# `+ B6 |/ j" ^
load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
) m& G( x, V" I4 z6 ~6 {A=data;                                        %将测量数据赋给A，此时A为N×2的数组
0 L% [7 p8 E9 g. ?. I0 U7 Y  Y0 D: bx=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
' n0 W, \" m$ a1 ?+ Yy=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
y=y';                                           %将列向量变成行向量
/ R/ f9 u, A4 m& F, L
+ q/ H- r- Z$ o%显示数据基本信息
fprintf('\n数据基本信息：\n')
5 F; o! g0 B9 q6 pfprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
5 u. N) C" r0 A6 Y& ^" Lfprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
! r* b2 Y- Y4 s* t$ tfprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
fprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
. m, f1 F" i2 v% S/ R- ]& lfprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
  F$ A. b. [& c6 z( b% G( jfprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
4 Y( x; Z8 i, w& k, M6 Dfprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
fprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
! Q% p: l: d* \9 C1 Efprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
fprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
! c/ t! W& ~( Q& f: b  
%显示原始数据曲线图（时域）
( H, I$ g7 Z5 Wsubplot(2,1,1);
plot(x,y)                                                                                %显示原始数据曲线图
axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
xlabel('时间 (s)');
: W2 {+ R) t) f) R8 lylabel('被测变量y');
1 Y0 _3 [2 J, |- `' ititle('原始信号(时域)');
grid on;
/ [) ^) P: d3 D' C+ ^
%傅立叶变换
3 t4 ?9 c9 i7 v4 Xy=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
Fs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
9 C( W* S! I$ e2 ?8 pN=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
* i; V3 O0 C2 [9 M' yz=fft(y);

%频谱分析
4 H# m# X( C5 m( M9 K! D: X2 @; Kf=(0:N-1)*Fs/N;
Mag=2*abs(z)/N;                                        %幅值，单位同被测变量y
Pyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
: W9 q# q6 {! P0 }: O1 q3 N( X
%显示频谱图(频域)
3 o2 G/ H5 N0 N& zsubplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
- \  m, a/ \$ k  L6 e2 \! L" X( h%                 |
%             将这里的Pyy改成Mag就是 幅值-频率图了
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)))])
  a+ F8 q( L% @2 pxlabel('频率 (Hz)')
/ G& j5 ^5 C% ~1 Tylabel('能量')
1 w2 Q0 p$ Q& C8 etitle('频谱图(频域)')
. L2 l- h5 ?( r8 \! k8 Pgrid on;

%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
fprintf('\n傅立叶变换结果：\n')
6 ~+ J. G' g' O: T  ?) }* f, Ofprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期
0 M  ]* a- j+ L9 i( C% t

yin123

离散信号MATLAB频谱分析程序