Series Filter Loop Error

4 次查看(过去 30 天)
S
S 2024-1-31
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)
  15 个评论
Mathieu NOE
Mathieu NOE 2024-2-13
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
Mathieu NOE 2024-2-13
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));

请先登录,再进行评论。

采纳的回答

Mathieu NOE
Mathieu NOE 2024-2-7
移动:Mathieu NOE 2024-2-13
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
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
  4 个评论
S
S 2024-2-26
@Mathieu NOE I just posted it as a new question in case you wanted me to repost!
Mathieu NOE
Mathieu NOE 2024-2-27
hello @S
ok , I'll have a look in your other post !

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Simulation, Tuning, and Visualization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by