While loop using up to input and fixing input definition
4 次查看(过去 30 天)
显示 更早的评论
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 个评论
Mathieu NOE
2024-3-18
- 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.
回答(1 个)
Mathieu NOE
2024-3-18
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 个评论
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!