Shift Filter outputs to align peaks

I am trying to add a time delay to each of my filters so that they are all overlayed on top of each other. They aren't all being shifted by the exact same amount as I want the center point (peak) of them all to align. I am unsure how I would either calculate this value or if there is a function which does this for you. I have found fftshift but no other shifting function. I am referring to figure 4! Thank you for your time!

clc

clearvars

close all

fs = 20e3;

numFilts = 32; %

filter_number = numFilts;

order = 4;

CenterFreqs = logspace(log10(50), log10(8000), numFilts);

figure

plot(CenterFreqs)

title('Center Frequencies')

% input signal definition

signal_freq = 100; % Hz

Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)

t = linspace(0,Nperiods/signal_freq,200*Nperiods); %

end

title('Impulse Response');

hold off

xlim([0 500]);

for ii = 1:numel(indf)

output(ii,:) = filter(IR_array{indf(ii)},1,input);

plot(t,output(ii,:))

LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each

end

output_sum = sum(output,1);

plot(t,output_sum)

LEGs{ii+1} = ['Sum'];

legend(LEGs{:})

legend('Show')

title('Filter output signals');

xlabel('Time (s)');

ylabel('Amplitude');

hold off

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

function IR = gammatone(order, centerFrequency, fs)

% Design a gammatone filter

earQ = 9.26449;

minBW = 24.7;

% erb = (centerFrequency / earQ) + minBW;

erb = ((centerFrequency/earQ)^order + (minBW)^order)^(1/order);% we use the generalized

[S]

@Mathieu NOE Tagging you in case!:)

[Mathieu NOE]

hello again @S

I am unsure why we have to do this - see my comment to the answer provided below

Tx

Answer 1

[Mathieu NOE]

hello @Chunru

sorry but I was a bit skeptical about your solution based on group delay

seems to me the "aligned" output signals are not that aligned ...

first I thought that maybe the group delay should be computed at the signal frequency (100 Hz) and not at the CF , but even with that mod I coudn't get the desired aligned plot

then also the sign of the time shift is incorrect as we need to add the delay to t and not retrieve it (the output of the grpdelay function is a positive number of samples, corresponding to a negative phase shift)

finally , I simply used the filter's phase value , interpolated at the signal frequency and I think now we can say the output signals are perfectly aligned

I simply removed the output sum from this time plot as I think it doesn't make sense here

fs = 20e3; 
numFilts = 32; % 
filter_number = numFilts;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
figure
plot(CenterFreqs)
title('Center Frequencies')

% input signal definition 
signal_freq = 100; % Hz
Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
t = linspace(0,Nperiods/signal_freq,200*Nperiods); % 
input = sin(2*pi*signal_freq*t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    IR = gammatone(order, CenterFreqs(ii), fs);
    [tmp,~] = freqz(IR,1,1024*2,fs);
    % scale the IR amplitude so that the max modulus is 0 dB
    a = 1/max(abs(tmp));
    %     % or if you prefer - 3dB
    %     g = 10^(-3 / 20); % 3 dB down from peak
    %     a = g/max(abs(tmp));
    IR_array{ii} = IR*a; % scale IR and store in cell array afterwards
    [h{ii},f] = freqz(IR_array{ii},1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(IR_array{ii})
end
title('Impulse Response');
hold off
xlim([0 500]);

figure
hold on
for ii = 1:filter_number
    plot(f,20*log10(abs(h{ii})))
end
title('Bode plot');
set(gca,'XScale','log');
xlabel('Freq(Hz)');
ylabel('Gain (dB)');
ylim([-100 6]);

figure
hold on
indf = find(CenterFreqs<signal_freq); % define up to which filter we need to work
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);
    %plot(t ,output(ii,:))         
    phase_cor(ii) = interp1(f,angle(h{ii}),signal_freq);
    t_shift(ii) = phase_cor(ii)/(2*pi*signal_freq);  % align
    plot(t + t_shift(ii),output(ii,:))            % align
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
output_sum = sum(output,1);
% plot(t,output_sum)
% LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals');
xlabel('Time (s)');
ylabel('Amplitude');
hold off

%% perform FFT of signal : 
[freq,fft_spectrum] = do_fft(t,output_sum);

figure
plot(freq,fft_spectrum,'-*')
xlim([0 1000]);
title('FFT of output sum signal')
ylabel('Amplitude');
xlabel('Frequency [Hz]')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)

time = time(:);
data = data(:);

dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length


%% use windowing or not at your conveniance
% no window
fft_spectrum = abs(fft(data))*2/nfft;

% % hanning window
% window = hanning(nfft);
% window = window(:);
% fft_spectrum = abs(fft(data.*window))*4/nfft;


% 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IR = gammatone(order, centerFrequency, fs)
% Design a gammatone filter
earQ = 9.26449;
minBW = 24.7;
%       erb = (centerFrequency / earQ) + minBW; 
erb = ((centerFrequency/earQ)^order + (minBW)^order)^(1/order);% we use the generalized
%       function (depending on order)

%       B = 1.019 .* 2 .* pi .* erb; % no, the 2*pi factor is implemented twice in your code
B = 1.019 * erb;
%       g = 10^(-3 / 20); % 3 dB down from peak % what is it for ? see main code above
f = centerFrequency;
tau = 1. / (2. .* pi .* B);
% gammatone filter IR
% t = (0:1/fs:10/f);
tmax = tau + 20/(2*pi*B);
t = (0:1/fs:tmax);

temp = (t - tau) > 0;
%       IR = (t.^(order - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2*pi*f*(t - tau)) .* temp;
IR = ((t - tau).^(order - 1)) .* exp(-2*pi*B*(t - tau)).* cos(2*pi*f*(t - tau)) .* temp;
end

Answer 2

[Chunru]

You can compute the group delay around the fc to align the output (approximately).

fs = 20e3; 
numFilts = 32; % 
filter_number = numFilts;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
figure
plot(CenterFreqs)
title('Center Frequencies')

% input signal definition 
signal_freq = 100; % Hz
Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
t = linspace(0,Nperiods/signal_freq,200*Nperiods); % 
input = sin(2*pi*signal_freq*t) + 0*rand(size(t));
figure
hold on
for ii = 1:filter_number
    IR = gammatone(order, CenterFreqs(ii), fs);
    
    % Get the group delay around center freq
    [gd,f] = grpdelay(IR, 1, 1024, fs);
    [~, ind] = min(abs(f-CenterFreqs(ii)));
    
    GroupDelay(ii) = gd(ind);   % in samples
    
    [tmp,~] = freqz(IR,1,1024*2,fs);
    % scale the IR amplitude so that the max modulus is 0 dB
    a = 1/max(abs(tmp));
    %     % or if you prefer - 3dB
    %     g = 10^(-3 / 20); % 3 dB down from peak
    %     a = g/max(abs(tmp));
    IR_array{ii} = IR*a; % scale IR and store in cell array afterwards
    [h{ii},f] = freqz(IR_array{ii},1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(IR_array{ii})
end

GroupDelay

GroupDelay = 1x32

1.0e+03 * 0.2382 0.2069 0.1890 0.1776 0.1402 0.1360 -0.9253 1.0474 1.0613 0.5075 0.4554 0.3999 0.3458 0.2968 0.2534 0.2158 0.1833 0.1558 0.1325 0.1124 0.0954 0.0809 0.0689 0.0583 0.0495 0.0419 0.0356 0.0305 0.0255 0.0217

title('Impulse Response');
hold off
xlim([0 500]);

figure
hold on
for ii = 1:filter_number
    plot(f,20*log10(abs(h{ii})))
end
title('Bode plot');
set(gca,'XScale','log');
xlabel('Freq(Hz)');
ylabel('Gain (dB)');
ylim([-100 6]);

figure
hold on
indf = find(CenterFreqs<signal_freq); % define up to which filter we need to work
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);
    plot(t ,output(ii,:))            
    %plot(t - GroupDelay(ii)/fs,output(ii,:))            % align
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
output_sum = sum(output,1);
plot(t,output_sum)
LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals');
xlabel('Time (s)');
ylabel('Amplitude');
hold off

figure
hold on
indf = find(CenterFreqs<signal_freq); % define up to which filter we need to work
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);
    %plot(t ,output(ii,:))            
    plot(t - GroupDelay(ii)/fs,output(ii,:))            % align
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
output_sum = sum(output,1);
plot(t,output_sum)
LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals (Aligned)');
xlabel('Time (s)');
ylabel('Amplitude');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IR = gammatone(order, centerFrequency, fs)
% Design a gammatone filter
earQ = 9.26449;
minBW = 24.7;
%       erb = (centerFrequency / earQ) + minBW; 
erb = ((centerFrequency/earQ)^order + (minBW)^order)^(1/order);% we use the generalized
%       function (depending on order)
%       B = 1.019 .* 2 .* pi .* erb; % no, the 3pi factor is implemented twice in your code
B = 1.019 * erb;
%       g = 10^(-3 / 20); % 3 dB down from peak % what is it for ? see main code above
f = centerFrequency;
tau = 1. / (2. .* pi .* B);
% gammatone filter IR
t = (0:1/fs:10/f);
temp = (t - tau) > 0;
%       IR = (t.^(order - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2*pi*f*(t - tau)) .* temp;
IR = ((t - tau).^(order - 1)) .* exp(-2*pi*B*(t - tau)).* cos(2*pi*f*(t - tau)) .* temp;
end

[S]

This looks great! Thank you! Could you explain what is happening here:

IR = gammatone(order, CenterFreqs(ii), fs);

% Get the group delay around center freq

[gd,f] = grpdelay(IR, 1, 1024, fs);

[~, ind] = min(abs(f-CenterFreqs(ii)));

GroupDelay(ii) = gd(ind); % in samples

Why are we taking the minimum of the absolute value? @Chunru

[Chunru]

grpdelay compute the Group delay as a function of frequency where is the phase response of the filter. So for different frequency, there is different group delay. For gamma tone filters, each filter can be approximated as a narrow band filter and the group delay can be considered as a constant approximately. From the group delay vs freq, we choose the group delay at CenterFreqs of the band. "min(abs(..))" is used to choose the closest freq to CenterFreqs.

[Mathieu NOE]

interesting , but I wonder why we have to do this.

I may not know everything about gammatone filters but this time alignment is something I haven't seen in publications - or I am unaware of something ?

also the "tau" time shift in the gammatone filters function is another point I wonder why it exist

tau = 1. / (2. .* pi .* B);
% gammatone filter IR
t = (0:1/fs:10/f);
temp = (t - tau) > 0;
IR = ((t - tau).^(order - 1)) .* exp(-2*pi*B*(t - tau)).* cos(2*pi*f*(t - tau)) .* temp;

this equation of the filter - taken from the attached pdf , dos not use any tau time shift in the IR computation.

anyway, this "tau" shift is then compensated by the newly introduced group delay time shift , so I wonder what we are doing here by shifting signals all over the place...

[Chunru]

Usually, we are not interested in the delays of individual filters over the filter bank. OP want to plot the filter bank output and try to align the peaks. Since different filter has different group delay, therefore the delay should be compensated to show a roughly aligned filter bank output.

[Chunru]

Hi @Mathieu NOE.

sorry but I was a bit skeptical about your solution based on group delay

seems to me the "aligned" output signals are not that aligned ...

==> Yes. It is an approximation and not ideally aligned. Since the filters are not linear phase filters (which has a constant delay over freq), signal waveform will not ideally aligned. The approximation works on the assumption that each subband has a "narrow" band frequency and the group delay can be approximated at the centre frequency (though we are aware that the group delay is not a constant over frequency in each band).

first I thought that maybe the group delay should be computed at the signal frequency (100 Hz) and not at the CF , but even with that mod I coudn't get the desired aligned plot

==> As explained above, we would like to find difference of the time delay of each subband. Each subband filter can be regarded as narrow band (approximately). The output waveform of each subband should be concentrated around CF. Therefore, group delay should aprroximately computed at CF.

then also the sign of the time shift is incorrect as we need to add the delay to t and not retrieve it (the output of the grpdelay function is a positive number of samples, corresponding to a negative phase shift)

==> I don't delay the signal. I just simple plot time vs signal by bring forward the time axis:

plot(t - GroupDelay(ii)/fs,output(ii,:))

finally , I simply used the filter's phase value , interpolated at the signal frequency and I think now we can say the output signals are perfectly aligned

I simply removed the output sum from this time plot as I think it doesn't make sense here

==> I think you are using phase delay instead of group delay. As approximation at CF, both can be valid estimates.

Group delay:

Phase delay:

Notice the -ve sign in above equations. I think you are not using time delay; instead you are using time advance. Therefore the sign change in your code.

t_shift(ii) = phase_cor(ii)/(2*pi*signal_freq);  % align
plot(t + t_shift(ii),output(ii,:))            % align

Hope this helps.

[S]

@Mathieu NOE Thank you! This looks great, I am doing this just for visual purposes to be able to show the differences between each filter output better so you are correct in that the sum is no longer needed. Just so I make sure I understand in this part:

phase_cor(ii) = interp1(f,angle(h{ii}),signal_freq);

t_shift(ii) = phase_cor(ii)/(2*pi*signal_freq); % align

plot(t + t_shift(ii),output(ii,:)) % align

interp1 is using interpolation to use the frequency and signal frequency to get the phase, and then we are using that phase to find the time shift?

Also it seems I can't thumbs up your answer! Can you submit it as an answer

[Mathieu NOE]

hello again

ok, I changed from comment to answer (tx in advance)

and yes , basically with angle(h{ii}) we have the filter phase response vs frequency f

we simply have to get the angle value corresponding to the signal frequency and that's why we use interp1 for this task

phase_cor(ii) = interp1(f,angle(h{ii}),signal_freq);

to convert a phase into a time delay it's simply a 2*pi*frequency factor (using the signal frequency here again)

t_shift(ii) = phase_cor(ii)/(2*pi*signal_freq);

[S]

Oh got it!! Thank you! @Mathieu NOE

So just looking through, and I was a bit confused on this part:

function [freq_vector,fft_spectrum] = do_fft(time,data)

time = time(:);

data = data(:);

dt = mean(diff(time));

Fs = 1/dt;

nfft = length(data); % maximise freq resolution => nfft equals signal length

I don't see where do_fft is used and why are we using data and not output? Also if I were to add zero padding to the segment right below:

fft_spectrum = abs(fft(data))*2/nfft;

Would I do that in the above line or this line?

fft_spectrum = abs(fft(data.*window))*4/nfft;

I was thinking of doing:

lpad = 2*length(output);

freq = 0:Fs/lpad:Fs/2;

plot(freq,output)

[Mathieu NOE]

hello again

from the code I posted above in my first answer, the last figure is about the fft spectrum of the sum of the filter's output

%% perform FFT of signal : 
[freq,fft_spectrum] = do_fft(t,output_sum);
figure
plot(freq,fft_spectrum,'-*')
xlim([0 1000]);
title('FFT of output sum signal')
ylabel('Amplitude');
xlabel('Frequency [Hz]')

here again the full code :

fs = 20e3; 
numFilts = 32; % 
filter_number = numFilts;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
figure
plot(CenterFreqs)
title('Center Frequencies')

% input signal definition 
signal_freq = 100; % Hz
Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
t = linspace(0,Nperiods/signal_freq,200*Nperiods); % 
input = sin(2*pi*signal_freq*t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    IR = gammatone(order, CenterFreqs(ii), fs);
    [tmp,~] = freqz(IR,1,1024*2,fs);
    % scale the IR amplitude so that the max modulus is 0 dB
    a = 1/max(abs(tmp));
    %     % or if you prefer - 3dB
    %     g = 10^(-3 / 20); % 3 dB down from peak
    %     a = g/max(abs(tmp));
    IR_array{ii} = IR*a; % scale IR and store in cell array afterwards
    [h{ii},f] = freqz(IR_array{ii},1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(IR_array{ii})
end
title('Impulse Response');
hold off
xlim([0 500]);

figure
hold on
for ii = 1:filter_number
    plot(f,20*log10(abs(h{ii})))
end
title('Bode plot');
set(gca,'XScale','log');
xlabel('Freq(Hz)');
ylabel('Gain (dB)');
ylim([-100 6]);

figure
hold on
indf = find(CenterFreqs<signal_freq); % define up to which filter we need to work
for ii = 1:numel(indf)
    output(ii,:) = filter(IR_array{indf(ii)},1,input);
    %plot(t ,output(ii,:))         
    phase_cor(ii) = interp1(f,angle(h{ii}),signal_freq);
    t_shift(ii) = phase_cor(ii)/(2*pi*signal_freq);  % align
    plot(t + t_shift(ii),output(ii,:))            % align
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end
output_sum = sum(output,1);
% plot(t,output_sum)
% LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
title('Filter output signals');
xlabel('Time (s)');
ylabel('Amplitude');
hold off

%% perform FFT of signal : 
[freq,fft_spectrum] = do_fft(t,output_sum);
figure
plot(freq,fft_spectrum,'-*')
xlim([0 1000]);
title('FFT of output sum signal')
ylabel('Amplitude');
xlabel('Frequency [Hz]')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)
time = time(:);
data = data(:);
dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length
%% use windowing or not at your conveniance
% no window
fft_spectrum = abs(fft(data))*2/nfft;
% % hanning window
% window = hanning(nfft);
% window = window(:);
% fft_spectrum = abs(fft(data.*window))*4/nfft;
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IR = gammatone(order, centerFrequency, fs)
% Design a gammatone filter
earQ = 9.26449;
minBW = 24.7;
%       erb = (centerFrequency / earQ) + minBW; 
erb = ((centerFrequency/earQ)^order + (minBW)^order)^(1/order);% we use the generalized
%       function (depending on order)
%       B = 1.019 .* 2 .* pi .* erb; % no, the 2*pi factor is implemented twice in your code
B = 1.019 * erb;
%       g = 10^(-3 / 20); % 3 dB down from peak % what is it for ? see main code above
f = centerFrequency;
tau = 1. / (2. .* pi .* B);
% gammatone filter IR
% t = (0:1/fs:10/f);
tmax = tau + 20/(2*pi*B);
t = (0:1/fs:tmax);
temp = (t - tau) > 0;
%       IR = (t.^(order - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2*pi*f*(t - tau)) .* temp;
IR = ((t - tau).^(order - 1)) .* exp(-2*pi*B*(t - tau)).* cos(2*pi*f*(t - tau)) .* temp;
end

[Mathieu NOE]

forgot to ask why the need for zero padding ?

do you want to increase the frequency resolution of the fft output ?

if this is what you want , it's simpler and more effective to increase the time duration of you signals

here you can increase again the number of periods :

Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)

[S]

Oh ok! So this way by increasing the number of periods we are using we will get more accurate results instead of adding zeros to the start and finish, is the idea right?

@Mathieu NOE

[Mathieu NOE]

yes, the fft frequency resolution get's better as your record length increases

df = nfft / Fs , so for a given sampling rate Fs, if you make a single fft computation (and not an averaged fft) , so nfft = nb of samples of your signal; a longer signal means a larger nfft value, so you see that df (your fft frequency resolution ) will be reduced (a better resolution)

you can see by yourself by choising Nperiods = 10, then 20, .. or higher values

[S]

That makes sense! Thank you! @Mathieu NOE