![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1341409/image.png)
I would like to calculate the frequency response function (FRF) for my input frequency and output response from accelerometer
78 次查看(过去 30 天)
显示 更早的评论
Hello. I collected the data from two single axial acceleroemters located at definite distance (channel C and D). My input excitation was given through a piezo at 8434 Hz (chennel A). From the recorded data, required details were loaded and channels were loaded and multiplied by corresponding sensitivity. (You can omit channel B-this the response from a load cell w.r.t input excitation).
%%Read file
file_1 = "load_file"
N1=file_1.RequestedLength;
time_interval=file_1.Tinterval;
sample_time=N1.*time_interval; %Frame size in seconds
Fs=(1/time_interval); %sampling rate
dt = (0:1:N1-1).*time_interval;
N=(N1); %nfft size
f_A=Fs*(0:N/2-1)./N;
Ny_f=Fs/2; %Nyguist frequency-Max Bandwidth
SL=N/2; %Spectral lines
res=Ny_f/SL;
sample_bins=1:1:N;
gain=81.82;
%loading channels and considering gain and sensitivity factors.
A1=detrend(file_1.A(1:N1).*81.82;); %piezo-after amplification - Detrending due to added offset to capyure data
B1=(file_1.B(1:N1).*1000./1.11); %load cell
C1=(file_1.C(1:N1).*1000./10.31); %accelerometer-D1
D1=(file_1.D(1:N1).*1000./10.37); %accelerometer-D2
Then I calculated the FFT for each channel.
%% FFT - Frquency Domain
% A-Piezo-After amplification
fft_S_A=fft(A1);
fft_S_A_onesided=fft_S_A(1:N/2); %Nyquist frequency-eliminating mirror part
S_A_meg=abs(fft_S_A_onesided);
% B-load cell
fft_S_B=fft(B1);
fft_S_B_onesided=fft_S_B(1:N/2); %Nyquist frequency-eliminating mirror part
S_B_meg=abs(fft_S_B_onesided);
% C-Accelerometer-D1
fft_S_C=fft(C1);
fft_S_C_onesided=fft_S_C(1:N/2); %Nyquist frequency-eliminating mirror part
S_C_meg=abs(fft_S_C_onesided);
% D-Accelerometer-D2
fft_S_D=fft(D1);
fft_S_D_onesided=fft_S_D(1:N/2); %Nyquist frequency-eliminating mirror part
S_D_meg=abs(fft_S_D_onesided);
Then autospectrum and cross spectrum are calculated, followed by transfer function for each individual acceleroemter with the input.
%% Manual FRF Calculation %%
%auto spectrum
Sxx=fft_S_A.*conj(fft_S_A); %Input-force/load
Syy1=fft_S_C.*conj(fft_S_C); %D1
Syy2=fft_S_D.*conj(fft_S_D); %D2
%Cross spectrum
Sxy1=conj(fft_S_A).*fft_S_C; %Input-receiver-Conj of i/p to o/p
Sxy2=conj(fft_S_A).*fft_S_D; %Input-receiver-Conj of i/p to o/p
Sxy3=conj(fft_S_C).*fft_S_D;
%TF-A,C (D1)
TF1=Sxy1./Sxx; % Transfer Function
TFreal1=real(TF1);
TFimag1=imag(TF1);
TFmag1=abs(TF1);
TFphase1=atan2d(imag(TF1),real(TF1)); %phase in degrees
COH1=abs(Sxy1.^2)./(Sxx.*Syy1);
[Cxy1,F1Coh] = mscohere(A1(1:N1),C1(1:N1),[],[],[],Fs,"onesided") ;
%TF-A,D (D2)
TF2=Sxy2./Sxx; % Transfer Function
TFreal2=real(TF2);
TFimag2=imag(TF2);
TFmag2=abs(TF2);
TFphase2=atan2d(imag(TF2),real(TF2)); %phase in degrees
COH2=((abs(Sxy2)).^2)./(Sxx.*Syy2);
[Cxy2,F2Coh] = mscohere(A1(1:N1),D1(1:N1),[],[],[],Fs,"onesided") ;
The FFT shows the correct excitation frequency in all the channels as expected becaused in forced vibration, we need to get the same harmonic response as input. But when i calculate the transfer function, it is heavily domainated by noise and hence the value are changing like anything. when i divide the output over input i must see the peak at the excitation frequency right? because others have lesser value than the highest peak. But noise is making the data to be random.
Can anyone please shed some light. I really need this help. Thanks. I will attache the code and plots. The input file are in this drive link Matlab file for input data . I put it in the drive because the file are big and unable to upload here. Both file are from the same test. When you plot FFT, you can see they are same. But when we plot FRF values are changing due to random noise. My notion is that for the same test, even if we collect multiple waveforms, the FRF should be same. But due to noise, FRF's are changing. Can anyone help?
0 个评论
采纳的回答
Mathieu NOE
2023-3-31
hello
why are you using a single frequency tone for measuring FRF's ?
normally we want to excite all frequencies (otherwie you notice that the coherence will drop at frequencies were either the input or the output is weak , or the output is too much contaminated with uncorrelated noise)
that a topic that modal analysis people do know
here a "good" example with appropriate wideband excitation signal
you notic that the coherence is close to one where the FRF has enough gain so that output noise is not an issue
if you would use that example with only a bandlimited excitation signal, you would have low coherence outside that freq range and your FRF estimation would be garbagge
NB that this code is what tfe and tfestimate and alike are doing in your Signal Processing Tbx
the data file is attached
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1341409/image.png)
data = load('beam_experiment.mat');
x = transpose(data.x); %input
y = transpose(data.y); %output
fs = data.fs; % sampling frequency
NFFT = 2048;
NOVERLAP = round(0.75*NFFT); % 75 percent overlap
%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
%% alternative with supplied sub function
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,hanning(NFFT),NOVERLAP);
% Txy = transfer function (complex), Cxy = coherence, F = freq vector
% Bode plots
figure(1),
subplot(3,1,1),plot(F,20*log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3,1,2),plot(F,180/pi*(angle(Txy)));
ylabel('Phase (°)');
subplot(3,1,3),plot(F,Cxy);
xlabel('Frequency (Hz)');
ylabel('Coherence');
%%% damping ratioes for modes
N = 2 ; % number of (dominant) modes to identify
[fn,dr] = modalfit(Txy,F,fs,N,'FitMethod','pp');
T = table((1:N)',fn,dr,'VariableNames',{'Mode','Frequency','Damping'})
%%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function and Coherence Estimate
% compute PSD and CSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x , y if length is less than the window length
x(nwind)=0;
y(nwind)=0;
n=nwind;
end
x = x(:); % Make sure x is a column vector
y = y(:); % Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);
Pyy = zeros(nfft,1);
Pxy = zeros(nfft,1);
for i=1:k
xw = window.*x(index);
yw = window.*y(index);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft);
Xx2 = abs(Xx).^2;
Yy2 = abs(Yy).^2;
Xy2 = Yy.*conj(Xx);
Pxx = Pxx + Xx2;
Pyy = Pyy + Yy2;
Pxy = Pxy + Xy2;
end
% 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);
Pyy = Pyy(select);
Pxy = Pxy(select);
else
select = 1:nfft;
end
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end
10 个评论
Mathieu NOE
2024-4-16
ok, let me explain it better
Txy is computed with a so called H estimator - see on internet the different H1,H2, Hv estimators used in experimental analysis. Those different estimators are used depending if there is noise contamination either on the input (x) , or on the output (y) or on both sides
most of the time we assume there is only noise (whatever the amount) on the output, so people tend to use the H1 estimator all the time
now, in a ideal world with noise free signals , a FRF can be obtained directly with the fft of x and y
X = fft(x) (complex vector)
Y = fft(y) (complex vector)
noise free case : H = Y/X , so it's a ratio of two frequency vectors , and to convert it in dB for a typical Bode plot, you do dB_gain = 20*log10(abs(H))
now, the "real life" situation, as explained above, is not " noise free signals " all the way, but you have to pick one the H estimator to reduce the impact of noise, and if as in many situations you may have the noise mostly on the measurement (output) channel , the you pick the H1 estimator , which is what is implemented in the matlab code above : Txy = Pxy ./ Pxx;
Pxy = Y*X' (crosspower of x and y) (X' is the conjugate of X)
Pxx = X*X' (autopower of x)
so Txy = (Y*X') / (X*X') - this is the H1 estimator and if there was no noise on Y, this could be simplified into Txy = Y / X
in short, whatever H estimator you pick, dB_gain = 20*log10(abs(H))
hope this helps
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!