While loop using up to input and fixing input definition

So I have 32 filters. I am putting an input through them all in parallel. I want to be able to add all of the outputs of these filters together until the filters reach the same frequency of the input then I want to discard the rest of the results the other filters have and not add them. To do this I am a bit confused as to how to ifx what I have. I want while the CenterFreqs value is less than the frequency of the input to add so I created the while loop at the bottom. I think the way I defined input may be the issue as I want to be able to put for example sin(6000) in the input line. Thank you for your time!!

clc
clearvars
close all
fs = 20e3; 
numFilts = 32; % 
filter_number = 5;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
figure
plot(CenterFreqs)
title('Center Frequncies')
% input signal definition 
Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
t = linspace(0,2*pi*Nperiods,200*Nperiods); % for 1 period (2*pi) we have 200 samples at fs = 20e3, so signal freq = 20e3/200 = 100 Hz
input = sin(t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    IR = gammatone(order, CenterFreqs(ii), fs);
    [tmp,f] = 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
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)');
figure
hold on
for ii = 1:filter_number
    output(ii,:) = filter(IR_array{ii},1,input);
    plot(t,output(ii,:))
    LEGs{ii} = ['Filter # ' num2str(ii)]; %assign legend name to each
end
legend(LEGs{:})
legend('Show')
hold off
figure
hold on
while input>CenterFreqs
    for ii = 1:filter_number
        output(ii,:) = filter(IR_array{ii},1,input);
        plot(t,output(ii,:))
        LEGs{ii} = ['Filter # ' num2str(ii)]; %assign legend name to each
    end
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 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]

@Mathieu NOE in case you can help!!:)

[Mathieu NOE]

no problem, but it's going late here , so this is something for tomorrow (if it can wait a few hours)...

[S]

No problem! Thank you!!

[Mathieu NOE]

ok , so this is it

to make your task easier , I wanted to explicitely define the input signal frequency ,

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

after that , what we need to do is to calculate ALL impulse responses (for all 32 center frequencies), because depending now on the signal frequency, I may need 5 or more (up to 32) computations of the IR filters outputs

that' why now : filter_number = numFilts;

you can try by yourself, if now you change from 100 to 200 Hz , the number of used outputs goes from 5 to 9

I added some lines to compute the sum of all outputs , added on the same plot as the individual outputs

hope it helps !!

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); % 
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,f] = 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 = indf
    output(ii,:) = filter(IR_array{ii},1,input);
    plot(t,output(ii,:))
    LEGs{ii} = ['Filter # ' num2str(ii)]; %assign legend name to each
end
output_sum = sum(output,1);
plot(t,output_sum)
LEGs{ii+1} = ['Sum'];
legend(LEGs{:})
legend('Show')
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]

@Mathieu NOE Thank you very much! So on the last figure what are the axis?

To do the opposite I just changed this line:

indf = find(CenterFreqs<signal_freq);

to

indf = find(CenterFreqs>signal_freq);

but for some reason that has much more noise, do you know why that could be happening?

[S]

Also if I wanted to find the frequency of the output could I just do:

plot(t,fft(output_sum))

but this does not give me a value, just two spikes. How can I get it to make the title the answer

Could you post your response as answer so I can thumbs up:)

[Mathieu NOE]

hello @S and welcome back !

• the last figure x axis is time (in seconds ) , the y axis is simply the amplitude of the time signals ; see the code updated in the answer section

• I also added a new code section and function for the FFT spectrum of the sum ouput (your second comment)

To answer your 1st comment : To do the opposite I just changed this line: indf = find(CenterFreqs<signal_freq) to indf = find(CenterFreqs>signal_freq); but for some reason that has much more noise, do you know why that could be happening?

Remember that you input signal is a tone + broadband noise (white noise , so it's energy is pread accross the entire spectrum from f = 0 to Fs/2)

in the first case (indf = find(CenterFreqs<signal_freq) , this signal will go through the first 5 filters , so the higher frequency content of the broadband noise will be attenuated by the filters (as they are band pass filters by definition)

in the other case (indf = find(CenterFreqs>signal_freq)) , your signal goes through the other remaining 27 filters that have significant gain at higher frequencies (now we have selected the filters with CF from 113 Hz to 8 kHz)

in this case , the white noise content in your input signal is not attenuated at the higher frequencies , so that's why you have more noise in the sum signal.

Answer 1

[Mathieu NOE]

Hello gain, so this is the code with the small changes mentionned above

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); % 
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,:))
    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


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

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

@Mathieu NOE Thank you so much! That explanation cleared up my question :)

If I were to change the second graph to have the filters on the horizontal axis instead of time, could I just do:

plot(t,output_sum)

with filter_number instead of t? Ultimately I was trying to make it into a 3D type plot to have time on the third axis. I was looking into this and was wondering whether you think a waterplot or 'mesh' would be better for this application?

To make the time on the third axis have the best looking results I am trying to add a time shift to each filter to try and line up their results but the only shifting I see results on here is to shift everything by the same amount, which defeats the purpose. Is there a function which shifts each output individually to line them all up? (I'm not sure if they are all timeshifted by the exact same amount from the previous filter)

[Mathieu NOE]

hello again

do you mean a waterfall plot ?

I am slightly confused, when you mention the second plot , which correspond to the filter's IR response vs time (and not the output_sum signal)

so I don't understand now why you say :

plot(t,output_sum) with filter_number instead of t?

[S]

@Mathieu NOE Yes! A waterfall plot in which the z axis is now time instead of the x.

In that segment of code I was trying to make it like how in the Center Frequencies graph we just have 1-32 on the x which represents the filters. This was the step between what we have and the waterfall plot I was thinking

[Mathieu NOE]

I am not 100% sure this is what you want, but there I plotted the IR data vertically side by side for the 32 filters

Is it what you wanted ?

also see that for filters above number 10 , the time vector is not long enough so the IR is truncated - we need to define the time vector in a better way

code :

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); % 
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
    
    % horizontal plot
%     plot(IR_array{ii}) 
    
    % vertical plot (side by side)
    xdata = IR_array{ii};
    scale_factor = 0.25/max(abs(xdata));
    xdata = ii+xdata*scale_factor; 
    ydata = (0:numel(xdata)-1);
    plot(xdata,ydata); % vertical plot
end
title('Impulse Response');
xlabel('Filter #');
ylabel('IR (samples)');

hold off

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,:))            
    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);
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]

so I change inside the gammatone function the upper limit of the time vector

instade of

t = (0:1/fs:10/f)

I changed that to

tmax = tau + 20/(2*pi*B);

t = (0:1/fs:tmax);

and now you get's IR duration that are related to the B factor (exponential decay rate in the IR equation)

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); % 
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
    
    % horizontal plot
%     plot(IR_array{ii}) 
    
    % vertical plot (side by side)
    xdata = IR_array{ii};
    scale_factor = 0.25/max(abs(xdata));
    xdata = ii+xdata*scale_factor; 
    ydata = (0:numel(xdata)-1);
    plot(xdata,ydata); % vertical plot
end
title('Impulse Response');
xlabel('Filter #');
ylabel('IR (samples)');

hold off

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,:))            
    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

[S]

So what I am trying to get is a 3D plot of the data with the variables of time, filter number, and output. Is there a better way do you think than a waterfall plot? @Mathieu NOE

Additionally I am now looking at how I could try to find the ERB of each filter but the only mention of it is here:

https://www.mathworks.com/help/audio/ref/gammatonefilterbank-system-object.html

towards the end ERB is mentioned. But not how it would be actually translated from each filter to its ERB

[Mathieu NOE]

hello again

ERB : do you think we are not generating the correct ERB in our code ? I used what is described in the pdf I sent you, which I think is also what you used for creating this code.

Waterfall plot : below, one way to create a kind of waterfall plot for the filter's outputs

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
    plot3(t + t_shift(ii),ii*ones(size(t)),output(ii,:))            % align
    LEGs{ii} = ['Filter # ' num2str(indf(ii))]; %assign legend name to each
end

view(-40,60); % view angles 
% the next two lines are simply ho display integer ticks on Y axis
ax = gca;
ax.YTick = unique(round(ax.YTick));

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('Output #');
zlabel('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

[S]

@Mathieu NOE Yes, so I understand that we are using ERB already! But I mean if we were to translate these curved filters to see the erb (rectangle shapes) to be able to see the Bn. Could I use the erb function within the gammatone function to plot this, or would I have to redefine it as it whole own thing is what I mean

Also this waterfall plot looks great! I don't know if this is maybe a silly question but how come the plots aren't touching each other? Is it because we are using just specific time points and not a continuos time equation? Or is it more due to the fact that the filters are individual?

[Mathieu NOE]

hello

sorry , here I don't understand what you mean with "see" the ERB and Bn (what is Bn by the way ? you mean B = 1.019 * erb ? ) ERB and B are just numbers, so what is to see here ?

An Efficient Implementation of Gammatone Filters (shef.ac.uk)

waterfall plot : yes , I played a bit with elevation and azimut angle to have the best rendering

Camera line of sight - MATLAB view - MathWorks France

this is done here :

view(-40,60); % view angles

[S]

So I meant transforming the filters into their equivalent rectangles and plotting them to see the comparison and see what their bandwidth would be.

Also for the bottom 3D plot you provided, I was thinking what if to better show the output, we choose just one timestep ( for ex. .04 sec) and then plotted what each filters output would be and then connect those outputs with a line. So by doing this time is set so it is no longer the variable on the axis but the variable for each line. And we keep filter output on the axis. Does that make sense? So instead of our legend being made up of filter numbers it would be made up of timestamps and the x axis would be filter number, and the y would be output.

@Mathieu NOE

[Mathieu NOE]

hello again

to the first part of your question, this is a modified code that allows you to plot the ERB vs CenterFreqs (as in the publication : [1] Slaney, Malcolm. "An Efficient Implementation of the Patterson-Holdsworth Auditory Filter Bank." Apple Computer Technical Report 35, 1993. : PattersonÕs Ear Model (purdue.edu) )

for this purpose , the gammatone function has now also erb in its outputs :

[IR,erb] = gammatone(order, centerFrequency, fs)

full code :

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); % 
input = sin(2*pi*signal_freq*t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    [IR,erb(ii)] = 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
plot(CenterFreqs,erb)
title('ERB vs CenterFreqs');
set(gca,'XScale','log');
xlabel('CenterFreqs(Hz)');
ylabel('ERB (Hz)');

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,erb] = 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]

hello again

to the second part of your request , I don't see what is this for ? what does this line represent or what is your idea behind that ?

Also for the bottom 3D plot you provided, I was thinking what if to better show the output, we choose just one timestep ( for ex. .04 sec) and then plotted what each filters output would be and then connect those outputs with a line.

[S]

@Mathieu NOE

Hi! So what I am trying to show is the values of each filter (filters now being the horizontal axis) and how it changes with the varying time inputs. So to do this I was thinking of changing the legend from filter numbers to timestamps to better be able to see this.