Indexing into a loop

So I am trying to index into this second loop in order to be able to for example at the top of the code be able to change filter_number and be able to see what the output of that filter is. Would I have to change how I wrote the loop? Thanks for your time!

close all
fs = 16e3;
numFilts = 32;
filter_number = 10;
%range = [50 8000];
CF1=linspace(50, 8000, numFilts+2) -50;
CF2=linspace(50, 8000, numFilts+2) +50;
for ii = 1:numel(CF1)-2
    bpfilt{ii} = designfilt( ...
        'bandpassfir', ...
        'FilterOrder',20, ...
        'CutoffFrequency1',CF1(ii+1), ...
        'CutoffFrequency2',CF2(ii+1), ...
        'SampleRate',fs);
    
end
[h{ii},f] = freqz(bpfilt{ii}.Coefficients,1,4*8192,fs);
figure
subplot(211)
plot(f,db(abs([h{:}])));
title('Magnitude')
ylabel('Magnitude (dB)')
xlabel('Frequency (Hz)')
subplot(212);
plot(f,180/pi*(angle([h{:}])));
title('Phase')
ylabel('Phase (degrees)')
xlabel('Frequency (Hz)')
impulse_input = 0*fs;
impulse_input(1) = 1;
%%
% Reference signal with some white noise to benchmarch the created filter performances
t = linspace(0,2*pi,200);
rng(13)                                        % Make it repeatable
x = sin(t) + 0.25*rand(size(t));  % Ref Signal with a noise
% Simulation of 1-D digital filter:  x_filtered = filter(b, a, x); 
a = 1;    
figure
hold on
CoLoR = rand(numel(bpfilt), 3);
for ii = 1:numel(bpfilt)  %THIS ONE!!!
    [x_filtered(ii,:),zf(:,ii)]=filter(bpfilt{1, ii}.Coefficients, a, x);
    plot(t,x_filtered(ii,:), 'LineWidth', 1.25, 'Color', CoLoR(ii,:))
    LEGs{ii} = ['Filter # ' num2str(ii)];
    legend(LEGs{:})
end
plot(t, x, 'k-', 'LineWidth', 2, 'DisplayName', 'Data')
xlabel('t')
ylabel('x(t) & x_{filtered} (t)')
grid on
legend('Show')
fprintf('Number of generated filters: %d \n', numel(bpfilt))

Answer 1

[Mathieu NOE]

hello

I did quite a few modifications / simplifications

there are several variables that seem to never be used , so I commented them

at the end there is only one single variable that drives how many filters are considered : numFilts

I supposed that you wanted to have all filter bode plots overlaid, so that's the first major modification (see 1st loop )

the second loop is simply based on the results obtained in the first loop , here also I had to modify the time vector and signal definition so it's frequency is clearly defined (and you can make it match (or not) one of the BP filter central frequency)

%% parameters 
fs = 20e3;  % 16e3 is not enough if you need to see the effect of a bandpass filter with a center frequency at 8000 Hz
numFilts = 5;
% filter_number = 10;
%range = [50 8000];
BW = 100; % filter bandwith
CentralFreq = linspace(1000, 8000, numFilts);
CF1=CentralFreq -BW/2;
CF2=CentralFreq +BW/2;

%% plots 
figure(1)
hold on

for ii = 1:numel(CentralFreq)
    bpfilt{ii} = designfilt( ...
        'bandpassfir', ...
        'FilterOrder',20, ...
        'CutoffFrequency1',CF1(ii), ...
        'CutoffFrequency2',CF2(ii), ...
        'SampleRate',fs);
    
    [h{ii},f] = freqz(bpfilt{ii}.Coefficients,1,1024,fs);
    
    subplot(211)
    plot(f,db(abs([h{:}])));
    title('Magnitude')
    ylabel('Magnitude (dB)')
    xlabel('Frequency (Hz)')
    subplot(212);
    plot(f,180/pi*(angle([h{:}])));
    title('Phase')
    ylabel('Phase (degrees)')
    xlabel('Frequency (Hz)')
    
end

hold off

% impulse_input = 0*fs;
% impulse_input(1) = 1;

%%
% Reference signal with some white noise to benchmarch the created filter performances
% t = linspace(0,2*pi,200);
rng(13)                                        % Make it repeatable

dt = 1/fs;
samples = 200;
freq = 4500; % signal frequency in Hz

t = (0:samples-1)*dt; % time vector according to sampling frequency fs; maybe you want to increase the number of samples 
x = sin(2*pi*freq*t) + 0.25*rand(size(t));  % Ref Signal with a noise


% Simulation of 1-D digital filter:  x_filtered = filter(b, a, x); 
a = 1;    
figure(2)
hold on
CoLoR = rand(numFilts, 3);
for ii = 1:numFilts  %THIS ONE!!!
    [x_filtered(ii,:),zf(:,ii)]=filter(bpfilt{1, ii}.Coefficients, a, x);
    plot(t,x_filtered(ii,:), 'LineWidth', 1.25, 'Color', CoLoR(ii,:))
    LEGs{ii} = ['Filter # ' num2str(ii)];
    legend(LEGs{:})
end
plot(t, x, 'k-', 'LineWidth', 2, 'DisplayName', 'Data')
xlabel('t')
ylabel('x(t) & x_{filtered} (t)')
grid on

legend('Show')
fprintf('Number of generated filters: %d \n', numFilts)

Number of generated filters: 5 

[S]

Thank you for your time! How come this filtered data looks so different from mine? Additionally to index into the second loop and see the output after a specific filter could I do:

for ii = 1:numFilts(filter_number)

[S]

Also why did you change the freq range is 50 - 8000 not good?

[Mathieu NOE]

hello again

sorry if I have messed up your code , but if we look back at your original code , there are a couple of things I don't understand in your way of processing the filters and use them , here for example in your first portion of the code (I simply moved the plot operation inside the first loop so we have the overlaid Bode traces

so this is the first portion of your original code , barely modified as explained above :

I simply reduced numFilts from 32 to 5 just to get a clearer plot in first instance

notice that you declared filter_number = 10 , which is not used after in the code . what was supposed this to be used for ?

fs = 16e3;
% numFilts = 32;
numFilts = 5;
% filter_number = 10;
%range = [50 8000];
CF1=linspace(50, 8000, numFilts+2) -50

CF1 = 1×7

0 1325 2650 3975 5300 6625 7950

CF2=linspace(50, 8000, numFilts+2) +50

CF2 = 1×7

100 1425 2750 4075 5400 6725 8050

figure
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,1024,fs);

        subplot(211)
        plot(f,db(abs([h{:}])));
        title('Magnitude')
        ylabel('Magnitude (dB)')
        xlabel('Frequency (Hz)')
        subplot(212);
        plot(f,180/pi*(angle([h{:}])));
        title('Phase')
        ylabel('Phase (degrees)')
        xlabel('Frequency (Hz)')

    
end

so , there is a few points that are not clear in my mind .

If numFilts = 5 , why then you add +2 to numFilts in these lines :

CF1=linspace(50, 8000, numFilts+2) -50

CF2=linspace(50, 8000, numFilts+2) +50

and then you retrieve 2 here

for ii = 1:numel(CF1)-2

and inside the for loop the current index is shifted by +1

'CutoffFrequency1',CF1(ii+1), ...

'CutoffFrequency2',CF2(ii+1), ...

you specifies more CF1 and CF2 values that you are using afterwards

I tried in my first answer to make a code that is cleaner regarding this aspect but maybe I failled to understand what your final goal is - and this is maybe where we should start again : what do you try to do ?

beside that , why the need to generate the Bode plot with such a huge amount of points (4*8192) ; try with amuch lower value like 1024 and see that the quality of the Bode plot is sufficient (as I did above)

later on we'll discuss about the second part of your code , but I need to better understand your project first

[S]

I appreciate your thoughts on this. So I ultimately use filter_number later to observe different filter numbers. I removed both that +2 and -2 since we last talked! Also I thought that the CF1(ii+1) was moving it through one CF1 value to the next one so the loop runs all the CF1 values. Is that not what it is doing? My final goal is cascading bandpass filters with differenct center frequencies in series. Also I decressed 4*8192 to 1024 like you suggested and it gives the following error

Error using horzcat

Dimensions of arrays being concatenated are not consistent.

Error in t121 (line 26)

plot(f,db(abs([h{:}])));

[Mathieu NOE]

what happens with your original code

CF1=linspace(50, 8000, numFilts+2) -50;
CF2=linspace(50, 8000, numFilts+2) +50;
for ii = 1:numel(CF1)-2
    bpfilt{ii} = designfilt( ...
        'bandpassfir', ...
        'FilterOrder',20, ...
        'CutoffFrequency1',CF1(ii+1), ...
        'CutoffFrequency2',CF2(ii+1), ...
        'SampleRate',fs);

is that the first and last value of CF1 and CF2 that you generated above are not used

so the simple question on my side is actually waht CF1 / CF2 values do you want ?

for the error you mention ( I don't have) , I will answer tomorrow as it's getting late here

[Mathieu NOE]

just to illustrate my previous comment , when I run the (almost) original code

close all
fs = 16e3;
% numFilts = 32;
numFilts = 5;
% filter_number = 10;
%range = [50 8000];
CF1=linspace(50, 8000, numFilts+2) -50

CF1 = 1×7

0 1325 2650 3975 5300 6625 7950

CF2=linspace(50, 8000, numFilts+2) +50

CF2 = 1×7

100 1425 2750 4075 5400 6725 8050

figure
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,1024,fs);

        subplot(211)
        plot(f,db(abs([h{:}])));
        title('Magnitude')
        ylabel('Magnitude (dB)')
        xlabel('Frequency (Hz)')
        subplot(212);
        plot(f,180/pi*(angle([h{:}])));
        title('Phase')
        ylabel('Phase (degrees)')
        xlabel('Frequency (Hz)')

    
end

you can see the plot does only show 5 curves , whereas CF1 and CF2 are both arrays of length = 7.

CF1 = 0 1325 2650 3975 5300 6625 7950

CF2 = 100 1425 2750 4075 5400 6725 8050

as I said the first and last values of CF1 / CF2 are not used as you can see there is no Bode plot corresponding to these frequencies.

also I don't have the error you mentionned

FYI, as you don't seem to use each individual h output , we can simplify your code and change these 2 lines

original code [h{ii},f] = freqz(bpfilt{ii}.Coefficients,1,1024,fs);

modified : [h,f] = freqz(bpfilt{ii}.Coefficients,1,1024,fs); % no need to index h

and also

original code plot(f,db(abs([h{:}])));

plot(f,180/pi*(angle([h{:}])));

can be modified : plot(f,db(abs(h)));

plot(f,180/pi*(angle(h)));

so all in all :

fs = 16e3;
% numFilts = 32;
numFilts = 5;
% filter_number = 10;
%range = [50 8000];
CF1=linspace(50, 8000, numFilts+2) -50

CF1 = 1×7

0 1325 2650 3975 5300 6625 7950

CF2=linspace(50, 8000, numFilts+2) +50

CF2 = 1×7

100 1425 2750 4075 5400 6725 8050

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

        subplot(211)
        plot(f,db(abs(h))); hold on
        title('Magnitude')
        ylabel('Magnitude (dB)')
        xlabel('Frequency (Hz)')
        subplot(212);
        plot(f,180/pi*(angle(h)));hold on
        title('Phase')
        ylabel('Phase (degrees)')
        xlabel('Frequency (Hz)')

    
end