Shift Filter outputs to align peaks

Question

S 2024-3-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098281-shift-filter-outputs-to-align-peaks

评论： S 2024-4-6

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

2 个评论
显示无隐藏无

S 2024-3-24

@Mathieu NOE Tagging you in case!:)

Mathieu NOE 2024-3-25

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 2024-3-26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098281-shift-filter-outputs-to-align-peaks#answer_1433271

移动：Mathieu NOE 2024-3-29

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Chunru 2024-3-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2098281-shift-filter-outputs-to-align-peaks#answer_1430246

在 MATLAB Online 中打开

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

13 个评论
显示 11更早的评论隐藏 11更早的评论

Chunru 2024-3-27

在 MATLAB Online 中打开

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.

Mathieu NOE 2024-4-2

在 MATLAB Online 中打开

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 2024-4-5

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 2024-4-6

That makes sense! Thank you! @Mathieu NOE

请先登录，再进行评论。

Shift Filter outputs to align peaks

2 个评论
显示无隐藏无

采纳的回答

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

13 个评论
显示 11更早的评论隐藏 11更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Shift Filter outputs to align peaks

2 个评论 显示 无隐藏 无

采纳的回答

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

更多回答（1 个）

13 个评论 显示 11更早的评论隐藏 11更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论