% m-file file name: fsd1.m
% Frequency Sampling Design Given Samples of the Amplitude Response.
% Written by Dr. James S. Kang, Professor
% Department of Electrical and Computer Engineering
% California State Polytechnic University, Pomona
% N = length (order) of the filter = Number of coefficients
% Av = amplitude vector = samples of the amplitude response = A(n)
% type = type (1,2,3,4) of the FIR filter (Enter 'type1' for type 1, etc.)
% A=[1 1 0 0 0 0 1]';
% fsd1(A,'type1')
function[h,H] = fsd1(Av,type)
N=size(Av,1);
disp('Number of Samples')
disp(N)
if (type=='type1')&(rem(N,2)==0)
    N=N+1;
    Av(N)=Av(2);
end
if (type=='type2')&(rem(N,2)==1)
    N=N+1;
    Av(N)=Av(2);
end
if (type=='type3')&(rem(N,2)==0)
    N=N+1;
    Av(N)=Av(2);
end
if (type=='type4')&(rem(N,2)==1)
    N=N+1;
    Av(N)=Av(2);
end
disp('Number of Samples (Adjusted)')
disp(N)
%
% Calculating the coefficients of the filter
%
if (type == 'type1')
   q = (N-1)/2+1;
   for k = 1:N
       h(k)=(1/N)*Av(1);
       for n = 2:q
           h(k)=h(k)+(1/N)*2*Av(n)*cos((2*pi/N)*(k-q)*(n-1));
       end
   end
end
if (type == 'type2')
   q = N/2;
   for k = 1:N
       h(k)=(1/N)*Av(1);
       for n = 2:q
           h(k)=h(k)+(1/N)*2*Av(n)*cos((2*pi/N)*(k-1-(N-1)/2)*(n-1));
       end
   end
end
if (type == 'type3')
   q = (N-1)/2+1;
   for k = 1:N
       h(k)=(1/N)*Av(1);
       for n = 2:q
           h(k)=h(k)+(1/N)*2*Av(n)*sin((2*pi/N)*((N-1)/2-(k-1))*(n-1));
       end
   end
end
if (type == 'type4')
   q = N/2;
   for k = 1:N
       h(k)=(1/N)*Av(N/2+1)*sin(pi*((N-1)/2-(k-1)));
       for n = 2:q
           h(k)=h(k)+(1/N)*2*Av(n)*sin((2*pi/N)*((N-1)/2-(k-1))*(n-1));
       end
   end
end
disp('Coefficients')
disp(h)
%
for k = 1:N
   t(k)=k-1;
end
figure
stem(h);grid;
title('Coefficients of the FIR Filter')
xlabel('Time')
ylabel('Amplitude')
%
% Calculating and plotting the frequency response.
%
M = 400;
for n = 1:M+1
   f(n) = (n-1)/M;
   q(n) = pi*f(n);
   HReal(n)=h(1);
   HImag(n)=0;
   for k = 2:N
      HReal(n) = HReal(n) + h(k)*cos(q(n)*(k-1));
      HImag(n) = HImag(n) - h(k)*sin(q(n)*(k-1));
   end
   mag(n) = sqrt(HReal(n)^2+HImag(n)^2);
   dB(n) = 20*log10(mag(n));
   if dB(n) < -200 dB(n) = -200; end
   pha(n) = atan2(HImag(n),HReal(n))/pi; 
end
figure
plot(f,mag);grid;
title('Magnitude Response')
xlabel('Frequency (1 = pi = fs/2)')
ylabel('Magnitude')
figure
plot(f,dB);grid;
title('Magnitude Response in dB')
xlabel('Frequency (1 = pi = fs/2)')
ylabel('Magnitude in dB')
figure
plot(f,pha);grid;
title('Phase Response')
xlabel('Frequency (1 = pi = fs/2)')
ylabel('Phase in Radians/pi')
%
% Calculating and plotting poles and zeros.
%
nozeros = 0;
k = 1;
while abs(h(k)) < 10^(-10) & k == nozeros + 1 & k < (N-1)/2
   nozeros = nozeros + 1;
   k = k+1;
end
stop = N - 2*nozeros;
for k = 1:stop
  h1(k) = h(k+nozeros);
end
dz=ROOTS(h1);
for k = 1:stop-1
  dp(k) = 0;
end
for m = 1:201
   x(m) = (m-101)/100;
   ty(m) = sqrt(1-x(m).*x(m));
   by(m) = - ty(m);
end
figure
plot(x,ty)
grid on
hold on
plot(x,by) 
plot(real(dp),imag(dp),'rx')
plot(real(dz),imag(dz),'bo')
axis equal
title('Pole-Zero Diagram')
xlabel('Real Part')
ylabel('Imaginary Part')
disp('Digital Pole Locations')
disp(dp)
disp('Digital Zero Locations')
disp(dz)