Half Power Method for Vibroacoustic Dataset
7 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Daniel Pauli
el 13 de Oct. de 2021
Comentada: Mathieu NOE
el 2 de Nov. de 2021
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 comentario
Star Strider
el 13 de Oct. de 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.
.
Respuesta aceptada
Daniel Pauli
el 18 de Oct. de 2021
5 comentarios
Mathieu NOE
el 2 de Nov. de 2021
hello Daniel
glad it's fine now for you
may I ask if you would accept my answer ?
tx
Más respuestas (7)
Mathieu NOE
el 13 de Oct. de 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 (°)');
0 comentarios
Daniel Pauli
el 14 de Oct. de 2021
1 comentario
Mathieu NOE
el 14 de Oct. de 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
Daniel Pauli
el 14 de Oct. de 2021
2 comentarios
Mathieu NOE
el 14 de Oct. de 2021
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
el 14 de Oct. de 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.
Daniel Pauli
el 14 de Oct. de 2021
1 comentario
Mathieu NOE
el 14 de Oct. de 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
Daniel Pauli
el 15 de Oct. de 2021
1 comentario
Mathieu NOE
el 15 de Oct. de 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
Daniel Pauli
el 15 de Oct. de 2021
1 comentario
Mathieu NOE
el 15 de Oct. de 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 !
Daniel Pauli
el 18 de Oct. de 2021
1 comentario
Mathieu NOE
el 18 de Oct. de 2021
hello Daniel
ok - thanks for the info
I was unsure if you still needed my support
all the best
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!