efilter.m

里面囊括了基于matlab滤波器设计的各种.m文件

% H = EFILTER(ALFA, STOPEDGE, PASSEDGE, ORDER)
%
% This program is to realize the Lowpass filter Algorithm of the paper
% 'Eigenfilters: A New Approach to Least-Squares FIR Filter Design and 
% Applications Including Nyquist Filters' IEEE TRANS. ON C&S VOL CAS-34,
% NO.1 JAN. 1987.
%
% Input arguments:
%   Alfa: The alfa that is in the range [0,1], controls the relative
%         accuracies of approximation in the pass and stopband.
%   Stopedge: The stopband edge of the lowpass filter
%   Passedge: The passband edge of the lowpass filter
%   Order : The order of the FIR lowpass filter
%
% Output H is the filter coefficients
%
% 

function h = efilter(alfa, stopedge, passedge, order);

s = stopedge;
p = passedge;

if mod(order,2) == 0
  M = order/2;
  for i = 1:M,
   for j = 1:M,
     if i == j 
       P(i,j) = (1-alfa)/pi*(1/i*(1.5*i*p-2*sin(i*p)+...
                 0.5*cos(i*p)*sin(i*p))) + alfa/pi*(1/(2*i))*...
                 (i*pi-cos(i*s)*sin(i*s)-i*s);
     else
       P(i,j) = (1-alfa)/pi*(p-1/j*sin(j*p)-1/i*sin(i*p)+...
                 0.5*(1/(i-j))*sin(i*p-j*p)+0.5*(1/(i+j))*...
                 sin(i*p+j*p)) - 0.5*alfa/pi*((1/(i-j))*...
                 sin(i*s-j*s)+(1/(i+j))*sin(i*s+j*s));   
     end
    end
  end

else
  
 M = (order+1)/2;
  for i = 1:M
    for j = 1:M
      if i == j
        P(i,j) = (1-alfa)/pi*(1/(i+0.5)*(1.5*(i+0.5)*p-2*sin((i+0.5)*p)+...
                  0.5*cos((i+0.5)*p)*sin((i+0.5)*p)))+alfa/pi*(1/(2*(i+0.5))*...
                  ((i+0.5)*pi-cos((i+0.5)*s)*sin((i+0.5)*s)-(i+0.5)*s));
      else
        P(i,j) = (1-alfa)/pi*(p-1/(j+0.5)*sin((j+0.5)*p)-1/(i+0.5)*sin((i+0.5)*p)+...
                  0.5*(1/(i-j))*sin(i*p-j*p)+0.5*(1/(i+j+1))*...
                  sin(i*p+j*p+p)) - 0.5*alfa/pi*((1/(i-j))*...
                  sin(i*s-j*s)+(1/(i+j+1))*sin(i*s+j*s+s));   
      end
    end
  end
end
[V,D] = eig(P);
for i = 1:M,
 ev(i) = D(i,i);
end
[eigv,index] = min(ev);
b = V(:,index);
b = b';
if mod(order,2) == 0
 h = [0.5*b(2:M)*antieye(M-1),b(1),0.5*b(2:M)];
else
 h = [0.5*b*antieye(M),0.5*b];
end
[H,w] = freqz(h,1,1000);
semilogy(w,abs(H));