filtfilt problem when trying to compensate for delay

For a school project i have to analyse AE data with matlab.

using the filtfilt function i want to zero phase the IIR filteroutput, but i get an error and i cant see why.

data = readtable('coupon 18 peak 1');    %load signal
signal = data.Value;                     %

Fs = 50000;                             % samplefreq
sig = signal * 100 ;                    
dt = 1/Fs;                              
[samples,channels] = size(sig);         

t = (0:dt:(samples-1)*dt);              
t = t';                                 

df = Fs/samples;
Hd = HdIIR;


%%%%%% Filter add %%%%%%


FilteredSignal = filter(Hd,sig);


FilteredSignalD = filtfilt(Hd,sig);  
%%% where sig is my signal and Hd is the IIR filter constructed by the filterdesigntool.

figure(2)
plot(t,sig,t,FilteredSignalD);
xlabel('tijd (s)')
ylabel('amplitude(mg)')
title('standaard signaal met gefilterd signaal zonder delay')

this is the part of the script where i try to filter my signal. I am using a IIR because for a FIR i get a minimum order above 400 (btw i managed to correct the delay for the FIR with grpdelay)

Answer 1

[Mathieu NOE]

hello

the Hd filter was not supplied so I suggest here some standard filter (butterworth lowpass or bandpass) or smoothing techniques using smoothdata (which acts as a zero phase lag lowpass filter)

you can easily change the code to fit your needs

I supposed that the reason for filtering was to smooth out the data , but this can be in contradiction to keep the transient trace accurate as possible

data = readtable('coupon 18 peak 1');    %load signal
signal = data.Value;                     %
Fs = 50000;                             % samplefreq
sig = signal * 100 ;                    
dt = 1/Fs;                              
[samples,channels] = size(sig);         
t = (0:dt:(samples-1)*dt);              
t = t';                                 
df = Fs/samples;
% Hd = HdIIR;
%%%%%% Filter add %%%%%%
f_low = 100;			
f_high = 5000;
N = 2;
% [B,A] = butter(N,[f_low f_high]*2/Fs, 'bandpass');	% bandpass
[B,A] = butter(N,f_high*2/Fs);	% lowpass

FilteredSignalD = filtfilt(B,A,sig);  

SmoothedSignalD = smoothdata(sig,'gaussian',25);

%%% where sig is my signal and Hd is the IIR filter constructed by the filterdesigntool.
figure(2)
plot(t,sig,t,FilteredSignalD,t,SmoothedSignalD);
xlabel('tijd (s)')
ylabel('amplitude(mg)')
title('standaard signaal met gefilterd signaal zonder delay')
legend('raw','IIR filter','smoothed');

[Peter Darrouj]

Thank you for commenting on my question. you helped me with your answer. believe it or not but the low freq is the noise i want to get rid of so i can analyse the peak and plot it in cumulative energy (i managed all of this, it is just the filter that i am kinda stuck on) :) thats why i am trying to construct the most efficient filter for me. I see u use a butterworth filter, i was thinking about an elliptic design so i get the order down even more. what syntax would you advise me to use and look into? i found out that there are alot of ways to create a filter in matlab. there is designfilt , design, fdesign.highpass

[Peter Darrouj]

think i got it. thanks alot again!

[B,A] = ellip(N,1,80,f_pass*2/Fs,"high"); %% used this line in the code

with N=8 and f_low = 600

i dont really know why the *2/Fs.

[Mathieu NOE]

yes

i was about to do the same suggestion - here a 6th order elliptical filter used in bandpass

the factor *2/Fs. is related to the Nyquist frequency (Fs/2) from which the digital filters coefficients are computed.

honestly , there is no huge difference with a N = 2 to 4 order Butterworth filter

data = readtable('coupon 18 peak 1');    %load signal
signal = data.Value;                     %
Fs = 50000;                             % samplefreq
sig = signal * 100 ;                    
dt = 1/Fs;                              
[samples,channels] = size(sig);         
t = (0:dt:(samples-1)*dt);              
t = t';                                 
df = Fs/samples;
% Hd = HdIIR;
%%%%%% Filter add %%%%%%
f_low = 3000;			
f_high = 15000;
% N = 2;
% [B,A] = butter(N,[f_low f_high]*2/Fs, 'bandpass');	% bandpass
% [B,A] = butter(N,f_high*2/Fs);	% lowpass
% [B,A] = butter(N,[f_low]*2/Fs, 'high');	% highpass
N = 6;
[B,A]=ellip(N,5,80,[f_low f_high]*2/Fs);    % Bandpass digital filter design         

FilteredSignalD = filtfilt(B,A,sig);  


%%% where sig is my signal and Hd is the IIR filter constructed by the filterdesigntool.
figure(2)
plot(t,sig,t,FilteredSignalD);
xlabel('tijd (s)')
ylabel('amplitude(mg)')
title('standaard signaal met gefilterd signaal zonder delay')
legend('raw','IIR filter');

[Mathieu NOE]

extra bonus

before designing a filter it's always good to have a look at the spectrogram of the signal to see what should be the filter's specs

here we can clearly see that the energy of the transient lies in the 3 to 15 kHz band - with the low freq noise clearly noticeable around 300 Hz

clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

data = readtable('coupon 18 peak 1');    %load signal
signal = data.Value;                     %
Fs = 50000; 

dt = 1/Fs;
[samples,channels] = size(signal);
time =  (0:samples-1)*dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.75;

% spectrogram dB scale
spectrogram_dB_scale = 60;  % 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;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;

%%%%%% legend structure %%%%%%%%
for ck = 1:channels
    leg_str{ck} = ['Channel ' num2str(ck) ];
end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % display 1 : time domain plot
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% figure(1),plot(time,signal);grid on
% title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
% xlabel('Time (s)');ylabel('Amplitude');
% legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : 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(2),plot(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),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+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Peter Darrouj]

thank you for your help. that is quite the script sir, i am not so educated in matlab so i will have to study it. i would like to understand how you created it, maybe you have a link to some in depth explanation.

could you explain what the overlap stands for. i see that it plays a role in the resolution of the spectrogram but how do u determine what to use? i dont find significant differences between 0.70 - 0.95 and how does it affect the computing time/power needed?

i believe the NFFT is for zero padding so the resolution is higher if there are not enough data points, right?

i undestand that i am asking too much in this threat so no hard feelings if you would like to stop the conversation here :).

[Mathieu NOE]

hehe

no problem at all

when you do the time / frequency analaysis with a spectrogram function , there are several factors that will play a role for frequency and time axis resolution

• frequency resolution df is driven by your sampling frequency Fs and the fft buffer length NFFT : df = Fs/NFFT. so large NFFT will increase the frequency resolution but for a gievn overlap , it will reduce the number of frequency spectra , so you loose time accuracy by the same amount
• you can increase your time accuracy if you increase your overlap, so the spectrogram looks better (smoother) with overlap compared to no overlap. but this increases the computation load as your asking for more fft spectra computation

NFFT & overlap : those are the two major parameters you can play with when you do short time fourier analysis

see for example :

37W_18839_1.pdf (tek.com)

Understanding FFT Overlap Processing Fundamentals | Tektronix