Multiple Bandpass Filter Creation

I am trying to make a cascade of bandpass filters. I want there to be 32 total filters covering the range of 50 to 8000. I am trying to do this efficiently by creating a variable for the upper and lower cutoff frequencies. How would I fix this so that the bpfilt cycles through all my CF1/CF2s in order? The goal would be to then put an input signal in and have it run through all the 32 filters in cascade. Thank you for your time!

close all

fs = 16e3;

range = [50 8000];

CF1=linspace(50, 8000, 32) -50;

CF2=linspace(50, 8000, 32) +50;

bpfilt = designfilt('bandpassfir', ...

'FilterOrder',20,'CutoffFrequency1',CF1, ...

'CutoffFrequency2',CF2,'SampleRate',fs);

freqz(bpfilt.Coefficients,1,[],fs)

Answer 1

[Paul]

fs = 16e3;
%range = [50 8000];
CF1=linspace(50, 8000, 32) -50;
CF2=linspace(50, 8000, 32) +50;

Can't use the first and last elements of CF1 and CF2 respectively. CF1(1) = 0, which is not a valid input for CutOffFrequency1. This issue could be dealt with by making the first filter a low pass. CF2(end) is greater than the Nyquist frequecy.

for ii = 1:numel(CF1)-2
    bpfilt{ii} = designfilt( ...
        'bandpassfir', ...
        'FilterOrder',20, ...
        'CutoffFrequency1',CF1(ii+1), ...
        'CutoffFrequency2',CF2(ii+1), ...
        'SampleRate',fs);
    
    [h{ii},f] = freqz(bpfilt{ii}.Coefficients,1,4*8192,fs);
    
end

Frequency responsee of each filter

figure
subplot(211)
plot(f,db(abs([h{:}])));
subplot(212);
plot(f,180/pi*(angle([h{:}])));

Frequency response of cascaded filters in series; this doesn't look like very useful.

Edit: In fact, the values on the magnitude plot on are so small that they are probably rounding error and are effectively zero (in absolute value). Are you sure you want the filters cascaded?

figure
subplot(211)
plot(f,db(abs(prod([h{:}],2))));
subplot(212);
plot(f,180/pi*(angle(prod([h{:}],2))));

Unfortunately, digitalFilter objects can't be mulitplied.

In theory, the coefficients of each filter can be combined as follows into a single filter

b = 1;
for ii = 1:numel(bpfilt)
    b = conv(b,bpfilt{ii}.Coefficients);
end

Compare the frequency response of the cascade to the cascade of the frequency responses

hcascade = freqz(b,1,f,fs);

figure
subplot(211)
hold on
plot(f,db(abs(prod([h{:}],2))));
plot(f,db(abs(hcascade)));
subplot(212);
hold on
plot(f,180/pi*(angle(prod([h{:}],2))));
plot(f,180/pi*angle(hcascade))

In practice, that doesn't look like such a good idea.

If going to pursue this path, it's probabably best to use a loop and filter the input in stages. That will work, in theory, for forward filtering. Not sure if that's appropriate for forward/backward filtering.

[S]

Thank you for your time! If I was to take an impulse/ signal and put it through the first image could I do it like this? Additionally where is the 8192 value from in the freqz line?

impulse_input = 0*fs;

impulse_input(1) = 1;

y=filter(impulse_input)

[Paul]

I don't know what you mean by "the first image."

Do you want the impulse response of each filter?

Or do you want the impulse response of the cascade of all the filters in series?

[S]

The second one! I want the impulse response of the cascade of all the filters in series. Then to be able to input a filter number and be able to see the output after that specfic filter in the cascade. Like

filter_number = 10;

[Paul]

Define the impulse input.

Write a loop to filter the input and use the output of the n-1th filter as the input to the nth filter.

filter only returns as many elements in the output as there are in the input, so make sure to include enough elements in the original impulse input to capture all of the elements of the output at the end of the cascade.

In the loop you can sequentially store the output of each successive filter for post-processing.

If still having problems after putting in some effort, feel free to post your code and explain what's going wrong or what you don't understand.

[S]

so by following: https://www.mathworks.com/help/matlab/ref/filter.html#d126e499602

I was trying to use the last option, but it says i do not have enough input arguments. Is this the correct way to get it to use the designed filter? Or is this not what you meant by use filter

impulse_input = 0*fs;

impulse_input(1) = 1;

for ii = 1:numel(CF1)-1

[y,zf]=filter(impulse_input)

end

[Paul]

impulse_input is the input to the filter, but you also have to tell filter what filter it needs to apply to the input. You can use the coefficients of each filter in the bpfilt cell array (assuming you used code that's posted above or something similar), or you can directly use each filter in the bpfilt cell array as the input to filter. For either case, see the Tips section of the doc page you linked.

As stated above, you'll need to design your loop so that the output of filter on one iteration is the input to filter on the next iteration.

If you want to continue the discussion, post your entire code that shows how your're building the array of filters and how you're trying to develop the impulse response of the cascade. It's difficult to help without seeing all that you're doing.

[Paul]

Are you sure you want the filters cascaded in series? Cascading all them makes the frequency response of the cascade effectively zero.

[Paul]

I suggest you continue this discussion with @Sulaymon Eshkabilov on your other question

[S]

Where is the 8192 coming from in this line?

[h{ii},f] = freqz(bpfilt{ii}.Coefficients,1,4*8192,fs);

[Paul]

In that usage of freqz, the third argument is the number of frequency points at which to evalute the response. I just picked a number. Probably started with 8192 to see what it looked like, and then reran after multiplying by 4 to see if the additional points revealed anything interesting. 8192 and 4*8192 are both powers of 2, which might make freqz run a bit more efficiently.

[S]

oh ok that makes sense! I though that was zero padding

So it is following this format correct?

[h,f] = freqz(___,n,fs)

If so why is there a 1 in our statement when the maximum allowed looks like its 3 items in the () based on this page:

https://www.mathworks.com/help/signal/ref/freqz.html

[Paul]

The convention used by the documentation for that underbar is sometimes hard to understand (I get confused myself sometimes). In this case, that underbar means that you can replace the underbar with any of the arguments to freqz that precede n in any of the syntax cases above it. So I was using the syntax

[h,f] = freqz(b,a,n,fs)

In this problem, each bpfilt is a FIR filter. For the bpfilt FIR filter, the denominator is 1, hence the entry of 1 in the second argument in the call to freqz.

[S]

OH! That makes sense! Thank you for explaining!