FFT plotting trouble: Concatenated inputs
显示 更早的评论
I am currently working on tap testing and am running into an issues with my FFT now that I am trying to Concatenate them to improve the "resolution". The plots max values fit what I expect however i need it to be a line and not have all the points underneath. I'm uncertian if my issue is in my measurements or in my code.
clear all
close all
clc
% each channel is notated with _n where n is the channel
% Hammer is on channel 1
% Accelerometer is on channel 2
% dataset_plot(11,'y')
% dataset_plot(12,'y')
% dataset_plot(13,'y')
%
% dataset_plot(21,'y')
dataset_plot(22,'y')
The goal is to have a single clear line so that I can create the transfer function to describe the system.
function dataset_plot(n,direction)
%% data Loading
out_2 = [];
out_1 = [];
for f_n = 1:5
fprintf('loading '+ string(n) + '_' + direction + '_' + string(f_n) + '_1.csv'+'\n')
raw_1 = table2array(readtable(string(n)+'_' + direction + '_' + string(f_n) + '_1.csv'));% read in hammer
raw_cut_1 = raw_1(3:height(raw_1),:); %remove header
if f_n == 1 %if its the first loop find signal length for window
window = hanning(length(raw_cut_1));
sample_rate = mean(diff(raw_cut_1(:,1))); %time between samples
Fs = 1/sample_rate; %Frequency in HZ
end
% Zero correction
Hammer_Zeroing = raw_cut_1(1:(length(raw_cut_1(:,2))*0.2)-1);%fix the offset
Hammer_adj = mean(Hammer_Zeroing);
raw_cut_1(:,2) = raw_cut_1(:,2) - (Hammer_adj);
%concatenate
out_1 = [out_1; raw_cut_1(:,2).*hanning(length(raw_cut_1))];
raw_2 = table2array(readtable(string(n)+'_' + direction + '_' + string(f_n) + '_2.csv'));% read in accel
raw_cut_2 = raw_2(3:height(raw_1),:); %remove header
%Zero correction
Accel_Zeroing = raw_cut_2(1:(length(raw_cut_2(:,2))*0.2)-1);%fix the offset
Accel_adj = mean(Accel_Zeroing);
raw_cut_2(:,2) = raw_cut_2(:,2) - (Accel_adj);
%concatenate
out_2 = [out_2; raw_cut_2(:,2).*hanning(length(raw_cut_2))];
clear raw_cut_2 raw_cut_1 raw_1 raw_2 padding Hammer_adj Hammer_Zeroing Accel_adj Accel_Zeroing
end
%% create X axis
Hammer = out_1;%(:,2);
Accel = out_2;%(:,2);
x_axis = linspace(0,sample_rate*length(Accel),length(Accel));
Accel = vertcat(Accel,ones(length(x_axis)-length(Accel),1).*mean(diff(Accel)));
meas_time_2 = sample_rate*length(x_axis);
%% run FFT functions
[H_f,H_P1] = FFT_ploting(smoothdata(Hammer,"gaussian",20),meas_time_2);
[A_f,A_P1] = FFT_ploting(smoothdata(Accel,"gaussian",20),meas_time_2);
%% fix units
H_P1 = H_P1./0.00225; %convert from V to N @ 2.25 mV/N
A_P1 = A_P1./0.0102; %convert from V to m/s^2 10.2 mV/(m/s²)
%% plot
close all
figure('Name','raw data check')
hold on
plot(Accel)
plot(smoothdata(Accel,"gaussian",10))
figure('Name','FFT Plots')
hold on
yyaxis left
plot(H_f,H_P1);
ylabel 'lbf'
yyaxis right
plot(A_f,A_P1);
ylabel 'm/s^2'
grid on
xlim([5 3000]); %accelerometer limit of 3000Hz
%we are focusing on 5 to 500 Hz
title('Accelerometer & Hammer FFT')
legend('Hammer FFT','Acceleromter FFT')
xlabel('Hz');
end
function [f,P1] = FFT_ploting(fun_data,meas_time)
% count = count + 1; %FFT itteration count
Fs = (meas_time/length(fun_data))^-1; %frequecy in Hz
T = 1/Fs; %meas time between samples (inverse)
L = length(fun_data);
t = (0:L-1)*T;
f = Fs*(0:(L-1)/2)/L;
fft_data = fft(fun_data,L);
P2 = (abs(fft_data/L));
P1 = P2(1:(L+1)/2);
P1(2:end) = 2*P1(2:end);
end
采纳的回答
hello
when you do modal analysis / FRF measurements with a hammer , you should split the record into individual hit input and output signals , then do the averaging of the FRFs (and coherence)
I may not directly give the code for your data , but the process is described here , based on the simulated response of a mass / spring / dash pot system. You can see the measured FRF pretty much overlays with the theoretical Bode plot of the system
%% Modal Analysis / hammer test / FRF acquisition and FRF computation tool
%% FFT processing parameters
Fs = 1024*4;
nfft = 1024*2;
channels = 1;
% dummy hammer test data - or real data if you have
Nhits = 8;
duration_hits = 1; % mean time duration between hits (in seconds)
% Mass / spring / dash pot model
M = 25; % kgs
fn = 100; % Hz (natural frequency)
K = M*(2*pi*fn)^2;
C = 0.05*sqrt(K*M);
% analog model acceleration vs force
num = [1 0 0];
den = [M C K];
% discretization
dt = 1/Fs;
samples_one_hit = duration_hits*Fs;
[numd,dend] = c2dm(num,den,1/Fs,'tustin');
% create some dummy hammer gaussian pulses with variable force amplitude
f0 = Fs/2.5; % pulse frequency
Var_Sq = 1; % pulse variance (squared)
t = dt*(0:samples_one_hit-1)';
t_pulse = max(t)/8;
offset = 0.5*duration_hits*(rand(Nhits,1)); % random time shifts
amplitude = 10*(1 + 1*(rand(Nhits,1))); % random amplitude
force_signal = [];
for k = 1:Nhits
pulse_signal = amplitude(k)*exp(-(((t-t_pulse-offset(k)).^2).*(f0^2))./Var_Sq);
force_signal = [force_signal;pulse_signal];
end
% now all samples are concatenated to make one signal with random delta t
% between pulses
samples = numel(force_signal);
time = dt*(0:samples-1)';
% apply model to force signal (input) to create acceleration signal (output)
accel_signal = filter(numd,dend,force_signal);
accel_signal = accel_signal + 0.001*randn(size(accel_signal));
% now let's trigger the data (using "zero crossing" on the force signal)
% this is basically the same way a FFT/network analyser operates
pre_trig_duration = 10e-3; % add some time to shift the pulse signal to the right
force_signal_level = 1;
[Zx] = find_zc(time,force_signal,force_signal_level);
Zx = Zx - pre_trig_duration;
%recommended max nfft
recommended_max_nfft = floor(min(Fs*diff(Zx)));
%% check valid segments (duration in samples must be above nfft
% pulses durations in samples :
pd = [Fs*diff(Zx); samples-Fs*Zx(end)];% add the last pulse duration :
data_valid = find(pd>=nfft);
figure(2),
subplot(2,1,1),plot(time,force_signal)
% warning message
if nfft>recommended_max_nfft
text(max(time)*0.1,max(force_signal)*0.9,['Selected nfft = ' num2str(nfft) ' is above recommended max value : ' num2str(recommended_max_nfft)]);
end
subplot(2,1,2),plot(time,accel_signal)
%% FFT processing
% windowing
windowx = ones(nfft,1); % no window on force signal
windowy = ones(nfft,1); % no window on response sensor signal
% windowy = exp(-3*(0:nfft-1)/nfft)'; % 0.1 exp window on response sensor signal
% initialisation
Pxx = zeros(nfft,1);
Pxy = zeros(nfft,channels);
Pyy = zeros(nfft,channels);
% check if data_valid is not empty (can happen if you choose nfft grossly
% oversized)
if isempty(data_valid)
error('No data or nfft too large - check data & nfft settings');
else % we have data and nfft is correctly defined
for ck = 1:numel(data_valid)
% do it only for valid segments
k = data_valid(ck);
% select time data
ind_start = find(time>=Zx(k),1,'first');
% check valid segments (duration in samples must be above nfft
ind = (ind_start :ind_start + nfft -1);
time_measure = time(ind);
time_measure = time_measure - time_measure(1);
force_signal_measure = force_signal(ind);
accel_signal_measure = accel_signal(ind);
figure(3), % just to show the individual gaussian pulses with different time offsets
subplot(2,1,1),plot(time_measure,force_signal_measure)
hold on
subplot(2,1,2),plot(time_measure,accel_signal_measure)
hold on
leg_str1{ck} = ['hit # ' num2str(k)];
leg_str2{ck} = ['record # ' num2str(k)];
% FFT processing
x = force_signal_measure;
y = accel_signal_measure;
xw = windowx.*x;
yw = (windowy*ones(1,channels)).*y;
Xx = fft(xw,nfft);
Yy = fft(yw,nfft,1);
Xx2 = abs(Xx).^2;
Xy2 = Yy.*(conj(Xx)*ones(1,channels));
Yy2 = abs(Yy).^2;
Pxx = Pxx + Xx2;
Pxy = Pxy + Xy2;
Pyy = Pyy + Yy2;
end
subplot(2,1,1),legend(leg_str1); % add legend of valid data
subplot(2,1,2),legend(leg_str2); % add legend of valid data
% Select first half
if ~any(any(imag([x y])~=0)) % if x and y are not complex
if rem(nfft,2) % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
Pxx = Pxx(select);
Pxy = Pxy(select,:);
Pyy = Pyy(select,:);
else
select = 1:nfft;
end
% set up output parameters
Txy = Pxy ./ (Pxx*ones(1,channels)); % transfer function estimate
Cxy = (abs(Pxy).^2)./((Pxx*ones(1,channels)).*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end
%% plot model vs identified FRF
flow = 10;
fhigh = Fs/2.56;
freq = logspace(log10(flow),log10(fhigh),300);
[gm,pm] = bode(num,den,2*pi*freq);
figure(4)
subplot(3,1,1),semilogx(freq,20*log10(gm),'b',f,20*log10(abs(Txy)),'r');
title(['FRF comparison / Valid records : ' num2str(data_valid')])
xlim([flow fhigh]);
ylabel('Magnitude (dB)');
legend('input model','measured');
subplot(3,1,2),semilogx(freq,pm,'b',f,180/pi*(angle(Txy)),'r');
xlim([flow fhigh]);
ylabel('Phase (°)');
legend('input model','measured');
subplot(3,1,3),semilogx(freq,ones(size(freq)),'b',f,Cxy,'r');
xlim([flow fhigh]);
xlabel('Frequency (Hz)');
ylabel('Coh');
legend('input model','measured');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)
% find time values of positive slope zero (or "threshold" value ) y
% crossing points
% put data in rows
x = x(:);
y = y(:);
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0); %define function: returns indices of +ZCs
ix=zci(y); %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Vibration Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
