2
$\begingroup$

I've been experimenting with the frequency sampling method for designing FIR filters. I created a low-pass filter that has a linear transition band. However, when I performed an IFFT on the frequency response samples, the resulting impulse response samples were complex and conjugate numbers, instead of a real value array. I am considering the possibility of IFFT round-off errors due to the small imaginary part represented as 0.000000000000000i or sometimes with an order of e-18 in MATLAB.

Since I am not experienced in filter design, I cannot rule out the possibility that the problem lies with the program rather than the round-off errors so I'm asking here. I'm wondering if anyone has any insights on why this might be happening and how I can fix it. Thanks in advance!

clear;
clc;

%{  
    FIR filter Type 1
    Filter periodic frequency response (needs to be sampled)

                     B                     fc-B
   1 |________________                       ________|
     |                \                     /        |
...  |                 \                   /         |   ...
     |                  \                 /          |
 __0_|___________________\_______________/___________|____ freq axis
     |                  B+WB        fc-(B+WB)        |
     0                                              fc
%}

clear;
clc;

fc = 1000;    % sampling freq [Hz]
T = 1/fc;
N = 300;
B = 200;      % monolater bandpass
df = 1/((N+1)*T);    % frequency sampling step
W = 0.1;

freq_samples = zeros(1,N+1);

% frequency sampling
for i=1:N+1
    if (i-1)*df <= B || (i-1)*df >= (fc - B)
        freq_samples(i) = 1;
    elseif (i-1)*df > B && (i-1)*df<=(B+W*B)
        freq_samples(i) = ((((i-1)*df)-B)/((B+W*B)-B))*(-1)+1;
    elseif (i-1)*df >= (fc-B-W*B) && (i-1)*df < (fc-B)
        freq_samples(i) = ((((i-1)*df)-(fc-B-W*B))/((fc-B)-(fc-B-W*B)));
    else
        freq_samples(i) = 0;
    end
end

%stem(freq_samples);
time_samples = ifft(freq_samples);
h = [time_samples(N/2+2:N+1) time_samples(1:N/2+1)];    % impulse response
fvtool(h,1);

N.B. I tried to call the IFFT() function with the 'symmetric' argument and I was able to fix the problem, however, I am interested in understanding why the problem arose and how to handle similar situations. I believe ignoring the imaginary part just because "it works" may not be the optimal solution and, therefore, would like to gain a better understanding of the issue.

EDIT These frequency response samples are for N = 20. enter image description here

EDIT 2 I used a MATLAB script to check for symmetry and found that the sampled frequency response is not symmetric.

format LONG E
test = [freq_samples(N/2+2:N+1) freq_samples(1:N/2+1)];
for i=1:N/2+1
    if test(i) ~= test(N+2-i)
        disp(['index: ' num2str(i) ', value: ']);
        disp(test(i))
        disp(['index: ' num2str(N+2-i) ', value: ']);
        disp(test(N+2-i))
        disp("-----------------------");
    end
end

Running this piece of code with N = 300, fc = 1000, B = 200, W = 0.1 I get this output:

index: 85, value: 
     3.654485049833624e-02

index: 217, value: 
     3.654485049833911e-02

-----------------------
index: 86, value: 
     2.026578073089695e-01

index: 216, value: 
     2.026578073089695e-01

-----------------------
index: 87, value: 
     3.687707641196027e-01

index: 215, value: 
     3.687707641196013e-01

-----------------------
index: 88, value: 
     5.348837209302303e-01

index: 214, value: 
     5.348837209302332e-01

-----------------------
index: 90, value: 
     8.671096345514968e-01

index: 212, value: 
     8.671096345514954e-01

-----------------------

So the symmetric index values are not strictly equal, but just very very close. Is this why IFFT() returns a complex impulse response?

$\endgroup$
14
  • 2
    $\begingroup$ It shouldn't be surprizing that a complex operation (and I mean "complex" in both manner of meaning) like the iFFT will result in a tiny imaginary part for a number that, if you had infinite precision, was meant to be purely real because of the finite numerical precision. As long as you're confident that the data going into the iFFT is Hermitian symmetric (that is: $$ H[N-k] = H^*[k] $$ ) then you know you should have a purely real result and you can and should just zero them imaginary mutherfuckers that ain't exactly zero. $\endgroup$ Commented May 2 at 21:52
  • 1
    $\begingroup$ I am good at MATLAB. Good enough to know what I don't like. That off-by-1 crap is the main thing. Frankly, I don't like your code. I don't really like an odd-sized FFT (301), but MATLAB should handle it. The way to insure that it's symmetric is first fill the lower half with the frequency spectrum you want. Then, except for DC, copy the lower half to the upper half, except mirror image it and complex conjugate it. But still because of numerical roundoff error, the imaginary parts of your result will not be exactly zero. It's because lotsa math and finite word size.. $\endgroup$ Commented May 2 at 23:35
  • 1
    $\begingroup$ In the context of analysis like what we do with MATLAB, there really isn't anything wrong with a non-integer power of 2. But if I'm using a FFT in an audio or music analysis/processing device, I will worry about MIPS and execution time, and I am most comfortable with $2^p$ FFT length. For fast-convolution, there is little to be gained by some other FFT length and sometimes the complexity gets worse. In the olden daze, MATLAB zero-padded up to the next power of 2, then after the FFT did sinc-based "bandlimited" interpolation between the bins. $\endgroup$ Commented May 4 at 0:54
  • 1
    $\begingroup$ One thing useful about odd-length FFT is that there is no Nyquist bin to worry about. There's DC and an equal number of positive and negative frequency bins. You have independent control of DC and the magnitude and phase of every frequency component. If there is a Nyquist bit (even-length FFT), then you can't independently have control over both phase and magnitude. It really has to be zero phase. $\endgroup$ Commented May 4 at 0:57
  • 1
    $\begingroup$ If you construct the input to the IFFT properly such that it is perfectly symmetric, MATLAB will recognize this and force the output to be real-valued (a complex-to-real IFFT is cheaper than a complex-to-complex one). I guess yours is not exactly symmetric due to rounding errors in your computation. $\endgroup$ Commented May 6 at 13:06

2 Answers 2

4
$\begingroup$

Just repeating what the comments say: The residual imaginary part of the impulse response is simply numerical noise. Matlab uses double precision floating point which has finite precicions

As simple example is

>> sqrt(2)*sqrt(2)-2
ans =
   4.4409e-16

The answer should be 0 but the limited precision introduces small errors in the calculations.

Note that they are not symmetrical.

Of course it's symmetrical It follows directly from the definition of the Discrete Fourier Transform that the signals are periodic with $N$ in both domains, i.e.

$$x[n] = x[n + mN],\quad X[k] = X[k + mN], \qquad \qquad \forall n, k, m \in \mathbb{Z}$$

Symmetric in this context means "circular symmetry" since $X[-k] = X[N-k]$. You can easily check this by flipping your frequency samples around $k = 0$, conjugating it and checking whether it's the same as the original.

freq_samples_flipped = conj([freq_samples(1) flip(freq_samples(2:end))]);
err =  freq_samples_flipped - freq_samples;
fprintf('Symmetry error = %6.2fdB\n', 10*log10(sum(abs(err).^2)/sum(abs(freq_samples.^2))));
$\endgroup$
2
$\begingroup$

Here is a dumb MATLAB program I wrote nearly 3 decades ago (and modified a little since with the 'var' thing for exist() and changing wavread() to audioread()) just to do a sanity check with my use of the fft() and fftshift() and such. When I'm doing an experiment to check out something, I often start with a copy of this file and then modify it to do something more specific.

    if ~exist('inputFile', 'var')
        inputFile = 'vibe.wav';
    end

    [inputBuffer, Fs] = audioread(inputFile);

    fileSize = length(inputBuffer);

    numSamples = 2.^(ceil(log2(fileSize))); % round up to nearest power of 2

    x = zeros(numSamples, 1);                   % zero pad if necessary

    x(1:fileSize) = inputBuffer(:,1);           % if multi-channel, use left channel only

    clear inputBuffer;                          % free this memory
    clear fileSize;

    t = linspace(0, (numSamples-1)/Fs, numSamples)';
    f = linspace(-Fs/2, Fs/2 - Fs/numSamples, numSamples)';

    X = fft(x);

    plot(t, x);
    xlabel('time (seconds)');
    ylabel('amplitude');
    title(['time-domain plot of ' inputFile]);
    sound(x, Fs);                                           % play the sound
    pause;




    % display both positive and negative frequency spectrum

    plot(f, real(fftshift(X)));
    xlabel('frequency (Hz)');
    ylabel('real part');
    title(['real part frequency-domain plot of ' inputFile]);
    pause;

    plot(f, imag(fftshift(X)));
    xlabel('frequency (Hz)');
    ylabel('imag part');
    title(['imag part frequency-domain plot of ' inputFile]);
    pause;

    plot(f, abs(fftshift(X)));                              % linear amplitude by linear freq plot
    xlabel('frequency (Hz)');
    ylabel('amplitude');
    title(['abs frequency-domain plot of ' inputFile]);
    pause;

    plot(f, 20*log10(abs(fftshift(X))+1.0e-10));            % dB by linear freq plot
    xlabel('frequency (Hz)');
    ylabel('amplitude (dB)');
    title(['dB frequency-domain plot of ' inputFile]);
    pause;





    % display only positive frequency spectrum for log frequency scale

    semilogx(f(numSamples/2+2:numSamples), 20*log10(abs(X(2:numSamples/2))));       % dB by log freq plot
    xlabel('frequency (Hz), log scale');
    ylabel('amplitude (dB)');
    title(['dB vs. log freq, frequency-domain plot of ' inputFile]);
    pause;

    semilogx(f(numSamples/2+2:numSamples), (180/pi)*angle(X(2:numSamples/2)));      % phase by log freq plot
    xlabel('frequency (Hz), log scale');
    ylabel('phase (degrees)');
    title(['phase vs. log freq, frequency-domain plot of ' inputFile]);
    pause;

    %
    %   this is an alternate method of unwrapping phase
    %
    %   phase = cumsum([angle(X(1)); angle( X(2:numSamples/2) ./ X(1:numSamples/2-1) ) ]);  
    %   semilogx(f(numSamples/2+2:numSamples), phase(2:numSamples/2));                  % unwrapped phase by log freq plot
    %

    semilogx(f(numSamples/2+2:numSamples), unwrap(angle(X(2:numSamples/2))));       % unwrapped phase by log freq plot
    xlabel('frequency (Hz), log scale');
    ylabel('unwrapped phase (radians)');
    title(['unwrapped phase vs. log freq, frequency-domain plot of ' inputFile]);


```
$\endgroup$

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.