er_fsample.m

dsp工具箱

function f = er_fsample(M,fpts,varargin)
% function f = er_fsample(M,fpts,varargin)
%   Computes Type-I FIR filter coefficients using frequency sampling approach.
%   Takes filter parameter in M (for filter length N = 2M+1, e.g., M = 15
%   for length 31 filter) and desired normalized bandedge frequency points in 
%   FPTS (max length = 2). These points correspond to the midpoint between 
%   pass/stop-bandedge frequency points.
%
%   If the argument in VARARGIN is specified 'smooth', a transition band of one
%   sample is added to the desired frequency response, with magnitude 0.5.
%
%   Author:  Evan Ruzanski, CU-Boulder, ECEN5632 MATLAB assignment, FA2004

% Check desired response
s = length(fpts);
if s == 1
    type = 0;
else
    type = 1;
end
% Set filter length
N = 2*M + 1;
% Calculate group delay
g = (N - 1)/2;
% Set exponential parameter
z = -j*g*2*pi/N;

if type == 0 % Lowpass, generate Hd values
    for k = 0:N-1
        w = 2*pi*k/N;
        if ((w >= fpts*pi) & (w <= (2 - fpts)*pi)) 
            Hd(k+1) = 0;
        else 
            Hd(k+1) = 1;
        end
        if (nargin == 3)
            if upper(varargin{1}) == 'SMOOTH'
                Hd(ceil(fpts*M)+1) = 0.5;
                Hd(N+1-ceil(fpts*M)) = 0.5;
            end
        end
        H(k+1) = Hd(k+1)*exp(z*k);
    end
end
if type == 1 % Bandpass, generate Hd values
    for k = 0:N-1
        w = 2*pi*k/N;
        if ((w >= fpts(1)*pi) & (w <= fpts(2)*pi))
            Hd(k+1) = 1;
        elseif ((w >= (2 - fpts(2))*pi) & (w <= (2 - fpts(1))*pi)) 
            Hd(k+1) = 1;
        else 
            Hd(k+1) = 0;
        end
    end
    if (nargin == 3)
            if upper(varargin{1}) == 'SMOOTH'
                Hd2 = findedge(Hd);
                Hd(Hd2) = 0.5;
            end
    end
    for k = 0:N-1
        H(k+1) = Hd(k+1)*exp(z*k);
    end
end

disp('Filter Parameters are: ')
b = ifft(H)

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

% Compute frequency response
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 -40 5]); grid on

%%%%%%%%%% DECLARE LOCAL FUNCTION %%%%%%%%%%
function a = findedge(b)
% FINDEDGE Finds one-to-zero transitions between sets of binary values
L = length(b);
a = [];
ctr = 1;
for n = 1:L
    if (n ~= L)
        if (b(n) == 1) & (b(n+1) == 0)
            a(ctr) = n + 1;
            ctr = ctr + 1;
        end
        if (b(n) == 0) & (b(n+1) == 1)
            a(ctr) = n;
            ctr = ctr + 1;
        end
    end
end