Frequency resolution using pwelch

Hi everyone,

i want to analyze the frequency domain of accelerometer data sets sampled at 500 kHz.

(figure: comparing different window-functions using pwelch with (in this case) windowlength=1024, nfft=1024, fs=500.000Hz, data length=100.000 samples)

When using the pwelch function in MATLAB i can use the windowlength (and a certain windowfunction of course) and the FFT-length as parameters affecting the frequency resolution. Looking at the figure above i'm facing two problems:

1) When changing these parameters in such way, that the frequency resolution is high i feel like it comes with high spectral leakage and therefore automatically becomes obsolete.

2) The system monitored with the accelerometer is a metal-to-metal contact, excited by a shock using a hammer. I am wondering why there is a peak at 0 Hz. In such systems i would expect high-frequency components rather than the spectrum seen in the figure.

Appreciate every help or idea on how to handle any of the two problems mentioned above!

Thanks in advance!

[Mathieu NOE]

hello

do you have the possibility to measure the hamer force signal ? i am simply trying to figureout why we would not rather use a transfer function estimate (with tfestimate) rather than just the output spectrum ?

seems to me (IMHO) a better way to characterize a system unles the input is very stable and repetitive. But I also wonderwhy your spectrum can show narrow peaks and wide peaks at the same time. So measuring the force would eliminate spurious signals (EMI) if we could average the TF estimate over multiple hits

[Oskar Kilgus]

Hi Mathieu, thanks for the reply!

I in fact have measurement data of the force but i´ve never heard of the tfestimate function and im not sure how to use it in this case. Can you help me with that? I feel like you have a clue about it.

Also I think the spectrum i attached shows narrow OR wide peaks depending on the window-function, which i tried to compare as seen in the legend.

[Mathieu NOE]

hello again

can you share your data (hammer + response time signals) may be for several hits so we can average them ?

[Oskar Kilgus]

I´ve got more than 10 Groups of Data-Sets, each with slightly different parameters to the model observerd.

In which way do you think averaging the data-sets with same parameters can contribute to a better solution? Could you explain the tfestimate approach to me?

Thanks in advance!

[Mathieu NOE]

I can make you a demo if you want to share some data

[Oskar Kilgus]

i attached a 2x2 cell array which contains:

{acceleration data of set 1, acceleration data of set 2;

force data of set 1, force data of set 2}.

Excited to see how the tfestimate works, thanks in advance!

[Mathieu NOE]

here we go

you have the choice to use tfestimate or my equivalent function mytfe_and_coh

now I have to admit the time signals gave me some headache...

I see on the force signal a min pulse folowed by some ringing , sometimes another hit appears- are you using a force sensor or are you computing the force indirectly from some acceleration signal ?

the acceleration response signal shows sometimes like delayed replications , but I assume this comes not from the input itself but is a exogeneous perturbation

here your data , limited to the first 50,000 samples (far enough)

set 1

set 2

to avoid those issues I computed the FRFs using only the first 5000 samples

the two individual FRF from each data set are quite appart in term of magnitude. This is not something I usually see when I do modal analysis with hammer testing.

so yes you can get an averaged curve, and again, my experience is that normally each hit should generate more or less the same FRF plot; averaging make simply the result better looking (smoother)

the bode plot below I limited the x axis to 100 kHz, as it gets even worse above.

Code used so far :

load('example.mat');

samples = 5000; % restrict data to first valid samples

fs = 500000;   % sampling frequency (Hz)
dt = 1/fs;
t = (0:samples-1)*dt;

fmax = 100e3; % must be < fs/2

    
% main loop

for set = 1:2
    
    x = example{2,set}; %input
    y = example{1,set}; %output
    
    % restrict to first valid samples 
    x = x(1:samples);
    y = y(1:samples);

    figure
    subplot(2,1,1),plot(t,x,'b')
    title(['Set #' num2str(set)]);
    xlabel('Time (s)')
    ylabel('Force')
    subplot(2,1,2),plot(t,y,'r')
    title(['Set #' num2str(set)]);
    xlabel('Time (s)')
    ylabel('Acceleration')

    % FFT parameters
    NFFT = samples;
    NOVERLAP = round(0*NFFT);  % 0 percent overlap

    %% solution 1 with tfestimate (requires Signal Processing Tbx)
    [Txy,F] = tfestimate(x,y,ones(NFFT,1),NOVERLAP,NFFT,fs);

    %% alternative with supplied sub function 
    % [Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,ones(NFFT,1),NOVERLAP); 
    % Txy = transfer function (complex), Cxy = coherence, F = freq vector

    % Bode plots
    ind = (F>fmax);
    F(ind) = [];
    Txy(ind) = [];
    
    figure,
    subplot(2,1,1),semilogy(F,abs(Txy));
    title(['FRF Set #' num2str(set)]);
    xlabel('Frequency (Hz)');
    ylabel('Mag');
    subplot(2,1,2),plot(F,180/pi*(angle(Txy)));
    xlabel('Frequency (Hz)');
    ylabel('Phase (°)');

    Txy_all(:,set) = Txy;
    leg{set} = ['set #' num2str(set)];

end

% averaged Txy
Txy_avg = mean(Txy_all,2);
leg{set+1} = ['average'];
figure,
subplot(2,1,1),semilogy(F,abs(Txy_all),F,abs(Txy_avg),'k');
title('FRFs');
xlabel('Frequency (Hz)');
ylabel('Mag');
legend(leg);
subplot(2,1,2),plot(F,180/pi*(angle(Txy_all)),F,180/pi*(angle(Txy_avg)),'k');
xlabel('Frequency (Hz)');
ylabel('Phase (°)');
legend(leg);




% %%% 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

[Oskar Kilgus]

first of all: thanks heaps for your effort!

I am using a force sensor directly and i guess the difference in magnitude might be related to the inconsistency of impacts because its just done with hammer and by hand.

Im curious why this way of analysing the data leads to such different results compared to the frequency spectrum in the original post above and also what the phase plot tells in this context?

[Mathieu NOE]

hello again

can you describe a bit more your physical set up ?

I am not sure to understand where is placed the force sensor in your case

when I do FRF measurement on a piece of structure , the force sensor is in the hammer and the response sensor is a accelerometer glued on the structure

[Oskar Kilgus]

its exactly like that. I do have a force sensor in the hammer and the acceleration is measured with a sensor glued on the structure.

Shouldnt be the frequency spectrum of the FRF and that of another FFT/pwelch analysis be similar?

[Mathieu NOE]

you say :

its exactly like that. I do have a force sensor in the hammer and the acceleration is measured with a sensor glued on the structure.

then I don't understand what is this record and how you obtained that data

[Oskar Kilgus]

I dont understand either. The problem is, that my task is just to analyse the data, i didnt accompany the tests themselves. I do know tho, that it should be one hit with a hammer and a force sensor to obtain force-data aswell as a acceleration sensor to obtain acceleration-data...

If we now first neglect the "artefacts" (repeated oscillations you marked in the plot above) and any confusions related to them, could you answer my questions from the last comment?

[Mathieu NOE]

your comment :

Shouldnt be the frequency spectrum of the FRF and that of another FFT/pwelch analysis be similar?

well we are comparing apples to oranges here

fft is the underlying "brick" for all frequency domain analysis of course

pwelch and alike are "single channel" fft operations that tells you about the spectral content of your signals. A hammer impact (pulse) should have a decaying fft amplitude as we move in the higher frequencies

FRF is a "dual channel" fft operations that computes cross and auto spectra from where you get the FRF estimate , but an FRF is not a spectrum as it has not relevance to a particuliar type of signal, it describes what is the magnitude / phase intrinsic characteristics of any device or structure from one input to one (or several) outputs.

FYI, there is lots of articles about modal analysis (with a hammer)

Modal Impact Postprocessing: Getting the Best FRF (siemens.com)

[Mathieu NOE]

see also

Modal Analysis of a Simulated System and a Wind Turbine Blade - MATLAB & Simulink - MathWorks France

[Mathieu NOE]

another good reference / reading :

excitation-techniques-dos-and-donts.pdf (modalshop.com)

[Oskar Kilgus]

Thank you so much Mathieur, that´ll do i suppose!

[Mathieu NOE]

my pleasure

do no hesitate to come back when you have new data to share !

Answer 1

[Gokul Nath S J]

Hi Oskar,

Changing the window length and FFT length in the pwelch function can impact the frequency resolution and spectral leakage of the resulting power spectral density (PSD) estimate. It is often a trade-off between these two factors, where increasing the window length improves the frequency resolution but also increases spectral leakage, and vice versa.

To minimize spectral leakage, it is important to choose an appropriate window function that can reduce the sidelobe levels of the FFT. Common window functions include the Hamming, Blackman, and Kaiser windows, among others. Choosing an appropriate window function and applying it properly can help you achieve a good balance between frequency resolution and spectral leakage.

Regarding the presence of a peak at 0 Hz in the spectrum, it is possible that this peak is due to the baseline, DC offset, or bias in the signal caused by some extraneous factors such as instrument noise, bias voltage, or other disturbances. This peak can also be related to the displacement or other non-linear properties of the system being monitored, especially if it is a mechanical or vibrational system.

To resolve this issue, you may want to consider subtracting the baseline or DC offset from the signal before performing the PSD estimation, or use high-pass filtering or other techniques to remove the low frequency components. You may also want to investigate the properties and behavior of the system being monitored, including its input-output relationship, the effects of boundary conditions, and other factors that may contribute to the observed low-frequency peak in the spectrum.

with regards,

Gokul Nath S J

[Oskar Kilgus]

Hey Gokul,

thanks for the reply.

I already tried different strategies to remove a DC offset, nothing worked properly to remove the 0 Hz peak. I also thought about simple high-pass or band-pass filtering the signal, but because its a mechanical system with possible low-frequency components of interest, i dont really want to filter anything below a certain frequency.

Maybe i didnt remove the DC offset the right way yet, what would you say is the best way to do so in case of accelerometer data?

Thanks in advance!

[Mathieu NOE]

you could use detrend to remove dc offset and maybe also a linear drift (then use detrend with 'linear' option)

[Oskar Kilgus]

I tried to remove dc offset through subtracting the mean value of a signal, which did not work properly. I also thought about using detrend (tried it too), but that feels rather false for a oscillating system, right?