离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
clear;
0 U2 V! i$ ^) w: K" V$ y0 N$ sclose all;

load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
" z2 C" U6 n) M0 TA=data;                                        %将测量数据赋给A，此时A为N×2的数组
- J% k4 ^1 Y1 e. }1 hx=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
: s/ S( l# A' A9 Z% ^6 R" G9 y( Jx=x';                                           %将列向量变成行向量
* ~' h* n$ R5 S! ny=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
y=y';                                           %将列向量变成行向量
4 b& y, N1 g9 P4 c" j
%显示数据基本信息
% @3 Q' a" s& s. i5 ~fprintf('\n数据基本信息：\n')
fprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
2 _9 }1 W5 f+ J, s+ C+ d' K" efprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
/ [, i  Q# y4 j# Kfprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
fprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
3 Q  o# Z) M6 Xfprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
fprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
7 R' B( r+ y' ~5 B  f% }2 a# Qfprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
4 }9 t0 f; I0 h7 ^$ c+ L  v; Vfprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
. a) G5 E% ~8 c* w1 t  u1 afprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
# X4 \7 a# B+ W  
%显示原始数据曲线图（时域）
subplot(2,1,1);
# c; P( ~  `, V6 kplot(x,y)                                                                                %显示原始数据曲线图
& n8 @, V1 ~5 `6 `0 vaxis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
xlabel('时间 (s)');
9 I1 K- S9 D4 s9 Y4 }ylabel('被测变量y');
title('原始信号(时域)');
grid on;

%傅立叶变换
. \( R9 Y! y" N# `  G2 @2 By=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
: H- M& ~0 o! `9 }7 jFs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
- U& ~- M2 P7 S% c: QN=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
, e, o8 R* s, ]& ]3 f: C: O- ]z=fft(y);

%频谱分析
f=(0:N-1)*Fs/N;
& m- e6 A3 u) M% {  IMag=2*abs(z)/N;                                        %幅值，单位同被测变量y
: l4 c1 z% Y' Z5 y3 D) UPyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式

%显示频谱图(频域)
subplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
%                 |
%             将这里的Pyy改成Mag就是 幅值-频率图了
1 o3 q/ O2 L, D4 kaxis([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)))])
  ]4 S1 E, F+ d) s# Y) n) U2 A# y- axlabel('频率 (Hz)')
. {  R* E' \, H' S  R2 p4 y  S8 Kylabel('能量')
title('频谱图(频域)')
grid on;

%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
6 J# f9 G8 }& Y  c/ @3 Jfprintf('\n傅立叶变换结果：\n')
fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期

yin123

离散信号MATLAB频谱分析程序