% 如果要处理很多数据,就要把数据分段,分别作变换后,平均。
%
% 比如,要做400Hz内的频响,分辨率为1600线,即需要800Hz的采样频率,4分钟的测量数据,共3200个测量点。
% 用函数来实现作频响。
clc
close
clear
t_time = 10;
Fs = 800;
[t,x_state,y_out] = sim('tt_1',[0,t_time]);
x = y_out(:,1); % excitation data
y = y_out(:,2); % response data
span = Fs/2; % span of fft
lines = 6400; % lines of fft
L_data = lines*2; % 每一段的数据长度
L_n = round(length(y)/L_data); % 总数据共分段数
for kk = 1:L_n
x_s = x(L_data*(kk-1)+1:min(L_data*kk,length(x)));
y_s = y(L_data*(kk-1)+1:min(L_data*kk,length(y)));
[f,H_s(kk,:)] = cal_data_FRF(x_s,y_s,span,lines);
end
for kk = 1:length(H_s(1,:))
H(kk) = mean(H_s(:,kk));
end
plot(f,H)
function [f,H] = cal_data_FRF(x,y,span,lines)
% calculate the frequency response function of the experiment data
% x, is the excitation signal
% y, is the response signal
% span, is the span of frequency of FFT ,for exzample 400,800,1600...
% lines, is the resolving of the FFT
Fs = span*2; % Sampling frequency
T = 1/Fs; % Sampling time
L = lines*2; % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of signal
X = fft(x,NFFT)/L; % fft translation of x
Y = fft(y,NFFT)/L; % fft translation of y
HH = (Y.*conj(X))./(X.*conj(X)); % FRF
f = Fs/2*linspace(0,1,NFFT/2); % x axis
H = abs(HH(1:NFFT/2)); % y axis
有时候做出的曲线不光滑,大家一般怎么处理的?
欢迎讨论
skynew2005@gmail.com |