I would like to calculate the frequency response function (FRF) for my input frequency and output response from accelerometer

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?

Answer 1

[Mathieu NOE]

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

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

[Yaser Arafath Gulam Dhasthagir]

@Mathieu NOE Thanks for the reply. Yes you are correct about modal analysis. I already did modal analysis and found the natural frequencies. I'm using the identified natural frequencies to compute FRF and compare between different cases. Thats why I'm using single frequency here. So when I do that, when computing TF the noise is interfering with the data because different captured files of the same test gives different FRF because of the uncorrelated noise. Hence I was wondering, how can i eliminate this noise to get consistent results with all the captured files of the same test.

[Mathieu NOE]

I am still not understanding what the single frequency computation is for ?

[Yaser Arafath Gulam Dhasthagir]

@Mathieu NOE its just for my research purpose. Trying out different solutions.

[Mathieu NOE]

ok, why not,simply I don't see the reason for it

never saw a FRF / modal analysis done one a single frequency ?

[Mathieu NOE]

Hello

Problem solved ?

would you mind accepting my answer ? thanks !

[Yaser Arafath Gulam Dhasthagir]

Hi, actually it did not but it is a good code for FRF. My aim is to see the change in amplitude for a single frequency using the similar transfer function method.

[Ben Halkon]

Nice @Mathieu NOE. But, since Txy is power, I believe plot(F,20*log10(abs(Txy))) should be 10log10 instead?

[Mathieu NOE]

@Ben Halkon

Txy is a ratio of powers, but not a power itself

Txy = Pxy ./ Pxx; % transfer function estimate

[Ben Halkon]

x^2/y^2 = (x/y)^2, not x/y.

[Mathieu NOE]

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

cross correlation - Real Data Complex Transfer Function using H0, H1, H2 Estimators - Signal Processing Stack Exchange

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