|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
%FFT变换,获得采样数据基本信息,时域图,频域图* |& [9 Y! ]: n) W" z0 g# a2 u
%这里的向量都用行向量,假设被测变量是速度,单位为m/s/ k7 s, L1 T1 n" t0 j2 g
clear;
0 U2 V! i$ ^) w: K" V$ y0 N$ sclose all;3 T5 p y0 N) W7 E- }" O9 F: m, Q
4 s2 [; O; i) B. f! k) ~) C# _
load data.txt %通过仪器测量的原始数据,存储为data.txt中,附件中有一个模版(该信号极不规则)
" z2 C" U6 n) M0 TA=data; %将测量数据赋给A,此时A为N×2的数组
- J% k4 ^1 Y1 e. }1 hx=A(:,1); %将A中的第一列赋值给x,形成时间序列
: s/ S( l# A' A9 Z% ^6 R" G9 y( Jx=x'; %将列向量变成行向量
* ~' h* n$ R5 S! ny=A(:,2); %将A中的第二列赋值给y,形成被测量序列% {% y& v. `( b+ ~! q- I
y=y'; %将列向量变成行向量
4 b& y, N1 g9 P4 c" j7 C4 G8 e7 f7 i- I# S* W# v0 k
%显示数据基本信息
% @3 Q' a" s& s. i5 ~fprintf('\n数据基本信息:\n') 8 f; i- h: F8 A
fprintf(' 采样点数 = %7.0f \n',length(x)) %输出采样数据个数
2 _9 }1 W5 f+ J, s+ C+ d' K" efprintf(' 采样时间 = %7.3f s\n',max(x)-min(x)) %输出采样耗时
/ [, i Q# y4 j# Kfprintf(' 采样频率 = %7.1f Hz\n',length(x)/(max(x)-min(x))) %输出采样频率5 {3 b7 d- u9 {+ i
fprintf(' 最小速度 = %7.3f m/s\n',min(y)) %输出本次采样被测量最小值
3 Q o# Z) M6 Xfprintf(' 平均速度 = %7.3f m/s\n',mean(y)) %输出本次采样被测量平均值1 L" e7 ]3 e1 }. q, F
fprintf(' 速度中值 = %7.3f m/s\n',median(y)) %输出本次采样被测量中值* A s4 V: j) A& z
fprintf(' 最大速度 = %7.3f m/s\n',max(y)) %输出本次采样被测量最大值
7 R' B( r+ y' ~5 B f% }2 a# Qfprintf(' 标准方差 = %7.3f \n',std(y)) %输出本次采样数据标准差
4 }9 t0 f; I0 h7 ^$ c+ L v; Vfprintf(' 协 方 差 = %7.3f \n',cov(y)) %输出本次采样数据协方差
. a) G5 E% ~8 c* w1 t u1 afprintf(' 自相关系数 = %7.3f \n\n',corrcoef(y)) %输出本次采样数据自相关系数
# X4 \7 a# B+ W ! {' k6 \9 q0 ], `/ A) h# r) [( [: d1 ?
%显示原始数据曲线图(时域)# k# B" J( t V! V+ m. H$ N
subplot(2,1,1);
# c; P( ~ `, V6 kplot(x,y) %显示原始数据曲线图
& n8 @, V1 ~5 `6 `0 vaxis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))]) %优化坐标,可有可无5 J9 D) z: A' U. |* n8 O
xlabel('时间 (s)');
9 I1 K- S9 D4 s9 Y4 }ylabel('被测变量y'); J$ T* l. N7 w' c. F* |% D
title('原始信号(时域)');# m/ ?; g' u) ` U' C* S8 k5 u8 f
grid on;% L6 F! S6 y. L# h
3 @8 W( O- E. x" P. S5 n
%傅立叶变换
. \( R9 Y! y" N# ` G2 @2 By=y-mean(y); %消去直流分量,使频谱更能体现有效信息
: H- M& ~0 o! `9 }7 jFs=2000; %得到原始数据data.txt时,仪器的采样频率。就是length(x)/(max(x)-min(x));
- U& ~- M2 P7 S% c: QN=10000; %data.txt中的被测量个数,即采样个数。其实就是length(y);
, e, o8 R* s, ]& ]3 f: C: O- ]z=fft(y);7 p8 @4 D% f6 \0 [2 H! g
0 r- | p) K! J) @4 U: S" e
%频谱分析3 q5 W3 z# p/ `9 b5 ~6 X
f=(0:N-1)*Fs/N;
& m- e6 A3 u) M% { IMag=2*abs(z)/N; %幅值,单位同被测变量y
: l4 c1 z% Y' Z5 y3 D) UPyy=Mag.^2; %能量;对实数系列X,有 X.*X=X.*conj(X)=abs(X).^2=X.^2,故这里有很多表达方式- Y$ ?+ q: t+ C7 {8 D% \: b( P
1 ?5 s/ k# U- w& q
%显示频谱图(频域)1 K; d ]" s$ q7 H* g' j+ f
subplot(2,1,2)3 D& `, k7 ^6 x6 ]6 @4 j0 t
plot(f(1:N/2),Pyy(1:N/2),'r') %显示频谱图. o7 _+ m' Q* C# P5 a% l) t, u
% |! ^* ^2 u- P$ ^9 b3 s9 V4 w" {
% 将这里的Pyy改成Mag就是 幅值-频率图了
1 o3 q/ O2 L, D4 kaxis([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)))])
]4 S1 E, F+ d) s# Y) n) U2 A# y- axlabel('频率 (Hz)')
. { R* E' \, H' S R2 p4 y S8 Kylabel('能量')3 B5 ^' P- h5 V4 A; M. W
title('频谱图(频域)')% Q. `( Y' g! _
grid on;6 m1 p& c2 @+ h2 [
) Q( @5 z. K: n
%返回最大能量对应的频率和周期值) P3 y' @) r% F, U, q2 y
[a b]=max(Pyy(1:N/2));
6 J# f9 G8 }& Y c/ @3 Jfprintf('\n傅立叶变换结果:\n') # o4 ^( T2 E; [/ Y
fprintf(' FFT_f = %1.3f Hz\n',f(b)) %输出最大值对应的频率; D% Y$ G% R& ~
fprintf(' FFT_T = %1.3f s\n',1/f(b)) %输出最大值对应的周期" ]4 s6 h3 T. M7 A% G9 E, w
|
|