Loop outputting deleted object?

I am trying to add a loop that plots each filters output, but am getting an error that I do not understand. Its the last loop trying to make figure 4. Thank you for your time!

close all
clear all
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)');

figure(3)
zer = -0.5; %Zero location
pol = 0.9*exp(j*2*pi*[-0.3 0.3]'); %complex pole location
zplane(zer,pol)

[b,a] = zp2tf(zer,pol,1); %Convert poles to transfer function form
fvtool(b,a)

This functionality is not available on remote platforms.

fvtool(b,a,'Analysis','polezero')
zplane(b,a) %zplane finds roots of numerator and denominator using the roots function
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

Error:

Error using te11

Invalid or deleted object.

[Dyuman Joshi]

There is no variable named te11 in the code above.

Share the full error message i.e. all of the red text, and provide a code (and any corresponding data required) that can be ran to reproduce the error, thus allowing us to provide suggestions/solutions accordingly.

Also, do not use builtin function names as names for variables (or scripts for that matter), the function in question is input. Rename the variable to something else.

[S]

Thats just the name of the file!

[Voss]

@S: Please show the complete error message, which includes the line number and the text of the line of code on which the error occurs.

[S]

>> t226

Error using t226

Invalid or deleted object.

>>

This is all it gives! @Voss

[S]

t226 is just the new file name I saved it again to see if it would fix

[S]

@Mathieu NOE

[Walter Roberson]

Start by doing

clear t226

at the command line, to get rid of the variable that is shadowing the file.

Note:

You should rarely have

clear all

in your code. Instead, explicitly initialize variables.

Answer 1

[Mathieu NOE]

hello and welcome back again @S

I see the added code below figure(4) is a copy paste from some older code , so this line is incorrect :

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

notice that the gammatone function we have today is not about applying the gammatone filter to an input signal but just to generate the IR of the filter

so what you should be doing instead is to apply the IR coefficients to your input signal

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

also I noticed that your input signal has only 200 samples (representing one period of a signal at 100 Hz, according to fs = 20e3 Hz) , but look again at the IR plots, you see we need a much longer input signal to reach a steady state output; here I modified slightly your code so we generate now 10 periods of input signal (you may just want to plot the last period - this can be done latter if you want)

so full code corrected : (NB I have commented what was in section figure (3) as this is not really the focus here)

clc
clearvars
close all

fs = 20e3; 
numFilts = 32; % 
filter_number = 5;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
% 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(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_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(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)');

% figure(3)
% zer = -0.5; %Zero location
% pol = 0.9*exp(j*2*pi*[-0.3 0.3]'); %complex pole location
% zplane(zer,pol)
% 
% [b,a] = zp2tf(zer,pol,1); %Convert poles to transfer function form
% % fvtool(b,a)
% fvtool(b,a,'Analysis','polezero')
% zplane(b,a) %zplane finds roots of numerator and denominator using the roots function

figure(4)
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

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

Thank you for this! It seems like the filter outputs are getting larger rather than smaller, what we want to see is that originally the output of filter 1 is the largest then as filter 1's output is filter 2's input the output of each filter should gradually be getting smaller. Would this issue be becuase if IR_array?

[Mathieu NOE]

hello

see slightly modified code below

we have here 5 center frequencies for the filters : 50.0000 58.8939 69.3699 81.7094 96.2438 Hz

if we use the same input signal at 100 Hz, we can compute what the filters gains are at 100 Hz :

g_100 = 0.0419 0.0758 0.1666 0.4380 0.9486

also I the code below, we can compute what are the output signal max amplitude (I took the last cycle to avoid the initial transient )

max_out = 0.0438 0.0790 0.1606 0.4300 0.9653

we can see here that the signal output amplitude is matching the gain of the filters at the same frequency 100 Hz (and the max amplitude is obtained for the 5th filter which is what we know already from the Bode plot)

code :

fs = 20e3; 
numFilts = 32; % 
filter_number = 5;
order = 4;
CenterFreqs = logspace(log10(50), log10(8000), numFilts);
% input signal definition 
Nperiods = 10; % we need more than 1 period of signal to reach the steady state output (look a the IR samples)
samples = 200*Nperiods; % % for 1 period (2*pi) we have 200 samples at fs = 20e3, so signal freq = 20e3/200 = 100 Hz
t = linspace(0,2*pi*Nperiods,samples); % 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(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_array{ii} = IR*a; % scale IR and store in cell array afterwards
    plot(IR_array{ii})
end
title('Impulse Response');
hold off

figure(2)
hold on
for ii = 1:filter_number
    [h{ii},f] = freqz(IR_array{ii},1,1024*2,fs); % now store h{ii} after amplitude correction
    plot(f,20*log10(abs(h{ii})))
    LEGs{ii} = ['Filter # ' num2str(ii)]; %assign legend name to each
   % filters gain at 100 hz
    g_100(ii) = interp1(f,abs(h{ii}),100);

end
title('Bode plot');
set(gca,'XScale','log');
xlabel('Freq(Hz)');
ylabel('Gain (dB)');
ylim([-100 6]);
legend(LEGs{:})
g_100     % filters gain at 100 hz


% figure(3)
% zer = -0.5; %Zero location
% pol = 0.9*exp(j*2*pi*[-0.3 0.3]'); %complex pole location
% zplane(zer,pol)
% 
% [b,a] = zp2tf(zer,pol,1); %Convert poles to transfer function form
% % fvtool(b,a)
% fvtool(b,a,'Analysis','polezero')
% zplane(b,a) %zplane finds roots of numerator and denominator using the roots function

figure(4)
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
    
    max_out(ii) = max(abs(output(ii,samples*(Nperiods-1)/Nperiods:samples))); % max amplitude of output signal - taken on the last period 
    
end
legend(LEGs{:})
% legend('Show')
hold off
max_out  % max amplitude of output signal - taken on the last period 


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

[Mathieu NOE]

there is something I don't understand in your previous comment : "as filter 1's output is filter 2's input "

for me the gammatone filters are used as a parallel filter bank, so each filter is feeded with the same input signal , but you are saying something different . why would filter 2's input be filter 1's output ?

[S]

You are completely right! Thank you!

[Mathieu NOE]

as always, my pleasure !

Answer 2

[Image Analyst]

>> test5

Error using input

Not enough input arguments.

Error in test5 (line 49)

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

Do not use "input" as the name of your variable in your argument list. It is the name of a built-in function and it requires input argument, just like the error said. Call it something else, like order instead of input.