function [impw]=jakes_filter(fd)
n=512;nn=2*n;
fs=0:fd/64:fd;
H=zeros(1,n);
for k=1:(n/8+1)
    jpsd(k)=1/((1-((fs(k))/fd)^2+eps)^0.5);
    if(jpsd(k))>1000
        jpsd(k)=1000;
    end
    H(k)=jpsd(k)^0.5;
end
for k=1:n
    H(n+k)=H(n+1-k);
end
[inv,time]=linear_fft(H,nn,fd/64);
imp=real(inv(450:577));
impw=imp.*hanning(128)';
energy=sum(impw.^2);
impw=impw/(energy^0.5);