% -------------------------------------------------------------------------
%     *** 2.161 Signal Processing - Continuous and Discrete  ***
%                      Tutorial MATLAB Function
%
% firfs - A simple Frequency-Sampling Linear-Phase FIR filter based on DFT %         interpolation.
%
% Source:   Class handout: Frequency Sampling FIR Filters 
%
% Usage :    h = firfs(samples)
%            where samples -  is a row vector of N equi-spaced, real values
%                    of the freq. response magnitude.
%                    The samples are interpreted as being equally spaced 
%                    around the top half of the unit circle at normalized
%                    (in terms of the Nyquist frequency f_N) frequencies 
%                    from 0 to 2(N-1)/(2N-1) x f_N, or at frequencies
%                    2n/(2N-1)xf_N for n = 0...N-1.
%                    Because the length is odd, the frequency response is not
%                    specified at f_N.
%                  h - is the output impulse response of length 2N-1 (odd).
%           The filter h is real, and has linear phase, i.e. has symmetric
%           coefficients obeying  h(k) = h(2N+1-k), k = 1,2,...,N+1.
%
% Version:  1.0
% Author:   D. Rowell   10/6/07
% -------------------------------------------------------------------------
% 
function h = firfs(samples)
%
% Find the length of the input array...
% The complete sample set on the unit circle will be of length (2N-1)
%
M = length(samples);
H_d = zeros(1,2*M-1);
%
% We want a causal filter, so the resulting impulse response will be shifted
% (2N-2)/2 = N-1 to the right, giving a phase shift of exp(-j*Omega*4 
% Move the samples into the upper and lower halves of H_d and add the
% linear phase shift term to each sample.
%
Phi = 2*pi*(M-1)/(2*M-1);
H_d(1) = samples(1);
for j = 2:M
   Phase        = exp(-i*(j-1)*Phi);
   H_d(j)       = samples(j)*Phase;
   H_d(2*M+1-j) = conj(H_d(j));
end
%
% Use the inverse DFT to define the impulse response.   We know the impulse
% response will be real.   (Note that MATLAB's fft routines can handle any
% length vector.)
%
h = (ifft(H_d));
%