Half Power Method for Vibroacoustic Dataset

25 views (last 30 days)
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
  1 Comment
Star Strider
Star Strider on 13 Oct 2021
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.
.

Sign in to comment.

Accepted Answer

Daniel Pauli
Daniel Pauli on 18 Oct 2021
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
  5 Comments
Mathieu NOE
Mathieu NOE on 2 Nov 2021
hello Daniel
glad it's fine now for you
may I ask if you would accept my answer ?
tx

Sign in to comment.

More Answers (7)

Mathieu NOE
Mathieu NOE on 13 Oct 2021
hello Daniel
the frequency vector is missing
the half power method is described here :
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 (°)');

Daniel Pauli
Daniel Pauli on 14 Oct 2021
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
  1 Comment
Mathieu NOE
Mathieu NOE on 14 Oct 2021
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

Sign in to comment.


Daniel Pauli
Daniel Pauli on 14 Oct 2021
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
  2 Comments
Mathieu NOE
Mathieu NOE on 14 Oct 2021
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.

Sign in to comment.


Daniel Pauli
Daniel Pauli on 14 Oct 2021
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
  1 Comment
Mathieu NOE
Mathieu NOE on 14 Oct 2021
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

Sign in to comment.


Daniel Pauli
Daniel Pauli on 15 Oct 2021
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
  1 Comment
Mathieu NOE
Mathieu NOE on 15 Oct 2021
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

Sign in to comment.


Daniel Pauli
Daniel Pauli on 15 Oct 2021
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
  1 Comment
Mathieu NOE
Mathieu NOE on 15 Oct 2021
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 !

Sign in to comment.


Daniel Pauli
Daniel Pauli on 18 Oct 2021
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
  1 Comment
Mathieu NOE
Mathieu NOE on 18 Oct 2021
hello Daniel
ok - thanks for the info
I was unsure if you still needed my support
all the best

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by