0
$\begingroup$

What steps are necessary to get the same impulse as before a FFT + IFFT.

At the moment I am trying out my code a bit and would like to do an FFT from an impulse response followed by an IFFT to generate the impulse response again. Actually I would like to change the frequency range a bit and look at the impulse response from that. However, I am already failing to get the same IR without changing the frequency range.

The impulse I generate looks like the original impulse, but the same impulse appears mirrored at the end.

What I tried before applying the IFFT:

  1. i set the result of the FFT from the half (to the Nyquist frequency) to 0.
  2. I mirrored the first half of the FFT result to the second half.

But the impulse response always has a mirrored impulse response at the end.

enter image description here

Perhaps the question is too trivial for many, however, I have not yet found the appropriate answer, perhaps not understood. Glad about any help.

Here also my Code:

void FFT::executeFFT(std::vector<double>* _data,std::vector<double>* _Intervall,std::string _Name)
    {

        if(_data !=nullptr|| real.empty() != true || imag.empty() != true){
        for (int i = 0; i < real.size(); i++)
        {
            real[i] = 0.0; imag[i] = 0.0;
        }
        
        int le, ip, m, i, nm1, j, t, le1, inv = 0;
        double ur, ui, tr, ti, tmpr, tmpi, nv2, k;
        if (_hin) inv = 1;
        else inv = -1;
        for (int i = 0; i <= Nmax - 2; i++)
        {
            real[i] = _data->at(i+1);
        }
        m = round((log(Nmax) / log(2)));
        nv2 = Nmax / 2;
        nm1 = Nmax - 1;
        j = 1;
        //Vertauschung der Eingangsfolge
        for (i = 1; i <= nm1; i++)
        {
            if (i < j) {
                tmpr = real[j - 1];
                tmpi = imag[j - 1];
                real[j - 1] = real[i - 1];
                imag[j - 1] = imag[i - 1];
                real[i - 1] = tmpr;
                imag[i - 1] = tmpi;
            }
            k = nv2;
            while (k < j)
            {
                j = round(j - k);
                k = k / 2.0;
            }
            j = round(j + k);
            
        }
        //Berrechnung der Butterflys
        for (int l = 1; l <= m; l = l + 1)
        {
            le = 1;
            for (t = 1; t <= l; t = t + 1)
            {
                le = le * 2;
            }
            le1 = round(le / 2);
            ur = 1;
            ui = 0;
            for (j = 1; j <= le1; j++)
            {
                i = j;
                while (i <= Nmax)
                {
                    ip = i + le1;
                    tr = real[ip - 1] * ur - imag[ip - 1] * ui;   //{ Berechnen des }
                    ti = imag[ip - 1] * ur + real[ip - 1] * ui;   //{ aktuellen     }
                    real[ip - 1] = real[i - 1] - tr;             //{ 'Butterflys'  }
                    imag[ip - 1] = imag[i - 1] - ti;
                    real[i - 1] = real[i - 1] + tr;
                    imag[i - 1] = imag[i - 1] + ti;
                    i = i + le;
                }
                ur = cos(M_PI * j / le1);
                ui = -inv * sin(M_PI * j / le1);
            }
        }
            if(!_hin){
                for (i = 0; i <= Nmax - 1; i++)
                {
                    real[i] = real[i] / (double)Nmax;
                    imag[i] = imag[i] / (double)Nmax;
                }
            }
        } else {std::cout<<"FFT Hat nen Nullptr";}
        //Betrag, Phase und Intervall
        double df = 0.0, f=0.0, fg = 48000.0;
        df = (fg / (float)Nmax);
        int index = 1;
   /* if(!_hin){
        //while (f < 22000){
        for(int index=0;index<real.size();index++){
            f = index * df;
            _Intervall->push_back(f);
            real[index] = Betrag(index,_hin);
            imag[index] = Phase(index);
            //index = index + 1;
            
        }
    }*/
       if (_hin) WriteFFT(_Name,&real,&imag,df);
    for (int u=1;u<real.size();u++)
    {
        _data->at(u)=real.at(u);
    }
    _data->at(0)=real.at(1);
}
  • Solution: *
void FFT::executeFFT(std::vector<double>* _real,std::vector<double>* _imag,std::vector<double>* _Intervall,std::string _Name)
    {
 if(_real !=nullptr|| real.empty() != true || imag.empty() != true){
            if(_hin){
                 for (int i = 0; i < real.size(); i++)
                 {
                     //real[i] = 0.0;
                     imag[i] = 0.0;
                     real[i]=_real->at(i);
                 }
            }else{ for (int i = 0; i < real.size(); i++)
                {
                    real[i] = _real->at(i);
                    imag[i] = _imag->at(i);
                }
            }
        int le, ip, m, i, nm1, j, t, le1, inv = 0;
        double ur, ui, tr, ti, tmpr, tmpi, nv2, k;
        if (_hin) inv = 1;
        else inv = -1;
        /*for (int i = 0; i <= Nmax - 2; i++)
        {
            real[i] = _data->at(i+1);
        }*/
        m = round((log(Nmax) / log(2)));
        nv2 = Nmax / 2;
        nm1 = Nmax - 1;
        j = 1;
        //Vertauschung der Eingangsfolge
        for (i = 1; i <= nm1; i++)
        {
            if (i < j) {
                tmpr = real[j - 1];
                tmpi = imag[j - 1];
                real[j - 1] = real[i - 1];
                imag[j - 1] = imag[i - 1];
                real[i - 1] = tmpr;
                imag[i - 1] = tmpi;
            }
            k = nv2;
            while (k < j)
            {
                j = round(j - k);
                k = k / 2.0;
            }
            j = round(j + k);
            
        }
        //Berrechnung der Butterflys
        for (int l = 1; l <= m; l = l + 1)
        {
            le = 1;
            for (t = 1; t <= l; t = t + 1)
            {
                le = le * 2;
            }
            le1 = round(le / 2);
            ur = 1;
            ui = 0;
            for (j = 1; j <= le1; j++)
            {
                i = j;
                while (i <= Nmax)
                {
                    ip = i + le1;
                    tr = real[ip - 1] * ur - imag[ip - 1] * ui;   //{ Berechnen des }
                    ti = imag[ip - 1] * ur + real[ip - 1] * ui;   //{ aktuellen     }
                    real[ip - 1] = real[i - 1] - tr;             //{ 'Butterflys'  }
                    imag[ip - 1] = imag[i - 1] - ti;
                    real[i - 1] = real[i - 1] + tr;
                    imag[i - 1] = imag[i - 1] + ti;
                    i = i + le;
                }
                ur = cos(M_PI * j / le1);
                ui = -inv * sin(M_PI * j / le1);
            }
        }
            if(!_hin){
                for (i = 0; i <= Nmax - 1; i++)
                {
                    real[i] = real[i] / (double)Nmax;
                    imag[i] = imag[i] / (double)Nmax;
                }
            }
        } else {std::cout<<"FFT Hat nen Nullptr";}
        //Betrag, Phase und Intervall
        double df = 0.0, f=0.0, fg = 48000.0;
        df = (fg / (float)Nmax);
        int index = 1;

    _real=&real;
    _imag=&imag;
    for (int u=1;u<real.size();u++) // u=0?
    {
        _real->at(u)=real.at(u);
        _imag->push_back(imag.at(u));
    }
Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ At first I would suggest to understand the theory about the periodic nature of DFT/IDFT. Trying random tweaks won't help. Afterwards, try to compare the output of each step with a third party FFT implementation (e.g. from MATLAB or Python). It's hard to tell what's going on by just looking at the output, without even knowing what your input impulse looked like. Also, I am confused by your code, especially the last lines. You copy the real part of the FFT result back to the data vector, but what happens with the imaginary part? $\endgroup$
    – Vito
    Nov 1, 2022 at 19:07
  • 1
    $\begingroup$ Did you try just taking the output of the FFT and feeding it to the IFFT should, in theory, return exactly the same result. If you have a math package that uses consistent scaling of the FFT and IFFT results, then xp = ifft(fft(x)) should hand you back an xp that's almost equal to x except for small differences due to rounding error within the FFT. $\endgroup$
    – TimWescott
    Nov 1, 2022 at 19:42
  • $\begingroup$ Ahh @TimWescott. Thats nice to know. Then I should have a look into my code again. $\endgroup$
    – Hannes Broschk
    Nov 1, 2022 at 19:54
  • $\begingroup$ And @Vito: Normally I have calculated the phase in the block I have excluded, but I don't need the imaginary part for the Ifft, do I? Therefore I had omitted this. $\endgroup$
    – Hannes Broschk
    Nov 1, 2022 at 19:54
  • $\begingroup$ Second half would be the complex conjugate of the mirror of the first half, IIRC - that means making the imaginary part its negative. $\endgroup$
    – user253751
    Nov 1, 2022 at 20:03

1 Answer 1

Reset to default
5
$\begingroup$

However, I am already failing to get the same IR without changing the frequency range.

Then there is something obviously wrong with your code. I suggest standard debugging procedures: Start with a single impulse, delayed impulse, double impulse, sine wave, dual sine wave etc. I.e. cases there the expected result can be easily calculated by hand and then find the spot where it goes wrong and and single step through it.

Once you have the simple test cases running you can write a comprehensive unit test.

Actually I would like to change the frequency range a bit and look at the impulse response from that

Direct manipulation of FFT coefficients in the frequency domain is mathematically fairly complicated. Unless you are comfortable with concepts like periodicity, circular convolution, spectral leakage, orthogonal bases, time domain aliasing, overlap-add/save etc. I would NOT go there and just stick with "standard" time domain filter algorithms like FIR or IIR.

but I don't need the imaginary part for the Ifft, do I?

You sure do. If you really don't know why, I strongly suggest spending a bit of quality time with text book or on-line class on the DFT before trying to write code. Most people's "intuition" on what the DFT does and how it works is plain wrong.

Improve this answer
$\endgroup$
1
  • $\begingroup$ You were right, @Hilmar. I thought I could set the imaginary part to 0 when returning to the IFFT. But I must have made a mistake there. If you don't set the imaginary part to 0 at the beginning, the code works and I get the desired result. $\endgroup$
    – Hannes Broschk
    Nov 2, 2022 at 13:31

Your Answer

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

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