离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
0 K; r1 x2 o0 C7 t$ Q& g' D( F; Rclear;
close all;

, ?& W  L( `: Q4 q! h' pload data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
A=data;                                        %将测量数据赋给A，此时A为N×2的数组
x=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
/ q; N. k9 e7 }' ^1 Sx=x';                                           %将列向量变成行向量
# K3 {+ j% E# Wy=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
y=y';                                           %将列向量变成行向量
* t5 x) q/ b* \* F% N2 o
%显示数据基本信息
& E9 ]& _0 [. j8 d" ?" T) w5 Vfprintf('\n数据基本信息：\n')
6 S) i3 w- q+ @/ sfprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
% ?# \; m9 ]$ r3 @4 U, C. C8 A% Rfprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
% y3 |* a  }2 S" u6 Y; X% K2 {# Q# rfprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
fprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
  Y/ W. A$ c9 rfprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
fprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
+ F, R: O5 J% Q1 Efprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
fprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
fprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
  
. F0 b9 o* t# W& H/ G%显示原始数据曲线图（时域）
! e% j" `+ B8 @- i7 T! {! W% esubplot(2,1,1);
plot(x,y)                                                                                %显示原始数据曲线图
axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
9 |: ]5 o% N9 z! [' f0 n' fxlabel('时间 (s)');
4 i  Y( g) Y' p0 E2 {3 q1 O! q2 v' dylabel('被测变量y');
% R' _2 ]/ \: a& Ctitle('原始信号(时域)');
" d. R1 s7 o' h4 _grid on;

%傅立叶变换
y=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
, r, ~- u3 p* E; S0 q4 dFs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
7 W- w( G# U8 uN=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
8 P2 _/ L4 H( m( ~6 @z=fft(y);
- F' A4 R: z1 {& j; e
%频谱分析
$ b- U+ N( F. s4 Y! o# ]f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N;                                        %幅值，单位同被测变量y
5 U% M# \  o' ]! I+ ePyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
% Q! X* F% s4 O7 ~. h3 {; d
  F4 i. n) z6 X! B%显示频谱图(频域)
4 \) F3 ~6 c2 h" G, @subplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
5 f; k9 X8 y% e) Z, u  \%                 |
. k+ Z, ~  e! k# }%             将这里的Pyy改成Mag就是 幅值-频率图了
$ U! V. x1 P/ ~9 u& v4 y$ x. Qaxis([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)))])
/ p7 Q( b# X  p$ [. X3 lxlabel('频率 (Hz)')
6 r; x" K# L" \( e* Vylabel('能量')
title('频谱图(频域)')
3 w& W: H, f- |  lgrid on;
0 R0 q; k) q! q$ E$ C0 N6 a/ H
% O' S8 l* I" [. y( q7 {; l%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
fprintf('\n傅立叶变换结果：\n')
fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期
0 S4 @; T+ p6 Z

yin123

离散信号MATLAB频谱分析程序