离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
9 h1 O8 _: P6 X, v; E6 Dclear;
close all;

load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
A=data;                                        %将测量数据赋给A，此时A为N×2的数组
/ u1 ~" U  u3 @; a. Vx=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
y=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
+ {- J" E1 ?& H* A/ J$ e. b4 Ry=y';                                           %将列向量变成行向量

. G+ D3 z+ j2 K* z; w$ {3 j* W5 E* ?%显示数据基本信息
5 Y* J/ A) g- k" h" B- |; Sfprintf('\n数据基本信息：\n')
/ R0 y" `, E; ufprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
fprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
% Z# z2 y6 x1 b4 \$ j4 bfprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
: b0 L# `2 l& j$ a/ ?5 ^& D1 r. Cfprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
9 Q; \7 ?0 `& H- D5 ^9 L  qfprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
- f- }9 s; q6 n; T; H6 cfprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
# V1 k" a, g4 R5 v, i$ Nfprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
fprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
# E! X( S) c" f; k* v0 V3 dfprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
  
: @8 q3 B8 i# o2 D. H8 z- o%显示原始数据曲线图（时域）
& h8 d( ]. W; \1 usubplot(2,1,1);
# {. X. R  e3 h: ]plot(x,y)                                                                                %显示原始数据曲线图
1 [) O( K7 c4 w. ]. y# vaxis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
" J+ W+ c/ N, V* C  uxlabel('时间 (s)');
, f) O# }- y& {ylabel('被测变量y');
title('原始信号(时域)');
  W9 x7 ]2 h0 A0 Z: m. L( \grid on;

- b1 _* t. s/ Y+ f: f+ _0 z) V4 |%傅立叶变换
y=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
& e0 k4 r4 s" f9 z, b+ zFs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
N=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
' j5 b% J3 `: B( _7 n) f4 P, Iz=fft(y);

1 Q3 p+ @) s+ Z( v: L- b9 ~5 A%频谱分析
* f" v7 i$ `4 Z0 ]1 [f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N;                                        %幅值，单位同被测变量y
+ z% U/ |+ ?; v4 A8 {Pyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式

. {5 m) H- v" n2 ~0 J%显示频谱图(频域)
subplot(2,1,2)
$ Y) `* H# G6 A3 M/ zplot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
& B% q' F) D/ W6 Q" m%                 |
%             将这里的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)')
: Y# Q! p& |: M3 i: y' D2 ?3 j5 wylabel('能量')
* \) P/ R& H8 t, F) o4 ftitle('频谱图(频域)')
grid on;

( A2 V0 s% y/ M$ R! O2 }, ^) r%返回最大能量对应的频率和周期值
[a b]=max(Pyy(1:N/2));
fprintf('\n傅立叶变换结果：\n')
fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
6 c( G: T6 U9 Efprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期
4 I$ f$ g& g: e

yin123

离散信号MATLAB频谱分析程序