马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
function myfrequencyplot
% all the files are save as t and s
clc;
clear;
load result_6000_30.mat;
length(t)
%%%-----------------------Program 1----------------------------------
%FFT变换,获得采样数据基本信息,时域图,频域图
%这里的向量都用行向量,假设被测变量是速度,单位为m/s
x = linspace(t(1),t(end),1e2);% 时间序列
x=x'; %将列向量变成行向量
y = interp1(t,s(:,4),x);
%将s中的第4列赋值给y,形成被测量序列
y=y'; %将列向量变成行向量
%显示数据基本信息
fprintf('\n数据基本信息:\n')
fprintf(' 采样点数 = %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)) %输出本次采样被测量最小值
fprintf(' 平均速度 = %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)) %输出本次采样数据标准差
fprintf(' 协 方 差 = %7.3f \n',cov(y)) %输出本次采样数据协方差
fprintf(' 自相关系数 = %7.3f \n\n',corrcoef(y)) %输出本次采样数据自相关系数
%傅立叶变换
y=y-mean(y); %消去直流分量,使频谱更能体现有效信息
Fs=length(x)/(max(x)-min(x)); %得到原始数据data.txt时,仪器的采样频率。其实就是length(x)/(max(x)-min(x));
N=length(y); %data.txt中的被测量个数,即采样个数。其实就是length(y);
z=fft(y);
%频谱分析
f=(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,故这里有很多表达方式
%显示频谱图(频域)
plot(f(1:N/2),Mag(1:N/2),'r') %显示频谱图
% 将这里的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)')
ylabel('幅值')
title('频谱图(频域)')
grid on;
%返回最大能量对应的频率和周期值
[a b]=max(Mag(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)) %输出最大值对应的周期 |