Effacer les filtres
Effacer les filtres

Shift Filter outputs to align peaks

10 vues (au cours des 30 derniers jours)
S
S le 24 Mar 2024
Commenté : S le 6 Avr 2024
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 commentaires
S
S le 24 Mar 2024
@Mathieu NOE Tagging you in case!:)
Mathieu NOE
Mathieu NOE le 25 Mar 2024
hello again @S
I am unsure why we have to do this - see my comment to the answer provided below
Tx

Connectez-vous pour commenter.

Réponse acceptée

Mathieu NOE
Mathieu NOE le 26 Mar 2024
Déplacé(e) : Mathieu NOE le 29 Mar 2024
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

Plus de réponses (1)

Chunru
Chunru le 25 Mar 2024
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 commentaires
Mathieu NOE
Mathieu NOE le 5 Avr 2024
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
S le 6 Avr 2024
That makes sense! Thank you! @Mathieu NOE

Connectez-vous pour commenter.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by