While loop using up to input and fixing input definition
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Mathieu NOE
el 18 de Mzo. de 2024
- 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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1644871/image.png)
- 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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1644876/image.png)
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1644881/image.png)
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.
Respuestas (1)
Mathieu NOE
el 18 de Mzo. de 2024
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 comentarios
Mathieu NOE
el 2 de Abr. de 2024
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.
Ver también
Categorías
Más información sobre Spectral Measurements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!