离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
* X: l- E8 z0 D: ?( L%这里的向量都用行向量，假设被测变量是速度，单位为m/s
clear;
5 H* d* O3 d; ^& a6 n% }% ~" gclose all;
; h! B5 F/ X& a, d6 W$ W% Q  s6 J
9 l' Z% q( T4 `; P* sload data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
  V9 ^( X- y3 M' C' \A=data;                                        %将测量数据赋给A，此时A为N×2的数组
x=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
/ U* Y! S4 [9 W# d- D! W+ [+ I4 `y=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
4 x8 y3 l- R+ m, d2 x5 wy=y';                                           %将列向量变成行向量
% l( @( v) |9 I, g& L
%显示数据基本信息
4 f% w0 G  G. E% @; [fprintf('\n数据基本信息：\n')
fprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
4 g$ a* C: c9 J0 ufprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
( ~- y) v9 M" Yfprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
; `* }0 h) [& u; Mfprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
fprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
fprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
9 D, S& A3 J0 S1 x  Mfprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
2 ?7 T! U* z# U% |fprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
fprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
* o  _; |! f1 S3 ^2 ^! }9 m: l9 t' f2 ffprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
6 N2 U( H* h$ ~' T& E& ]- {6 \  
%显示原始数据曲线图（时域）
3 j. d4 s8 @% e. i  g* f% {subplot(2,1,1);
; C; |# m# b2 k: d( c5 Pplot(x,y)                                                                                %显示原始数据曲线图
axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
* w( P5 U5 |3 D: r5 xxlabel('时间 (s)');
- X! t8 P$ L2 @ylabel('被测变量y');
% ~( v& e2 u) e$ p- W& Otitle('原始信号(时域)');
' s" e+ U" K, [grid on;

4 Y( N4 }4 L) w" N0 V& E%傅立叶变换
y=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
* ]  ~8 C8 t. }- qFs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
: W8 T/ z/ }/ Q2 @* i, a! a, a9 TN=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
% c, n" T- T, k) Nz=fft(y);

- J" E* v; F* N# r" f5 Q0 G/ V, h%频谱分析
f=(0:N-1)*Fs/N;
7 ]2 _4 j4 y% `6 U% X' _! X, b4 ]7 xMag=2*abs(z)/N;                                        %幅值，单位同被测变量y
+ W- h3 g7 ^" L' x2 P+ g9 {& SPyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
/ k- \, ]* _3 Z$ `- t! W
, g/ |/ o4 `2 {5 W%显示频谱图(频域)
5 `' y5 ]# h7 V; R) s7 |5 Jsubplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
%                 |
%             将这里的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)))])
xlabel('频率 (Hz)')
, O! V; y6 |* q. fylabel('能量')
title('频谱图(频域)')
grid on;

  E9 U9 H- v. c3 c%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
- G% P) I# e9 y. r  p! Z/ s! ffprintf('\n傅立叶变换结果：\n')
( i" p3 w; T) R8 ]: D; `6 e6 G6 D) @fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期

yin123

离散信号MATLAB频谱分析程序