MATLAB Answers

Multiplying from a loop

1 view (last 30 days)
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

Accepted Answer

William Rose
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 Comments
Michael Akiva
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.

Sign in to comment.

More Answers (1)

William Rose
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.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!