Shaker Testting data - Plot sweep sinusoid in frequency domain

Hi,

I'm using data from vibration test in lab. The shaker input was a sweep sinus 2400Hz/minute, 1G (9.81m/s2) 40Hz to 30kHz. My objetive is to check using matlab script the true acceleration in all range of frequencies. Like the next picture.

Sample rate : 200kHz

Lenght of the signal : 171340508 points

Time record: 856.7 seconds

First I used fft function using 512 points and the amplitude value is almost correct. In the next try I used 20000000 points and the amplitude the amplitude decreases drastically as expected.

I think that I need to use a function, dividing the signal and defining windowing and overlap.

anyone who can help me?

[Mathieu NOE]

hello André

seems to me what you need is to extract the frequency and amplitude of the output signal of your system (the amplitude at the input is fixed at 1 G)

you can either simply compute the amplitude and frequency of a chirp signal in time domain or in frequency domain (you do a spectrogram)

IMHO the benefit of the time domain approach is that you can extract the amplitude and frequency at each cycle of your signal , whereas the fft approach would require to split the signal in smaller chunks and you will have to do compromises between frequency resolution (which is better if you consider longer buffer size) and time resolution (which is better with smaller buffer size)

Answer 1

[Mathieu NOE]

hello again

see demo below for a time domain approach

% dummy data
n = 1000;
t = (0:n-1)/n;
y = (1+ sin(pi*t/max(t))).*sin(2*pi*10*(t+1*t.^3)); % frequency and amplitude varying signal

threshold = 0; % zero crossing level

t0_pos1 = find_zc(t,y,threshold);
period = diff(t0_pos1); % delta time of crossing points
freq = 1./period; % signal frequency

% signal amplitude
for k = 1:numel(t0_pos1)-1
    ind = (t>=t0_pos1(k) & t<t0_pos1(k+1));
    amp(k) = max(abs(y(ind)));
end

t0 = t0_pos1(1:end-1) + 0.5*diff(t0_pos1); % define t index at the center of the buffer (between ZC events)

figure(1)
subplot(3,1,1),plot(t,y,'b',t0_pos1,threshold*ones(size(t0_pos1)),'.r','linewidth',2,'markersize',25);grid on
xlim([t(1) t(end)]);
title('Signal plot');
legend('signal','crossing points');
xlabel('time (s)')
ylabel('Amplitude')
subplot(3,1,2),plot(t0,freq,'b.-','linewidth',2,'markersize',12);grid on
xlim([t(1) t(end)]);
title('Signal frequency ');
xlabel('time (s)')
ylabel('Hz')
subplot(3,1,3),plot(t0,amp,'b.-','linewidth',2,'markersize',12);grid on
xlim([t(1) t(end)]);
title('Signal amplitude ');
xlabel('time (s)')
ylabel('EU')

% now we can plot the amplitude vs frequency graph
figure(2)
plot(freq,amp,'b.-','linewidth',2,'markersize',12);grid on
title('Signal amplitude vs frequency ');
xlabel('frequency (HZ)')
ylabel('Amplitude')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Zx] = find_zc(x,y,threshold)

% put data in rows
x = x(:);
y = y(:);

% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
Zx = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

[Mathieu NOE]

and this is a fft based solution (basically make a spectrogram of the signal and do this

sp_peak_hold = max(sg,[],2); % to get a "peak hold" amplitude spectrum vs frequency

to get a "peak hold" amplitude spectrum vs frequency)

code

% log sweep demo
f1 = 50;        % start freq
f2 = 250;        % stop freq
Fs = 2e3;       % sampling frequency
duration = 30;  % s

%%%%%%%%%%%
dt = 1/Fs;
samples = floor(duration*Fs)+1;

t = (0:dt:(samples-1)*dt);
beta = log10(f2-f1)/t(end);
f = f1+10.^(beta.*t);
omega = 2*pi*f;
angle_increment = omega.*dt;
angle = cumtrapz(angle_increment);  % angle is the time integral of omega.

signal = sin(angle); % constant amplitude (1) log chirp
% signal = (1+ sin(pi*t/max(t))).*sin(angle); % varying amplitude log chirp


%%%%%%%%%%%

% spectrogram demo
NFFT = 200;    % 
overlap = 0.75;
w = hamming(NFFT);

[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),overlap*NFFT);
% FFT normalisation 
sg = abs(sg)*4/NFFT;     % NB : X=fft(x.*hanning(N))*4/N; % hanning window amplitude correction

% plots sg
figure(1)
imagesc(tsg,fsg,sg);
colormap('jet')
axis('xy');
colorbar('vert');
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

% extract amplitude / frequency data
sp_peak_hold = max(sg,[],2); % to get a "peak hold" amplitude spectrum vs frequency

% now we can plot the amplitude vs frequency graph
figure(2)
plot(fsg,sp_peak_hold,'b.-','linewidth',2,'markersize',12);grid on
ylim([0 2]);
title('Signal amplitude vs frequency ');
xlabel('frequency (HZ)')
ylabel('Amplitude')

[André]

Hi Mathieu

Thanks for the Answer. Using spectogram and plot amplitude peak / frequency solves my problem.

I'm using the same script in diferent data (DATA 2) and I'm having another problem about NFFT size.

Data 2

Sample rate = 10kHz

Sweep sinus 12,22Hz/second, 1G (9.81m/s2) 20Hz to 3.3kHz

Time record: 270 seconds

In the first data (data1) the amplitude is the same regardless of NFFT size.

In the second data (data 2) the amplitude values change dependent on the size of the NFFT. In this case the results are more plausible considering NFFT=1024. There is any criteria to define the NFFT size ?

[Mathieu NOE]

hello

this is the criteria I have found correct my my own needs :

the frequency bin width (df = Fs/Nfft) must be large enough so that during the corresponding time duration of the buffer (tw = Nfft/Fs) the signal variation of frequency must be below df

if we know the sweep_rate SR (in Hz / s) this becomes :

Fs/Nfft > SR * Nfft/Fs

rearranging :

Nfft < Fs / sqrt(SR)

in your case for Data 2 this is

Nfft < 10,000 / sqrt(12.22) = 2860

so I would suggest Nmax = 2056 here

a higher overlap factor can marginaly improve the result , but it's a not a key factor here.

[André]

Hi,

I got it. Thanks for the support.

Have a nice day.