|
|
MATLAB頻譜分析程序
%FFT變換,獲得采樣數(shù)據(jù)基本信息，時(shí)域圖，頻域圖
%這里的向量都用行向量，假設(shè)被測(cè)變量是速度，單位為m/s
clear;
close all;
load data.txt %通過(guò)儀器測(cè)量的原始數(shù)據(jù)，存儲(chǔ)為data.txt中，附件中有一個(gè)模版(該信號(hào)極不規(guī)則)
A=data; %將測(cè)量數(shù)據(jù)賦給A，此時(shí)A為N×2的數(shù)組
x=A(:,1); %將A中的第一列賦值給x，形成時(shí)間序列
x=x'; %將列向量變成行向量
y=A(:,2); %將A中的第二列賦值給y，形成被測(cè)量序列
y=y'; %將列向量變成行向量
%顯示數(shù)據(jù)基本信息
fprintf('\n數(shù)據(jù)基本信息：\n')
fprintf(' 采樣點(diǎn)數(shù) = %7.0f \n',length(x)) %輸出采樣數(shù)據(jù)個(gè)數(shù)
fprintf(' 采樣時(shí)間 = %7.3f s\n',max(x)-min(x)) %輸出采樣耗時(shí)
fprintf(' 采樣頻率 = %7.1f Hz\n',length(x)/(max(x)-min(x))) %輸出采樣頻率
fprintf(' 最小速度 = %7.3f m/s\n',min(y)) %輸出本次采樣被測(cè)量最小值
fprintf(' 平均速度 = %7.3f m/s\n',mean(y)) %輸出本次采樣被測(cè)量平均值
fprintf(' 速度中值 = %7.3f m/s\n',median(y)) %輸出本次采樣被測(cè)量中值
fprintf(' 最大速度 = %7.3f m/s\n',max(y)) %輸出本次采樣被測(cè)量最大值
fprintf(' 標(biāo)準(zhǔn)方差 = %7.3f \n',std(y)) %輸出本次采樣數(shù)據(jù)標(biāo)準(zhǔn)差
fprintf(' 協(xié) 方 差 = %7.3f \n',cov(y)) %輸出本次采樣數(shù)據(jù)協(xié)方差
fprintf(' 自相關(guān)系數(shù) = %7.3f \n\n',corrcoef(y)) %輸出本次采樣數(shù)據(jù)自相關(guān)系數(shù)
%顯示原始數(shù)據(jù)曲線(xiàn)圖（時(shí)域）
subplot(2,1,1);
plot(x,y) %顯示原始數(shù)據(jù)曲線(xiàn)圖
axis([min(x) max(x) 1.1*floor(min(y)) 1.1*ceil(max(y))]) %優(yōu)化坐標(biāo)，可有可無(wú)
xlabel('時(shí)間 (s)');
ylabel('被測(cè)變量y');
title('原始信號(hào)(時(shí)域)');
grid on;
%傅立葉變換
y=y-mean(y); %消去直流分量，使頻譜更能體現(xiàn)有效信息
Fs=2000; %得到原始數(shù)據(jù)data.txt時(shí)，儀器的采樣頻率。就是length(x)/(max(x)-min(x));
N=10000; %data.txt中的被測(cè)量個(gè)數(shù)，即采樣個(gè)數(shù)。其實(shí)就是length(y);
z=fft(y);
%頻譜分析
f=(0:N-1)*Fs/N;
Mag=2*abs(z)/N; %幅值，單位同被測(cè)變量y
Pyy=Mag.^2; %能量；對(duì)實(shí)數(shù)系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2，故這里有很多表達(dá)方式
%顯示頻譜圖(頻域)
subplot(2,1,2)
plot(f(1:N/2),Pyy(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;
%返回最大能量對(duì)應(yīng)的頻率和周期值
[a b]=max(Pyy(1:N/2));
fprintf('\n傅立葉變換結(jié)果：\n')
fprintf(' FFT_f = %1.3f Hz\n',f(b)) %輸出最大值對(duì)應(yīng)的頻率
fprintf(' FFT_T = %1.3f s\n',1/f(b)) %輸出最大值對(duì)應(yīng)的周期 |
|
回復(fù)此樓
