# Multiplying from a loop

1 view (last 30 days)
Michael Akiva on 22 May 2021
Commented: Michael Akiva on 23 May 2021
Hi,
So i made a notch filter in my class and i am having a bit of a problem down here.
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

William Rose on 22 May 2021
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
##### 2 CommentsShowHide 1 older comment
Michael Akiva on 23 May 2021
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.

William Rose on 22 May 2021
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.