Multiplying from a loop

Hi,

So i made a notch filter in my class and i am having a bit of a problem down here.

the problem is, mistakenly i made 3 filters instead of 1.

Adding the plot i got from this code -

so, instead of having 3 graph in the same plot, i need to multiply them and recieve 1 filter.

this is my code -

clear ;
r=0.9;
for n=[1:3] 
    z=[exp(-i*n*pi/4)  exp(i*n*pi/4) ];
    p=r*z;
    phi=pi*(0:500)/500;
    b=poly(z);
    a=poly(p);
    Freq = freqz(b,a,phi);
    H = 20*log10(r*abs(Freq));
    idx = find(abs(phi-(n*pi/4-0.05))<0.2);
    idx = mean(idx);
    while  (H(idx)>-0.9 || H(idx)<-1)
            r = r+0.001;
            p=r*z;
            b=poly(z);
            a=poly(p);
            H = 20*log10(r*abs(freqz(b,a,phi)));
    end
    plot(phi,H);
    hold on
end
axis([0.4 3 -2 0])
hold off

Answer 1

[William Rose]

@Michael Akiva,

Get rid of the loop

for n=[1:3] 
    ... 
end

but keep the code inside the loop, and define all 6 zeros at once:

%z=[exp(-i*n*pi/4)  exp(i*n*pi/4) ];
z=[exp(-i*pi/4),exp(i*pi/4),exp(-i*2*pi/4),exp(i*2*pi/4),exp(-i*3*pi/4),exp(i*3*pi/4)];

I do not totally understand what you are doing with

idx=mean(idx)

I think you are finding a frequency near the frequency of each pole pair, and adjusting the radius of that pole pair, in order to make the transfer function have a value between -0.9 and -1.0 dB at that frequency. But I don't understand how you are choosing the frequency f at which to enforce 0.9<H(f)<1.0.

Therefore I commented out that part of the code, in order to see the transfer function without the adjustment of r. Please comment on the rationale for the radius adjustment, and the choice of frequencies at which to do the adjustment. Figure below, then the code.

%AkivasFilterWR
%Notch filter code by Michal Akiva, 2021-05-22
%Modified by WCR 20210522
clear ;
r=0.9;
z=[exp(-i*pi/4),exp(i*pi/4),exp(-i*2*pi/4),exp(i*2*pi/4),...
    exp(-i*3*pi/4),exp(i*3*pi/4)];
p=r*z;
phi=pi*(0:500)/500;
b=poly(z);
a=poly(p);
Freq = freqz(b,a,phi);
H = 20*log10(r*abs(Freq));
%idx = find(abs(phi-(n*pi/4-0.05))<0.2);
%idx = mean(idx);
%while  (H(idx)>-0.9 || H(idx)<-1)
%        r = r+0.001;
%        p=r*z;
%        b=poly(z);
%        a=poly(p);
%        H = 20*log10(r*abs(freqz(b,a,phi)));
%end
plot(phi,H,'b.-');
hold on
axis([0 3.1 -2 2])
hold off

[William Rose]

@Michael Akiva

The code you provided in the original post ran without error on the first try. That is a very nice way to start!

[Michael Akiva]

Hi,

first of all, thank you for the answer.

1. I defiened all zeroes as you suggested and it worked! (also got rid of the loop)
2. the reason i used the while function is to find the radius that gives me -1dB. so i ran on all radiuses to find that radius.

Answer 2

[William Rose]

@Michael Akiva,

Here is code that adjusts r at the frequencies near each pole pair, until -0.9<=H(w)<=-1.0. Since there are three pole pairs, there are three values for r: one for each pole pair. The plot below shows Hinit (=H before adjusting the three r's), and final H, after adjusting three r's. The later r's are affected by the r's already adjusted. If you change the order of adjusting the r's (3,2,1 versus 1,2,3), you get slightly different values for the three r's, but the transfer function is not noticeably different. Code to make the figure is attached.

In your original code, you have

H = 20*log10(r*abs(Freq));

and

H = 20*log10(r*abs(freqz(b,a,phi)));

Why do you multiply abs(Freq) by r, and abs(freqz()) by r? That seems unjustified - in the sense that, if you filter a signal with

y=filter(b,a,x)

then the transfer function from x to y will not inclue the factor r that you are including in H. Also, since there are three r's (one for ech pole pair), it is not obvious which r shoudl be used. Therefore I have removed the factor of r from Hinit and from H.