MATLAB: Obtaining displacement from acceleration with omega arithmetic

fftomega arithmetic

Hi everyone,

I have read on different posts, that another method to obtain the displacement from the acceleration is by using omega arithmetic. I have tried different techniques, however I get very confusing information / results. I am getting the accelerometer data from a tri-axis accelerometer, for which I use a sampling frequency of 400Hz. I have written the following code:

Fs = 400; %Sampling frequency from my Accelerometer
Fn = Fs/2; %Nyquist Frequency
dt = 1/Fs; 
N = numel(Time); % Amount of Samples
%First I am doing a High-Pass-Filter on my data, using the Buttord, butter, zp2sos and filtfilt functions
Wp = [5.0/400 20.0/400];    % Passband (Hz)
Ws = [0.1/400 40.0/400];    % Stopband (Hz)
Rp = 1;                     % Passband Riple (dB)
Rs = 10;                    % Stopband Riple (dB)
[n,Wn]  = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn);
[sos,g] = zp2sos(z,p,k);
AccelVek_filtered = filtfilt(sos,g,AccelVek); %AccelVek is three dimensional (Acc_x, Acc_y and Acc_z in each column for the given Time)
% Next doing the FFT
NFFT = 2^nextpow2(N);
AFFT = fft(AccelVek_filtered, NFFT)*dt; %dAccelVek_filtered is simply the acceleration filtered using filtfilt
f_fft = Fs*(0:NFFT-1)/NFFT; % Frequency
omegai = 2*pi*Fs*sqrt(-1); % Omega to use in the Frequency domain
% dividing by omega*i^2 
for i = 1:NFFT
    if f_fft(i) > 10 && f_fft(i) < Fn
        DFFT(i,1) = AFFT(i,1)/(omegai^2);
        DFFT(i,2) = AFFT(i,2)/(omegai^2);
        DFFT(i,3) = AFFT(i,3)/(omegai^2);
    else
        DFFT(i,:) = 0;
    end
    
end
DisVek = ifft(DFFT, N)*Fs; %Transforming it back using ifft

I know it is a pretty long code, but my results don't seem to make sense. Could someone please have a look and tell me if i made a mistake somewhere?

Thank you 🙂

Sam

Best Answer

Hello Samuel!

If you want to perform a single integration in the frequency domain, you have to multiply each spectral line by

1 ./ (2i * pi * f_fft),

where f_fft is a frequency of the spectral line rather than the sampling frequency.

For a full code, please refer to this submission: autoFFT

EDIT: Moreover, there is no need for low-pass filtering if you do not want to apply, e.g. frequency zoom. Note that each filter can alter magnitudes and phases of components in your signal and it in turn adds some degree of uncertainty to your experimental data. If you measure vibrations, it is usually better to perform only analog anti-aliasing filtering (your analyser should do it) and digital or analog high-pass filtering. If you do not need higher frequencies, you can discard them after the Fourier transform.

Related Solutions

MATLAB: 50 Hz interference and noise reduction from ECG.

hello

If you want only to remove the 50 Hz tone , use a notch filter as shown in example below. I used your data and guessed the sampling frequency from the fact that the peak was at 50 Hz.

So sampling at 10 kHz is pretty high for an ECG signal, you could easily go down by factor 10.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
signal = importdata('Sinal1.txt');
dt = 1e-4; % 0.1 milli seconds
samples = length(signal);
Fs = 1/dt; % sampling frequency (Hz)
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = Fs;    % 
Overlap = 0.75;
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 50;       % notch freq
wc = 2*pi*fc;
Q = 7;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)

[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid on
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
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 overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    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;
        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

MATLAB: Efficient pass filters for EEG signals

The lowpass function (and those like it) have 2 outputs, the second being a digital filter object. Save the digital filter object and use it with filtfilt to use it with other signals.

That is much more efficient than designing it each time you use it, which is what it appears you are doing now. So, design it once and use it as often as necessary afterwards.

Best Answer

Related Solutions

MATLAB: 50 Hz interference and noise reduction from ECG.

MATLAB: Efficient pass filters for EEG signals

Related Question