离散信号MATLAB频谱分析程序

baqiao

%FFT变换,获得采样数据基本信息，时域图，频域图
%这里的向量都用行向量，假设被测变量是速度，单位为m/s
clear;
; b* S0 M- S8 O& ]- _: m. zclose all;
* a7 J/ }! |6 Q- H* \- n7 U; f
load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
A=data;                                        %将测量数据赋给A，此时A为N×2的数组
( x6 a% q# M' Ex=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
8 V" h; X6 K/ l( Jx=x';                                           %将列向量变成行向量
' C7 M% s% S/ B0 p9 B1 d2 oy=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
( z! e+ u9 i- r5 qy=y';                                           %将列向量变成行向量
2 T* j) ~0 ~" R! y2 R+ m
%显示数据基本信息
) X4 m! q4 |# R  ~& K* C. jfprintf('\n数据基本信息：\n')
9 F. D0 n- G( Xfprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
fprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
4 \9 j& E+ n; afprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
fprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
fprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
$ W+ @0 Q7 j# D4 b4 ffprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
: J" \9 f5 w2 L' A& Y3 `fprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
fprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
* ]5 ^1 @/ S/ Y. ifprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
  
3 W: A6 @! {9 E4 p  Y* }) {7 G% i%显示原始数据曲线图（时域）
+ P3 s( q/ i+ S, z9 C. msubplot(2,1,1);
plot(x,y)                                                                                %显示原始数据曲线图
8 C8 B  {1 n( I! g- }9 k. vaxis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
' k2 R6 ]/ _  t5 }- p1 K9 s$ \/ \3 \xlabel('时间 (s)');
. F! [2 ?5 L5 |# A: C5 Aylabel('被测变量y');
% M4 H. M" E' `title('原始信号(时域)');
, [4 k0 c$ Q; x* v/ D& L/ Vgrid on;

+ w& d" X( k5 t. x* P" H4 P%傅立叶变换
, N8 W) w' D& |! `6 oy=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
Fs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
- N. s* A3 U* o) a/ D4 @N=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
z=fft(y);

%频谱分析
f=(0:N-1)*Fs/N;
6 s9 \! y' J1 ]& D" BMag=2*abs(z)/N;                                        %幅值，单位同被测变量y
7 m3 Y! j0 K1 \: r+ vPyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式
& v6 u8 Q" Q* M
%显示频谱图(频域)
subplot(2,1,2)
) `7 g% D: N6 j& ]  N- Oplot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
6 u, h% @3 b" q4 P; M, l: c%                 |
%             将这里的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)')
! W5 H8 x9 u, B8 Q8 x# ?8 nylabel('能量')
3 b2 b! l' V' T/ v2 l0 i( ^title('频谱图(频域)')
grid on;
  m, C9 C3 _' n% u4 ?
%返回最大能量对应的频率和周期值
[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))          %输出最大值对应的周期
; o7 t# |% [, V

yin123

离散信号MATLAB频谱分析程序