Series Filter Loop Error

I am trying to make a loop that goes through different filters in cascade, so the output of one is the input to the next. I want to use the built in matlab gammatone function to do this. I have been stuck here for a bit. Thank you for your time!

fs = 20e3; 
numFilts = 32;
BW = 100; %Filter Bandwidth
filter_number = 5;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
CF1 = CenterFreqs - BW/2; %Lower cutoff frequency
CF2 = CenterFreqs + BW/2; %Upper cutoff frequency
t = linspace(0,2*pi,200);
input = sin(t) + 0.25*rand(size(t));
for ii = 1:filter_number
    output = gammatone(input, CenterFreqs, fs)
end

Unrecognized function or variable 'gammatone'.

[h{ii},f] = freqz(gammatone{ii}.Coefficients,1,4*8192,fs);

Error: (I do not have a line 32?, also I did not make gammatone so I don't think there is anything wrong with it)

Not enough input arguments.

Error in gammatone (line 32)

A = 10^((q_factor) / 20);

Error in t131 (line 13)

output = gammatone(input, CenterFreqs, fs)

[Mathieu NOE]

hello and welcome back

I see you are still working hard on the same topic....

the error message tells you there is no know function gammatone in your path or working directory

did you create one / downloaded one from Fex ?

there is a similar function : 'gammatoneFilterBank' , from the Audio Toolbox.

type this to double check :

which gammatone -all

[S]

yes, it is from a toolbox:

> help gammatone

% Function for gammatone filter

function output = gammatone(input, centerFrequency, fs)

% Design a gammatone filter

earQ = 9.26449;

minBW = 24.7;

erb = (centerFrequency / earQ) + minBW;

B = 1.019 .* 2 .* pi .* erb;

g = 10^(-3 / 20); % 3 dB down from peak

f = centerFrequency;

tau = 1. / (2. .* pi .* B);

% Apply the gammatone filter

t = 0:1/fs:(length(input)-1)/fs;

temp = (t - tau) > 0;

output = (t.^(4 - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2 .* pi .* f .* (t - tau)) .* temp;

end

%plot freq eresponse and impulse response try at 250 hz make sure that freq

%shows up at that signal freq, initialize all arrays and zero them out

>> which gammatone -all

/Users/s/Desktop/E/gammatone.m

[Mathieu NOE]

well , seems to me there is more than just one problem

according to your first post (error message : Unrecognized function or variable 'gammatone'.) it was only a problem that the gammatone function seemed to be unknow from the matlab path

when I use your gammatone code above and create that function , it works in terms of calling this function and getting some output, but what I see is not what I expected (or there is missing code )

also you must correct the line inside the for loop as CenterFreqs must be indexed

for ii = 1:filter_number

output = gammatone(input, CenterFreqs(ii), fs)

end

output

For me a gammatone should be some king of bandpass filter , which is obviously not what the function does here

if we look at the ouput , we have this , and most funny the output is even not function of the input (except i has same size)

where does this function come from and what was the supposed purpose of it ?

I suspect it has to do with generating the impulse response of the gammatone filter , but for me the function (as provided) is incomplete

Gammatone filter - Wikipedia

are you trying to combine input and output to generate the impulse response of the gammatone filter ? then the "true" filtering must happen latter in the code

fs = 20e3; 
numFilts = 32; % why this ? not used
BW = 100; %Filter Bandwidth
filter_number = 5;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
CF1 = CenterFreqs - BW/2; %Lower cutoff frequency
CF2 = CenterFreqs + BW/2; %Upper cutoff frequency
t = linspace(0,2*pi,200);
input = sin(t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    output = gammatone(input, CenterFreqs(ii), fs);
    plot(output)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function output = gammatone(input, centerFrequency, fs)
      % Design a gammatone filter
      earQ = 9.26449;
      minBW = 24.7;
      erb = (centerFrequency / earQ) + minBW;
      B = 1.019 .* 2 .* pi .* erb;
      g = 10^(-3 / 20); % 3 dB down from peak
      f = centerFrequency;
      tau = 1. / (2. .* pi .* B);
  
      % Apply the gammatone filter
      t = 0:1/fs:(length(input)-1)/fs;
      temp = (t - tau) > 0;
      output = (t.^(4 - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2 .* pi .* f .* (t - tau)) .* temp;
  end

[Mathieu NOE]

FYI there is this Fex submissions for computing gammatone filterbanks

Gammatone filterbank - File Exchange - MATLAB Central (mathworks.com)

see also

Gammatone-based (auditory) spectrograms - File Exchange - MATLAB Central (mathworks.com)

still wondering .... I coud not find any gammatone function in any matlab toolbox - where does it come from ?

[S]

numFilts is used in the CenterFreqs line! Also when I copy that exact code you gave at the top I just get a blank graph output. I did go through that file exchange, that must be where I got the function from. How come I am not getting that output graph? I am still getting:

Not enough input arguments.

Error in gammatone (line 32)

A = 10^((q_factor) / 20);

Error in t131 (line 37)

output = gammatone(input, CenterFreqs(ii), fs);

[Mathieu NOE]

you must be calling another gammatone function

if you look again to the code I posted above , the gammatone function (from your own post above) has only 13 lines

so your error message at line 32 inside YOUR gammatone function clearly indicates we are not using the same function !

Error in gammatone (line 32)

A = 10^((q_factor) / 20);

NB also that there is no such code ( A = 10^((q_factor) / 20); ) below

can you please share your gammatone function ? It must have more than 32 lines as we know from now.

again , the code I posted above (according to your previous inputs) :

I highly suspect this is not the right gammatone function

fs = 20e3; 
numFilts = 32; % 
BW = 100; %Filter Bandwidth
filter_number = 5;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
CF1 = CenterFreqs - BW/2; %Lower cutoff frequency
CF2 = CenterFreqs + BW/2; %Upper cutoff frequency
t = linspace(0,2*pi,200);
input = sin(t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    output = gammatone(input, CenterFreqs(ii), fs);
    plot(output)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function output = gammatone(input, centerFrequency, fs)
      % Design a gammatone filter
      earQ = 9.26449;
      minBW = 24.7;
      erb = (centerFrequency / earQ) + minBW;
      B = 1.019 .* 2 .* pi .* erb;
      g = 10^(-3 / 20); % 3 dB down from peak
      f = centerFrequency;
      tau = 1. / (2. .* pi .* B);
  
      % Apply the gammatone filter
      t = 0:1/fs:(length(input)-1)/fs;
      temp = (t - tau) > 0;
      output = (t.^(4 - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2 .* pi .* f .* (t - tau)) .* temp;
  end

[S]

I did not notice that!! You're right! I found a 32 line file with the same title but different function name! I now resaved the function I shared above under that file and am able to get the same image as you but also am getting:

Brace indexing into the result of a function call is not supported. Assign the result of 'gammatone'

to a variable first, then brace index into it.

Error in t131 (line 22)

[h{ii},f] = freqz(gammatone{ii}.Coefficients,1,4*8192,fs);

[S]

also, what would the labeling on the axes be for the image?

[S]

Also I for some reason am not seeing an accept this answer/ thumbs up option on your answer?

[Mathieu NOE]

To your first comment above

I think youy are mixing some code from different horizons , I believe I have already seen this particular line (or something very similar) in one of your other post (same topic indeed) :

[h{ii},f] = freqz(gammatone{ii}.Coefficients,1,4*8192,fs);

Notive here gammatone is not a structure, it's a function that simply returns output (array) as we have declared it :

function output = gammatone(input, centerFrequency, fs)

so there is not gammatone{ii}.Coefficients field !! and you cannot do a freqz action on a non existent field

maybe you can do this (output is assumed to be the FIR coefficients as in the other post) , but I doubt that the result is matching your expectations , as the Bode plot does not look like a bank of gammatone filters

fs = 20e3; 
numFilts = 32; % 
BW = 100; %Filter Bandwidth
filter_number = 5;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
CF1 = CenterFreqs - BW/2; %Lower cutoff frequency
CF2 = CenterFreqs + BW/2; %Upper cutoff frequency
t = linspace(0,2*pi,200);
input = sin(t) + 0.25*rand(size(t));
figure
hold on
for ii = 1:filter_number
    output = gammatone(input, CenterFreqs(ii), fs);
%     plot(output)
    [h{ii},f] = freqz(output,1,1024,fs);
    plot(f,20*log10(abs(h{ii})))
end
set(gca,'XScale','log');
xlabel('Freq(Hz)');
ylabel('Gain (dB)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 function output = gammatone(input, centerFrequency, fs)
      % Design a gammatone filter
      earQ = 9.26449;
      minBW = 24.7;
      erb = (centerFrequency / earQ) + minBW;
      B = 1.019 .* 2 .* pi .* erb;
      g = 10^(-3 / 20); % 3 dB down from peak
      f = centerFrequency;
      tau = 1. / (2. .* pi .* B);
  
      % Apply the gammatone filter
      t = 0:1/fs:(length(input)-1)/fs;
      temp = (t - tau) > 0;
      output = (t.^(4 - 1)) .* exp(-2 .* pi .* B .* (t - tau)) .* g .* cos(2 .* pi .* f .* (t - tau)) .* temp;
  end

[Mathieu NOE]

hello again

how is it going ?

probelm (finally) solved ?

all the best

[S]

Thank you so much @Mathieu NOE that was very helpful! I really appreciate your help:)

In the final line IR = why did you remove the g from cos?

[S]

Can you leave an 'answer this question' so I can accept your answer, for some reason I am not seeing that option right not

[Mathieu NOE]

I will transform my comment above into a answer so you can accept it (and tx for that ! )

for the g factor , we don't need it inside the gammatone function itself as in the main code we are doing it in these lines :

(up to you to choose a 0 or -3 dB peak amplitude of the filters)

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

[Mathieu NOE]

FYI, the lines of code above could also be written this way if you prefer

    % scale the IR amplitude to match a required Bode plot magnitude
    g = 1; % the max modulus is 0 dB
%     g = 10^(-3 / 20);    % or if you prefer - 3dB
     a = g/max(abs(tmp));

Answer 1

[Mathieu NOE]

hello again

so I wanted to learn a bit about gammatone filters , so I read this , which I suppose i also from where you started your code

PattersonÕs Ear Model (purdue.edu)

there are a few modifications I made in the code , see the original lines commented are replaced by new code

so the function output is the IR of the filter , we can use it latter on (not yet coded) to filter any time signal

IMHO, now the IR and Bode plots are correct

here we plot only the first 5 filters , even though the CF are computed for 32 filters

fs = 20e3; 
numFilts = 32; % 
% BW = 100; %Filter Bandwidth
filter_number = 5;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
% CF1 = CenterFreqs - BW/2; %Lower cutoff frequency
% CF2 = CenterFreqs + BW/2; %Upper cutoff frequency
% t = linspace(0,2*pi,200);
% input = sin(t) + 0.25*rand(size(t));
figure(1)
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 = IR*a;
    [h{ii},f] = freqz(IR,1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(IR)
end
title('Impulse Response');
hold off

figure(2)
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)');

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

Hi! If I wanted to use the input and have it go through all the filters up to filter number, would I add a loop at the bottom? Should I post this as a new question? Thanks again!:)

[S]

when I add:

figure(4)

hold on

for ii = 1:filter_number

output(ii,:) = gammatone(input, CenterFreqs(ii), fs);

plot(output(ii,:))

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

legend(LEGs{:})

legend('Show')

end

hold off

I get:

Error using t226

Invalid or deleted object.

[S]

@Mathieu NOE I just posted it as a new question in case you wanted me to repost!

[Mathieu NOE]

hello @S

ok , I'll have a look in your other post !