我正在尝试实现任意 IIR 滤波器设计。我正在将所有算法分解成可以测试的较小部分。

我目前正在根据任意响应设计 FIR 滤波器。

Matlab的主要功能在这里：

function h = lslevin (N ,om, D, W)
% Complex Least Squares FIR filter design using Levinson' s algorithm
%
% h filter impulse response
% N filter length
% om frequency grid (0 <= om <= pi)
% D complex desired frequenccy response on the grid om
% W positive weighting function on the grid om
om = om(:);
D = D(:); 
W = W(:); 
L = length(om);
DR = real(D) ;
DI = imag(D);
 a = zeros(N,1); 
 b = a;
% Set up vectors for quadratic objective function
% (avoid building matrices)
dvec = D;
evec = ones(L, 1); 
e1 = exp(j*om);
for i=1:N,
    a(i) = W.'*real (evec);
    b(i) = W.'*real (dvec);
    evec = evec.*e1; 
    dvec = dvec.*e1;
end
a=a/L;
b=b/L;
% Compute weighted l2 solution
h = levin(a, b) ;

和莱文功能：

function x = levin(a ,b)
% function x = levin(a ,b)
% solves system of complex linear equations toeplitz (a)*x = b
% using Levinson's algorithm
% a ... first row of positive definite Hermitian Toeplitz matrix
% b ... right hand side vector
%
% Author: Mathias C. Lang, Vienna University of Technology, AUSTRIA
% 9-97
a = a(:);
b = b(:); 
n = length(a);
t = 1; 
alpha = a(1); 
x = b(1) / a(1);
for i = 1:n-1,
    k = -(a(i+1 :-1: 2)' * t)/alpha;
    t = [t;0] + k*flipud( [conj(t);0]);
    alpha = alpha*( 1 - abs(k)^2);
    k = (b(i+1) - a(i+1:-1: 2)' * x) /alpha;
    x = [x;0] + k*flipud(conj(t));
end

可以在这里完成一个示例函数（它应该是一个简单的带通）：

om1=pi*[linspace(0,0.23,230),linspace(0.3,0.5,200),linspace(0.57,1,430)];
D1=[zeros(1,230), exp(-j*om1(231:430) *2), zeros(1,430)]; % zeros are bandstop others bandpass
W1=[10* ones(1,230), ones(1,200), 10 *ones(1,430) ];%weight function
h = lslevin(61,om1, D1, W1) ;
wavwrite(h, 44100,  'e:\tmp2.wav'); % write the impulse response

但是当我计算脉冲响应时，我有类似的东西：

这与原始设计相差甚远

你能帮我找出代码中的错误吗？

多谢

杰夫

EDIT1：这是来自下面给出的答案的 IR 的最终 FFT：