2
$\begingroup$

i'm doing overlap save method in frequency domain. i make hankel function to digital filter using frequency sampling method(to make arbitrary magnitude phase filter). and i do ifft and zero padding fft

i heard that if the filter is non-causal, zero padding fft have gibb's phenomenon. my filter have same problem. i'm doing shift FIR filter to make causal filter but the problem is not solved.... please help me



    clear all;
    close all;

    %% parameter
    
    fs = 2000;
    ts = 1/fs;
    c = 340;
    R = 1;
    
    %% hankel filter     
    Lp = 50;
    N = 2*Lp+1;                     %filter length
    resol = 0:N-1;                  % frequency resolution
    f = resol/N*fs;
    k = 2*pi*f/c;
    M=5;
    R = 2;
    nu = 2;                   %hankel order
    
    
    % frequency sampling method

    right = besselh(nu,k(1:(N-1)/2 +1)*R); %Hankel function 0~pi
    right(1:5) = 0;                     %eliminate infinite part
    left = fliplr(right);              
    left = conj(left);                  
    left(end) = [];
    
    H = [right left];               % symmetric freuqnecy response
    h = ifft(H);                % impulse response   
    
    nn = 80;                        %shift index 
    h1 = [h(nn : end) h(1:nn-1)];   % shift impulse response to make causal filter
    H1 = fft(h1,2*N);               %zero padding fft

    
    subplot(2,1,1)
    plot(h1)
    subplot(2,1,2)
    plot(abs(H1))

$\endgroup$
2
  • $\begingroup$ Your H is not conjugate symmetric so that h is not real, please check your code. To eliminate the ininite part, why do you set right(1:5) to zero? $H_\nu(z)$ goes to infinity only when $z=0$ due to the property of Neumann function. $\endgroup$
    – ZR Han
    Feb 7, 2022 at 4:59
  • $\begingroup$ i edit code Thank you. 1:5 is when frequency resolution is high first some part have very high value. So i do 1:5 $\endgroup$
    – gg h
    Feb 7, 2022 at 5:01

1 Answer 1

2
$\begingroup$

The Gibbs phenomenon is caused by frequency sampling method. I recommend you to use window method instead. It is fast, convenient and robust, although not optimal. BTW, the frequency range should be 0~fs/2, I've fixed it in the following code

clear;
close all;

%% parameter
fs = 2000;
ts = 1/fs;
c = 340;

%% hankel filter     
Lp = 50;
N = 2*Lp+1;                     % filter length
NN = 2^nextpow2(N * 100);       % NN >> N
resol = 0:NN-1;                  % frequency resolution
f = resol/NN*fs/2;              % you miss /2 here
k = 2*pi*f/c;
M = 5;
R = 2;
nu = 2;                   %hankel order

% window method
right = besselh(nu,k(1:(NN-1)/2 +1)*R); %Hankel function 0~pi
right(1:200) = 1;                     %eliminate infinite part
left = fliplr(right);              
left = conj(left);                  
left(end) = [];

H = [right left];               % symmetric freuqnecy response
h = ifft(H);                % impulse response   

h1 = ifftshift(h);   % shift impulse response to make causal filter
win = hann(N).';
h2 = h1(NN/2-N/2:NN/2+N/2-1) .* win; % windowed impulse response
H2 = fft(h2,2*N);               %zero padding fft

figure;
subplot(2,1,1)
plot(h2)
subplot(2,1,2)
plot(abs(H2))

Frequency sampling method:

Frequency Sampling

Window method:

Window

$\endgroup$
7
  • $\begingroup$ OMG............Very Very Very Thank you so much $\endgroup$
    – gg h
    Feb 7, 2022 at 5:37
  • $\begingroup$ if error is large, is this normal? $\endgroup$
    – gg h
    Feb 7, 2022 at 6:00
  • $\begingroup$ @ggh I don't see significant error between these two methods. $\endgroup$
    – ZR Han
    Feb 7, 2022 at 6:10
  • $\begingroup$ there is some wrong part. i fixed. Thank you very much $\endgroup$
    – gg h
    Feb 7, 2022 at 6:33
  • $\begingroup$ Could i ask one more question?? How can i make non-causal filter? that didn't use ifftshift command?? $\endgroup$
    – gg h
    Feb 8, 2022 at 7:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.