Half Power Method for Vibroacoustic Dataset

Hello MatLab users,

I have recorded several data sets using laser double vibrometry to determine the damping of individual resonance frequencies.

It was recommended to use the -3DB (half power rule) rule to evaluate them.

Unfortunately, I could not find an answer by searching the WWW and this forum.

I hope you can help me with my problem.

I have added the data set in the attachment. We collect this Dataset with "ITA Toolbox".

Best regards

Daniel

[Star Strider]

The file appears to be a complex Fourier transform, however there is no associated frequency vector, so it is not possible to analyse those data with respect to frequency.

.

Answer 1

[Daniel Pauli]

Hi Mathieu,

sorry to bother you again.

I think if we apply the smoothing as you suggested, too much data could be lost and points where the -3db rule could be applied also. Is it possible to reduce the smoothing a bit?

Furthermore, when applying the DB rule, only the curve that lies highest is marked. All other peaks are not taken into account.

How can we fix that?

KR

Daniel

[Mathieu NOE]

hello Daniel

I will come back to you later this week or beginning next week as I am traveliling this week for customer support

[Daniel Pauli]

Hi Mathie,

thank u very much

have a good trip.

KR

Daniel

[Mathieu NOE]

hello Daniel

i'm back !!

so , as far as I remember I had still an open point about having also the exciation channel recorded simultaneously with the response

do you have that ?

[Daniel Pauli]

Hello Mattiau,

I'm sorry I haven't been able to get back to you.

We have all the information we need.

Thank you again!

Best regards

Daniel

[Mathieu NOE]

hello Daniel

glad it's fine now for you

may I ask if you would accept my answer ?

tx

Answer 2

[Mathieu NOE]

hello Daniel

the frequency vector is missing

the half power method is described here :

Microsoft Word - 41FC527E-54FA-00270E.doc (vibrationdata.com)

also available in the attachement

I am a bit dissapointed when looking at your data

the phase plot shows so much phase rotation I wonder how the measurement was performed

normally we look for a clearly dominant (isolated) peak (showing also a phase rotation of only 180 °)

but here I don't know what to do as we have profusion of peaks (and ripples)

is that a complex spectrum or a FRF (what I supposed in first place) ?

so far what I could do, before you tell me which peak you see / want to analyse

frfc = importdata('Dataset.txt'); % complex FRF
frfc = str2double(frfc);
frf_dB = 20*log10(abs(frfc));
frf_phase = angle(frfc)*180/pi;


freq = (0:1:length(frf_dB)-1);

figure(1)
subplot(211),semilogx(freq,frf_dB,'b');
ylabel('Mag (dB)');
subplot(212),semilogx(freq,frf_phase,'b');
xlabel('Frequency');
ylabel('Phase (°)');

Answer 3

[Daniel Pauli]

Hello together,

thank you first of all for your answers.

I would like to describe our experimental set-up to you. To investigate the vibrating body, we used a laser double vibrometer, an accelerometer and a force transducer. We have analogue data sets for the other sensors. Our task is to find similar peaks of the same frequency between the different sensors with the half power / -3DB method.

Here is a plot from a messurement

In the appendix you will find data sets called Wavenumber & Time. We think (in our case) the wavenumber must be our frequenz (in kHz) abd the so called time should be the dB (Amplitude).

Thanks so lot for your help!

KR Daniel

[Mathieu NOE]

hello Daniel

so Time.txt seems to have the time response signal from the vibrometer (excitation is a constant amplitude chirp signal as I guess from a spectrogram plot)

I reconstructed the averaged fft spectrum and I have so far a smoother curve compared to the one you posted firts

I think that from there you can easily apply the - 3dB method - just follow the explanations in the publication I attached above.

all the best

you may have to correct for the sampling frequency , which I also guessed from your plots

clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal = load('Time.txt');
[samples,channels] = size(signal);

selected_channel = 1; % select your data channel here
signal = signal(:,selected_channel); 

Fs = 7500; % to be checked 
dt = 1/Fs;

time = (0:samples-1)*dt;


%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
samples = length(signal);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 256;    % fft length
OVERLAP = 0.75; % between 0 and 1 max ;  buffer overlap = OVERLAP*NFFT

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(1),plot(freq,sensor_spectrum_dB,'b');grid on
title(['Averaged FFT Spectrum  / Fs = ' num2str(0.1*round(10*Fs)) ' Hz / Delta f = ' num2str(0.1*round(10*(freq(2)-freq(1)))) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  

% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(fsg);
    sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
    my_title = ('Spectrogram (dB (A))');
else
    my_title = ('Spectrogram (dB (L))');
end

% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(0.1*round(10*Fs)) ' Hz / Delta f = ' num2str(0.1*round(10*(fsg(2)-fsg(1)))) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');




function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end


function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Answer 4

[Daniel Pauli]

Hi Mathieu,

Thanks again for your help.

If I understand correctly, you have developed a script that smoothes the signals. That helps me a lot and makes it clearer. Now I have visualised my 3 channels and have arrived at the following representation:

The two upper curves behave similarly. Is it possible to determine the matching peaks of the two curves? Then the -3dB method would have to be applied. Which commands are necessary for this? Is it possible to filter out these frequencies and apply the -3DB rule there?

Best regards

Daniel

[Mathieu NOE]

hello daniel

my code is simply doing FFT spectral averaging .

now I don't know which peak are we talking about ,; I can see several peaks but IMHO only the very first one on the left is common to the three channels . I assume you have measured your structuralresponseon 3 locations , and the first frequency is corresponding to the first eigenmode of the structure so it will be present on all channels.

can you share the 3 channels data ?

I will try to work this out

still some recommendation : it would be better not to work from only the response signals , but I would like to have a copy of the excitation signal as well.

also if you are only trageting the low frequencies , do we need the chirp signal to go so high in frequency ?

all the best

[Mathieu NOE]

to clarify : it would be better not to work from only the response signals , but I would like to have a copy of the excitation signal as well , because then we can compute a transfer function (using tfestimate) which should be cleaner than just the response spectrum

to the -3dB rule : it works well if your peak zone is clean; if your data show ripples and distorsion, this will not work very well , but you can use a curve fitting technique to get the damping factor in a more robust manner.

Answer 5

[Daniel Pauli]

Hi Mathieu,

In the appendix you will find data sets for the channels.

I have added the files twice. Once as a .mat file and as a .txt file.

Thanks for the tip with the curve adjustment.

Our main question is which commands to execute (after editing the data/curves) to apply the -3 dB method?

Ideally then for all our data.

Best regards

Daniel

[Mathieu NOE]

hello Daniel

homework finished !!!

I have implemented the method on my spectrum output. see the code and results

clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal = load('Time.txt');
[samples,channels] = size(signal);

selected_channel = 1; % select your data channel here
signal = signal(:,selected_channel); 

Fs = 7500; % to be checked 
dt = 1/Fs;

time = (0:samples-1)*dt;


%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 4;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
samples = length(signal);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 64;    % fft length
OVERLAP = 0.75; % between 0 and 1 max ;  buffer overlap = OVERLAP*NFFT

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

ind = find(freq>200 & freq <800);
freq = freq(ind);
sensor_spectrum_dB = sensor_spectrum_dB(ind);


% find both - 3dB crossing points
threshold = max(sensor_spectrum_dB) - 3; % your value here

[freq_low,~,freq_high,~]= crossing_V7(sensor_spectrum_dB,freq,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 

figure(1)
plot(freq,sensor_spectrum_dB,'b',freq_low,threshold,'+r',freq_high,threshold,'+g','linewidth',2,'markersize',12);grid on
legend('signal','freq low -3 dB crossing point','freq high -3 dB crossing point');
title(['Averaged FFT Spectrum  / Fs = ' num2str(0.1*round(10*Fs)) ' Hz / Delta f = ' num2str(0.1*round(10*(freq(2)-freq(1)))) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);

% delta f and fn computation
df = freq_high - freq_low;
% in a perfect resonating system , the central freq is the geometric average
% of the two -3 dB points
fn = sqrt(freq_high*freq_low); % magic ! we have fn = 381 Hz

% and finally : 
Q = fn/df
dzeta = 1/(2*Q)

%% tada !!!!!!!!!%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal  (example sinus amplitude 1   = 0 dB after fft).
% Linear averaging
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    s_tmp((1:samples),:) = signal;
    signal = s_tmp;
    samples = nfft;
end

% window : hanning
window = hanning(nfft);
window = window(:);

%    compute fft with overlap 
 offset = fix((1-Overlap)*nfft);
 spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
%     % for info is equivalent to : 
%     noverlap = Overlap*nfft;
%     spectnum = fix((samples-noverlap)/(nfft-noverlap));	% Number of windows

    % main loop
    fft_spectrum = 0;
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
        fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft);     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end
    fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Answer 6

[Daniel Pauli]

Good morning Mathieu,

I am very impressed with the effort you have put in for me. It works just as I imagined it would. But I still have one question. Is it also possible that it documents the values directly for all 3 channels and optionally visualises them? The command this by means of : does not work.

Greetings and many thanks

Daniel

[Mathieu NOE]

hello daniel

I will see that a bit later, I have a travel to customer to do next week and things to do

In the mean time , could it be possible to send a record of the exciation signal as well as the 3 responses ?

I would try a different approach but I need the exciation signal too (synchronous acquisition with the 3 responses signals)

tx

Answer 7

[Daniel Pauli]

Hi Mathieu,

this is no problem.

The oscillating body we measured was excited by a small hammer. The force of the hammer was measured. Do you mean that?

This Signal ist always Channel 2.

KR

Daniel

[Mathieu NOE]

hello Daniel

I am surprised

a hammer test generates a decaying pulse , not a chirp...

the data you provided so far are response to chirp excitation - this can be clearly seen in a spectrogram

example channel 1

do you have the record of the excitation signal as well ?

if you do hammer testing, I works for me too, usually I do multiple hits to average the results. In this case also, please provide the hammer signal synchronous with structure response

thanks !

Answer 8

[Daniel Pauli]

Hi Mathieu,

we know that our hammer brings a pulse with it.

It is interesting for us to see how the damping changes over time due to the impulse hammer. For this we needed the -3DB values. Thanks to you, we are now in a position to display them (many thanks again for this). Now we want to look at the totality of the measurements. To understand: we have measured a new food item every day with our technical setup and would now like to determine the differences in the damping values and correlate these with e.g. ageing processes.

Best regards

Daniel

[Mathieu NOE]

hello Daniel

ok - thanks for the info

I was unsure if you still needed my support

all the best