function f = er_fwindow(M,wc,option)
% function f = er_fwindow(M,wc,option)
%   Computes FIR filter coefficients using frequency sampling approach.
%   Takes filter parameter in M (for filter length N = 2M + 1) and desired 
%   normalized bandedge frequency points in WC (ordered low to high, max 
%   length = 2).  These points correspond to the midpoint between respective
%   pass/stop-bandedge frequency points.  Window type is specified by OPTION:
%
%       0 = rectangular window
%       1 = raised cosine window
%
%   Author:  Evan Ruzanski, CU-Boulder, ECEN5632 MATLAB assignment, FA2004

% Set indices
n = -M:M;
N = length(n);

% Check desired response
s = length(wc);
if s == 1
    type = 0;
else
    type = 1;
end

% Make filter kernel
if type == 0
    num = wc*sinc(wc*n);
elseif type == 1 
    num = wc(2)*sinc(wc(2)*n) - wc(1)*sinc(wc(1)*n);
end

% Make window sequence
if option == 0
    % Rectangular window
    w = ones(1,2*M+1);
elseif option == 1
    % Raised cosine window
    w = 0.5 + 0.5*cos(2*pi*n/(2*M));
end

% Find filter coefficients
disp('Filter coefficients are: ')
b = num.*w

% Pack parameters into struct
f = struct('tf_complete',{b,1});

% Optional plotting
% Compute frequency response
% [FF,w] = freqz(f,1,512);
nnumd = max(size(N));
fftn = 1024; % Take default 2^10 pt. FFT
ss = 2; % For half of unit circle
omega = (0:fftn-1)'*2*pi/fftn/ss; % Set frequency vector
nfft = lcm(fftn,N);
FF = (fft([b zeros(1,ss*nfft - nnumd)])); % Perform FFT
FF = FF(1+(0:fftn-1)*nfft/fftn);
k = 0:N-1;
subplot(211)
stem(k,real(b),'filled');
xlabel('Time index n'); ylabel('Amplitude')
grid on
subplot(212)
plot(omega/pi,20*log10(abs(FF)));
xlabel('\omega/\pi'); ylabel('Gain, dB');
axis([0 1 -50 5]); grid on

%%%%%%%%%% Declare local functions %%%%%%%%%%
function y = sinc(x)
% Y = SINC(X) Creates sinc function
y = ones(size(x));
i = find(x);
y(i) = sin(pi*x(i))./(pi*x(i));