离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
clear;
close all;
, o, L; h  z! i' y) p
load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
A=data;                                        %将测量数据赋给A，此时A为N×2的数组
5 J! Y. k, D8 _6 O) p2 V. ix=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
: k' `2 w- L8 `8 A+ A& [y=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
) q( k+ r- w. O0 j; cy=y';                                           %将列向量变成行向量

6 y- [; z. Q$ t8 w; ?* m  i%显示数据基本信息
% i7 M/ \% j% c  z9 m( b: Vfprintf('\n数据基本信息：\n')
" ^, [! @8 a; P4 M# Dfprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
fprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
fprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
fprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
- D: ^2 s* M8 q+ X! Q6 Afprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
fprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
fprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
6 w0 P2 W6 b5 f3 t. p1 \fprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
7 ]7 L# h. i9 P5 y5 t, K! M3 dfprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
- P( y) I& E: o' v4 R7 w4 Z: R  
% o; f* B) J; \%显示原始数据曲线图（时域）
subplot(2,1,1);
, z9 n# D5 ?- G. `/ }/ uplot(x,y)                                                                                %显示原始数据曲线图
axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
9 A5 K5 i6 V* n, y: V4 lxlabel('时间 (s)');
ylabel('被测变量y');
/ b8 Z3 l3 p' b9 H7 G' O8 Ititle('原始信号(时域)');
2 g6 a* E2 W, J- R7 ?grid on;
" G! n3 O2 ~1 ]: ~2 f! Q% r" @
1 o) E7 Y. ~+ Y%傅立叶变换
+ \  I; D4 S: h+ D! v1 }/ Oy=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
Fs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
N=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
z=fft(y);

%频谱分析
2 X# g* }& ^" I" o. D4 Wf=(0:N-1)*Fs/N;
Mag=2*abs(z)/N;                                        %幅值，单位同被测变量y
+ ~1 W5 y" s3 R; D  sPyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
3 Q0 ^: K% C- ?9 \8 m
0 {3 M' j% {+ u! f5 ~# |%显示频谱图(频域)
subplot(2,1,2)
8 P5 ?4 A2 V* R- J" Yplot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
%                 |
5 w! s) A* U( w/ P1 k%             将这里的Pyy改成Mag就是 幅值-频率图了
6 k, f+ a2 E( L4 q5 D7 {4 b* }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)))])
; ^" N+ Y1 v% f0 u/ Pxlabel('频率 (Hz)')
ylabel('能量')
title('频谱图(频域)')
grid on;
9 s5 ^5 D, T$ I' j  r
1 L2 [5 j5 N3 _$ x%返回最大能量对应的频率和周期值
. ]8 N0 I8 o5 M- f. Y[a b]=max(Pyy(1:N/2));
% |/ ]8 z2 ]$ M8 G6 I% G3 hfprintf('\n傅立叶变换结果：\n')
fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期

yin123

离散信号MATLAB频谱分析程序