Skip to content

Instantly share code, notes, and snippets.

@mattdsp
Created March 11, 2023 11:19
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mattdsp/dde270c64ace9d7bcad30333c6b9c90a to your computer and use it in GitHub Desktop.
Save mattdsp/dde270c64ace9d7bcad30333c6b9c90a to your computer and use it in GitHub Desktop.
Frequency sampling design of linear phase FIR filters
function h = fsamp(mag,N)
% function h = fsamp(mag,N)
% frequency sampling design of linear phase FIR filter
%
% mag ... desired magnitude response on an equidistant grid
% between 0 and Nyquist
% N ... filter length (N <= 2*length(mag) - 1 )
%
% for even N the frequency grid is defined by (M = length(mag)):
%
% f_k = k/(M-1), k=0,1,...,M-1 (includes Nyquist = 1)
%
% for odd N the grid is f_k = 2/(2*M-1), k=0,1,...,M-1 (Nyquist is not included!)
%
% M. Lang, mattsdspblog@gmail.com
M = length(mag);
Neven = (rem(N,2) == 0);
if Neven,
n = 2*(M-1); % FFT length
if ( N > n ),
error('fsamp: filter length N must not be greater than 2*(length(mag)-1).');
end
w = (0:M-1)'/(M-1)*pi; % frequency grid
D = mag(:).*exp(-1i*w*(n-1)/2); % desired complex frequency response
D = [D;conj(D(M-1:-1:2))]; % extend to 2*Nyquist
else
n = 2*M-1;
if ( N > n ),
error('fsamp: filter length N must not be greater than 2*(length(mag)-1).');
end
w = (0:M-1)'*2*pi/(2*M-1);
D = mag(:).*exp(-1i*w*(n-1)/2);
D = [D;conj(D(M:-1:2))];
end
% impulse response via IFFT
h = real(ifft(D));
% windowing (Hamming)
I = (n-N)/2+1;
ww = hamming(N);
h = h(I:I+N-1).*ww(:);
return
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment