While loop using up to input and fixing input definition

Question

S 2024-3-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2090171-while-loop-using-up-to-input-and-fixing-input-definition

评论： S 2024-4-4

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

Mathieu NOE 2024-3-5

编辑：Walter Roberson 2024-3-5

在 MATLAB Online 中打开

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 2024-3-16

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 2024-3-18

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 2024-3-18

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2090171-while-loop-using-up-to-input-and-fixing-input-definition#answer_1426866

在 MATLAB Online 中打开

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

13 个评论
显示 11更早的评论隐藏 11更早的评论

Mathieu NOE 2024-3-25

在 MATLAB Online 中打开

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 2024-3-25

在 MATLAB Online 中打开

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

Mathieu NOE 2024-3-26

在 MATLAB Online 中打开

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

Mathieu NOE 2024-4-2

在 MATLAB Online 中打开

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 2024-4-2

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 2024-4-4

@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.

请先登录，再进行评论。

While loop using up to input and fixing input definition

7 个评论
显示 5更早的评论隐藏 5更早的评论

回答（1 个）

13 个评论
显示 11更早的评论隐藏 11更早的评论

另请参阅

类别

标签

Community Treasure Hunt

While loop using up to input and fixing input definition

7 个评论 显示 5更早的评论隐藏 5更早的评论

回答（1 个）

13 个评论 显示 11更早的评论隐藏 11更早的评论

另请参阅

类别

标签

Community Treasure Hunt

7 个评论
显示 5更早的评论隐藏 5更早的评论

13 个评论
显示 11更早的评论隐藏 11更早的评论