function [HD,line, hd, h, H, H_ABS, Hdb, PHASE]=FINAL(WINDOW,FLTR,N,W1,W2)
%
% Code for filter final project
%
% by Nasser Abbasi
%
% Northeastern UNIV.  (Network Student)
%
% Digital Signal Processing
%
% input 
% ----- 
% WINDOW = this is a number , has value 1,2,3,4 one for each type
% of window, 1=rectangular,2=hanning,3=hamming,4=blackman
% the user when calling function PRJ1 will supply this number
%
% FLTR= this is number, has values 1,2,3,4 , 1= BANDPASS, 2=BAND REJECT
%       3= HIGH-PASS , 4= TEST CASE, uses example of equation
%       8.1.18 page 583 in book, i need this to make sure all the
%       time my work is consistant by verifying this output with
%       book
%
% N = window size, this program assumes N is ODD.
%
% W1,W2 are the band-pass edge frequncies in RAD/SAMPLE
%
% output
% ------
% LINE= the line scaling from 0..Pi scaled such that resolution is
% pi/256, used to plot H(W) and PHASE
%
% ABS_H(W) = the MAGNITUDE of H(W) vector
% H = H(W) 
% PHASE = the phase angle vector 
% Hdb(W) = the mangnitued of H(w) in DB (frequency responce)
% hd = the desired h(n)
% h = the final h = h_window(n)*h_desired(n)
% HD = the desired H(w)  (this is given in the handout, i just plot
%      for sake of completion
%

%
% MAIN LOGIC OF PROGRAM
%
%  Depnding on Window, build the window signal, call it hw(n)
% 
%  build the desired h(n) based on the analytical solution allready
%  done by hand (see report), and cal it h(n)
%
%  establish h(n) =h(n) * hw(n)
%
%  find H(w) by taking the DFT of h(n)
%
%  find the magnitue of H(w) and phase of H(W)
%
%  IN ADDITION: i find the mangnitued of H(w) in DB (frequency repsonce)
%  altough this was not rasked for in the project handout
%
%  plot the graphs, and determine the pass-band ripple and stop band
%  attenuations from the graph (see report)
%
%  

%
%Choose DFT points to bee 512, since we want 256 resolution
%and H(W) is of periord 2 PI, so 256 is what we want for 0..Pi plots
% 
DFT=512;

if WINDOW==1
   %rectangular
   fprintf('building the rectangular widnow ...\n\n');

   %
   %build the rectangular window
   %
   for n=1:1:N
       hw(n)=1;
   end

elseif WINDOW==2
   %hanning
   fprintf('building the hanning window ...\n\n');

   %
   % Notice slight difference than in book, i divid by N not
   % N-1 since MATLAB starts arrays from 1 not 0
   % while book assumed h(n) to start from 0
   %
   for n=0:1:N-1
       ARG= 2*pi*n/N;
       hw(n+1)= .5*(1-cos(ARG));
   end
   
elseif WINDOW==3
   %hamming
   fprintf('building the hamming window ...\n\n');

   for n=0:1:N-1
       ARG= 2*pi*n/N;
       hw(n+1)= .54- .46*cos(ARG);
   end
   
elseif WINDOW==4
   %blackman
   fprintf('building the blackman window ...\n\n');

   for n=0:1:N-1
       ARG=  2*pi*n/N;
       ARG1= 2*ARG;  
       hw(n+1)= .42 - .5*cos(ARG) + .08*cos(ARG1);
   end
end

fprintf('Window signal hw(n) has been constructed ... \n\n');

%
% For each filter , build its desired h(n) (the nalytical derivation
% was dont by hand, see report
%

%
% evaluate the TAO argument once, it will be same and used by all
% filter below
% 
TAO = (N-1)/2

if FLTR==1
   % 
   %   B A N D       P A S S
   %
   %
 
   %
   % As extra, I evaluate the Desired H(W) (which is given in handout)
   % so that I may plot it if I want, although this is not need for
   % the main work to follow
   %
   k=0;
   for w=0:pi/256:pi
       k=k+1;
       if (w>0 & w<W1) | (w>W2 & w<pi)
          HD(k)=0;
       else
          HD(k)= exp(-j*w*TAO);
       end
   end
                
   %   
   % build the desired h, this is the inverse fourie transform from the
   % desired H(W) given in the problem, desired h has allready been
   % calcluated analytically by hand befor as: (see report for derivation)
   %
   %
   %        1  (  SIN W2(n-TAO)     SIN W1(n- TAO)     )
   % h(n)= --- ( --------------  - ------------------  )
   %  d    PI  (  n - TAO            n- TAO            )
   %
   
   fprintf('Building the desiered h(n) signal for BAND PASS ..\n\n');
   
   for n=0:1:N-1

       %
       %
       % Notice to SKIP OVER THE n=(N-1)/2  CASE
       %
       
       if(n ~= TAO)
          ARG= n-TAO;
          temp1= sin(W1*ARG);
          temp2= sin(W2*ARG);
          hd(n+1)= (temp2 - temp1) /(ARG*pi);
       else
          %
          % See report for how this evaluation of hd(n) when n=TAO
          % was done
          %
          hd(n+1)= W2/pi - W1/pi;
       end
   end
   
   fprintf('the desired hd(n) for BAND-PASS has been constructed ... \n\n');
elseif FLTR==2
   %
   %  B A N D -  R E J E C T
   %

   %
   % As extra, I evaluate the Desired H(W) (which is given in handout)
   % so that I may plot it if I want, although this is not need for
   % the main work to follow
   %

   k=0;
   for w=0:pi/256:pi
       k=k+1;
       if (w>W1 & w<W2)
          HD(k)=0;
       else
          HD(k)= exp(-j*w*TAO);
       end
   end
   hd_temp=ifft(HD,512);
   
   %build the desired h, this is the inverse fourie transform from the
   %desired H(W) given in the problem, desired h has allready been
   %calcluated analytically by hand befor as: (see report for derivation)
   %
   %
   %        1  (  SIN W1(n-TAO)     SIN W2(n- TAO)     )
   % h(n)= --- ( --------------  - ------------------  )
   %  d    PI  (  n - TAO            n- TAO            )
   %
   fprintf('Building the desiered h(n) signal for BAND REJECT ..\n\n');
   
   for n=0:1:N-1
       %
       %
       if(n ~= TAO)
          TAO
          ARG= n-TAO;
          temp1= sin(W1*ARG);
          temp2= sin(W2*ARG);
          hd(n+1)= (temp1 - temp2) /(ARG*pi);
       else
          %
          %                                                /       jwn
          % This was calculated by evaluating the integral | H(w) e
          %                                               /
          % when n=TAO
          %
          hd(n+1)= W1/pi - W2/pi +1;
       end
   end
   fprintf('the desired hd(n) for BAND-REJECT has been constructed ... \n\n');
 
elseif FLTR == 3
   %
   %  H I G H   P A S S 
   %

   %
   % As extra, I evaluate the Desired H(W) (which is given in handout)
   % so that I may plot it if I want, although this is not need for
   % the main work to follow
   %
   
   k=0;
   for w=0:pi/256:pi
       k=k+1;
       if (w>0 & w<W2) 
          HD(k)=0;
       else
          HD(k)= exp(-j*w*TAO);
       end
   end
             
   %build the desired h, this is the inverse fourie transform from the
   %desired H(W) given in the problem, desired h has allready been
   %calcluated analytically by hand befor as: (see report for derivation)
   %
   %
   %        1  (     SIN W2(n-TAO)  )
   % h(n)= --- ( -  --------------  )
   %  d    PI  (       n - TAO      )
   %
   fprintf('Building the desiered h(n) signal for HIGH-PASS..\n\n');
   
   for n=0:1:N-1
       %
       %
       %SKIP OVER THE n=(N-1)/2  CASE
       %
       if(n ~= TAO)
          ARG= n-TAO;
          temp2= sin(W2*ARG);
          hd(n+1)= - temp2 /(ARG*pi);
       else
          %
          % See report to see why this is 1 at mid-point
          %
          hd(n+1)= 1;
       end
   end
   
   fprintf('the desired hd(n) for HIGH-PASS has been constructed ... \n\n');
   
elseif  FLTR==4
   %
   % T E S T   C A S E   F R O M    B O O K
   %
   %
   %
   % This case is not part of the project, I have here to check
   % my code works ok, by testing a solved problem from text book
   % and comparing the outputs
   %

   %
   %build the desired h, this is given in 8.1.20 in book
   %
   %
   %        1  (     SIN Wc(n-TAO)  )
   % h(n)= --- (    --------------  )
   %  d    PI  (       n - TAO      )
   %
   Wc= 0.2*pi;  % from graph in book , page 585
   
   fprintf('Building the desiered h(n) signal for TEST case..\n\n');
   
   for n=0:1:N-1
       %
       %
       %SKIP OVER THE n=(N-1)/2  CASE
       %
       if(n ~= TAO)
          ARG= n-TAO;
          temp2= sin(Wc*ARG);
          hd(n+1)=  temp2 /(ARG*pi);
       else
          hd(n+1)= Wc/pi;
       end
   end
   
   fprintf('the desired hd(n) for TEST has been constructed ... \n\n');
end
  

   %
   % multiply both signals in time domain, the desired h, and the h window
   %
   h=hd(1:N) .* hw;

   %
   % find the discrete fourier transform of h
   %
   
    fprintf('Doing the FFT ... \n\n');

    H=fft(h,DFT);

    fprintf('FFT calculations complete, H(w) calculated.. \n\n');
 
   %
   % Take half of the H(W) , since we want 1/256 resolution
   %   
    H=H(1:DFT/2);   % because symmetry of H(w) take only half points
   
   %
   % construct the line scaling (x-axis) for the plots
   % we want pi/256 resolution (that is why also we used 512 FFT)
   %
   line=linspace(0,pi,DFT/2);    
   
   %
   %
   %  P H A S E   C A L C U L A T I O N S
   %
   %
   
   fprintf(' finding the PHASE ... \n\n');

   for n=1:1:DFT/2

       %
       % Get the x,y coorinates of the complex number H(W)
       %
       
       x=real(H(n));
       y=imag(H(n));
       
       %        
       %  adjust according to quardrent
       %
       if x >= 0
          if y >= 0
             % first quadrent
             if x==0
                PHASE(n) = 0;
             else
                ANGLE_R= atan(abs(y)/abs(x));
                PHASE(n) = ANGLE_R;
             end
          else
             % fourth quadrent
             if x==0
                PHASE(n) = -pi/2;
             else
               ANGLE_R= atan(abs(y)/abs(x));
               PHASE(n)= -ANGLE_R;
             end
          end;
       else
          if y >= 0
             % second quadrent
             if y==0
                PHASE(n)= pi;  
             else
                ANGLE_R= atan(abs(x)/abs(y)); 
                PHASE(n)= (pi/2)+ANGLE_R;
             end 
          else
             ANGLE_R= atan(abs(x)/abs(y));
             PHASE(n)= -( (pi/2)+ANGLE_R);
          end
       end
   end
   
   fprintf('Phase calculations completed .. \n\n');
   
   %
   % Now Get the magnitude of H
   %

   H_ABS=abs(H);
   H_MAX=max(H_ABS);
   
   %
   % N O R M A L I Z E
   %
   % Note: we normalize magnitude
   %
   
   H_ABS=H_ABS/H_MAX;
   
   %
   % find frequency response
   % This is extra work, I do it because I want to plot db of |H(W)|
   %
   fprintf('Finding the Frequency responce in db.. \n\n');
   
   for n=1:1:DFT/2
       Hdb(n)= 20*log10(H_ABS(n));
   end
   
   %
   % Finally, Normalize the h(n)
   %
   h=h/max(h);  


   %
   %   E N D    O F   P R O J E C T   C O D E
   %      
end