H1 H2 Hv estimators all the same

Hi, I'm trying to analyse a response to impact of a fixed structure. I have a force signal recorded by a load cell on the hammer and an acceleration measured by the accelerometer on my structure.

I would like to compute the FRF using all the estimators H1 H2 Hv and compare them, but i get the exact same result in all cases, this seem impossible since it means that no noise has been recorded. Also i find quite strange that the coherence is always 1.

This is what I have written so far, unfortunately I can't see my mistakes. If anybody has any idea on why this is happenig it would be of great help.

Thanks in advance

clc

clear

close all

load data.mat

force = force - mean(force);

acc = (acc - mean (acc))*9.81;

% sampling frequency

fs = 6600;

N = 30000;

freq = (0:N/2)'*fs/N;

% time domain

figure;

subplot(2,1,1)

plot(t,force)

title('Force - time domain')

ylabel('F [N]')

xlabel('Time [s]')

grid minor

subplot(2,1,2)

plot(t,acc)

title('Acceleration - time domain')

ylabel('Acc. [ms^{-2}]')

xlabel('Time [s]')

grid minor

ACC = fft(acc);

FORCE = fft(force);

phase = angle(ACC./FORCE);

ACC = ACC/length(ACC);

FORCE = FORCE/length(FORCE);

Sxx = conj(FORCE).*FORCE;

Syy = conj(ACC).*ACC;

Sxy = conj(FORCE).*ACC;

Syx = conj(ACC).*FORCE;

% only the positivie part

Gxx=2*Sxx(1:length(ACC)/2+1);

Gyy=2*Syy(1:length(ACC)/2+1);

Gxy=2*Sxy(1:length(ACC)/2+1);

Gyx=2*Syx(1:length(ACC)/2+1);

ACC=2*ACC(1:length(ACC)/2+1);

FORCE=2*FORCE(1:length(FORCE)/2+1);

% estimators FRF

H1 = Gxy./Gxx;

H2 = Gyy./Gyx;

Hv =(Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx); % SHIN-HAMMOND

Coe =(abs(Gxy).^2)./(Gxx.*Gyy);

H = ACC./FORCE;

figure;

subplot(2,1,1);

plot(freq,20*log10(abs(Hv)),'linewidth',1);

xlabel('Frequency [Hz]')

ylabel('|H| (dB)')

grid on;

subplot(2,1,2);

plot(freq,phase(1:N/2+1)./pi)

% semilogy(freq,Coe,'linewidth',1); grid on;

xlabel('Frequency [Hz]')

ylabel('Phase');

grid on;

figure;

plot(freq,20*log10(abs(H1)))

hold on

plot(freq,20*log10(abs(H2)))

plot(freq,20*log10(abs(Hv)))

plot(freq,20*log10(abs(H)))

hold off

grid minor

xlabel('Frequency [Hz]')

ylabel('Amplitude [dB]')

legend('H1','H2','Hv','H')

[Paul]

Hi Riccardo,

What is the source of the equations for Sxy and Syx?

[Riccardo]

Hi Paul,

It's the cross psd you can read about it here https://community.sw.siemens.com/s/article/what-is-a-frequency-response-function-frf

Answer 1

[William Rose]

@Riccardo,

[Edit one sentence fo rclarity. No changes to code.]

When you get coherence=1, it is an indicator that you have not chopped the data into segments. You need to chop the time-domain data into segments and estimate Sxx, Sxy, Syy for each segment, then average the results from the different segments. Estimating the tranfer function with one long segment is kind of like estimating the slope of a linear relationship with just two points: the correlation will be 1.00.

Because the segments are shorter than the full length sequence, this approach produces esitmates with less frequency resolution. (The resolution in Hz, i.e. the spacing between elements of the spectrum, equals the reciprocal of the segment length, in seconds.) But you gain a measure of confidence in the spectral estimates, because the spectrum estimate at every frequency is based on multiple measurements (i.e. multiple segments).

For this data, the segmenting approach probably won't give good results, because there's only one hammer strike. Therefore most of the segments will contain mainly noise, with little to no signal, in both channels. If you could record N hammer strikes at regular intervals (long enough intervals for the response to decay to zero, or close to it) then you could adjust the segments to have one strike per segment. Then the spectrum estimates would be the average of the responses to N hammer strikes.

I'll show how you can do the segmenting and averaging with Matlab's cpsd(). But I don't expect nice results, for the reasons above.

You could also just use tfestimate(), which gives the same results as shown below.

Note: Sxy=cpsd(x,y,...)=E(X.Y*) where E() denotes expected value or mean, X=fft(x), Y=fft(y), Y*=conj(Y). Some books and articles, including Wikipedia, define Sxy=E(X*.Y), which is the complex conjugate of Matlab's Sxy, computed with cpsd(x,y). So beware.

load data.mat
force = force - mean(force);        % force=x (N)
acc = (acc - mean (acc))*9.81;      % accel=y (m/s^2)
fs = 6600;                          % sampling frequency (Hz)
Ntot=length(acc);
%% Compute spectra
nseg=29;  % number of half-overlapped segments
nfft=floor(Ntot/((nseg+1)/2));  % window length and number in each FFT
nover=floor(nfft/2);
% fprintf('Ntot=%d, nfft=%d, noverlap=%d.\n',Ntot,nfft,nover)
[Sxx,f]=cpsd(force,force,nfft,nover,nfft,fs);   % E(X.X*)
Syx=cpsd(acc,force,nfft,nover,nfft,fs);         % E(Y.X*)
Sxy=cpsd(force,acc,nfft,nover,nfft,fs);         % E(X.Y*)
Syy=cpsd(acc,acc,nfft,nover,nfft,fs);           % E(Y.Y*)
H1=Syx./Sxx;
H2=Syy./Sxy;
Coh=Syx.*Sxy./(Sxx.*Syy);

%% Plot time domain data
t=(0:Ntot-1)/fs;
figure; subplot(211); plot(t,force); grid on
title('Force'); ylabel('F (N)')
subplot(212); plot(t,acc); grid on
title('Acceleration'); ylabel('Acc. (m/s^2)'); xlabel('Time (s)')

%% Plot transfer function estimates
figure; subplot(311)
plot(f(1:nfft/2),abs(H1(1:nfft/2)),'-r',f(1:nfft/2),abs(H2(1:nfft/2)),'-b'); 
grid on; title('|H|'); legend('H1','H2')
ylabel('kg^{-1}'); xlabel('Frequency (Hz)')
subplot(312)
plot(f(1:nfft/2),angle(H1(1:nfft/2)),'-r',f(1:nfft/2),angle(H2(1:nfft/2)),'-b'); 
grid on; title('angle(H)'); legend('H1','H2')
ylabel('radians'); xlabel('Frequency (Hz)')
subplot(313)
plot(f(1:nfft/2),Coh(1:nfft/2),'g'); grid on
title('Coherence'); ylabel('\gamma^2'); xlabel('Frequency (Hz)')

The H1 and H2 estimates are so similar that H2 covers up H1 in the plots above. The coherence is very close to unity at most frequencies.

Answer 2

[Paul]

Hi Riccardo,

I'm only marginally familiar with the H1/H2 methods of estimation. Keeping that in mind ...

load data.mat
force = force - mean(force);
acc = (acc - mean (acc))*9.81;
% sampling frequency
fs = 6600;
N = 30000;
freq = (0:N/2)'*fs/N;

ACC = fft(acc);
FORCE = fft(force);
phase = angle(ACC./FORCE);
ACC = ACC/length(ACC);
FORCE = FORCE/length(FORCE);
Sxx = conj(FORCE).*FORCE;
Syy = conj(ACC).*ACC;
Sxy = conj(FORCE).*ACC;
Syx = conj(ACC).*FORCE;
% only the positivie part
Gxx=2*Sxx(1:length(ACC)/2+1);
Gyy=2*Syy(1:length(ACC)/2+1);
Gxy=2*Sxy(1:length(ACC)/2+1);
Gyx=2*Syx(1:length(ACC)/2+1);
ACC=2*ACC(1:length(ACC)/2+1);
FORCE=2*FORCE(1:length(FORCE)/2+1);

Shouldn't the way that H1 and H2 are computed (I didn't go through the Hv calculation) result in both of them being the same as H?

% estimators FRF
H1 = Gxy./Gxx; % = Sxy/Sxx = conj(FORCE).*ACC./(conj(FORCE).*FORCE) = ACC./FORCE
H2 = Gyy./Gyx; % = Syy/Syx = conj(ACC).*ACC./(conj(ACC).*FORCE) = ACC./FORCE
Hv =(Gyy-Gxx + sqrt((Gxx-Gyy).^2 + 4*abs(Gxy).^2))./(2*Gyx); % SHIN-HAMMOND
Coe =(abs(Gxy).^2)./(Gxx.*Gyy);
H = ACC./FORCE;

If we use tfestimate w/ a rectangular window with no overlap

[txy1,f] = tfestimate(force,acc,ones(N,1),N-1,[],fs,'Estimator','H1');
[txy2,f] = tfestimate(force,acc,ones(N,1),N-1,[],fs,'Estimator','H2');

then txy1 and txy2 are essentially the same

figure
plot(f,abs(txy1-txy2))

●

If we use the default window and overlap we do see a bit of difference between 1500-2000 Hz

[txy1,f] = tfestimate(force,acc,[],[],[],fs,'Estimator','H1');
[txy2,f] = tfestimate(force,acc,[],[],[],fs,'Estimator','H2');
figure
plot(f,abs(txy1-txy2))

figure
plot(f,db(abs(txy1)),f,db(abs(txy2))),grid

@William Rose has some Answers on this topic, I believe. Maybe they will help if you can find them.