EDA365电子论坛网

标题: 离散信号MATLAB频谱分析程序 [打印本页]

---

作者: baqiao    时间: 2020-5-9 10:25
标题: 离散信号MATLAB频谱分析程序
%FFT变换,获得采样数据基本信息，时域图，频域图
" @3 c$ g% @, O5 U. w: @%这里的向量都用行向量，假设被测变量是速度，单位为m/s
) J& G5 i6 a# U+ t; Rclear;
close all;
$ Y1 j5 p3 X* N8 N# J; ]
) E* J/ \5 y$ _+ B' G/ k. @load data.txt              %通过仪器测量的原始数据，存储为data.txt中，附件中有一个模版(该信号极不规则)
A=data;                                        %将测量数据赋给A，此时A为N×2的数组
5 D1 W! j; Q! Y" px=A(:,1);                                     %将A中的第一列赋值给x，形成时间序列
x=x';                                           %将列向量变成行向量
, o* u* |2 a- I6 T( n4 O. p$ a1 L. Py=A(:,2);                                     %将A中的第二列赋值给y，形成被测量序列
5 w3 O5 Z, u' A2 K* L* Y) s* vy=y';                                           %将列向量变成行向量

, e" m' _9 b/ s%显示数据基本信息
2 A, T) N+ t* a) vfprintf('\n数据基本信息：\n')
" J7 S2 Y4 K# h* w4 v! l% `& Hfprintf('        采样点数 = %7.0f \n',length(x))                         %输出采样数据个数
: C6 W1 n5 A( p4 x5 Ifprintf('        采样时间 = %7.3f s\n',max(x)-min(x))                    %输出采样耗时
) ^' o. s4 n. c" M- `) Ofprintf('        采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x)))   %输出采样频率
fprintf('        最小速度 = %7.3f m/s\n',min(y))                         %输出本次采样被测量最小值
$ x2 i- W- f0 F- H) Ffprintf('        平均速度 = %7.3f m/s\n',mean(y))                      %输出本次采样被测量平均值
3 g2 d* M- j3 N2 C) X2 G3 y, u) A* efprintf('        速度中值 = %7.3f m/s\n',median(y))                   %输出本次采样被测量中值
fprintf('        最大速度 = %7.3f m/s\n',max(y))                          %输出本次采样被测量最大值
( N* |/ ~# x) V0 C( p# ^8 ufprintf('        标准方差 = %7.3f \n',std(y))                               %输出本次采样数据标准差
8 P6 W1 ^' B' X* S; ^6 ?/ n% Vfprintf('       协 方 差 = %7.3f \n',cov(y))                                %输出本次采样数据协方差
fprintf('     自相关系数 = %7.3f \n\n',corrcoef(y))                       %输出本次采样数据自相关系数
  
! p' A4 E2 X7 _# c3 _: ~5 h%显示原始数据曲线图（时域）
subplot(2,1,1);
4 g/ J: N7 F! ~  G1 Lplot(x,y)                                                                                %显示原始数据曲线图
" Z7 r* w8 F% s  f% V% qaxis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))])             %优化坐标，可有可无
7 |# N, ]* e/ t- G" Z# S7 Nxlabel('时间 (s)');
( {9 i4 D) }# T9 B# Aylabel('被测变量y');
title('原始信号(时域)');
$ O. b' @; Z3 c: Igrid on;

%傅立叶变换
; d0 X* w& |0 zy=y-mean(y);                                               %消去直流分量，使频谱更能体现有效信息
( c  Q1 \+ L; U1 U$ f  VFs=2000;                %得到原始数据data.txt时，仪器的采样频率。就是length(x)/(max(x)-min(x));     
$ t+ g; @5 N: q, XN=10000;                                                 %data.txt中的被测量个数，即采样个数。其实就是length(y);
9 s7 I$ l) s9 a& X3 rz=fft(y);

7 C- P. O# z5 K% H0 t%频谱分析
1 ~" p7 S- r/ i$ _' r+ Of=(0:N-1)*Fs/N;
Mag=2*abs(z)/N;                                        %幅值，单位同被测变量y
Pyy=Mag.^2;          %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故这里有很多表达方式

%显示频谱图(频域)
subplot(2,1,2)
plot(f(1:N/2),Pyy(1:N/2),'r')                         %显示频谱图
%                 |
%             将这里的Pyy改成Mag就是 幅值-频率图了
8 c/ Z$ B6 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)))])
+ E0 j0 N$ h: G9 A9 G# Sxlabel('频率 (Hz)')
ylabel('能量')
# Y& A, a* O1 s% L* w5 c4 xtitle('频谱图(频域)')
* B( U' [/ X$ x( o# M, s3 q( Igrid on;

%返回最大能量对应的频率和周期值
+ L( |; Y& B2 T8 C& f  d[a b]=max(Pyy(1:N/2));
# d. V4 I' a1 S! Pfprintf('\n傅立叶变换结果：\n')
fprintf('           FFT_f = %1.3f Hz\n',f(b))             %输出最大值对应的频率
fprintf('           FFT_T = %1.3f s\n',1/f(b))          %输出最大值对应的周期
9 ]6 r6 a% s" d4 ?  S* ^

---

作者: yin123    时间: 2020-5-9 13:23
离散信号MATLAB频谱分析程序

---

欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)

Powered by Discuz! X3.2