Old 09-08-2016, 03:11 AM   #1
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default JSFX real FFT (v5.25+)

I thought I would post some examples of the new fft_real() / ifft_real() JSFX functions.

Here is a pulse tone generator, that uses the FFT to convert a naive waveform into its Fourier coefficients (in @init). It then uses the inverse FFT to generate an alias-free waveform (in @slider). It seems to be about twice as fast as when doing the same using the complex FFT, while using only half the memory, plus you don't have to worry about mirroring the negative frequency bins.

Code:
desc:Real FFT tone generator
slider1:440<20,20000,1>Frequency (Hz)
slider2:-12.0<-120.0,6.0,0.1>Gain (dB)

@init

size = 1024;

buf = 0;
coef = buf + size;

// Render naive waveform.
scale = 0.5/size;
i = 0;
loop(size,
  naive_saw = 2*i / (size - 1) - 1;
  coef[i] = naive_saw * scale;
  i += 1;
);

// Real FFT i.e. real to complex.
fft_real(coef, size);
fft_permute(coef, size/2);

// Clear DC (optional).
coef[0] = 0;

// Clear Nyquist.
coef[1] = 0;

@slider

dt = slider1 / srate;

dt < 0.5 ? (
  // Copy coefficients up to Nyquist.
  n = min(2 * ceil(0.5 / dt), size);
  memcpy(buf, coef, n);

  // Clear coefficients beyond Nyquist.
  memset(buf + n, 0, size - n);

  // Inverse real FFT i.e. complex to real.
  fft_ipermute(buf, size/2);
  ifft_real(buf, size);
);

gain = exp(log(10)/20 * slider2);

@sample

dt >= 0.5 ? out = 0 : (
  // Linear interpolation.
  lerp = size * t;

  t += dt;
  t >= 1 ? t -= 1;

  i = lerp|0;
  lerp -= i;

  out = (1 - lerp) * buf[i] + lerp * buf[(i + 1) & (size - 1)];
  out *= gain;
);

spl0 += out;
spl1 += out;

Last edited by Tale; 03-09-2017 at 04:18 AM. Reason: Replaced weird 30% pulse with saw
Tale is offline   Reply With Quote
Old 09-08-2016, 04:35 AM   #2
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Here is another example showing how to do convolution using the real FFT:

Code:
desc:FIR filter using real FFT convolution

@init

size = 1024;

in = 0;
ir = in + size;
out = ir + size;

cutoff = min(1000 / srate, 0.5);
i = 0;
loop(size/2,
  // Low-pass sinc filter.
  x = 2*$pi * (i - (size/2 - 1)/2);
  ir[i] = x != 0 ? sin(cutoff * x) / x : cutoff;

  // Blackman-Harris window.
  x = 2*$pi * (i + 1) / (size/2 + 1);
  ir[i] *= 0.35875 - 0.48829 * cos(x) + 0.14128 * cos(2*x) - 0.01168 * cos(3*x);

  i += 1;
);

memset(ir + size/2, 0, size/2);
fft_real(ir, size);
ir[1] = 0; // Correct DC

i = 0;

@sample

in[i] = 0.5*(spl0 + spl1);

(i += 1) >= size/2 ? (
  memset(in + size/2, 0, size/2);
  fft_real(in, size);
  in[1] = 0; // Correct DC

  convolve_c(in, ir, size/2);
  ifft_real(in, size);

  // Overlapp-add.
  memcpy(out, out + size/2, size/2);
  memset(out + size/2, 0, size/2);
  i = 0;
  loop(size,
    out[i] += in[i] * 0.5/size;
    i += 1;
  );
  i = 0;
);

spl0 = spl1 = out[i];
Tale is offline   Reply With Quote
Old 09-08-2016, 11:28 AM   #3
Justin
Administrator
 
Justin's Avatar
 
Join Date: Jan 2005
Location: NYC
Posts: 15,721
Default

Awesome, thanks! And thanks for the real FFT implementation too
Justin is offline   Reply With Quote
Old 09-08-2016, 11:54 AM   #4
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Sounds great !

I am not an expert on that issue at all, but from some discussion about FFT in JSFX I learned that there are some standard tasks that seemingly nearl always need to be done (e.g. managing the FFT window with respect to the sample block size). This is said to be rather time consuming when done in loops in the JSFX.

I wonder if there could be a chance to do native code support functions that allow for speeding up certain standard tasks that are related to FFT handling.

Thanks for listening,
-Michael
mschnell is offline   Reply With Quote
Old 09-08-2016, 01:08 PM   #5
clepsydrae
Human being with feelings
 
clepsydrae's Avatar
 
Join Date: Nov 2011
Posts: 3,409
Default

Thanks Tale, much obliged.
clepsydrae is offline   Reply With Quote
Old 09-10-2016, 03:07 AM   #6
IXix
Human being with feelings
 
Join Date: Jan 2007
Location: mcr:uk
Posts: 3,889
Default

Thanks Tale!

(I think that "Thanks Tale! " is my most posted message)
IXix is offline   Reply With Quote
Old 09-26-2016, 01:57 PM   #7
Smashed Transistors
Human being with feelings
 
Smashed Transistors's Avatar
 
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
Default

Big Thanks Tale
__________________
JSFX plugins and synths. See you here and there: SoundCloud, Youtube, Google Play...
Smashed Transistors is offline   Reply With Quote
Old 11-15-2016, 04:01 PM   #8
Smashed Transistors
Human being with feelings
 
Smashed Transistors's Avatar
 
Join Date: Jul 2014
Location: Là bas les huîtres (FR)
Posts: 424
Default

I'm playing with the fft_real for some mipmapping, it is really great !!
__________________
JSFX plugins and synths. See you here and there: SoundCloud, Youtube, Google Play...
Smashed Transistors is offline   Reply With Quote
Old 11-16-2016, 02:17 AM   #9
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Cool!

BTW, maybe someone could move this thread to the JSFX forum, now that real FFT support is "official"?
Tale is offline   Reply With Quote
Old 11-17-2016, 07:10 AM   #10
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

I analyzed the "FFT SPlitter" JSFX and did some testing.

I understand that the (complex) FFT can convert two (real) signals at the same time and provided that the action done to the spectrum are linear (e.g. a multiplication with another spectrum to do a convolution), the inverse FFT will provide the action (e.g. convolution) applied to both signals.

Is this correct ?

Is it correct to assume that the "real FFT" only handles a single signal and runs faster with the same window size ?

-Michael
mschnell is offline   Reply With Quote
Old 11-17-2016, 09:18 AM   #11
eugen2777
Human being with feelings
 
eugen2777's Avatar
 
Join Date: Aug 2012
Posts: 271
Default

Tale, Thank you! I do not know exactly how it works mathematically, but I used it once in my script. It works twice as fast (I have not used a window)
__________________
ReaScripts
eugen2777 is offline   Reply With Quote
Old 11-17-2016, 10:11 AM   #12
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by eugen2777 View Post
(I have not used a window)
??? FFT always works in Windows (= chunks) of 2^n samples.

-Michael
mschnell is offline   Reply With Quote
Old 11-17-2016, 11:43 AM   #13
eugen2777
Human being with feelings
 
eugen2777's Avatar
 
Join Date: Aug 2012
Posts: 271
Default

I used a rectangular window. This works fine for my needs. You can check it. Link
I made a large window size and use the middle part
Let me explain why I did - it works much faster without much loss
============
Of course, we get some artifacts, but with a large window size are not significant. I wrote about it in the comments
__________________
ReaScripts

Last edited by eugen2777; 11-17-2016 at 12:01 PM.
eugen2777 is offline   Reply With Quote
Old 11-17-2016, 01:25 PM   #14
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

You seem to use "block_size = 1024*16". So FFT-Windowsize = 16384.

-Michael
mschnell is offline   Reply With Quote
Old 11-17-2016, 01:41 PM   #15
eugen2777
Human being with feelings
 
eugen2777's Avatar
 
Join Date: Aug 2012
Posts: 271
Default

Oh sure. I just did not use window function. Rectangular suits me here.
I did a lot of tests, and this works well at large window sizes.
I select middle part of the window to avoid the large artifacts.
Then I connect these areas. I can show you how it sounds
===
This function creates a file from an array of numbers link-lua
You can use it to test filter.
Also, there is eel version link-eel
__________________
ReaScripts

Last edited by eugen2777; 11-17-2016 at 02:06 PM.
eugen2777 is offline   Reply With Quote
Old 11-17-2016, 02:32 PM   #16
DJ Saint-Hubert
Human being with feelings
 
Join Date: Jan 2013
Posts: 257
Default

sorry I have to say this but when did sh*t literally get real???

also I'm messing with FFT to do multiplication, am curious if this gives better results when used for that purpose or whatever?
DJ Saint-Hubert is offline   Reply With Quote
Old 11-17-2016, 02:43 PM   #17
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by mschnell View Post
I understand that the (complex) FFT can convert two (real) signals at the same time and provided that the action done to the spectrum are linear (e.g. a multiplication with another spectrum to do a convolution), the inverse FFT will provide the action (e.g. convolution) applied to both signals.
Well, the regular (complex) FFT transforms a single complex signal from the time domain to the frequency domain. Indeed you can use the complex FFT to transform two real signals at the same time. In fact, that is how the real FFT works internally.

Quote:
Originally Posted by mschnell View Post
Is it correct to assume that the "real FFT" only handles a single signal and runs faster with the same window size ?
Exactly. Note also that the real FFT returns only half the complex frequency bins (only the positive ones). But because the negative frequencies are the complex conjugates of the positive frequencies, nothing is lost (and thus the inverse real FFT can perfectly reconstruct the original signal).
Tale is offline   Reply With Quote
Old 11-17-2016, 02:46 PM   #18
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by eugen2777 View Post
Tale, Thank you! I do not know exactly how it works mathematically, but I used it once in my script. It works twice as fast (I have not used a window)
You are welcome. BTW, if you are interested in how it works, then lookup discrete Fourier transform (DFT). FFT is simply an optimized version of this.
Tale is offline   Reply With Quote
Old 11-17-2016, 03:24 PM   #19
eugen2777
Human being with feelings
 
eugen2777's Avatar
 
Join Date: Aug 2012
Posts: 271
Default

Oh, well, I can use it, I understand how it works, but not inside(mathematics).
In any case, thank You. I look like a fool, but I'm just learning at the age of 33
__________________
ReaScripts
eugen2777 is offline   Reply With Quote
Old 11-18-2016, 01:58 AM   #20
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by eugen2777 View Post
I look like a fool, but I'm just learning at the age of 33
I think I also started with DSP at 33, so that's a good age for it I guess.
Tale is offline   Reply With Quote
Old 11-18-2016, 05:29 AM   #21
eugen2777
Human being with feelings
 
eugen2777's Avatar
 
Join Date: Aug 2012
Posts: 271
Default

Quote:
I think I also started with DSP at 33, so that's a good age for it I guess.
You just give me hope.
I knew mathematics(including mathematical analysis), but now everything is forgotten, but it all can be restored
__________________
ReaScripts
eugen2777 is offline   Reply With Quote
Old 11-18-2016, 05:37 AM   #22
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by Tale View Post
I think I also started with DSP at 33, so that's a good age for it I guess.
I hate you for saying this. I am 63 and eager to learn ...

-Michael
mschnell is offline   Reply With Quote
Old 11-18-2016, 05:45 AM   #23
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by Tale View Post
Exactly. Note also that the real FFT returns only half the complex frequency bins (only the positive ones). But because the negative frequencies are the complex conjugates of the positive frequencies, nothing is lost (and thus the inverse real FFT can perfectly reconstruct the original signal).
Great !!!
Thanks a lot !!

I'll try to modify the example JSFX code I am working on (that had been derived from FFTSplitter) to use the Real FFT calls.

I added a linear (in frequency) graph, so I can easily see what happens.

With "negative frequencies" I suppose you mean those that are in the buffer above index "fftsize" (or 2*fftsize, as complex values take two places). I understand that the negative values are shown there (as in DFT the spectrum is a periodic function).

I fear I will be back with more questions .

-Michael
mschnell is offline   Reply With Quote
Old 11-18-2016, 11:23 AM   #24
eugen2777
Human being with feelings
 
eugen2777's Avatar
 
Join Date: Aug 2012
Posts: 271
Default

Quote:
Indeed you can use the complex FFT to transform two real signals at the same time. In fact, that is how the real FFT works internally.
Ops. I used it. "Stereo FFT" - re=left channel, im=right channel. I found it here(second answer)
http://stackoverflow.com/questions/1...nal-seperately
Maybe it will useful for someone else
__________________
ReaScripts
eugen2777 is offline   Reply With Quote
Old 11-19-2016, 03:35 AM   #25
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by mschnell View Post
I fear I will be back with more questions .
Hi, Tale, as I have your attention right now two questions that came up when I did tests with the FFT (FFT -> do some stuff with the spectrum -> reverse FFT, testing e.g. with a "band limited" rectangle signal as an input):

1)
copied from the FFTSplitter I have this code snippet:
Code:
    i = 2;
    sumRe = sumIm = 0.0;
    loop(fftsize,
      sumRe += result1[i];
      sumIm += result1[i+1];
      i += 2;
    );
    result1[0] = -sumRe;
    result1[1] = -sumIm;
that is executed before the reverse FFT.

Seemingly this should set the DC offset (Frequency=0) to zero, by compensating the sum of the other frequencies.

I found that with my "good" experiments, this does (close to) nothing.

With certain "bad" (not really doing what I want them to) experiments (e.g. trying to shift the phase or the frequency) the result seems to be better when I drop this correction.

What is the said code really supposed to accomplish ?



2)
With the spectrum values I did this:
Code:
      loop(fftsize-1,  
        rea = result1[i];
        ima = result1[i+1];
        absol = sqrt(rea*rea+ima*ima);
        phase = atan2(ima, rea);
        tanph = tan(phase);

        phase += 0*3.14159;

        rea = absol*cos(phase);
        ima = absol*sin(phase);
        result1[i] = rea;
        result1[i+1] = ima;
        i += 2;
      );
With "phase += 0*3.14159" of course nothing is modified. With "phase += 1*3.14159" the signal is inverted, etc (as expected). But with non integer values I do not get a phase shift but a rather chaotic oscillogram, e.g. all resulting samples negative.



3)
I tried to do a frequency shift, e.g. doubling the frequency by moving all spectrum values up to the double index for the lower half and appropriately shifting down (or recreating from the lower = same result) the upper half, filling in zeroes to the gaps.

Here again the output signal looks rather funny. I do get the correct wave form, but it does not stay in a position (the virtual oscilloscope is triggered by the original wave form) but the phase and the volume seems to be constantly modifying. The sound is "chopped", supposedly according to the frame frequency.

Any analyzing, tips and tricks ?

Thanks a lot for any answer,
-Michael

Last edited by mschnell; 11-19-2016 at 03:40 AM.
mschnell is offline   Reply With Quote
Old 11-20-2016, 03:56 AM   #26
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by mschnell View Post
Seemingly this should set the DC offset (Frequency=0) to zero, by compensating the sum of the other frequencies.
In the time domain the DC is the average of the signal, and you can remove it by subtracting the average from the signal. However, in the frequency domain it is simply the first bin, so you can simply set it to zero. See also my real FFT tone generator example above.

Note that the Nyquist real value is stored in the imagnary part of the DC bin, i.e. there is no imaginary part for DC (or Nyquist).

Quote:
Originally Posted by mschnell View Post
Any analyzing, tips and tricks ?
Maybe you forgot to do [i]fft_permute()?
Tale is offline   Reply With Quote
Old 11-20-2016, 04:46 AM   #27
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by Tale View Post
Note that the Nyquist real value is stored in the imagnary part of the DC bin, i.e. there is no imaginary part for DC (or Nyquist).
Obviously there is no phase to frequency 0.
I don't know yet what the "Nyquist real value" is.
I'll try to set the real and/or imaginary part of the DC bin to Zero and see what happens...
Quote:
Originally Posted by Tale View Post
Maybe you forgot to do [i]fft_permute()?
No it's in place (as I started by copying FFTSplitter)

So you suggest that (e.g)adding a constant to phase in fact should create a correctly shifted waveform even for non-integer multiple of pi ? I just tested that with my code trying to shift the phase by 0.25*pi provides silence .

EDIT
Maybe this even is correct as phase shifting a random signal other than n*pi does not really make sense.

Time-shifting would make sense. I understand that this is a convolution with a delta signal and hence multiplying the spectrum with a complex cos (or similar)...

Supposedly my initial idea of introducing a frequency shifting is incorrect as well.


Thanks for you thoughts !
-Michael

Here is the complete code, in case you are interested:
Code:
desc: FFT Test
//tags: processing FFT routing
//author: Schwa

slider1:3<0,4,1{64,128,256,1024,2048,4096,8192}>FFT Size
//slider2:5000<0,20000,1>Split Frequency (Hz)
//slider3:0<0,4,1{1/2,3/4,5/6,7/8}>Low Band Destination
//slider4:1<0,4,1{1/2,3/4,5/6,7/8}>High Band Destination
// slider5:0<0,1,0.1>Debug

slider2:1<0,2,1{lin,log}>Graph
slider3:1<0,2,1{full,half}>Graph
slider4:0<0,5,1{no,lowpass,pitch,both, phase_A, Phase_B}>Filter
slider5:5000<0,20000,1>Low Pass  Frequency (Hz)
slider6:1<0,2,1{no,yes}>reconstruct upper half

@init 
  maxmax = 1000000000;


  fftsize = -1;
  pdc_bot_ch = 0;
  pdc_top_ch = 8;

  // array pointers
  window=0;
  buf1     = 16384;
  buf2     = buf1*2;
  buf3     = buf1*3;
  buf4     = buf1*4;
  result1  = buf1*5;
  result2  = buf1*7;
  
  buf_21   = buf1*8;
  buf_22   = buf1*9;
  
  bufg0    = buf1*10;
  bufg1    = buf1*11;
  
  do_display = 0;

@slider
  fftsize2 = fftsize*2; 
  pos = 0;
  memset(buf1,    0, fftsize2);
  memset(buf2,    0, fftsize2);
  memset(result1, 0, fftsize2);
  memset(result2, 0, fftsize2);

  sliderfft = (2^(slider1+7))|0;
  fftsize != sliderfft ? (
    fftsize2 = fftsize*2;
    memset(window, 0, fftsize); 
    fftsize = sliderfft;
    w = 2.0*3.14159/fftsize;
    i = 0;
    loop(fftsize/2,
      window[i] = 0.42-0.50*cos(i*w)+0.08*cos(2.0*i*w);
      i += 1;
    ); 
    pdc_delay = fftsize;
  );  



@sample

  pos >= fftsize ? (
    //memcpy(buf1, buf2+fftsize, fftsize);
    tmp    = buf1;
    buf1   = buf2;
    buf2   = tmp;
 
    tmp    = buf_21;
    buf_21 = buf_22;
    buf_22 = tmp;
 
 
    tmp = result1;
    result1 = result2;
    result2 = tmp;    
    
    fft(buf1, fftsize);
    fft_permute(buf1, fftsize);
 
//    fft(buf_21, fftsize);
//    fft_permute(buf_21, fftsize);

    memcpy(result1+2, buf1+2, fftsize2-2);
    
    !do_display ? (
      memcpy(bufg0, result1, fftsize2); 
    );  

  
//Test1
    !(slider4 & 4) && slider4 & 1 ? (
      splitfreq = slider5;
      splitfreq = 2*((fftsize * splitfreq / srate)|0);
      memset(result1+splitfreq+2, 0, 2*(fftsize-splitfreq-1));
//      memset(result1+fftsize, 0, fftsize);
//      memset(result1+2, 0, fftsize);
      i = 0;
      ii = fftsize;
      loop(i,     
        result1[fftsize2-i-2] = result[i];
        result1[fftsize2-i-1] = -result[i+1];
        i += 2;
      );
    );  
    
// - Test 1



// Test 2
    !(slider4 & 4) && slider4 & 2 ? (
      buf3[0] = 0;
      buf3[1] = 0;
      i = 2;
      loop(fftsize/4,     
        buf3[2*i+0] =  result1[i];
        buf3[2*i+1] =  result1[i+1];
        buf3[2*i+2] =  0;
        buf3[2*i+3] =  0;
      
//        buf3[fftsize2-2*i-4] =  0;
//        buf3[fftsize2-2*i-3] =  0; 
        buf3[fftsize2-2*i-4] =  result1[fftsize2-i-2];
        buf3[fftsize2-2*i-3] =  result1[fftsize2-i-1]; 
        buf3[fftsize2-2*i-2] =  0;
        buf3[fftsize2-2*i-1] =  0;
      
        i += 2
      );
      memcpy(result1+2, buf3+2, fftsize2-2);
    );  

// - Test2

// Test 3
    slider4 & 4 ? (
      i = 0;
      loop(fftsize-1,  
        rea = result1[i];
        ima = result1[i+1];
        slider4 & 1 ? (
          absol = sqrt(rea*rea+ima*ima);
          phase = atan2(ima, rea);
          tanph = tan(phase);
/*          
          phpo  = rea > 0;          
          ret = absol / sqrt(1+tanph*tanph);
          imt = rea * tanph;          
          phpo > 0 ? (
            rea = ret;
           ) : ( 
            rea = -ret;
          );               
          ima = imt;                    
*/        
          phase += 0*3.14159;

          rea = absol*cos(phase);
          ima = absol*sin(phase);
         ) : (
          tmpa = rea;
          rea = -ima;
          ima = -tmpa;
        );
        result1[i] = rea;
        result1[i+1] = ima;
        i += 2;
      );
    );  
// -Test 3


    slider6 ? (
      ii = fftsize;
      i = 0;
      loop(i,     
        result1[fftsize2-i-2] = result[i];
        result1[fftsize2-i-1] = -result[i+1];
        i += 2;
      );
    );  


    !do_display ? (
      memcpy(bufg1, result1, fftsize2); 
      do_display = 1;
    );  

    ii = fftsize;                       /* wozu soll das gut sein ? */
//    ii = 0;
    i = 2;
    sumRe = sumIm = 0.0;
    loop(ii,
      sumRe += result1[i];
      sumIm += result1[i+1];
      i += 2;
    );
    result1[0] = -sumRe;
    result1[1] = -sumIm;

    fft_ipermute(result1, fftsize);
    ifft(result1, fftsize);

    pos=0;
  );

  w1 = window[pos/2];
  w2 = window[(fftsize-pos)/2-1];
  sw = (w1+w2)*fftsize;

  buf1[pos] = w1*spl0;
//  buf1[pos+1] = w1*spl1;
  buf1[pos+1] = w1*spl0;
  buf2[fftsize+pos] = w2*spl0;
//  buf2[fftsize+pos+1] = w2*spl1;
  buf2[fftsize+pos+1] = w2*spl0;

  spl0r = (result1[pos]+result2[fftsize+pos])/sw;
//  spl1r = (result1[pos+1]+result2[fftsize+pos+1])/sw;    // ??? imaginaerteil ??? 

  
  spl4 = spl0r;
  spl5 = spl0r;

//  spl4 = spl0;
//  spl5 = spl1;

  pos += 2;


@gfx 640 400


  function f0 (x) (
    x2 = x * 2;
    re = bufg0[x2];
    im = bufg0[x2+1];
    y = re*re+im*im;
    y = sqrt(y);
    slider2 ? (
      y = log(y);
     ) : ( 
      y > maxmax ? y = maxmax;
    );  
//    x/2 + fftsize/4;
  );

  function f1 (x) (
    x2 = x * 2;
    re = bufg1[x2];
    im = bufg1[x2+1];
    y = re*re+im*im;
    y = sqrt(y);
    slider2 ? (
      y = log(y);
     ) : ( 
      y > maxmax ? y = maxmax;
    );    
  );

do_display ? (
  xgmin = 2;

  do_display = 0;
  gfx_r=gfx_g=gfx_b=0; gfx_a=1;
  gfx_x=gfx_y=0;
  gfx_rectto(gfx_w,gfx_h);

  q1 = gfx_w;
  q2 = gfx_h;


//gfx_line();
  gfx_r=gfx_g=gfx_b=0;
  gfx_g= 1;
  gfx_y = 0;
  gfx_x = 0;
  x = xgmin;
  xw = fftsize;
//  xw = fftsize / 2;
  xw = fftsize / (slider3+1);
  yymax = 0;
  yymintest = 999999999;
  
  
  while (x<xw) (
    y = f0(x);
    y > yymax ? (
      yymax = y;
    );
    y < yymintest ? (
      yymintest = y;
    );
    x = x+1; 
  );
  
  
   
//  yymax = 40;
  
  yy = gfx_h / yymax;
  
  
  dg =  xw / gfx_w;
  dg < 1 ? dg = 1;
  x  = xgmin;
  xg = 0;
  xd = dg;
  ymax = 0;
  while (x<xw) (
    y = f0(x);
    y > ymax ? (
      ymax = y;
    ); 
    x > xd ? ( 
      xd = xd + dg;
      gfx_lineto(xg, gfx_h-ymax*yy, 1);
//      gfx_lineto(xg, 50, 1);
      ymax = 0;
      xg = xg + 1;
    );    
    x = x + 1;
  );


//  yymax = 40;
  
  yy = gfx_h / yymax;
  gfx_x=gfx_y=0;


  gfx_b = gfx_g = 0;
  gfx_r = 1;
  x  = xgmin;
  xg = 0;
  xd = dg;
  ymax = 0;

  while (x<xw) (
    y = f1(x);
    y > ymax ? (
      ymax = y;
    ); 
    x > xd ? ( 
      xd = xd + dg;
      gfx_lineto(xg, gfx_h-ymax*yy /*- 20*/, 1);
//      gfx_lineto(xg, 50, 1);
      ymax = 0;
      xg = xg + 1;
    );    
    x = x + 1;
  );
);

Last edited by mschnell; 11-20-2016 at 05:14 AM.
mschnell is offline   Reply With Quote
Old 11-20-2016, 05:40 AM   #28
octarone
Human being with feelings
 
octarone's Avatar
 
Join Date: Apr 2016
Posts: 30
Default

I just stumbled on this since I'm not very active around here... but this is absolutely fantastic! Now I can make a proper mono version of a JS I'm working on.

I understand it's obviously missing from ReaJS since it was only recently added to REAPER, I hope you add fft_real there soon since it is very useful. (will next update of ReaJS be early next year like it was this year? I hope not too much longer )


@mschnell: The "Nyquist frequency" is simply the last bin, but it is stored in the imaginary part of the 'DC' bin (which is the first bin, zeroth) because there's no room for it at the end, it's just a gimmick, not mathematically correct.

This makes some calculations with loops a bit more involved because you have to treat the Nyquist separately, and its position (as second and imaginary) makes no mathematical sense. But this is the only way to store the fft_real with the same number of outputs as inputs so that it can be used in place in a buffer. Some other FFT implementations/libraries do a similar thing.

The DC component has no phase and thus no imaginary component, and the Nyquist has no imaginary component (always 0). Think of a frequency bin as a 2D vector, with x and y coordinates. The phase is the angle, while the magnitude is... well the magnitude of the vector.

Why does the Nyquist have no imaginary component? Think of it like this. The Nyquist frequency is always 2 samples long: a high value followed by a low value (to make an oscillation). This can be either +/- or -/+ (high->low or low->high), so it does have phase, but it is restricted to two phases only: zero degrees, and pi (180 degrees) which inverts its polarity.

Both of those phases, on the vector, are a straight line on 'x', with 'y' being zero. It's just that the second case has a negative 'x'. Magnitude is of course calculated as sqrt(x*x+y*y) and thus always positive. So it does have phase, but only either 0 or pi, which still has no imaginary component (y), so it is not stored.


Lastly, I'm telling you about the vectors because the way you add phases is hopelessly tedious and slow. You can simplify it by just rotating the vector.

Remember: the phase is the angle of the vector (compared to the x axis, with positive angle being counter-clockwise rotation). So, if you want to add a phase X, just rotate the vector by X counter-clockwise.

No need for any trigonometry when rotating, you only need sin/cos for the angle, but you can precalculate that, because it's the same for every bin right? So only one sin/cos is needed before the loop.

(I can give example but 2D vector rotation should be easy to just google with proper diagrams and better explanations than mine )
__________________
My JSFX & other things at my simple site, if you are interested: http://sites.google.com/site/octarone
octarone is offline   Reply With Quote
Old 11-20-2016, 03:49 PM   #29
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Thanks a lot for the explanations.

So the "Nyquist Frequency" has only to be considered with real FFT. Right now I am just doing the normal complex stuff, just to get some understanding of what is happening.

But still I don't understand the point with the loop to set the real and imaginary part of the DC bin (with the complex FFT). Seemingly setting both ot zero or leaving them as they are does not harm unless the operation on the spectrum produces something unusable, anyhow.

This also is the motivation for the tests with certain operations on the spectrum. I did not consider any optimization right now.

An application I have in mind is "spectral subtraction". I will try to implement this once I feel I have enough understanding on the algorithms involved.

-Michael

Last edited by mschnell; 11-20-2016 at 03:55 PM.
mschnell is offline   Reply With Quote
Old 11-20-2016, 04:34 PM   #30
octarone
Human being with feelings
 
octarone's Avatar
 
Join Date: Apr 2016
Posts: 30
Default

The Nyquist always exists, but with normal complex FFT, it's the last bin -- in its proper place.

The "Nyquist" frequency is simply the highest possible frequency for your sample rate -- i.e half your sample rate. Such a waveform is constructed by alternating two sample points in oscillation; you can't go "faster" than this (because of your sampling rate).

But, with normal complex FFT, it does have imaginary component, because the input is complex too. Even the DC has an imaginary component (which corresponds to the DC of the imaginary input). The reason it's stored weirdly with the real fft is again, because there's "no space" in the buffer otherwise. If the input is real, then both the DC and Nyquist have imaginary component set to 0, so it's not stored, so Nyquist's real value (only value that is stored) was placed where the DC had its imaginary component.

As for your question about setting the DC to zero, yes I agree.

I don't understand that code either. If you want to get rid of DC just set it to zero, but IMHO if you do "spectral processing" (on arbitrary input), you should treat DC properly for its magnitude (no phase though), with low overlap it may make a difference.

You can set Nyquist to 0 if you don't want to bother processing it, it's probably inaudible, but I'm doing it "properly" because I'm a perfectionist.

Anyway the Nyquist (last bin, highest frequency) being put in the first frequency's imaginary part is a "trick" used by real fft, it doesn't actually make any mathematical sense. So don't think too hard about it. (because mathematically, FFT is a complex-input, complex-output process, so real FFT itself is a 'trick')
__________________
My JSFX & other things at my simple site, if you are interested: http://sites.google.com/site/octarone
octarone is offline   Reply With Quote
Old 11-21-2016, 08:57 AM   #31
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

I think I understand. Thanks a lot for the explanation !

Now trying to find out how a frequency shift (e.g. doubling it) can be done by tweaking the spectrum. Supposedly my initial test handles phase in a wrong way (i.e: not at all)

Thanks again,
-Michael
mschnell is offline   Reply With Quote
Old 11-23-2016, 02:04 PM   #32
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Trying to wrap my head around the Real FFT, I enhanced my example program to be able to select either normal or Real FFT.

Up till now I did not get the real FFT to work as it should with my simple example doing FFT, showing a graph of the spectrum, and doing reverse FFT, that then should reproduce the signal.


Please let me know if I am correct with my assumptions:

I understand that with normal (complex) FFT, we have two input "channels" (officially a complex functions consisting of real and imaginary part, but as FFT is linear it also can be two real functions that independently get regenerated by the inverse FFT).

I understand that FFT converts a real function in a symmetric function, so the upper half of the spectrum is "dummy", which is perfectly according to the Nyquist theorem that says that frequencies above half the sample rate are forbidden.

I understand that FFT always converts a (presumed) periodic function in a periodic function. To tame that effect, we need to use the "Window" function that dims the signal at both ends.

That is why my code (derived form the "FFT-Splitter") uses two overlapping buffers, filled in a way that a new one start at the half of the previous one.

Now I am wondering why the "Fir Filter" example given above does not use overlapping buffers. IMHO this would result in certain sample not being decently used in the process.


I understand that the Real FFT does not use complex numbers in time domain. So only half of the buffer is used for samples. That is why the upper half is filled with zeroes.

Now I wonder why the fft_real function does not fill in the zeroes itself, as this action seemingly is generally necessary.

I understand that the spectrum created by fft_real should be the same as the lower half of that crated by a complex fft, with the imaginary part of the time domain function set to zero.

Unfortunately in my example code it is similar but not equal: with a "saw" example waveform, the peaks are not sharp any more, they get lower more quickly, and they wobble slightly.

I understand that doing a reverse real FFT from a spectrum generated by a real FFT should reproduce the real waveform.

In my example the resulting wave form is a lot smaller, nearly a sine, and the phase is not fixed.

Any help ?

-Michael

Last edited by mschnell; 11-23-2016 at 11:01 PM.
mschnell is offline   Reply With Quote
Old 11-24-2016, 02:19 AM   #33
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by mschnell View Post
I understand that with normal (complex) FFT, we have two input "channels" (officially a complex functions consisting of real and imaginary part, but as FFT is linear it also can be two real functions that independently get regenerated by the inverse FFT).
Correct, although you could run into trouble if the two real signals are orders of magnitude apart (I guess because there is a relation bewteen the real and imaginary parts of a complex number, but there is no relation between the two real signals).

Quote:
Originally Posted by mschnell View Post
I understand that FFT converts a real function in a symmetric function, so the upper half of the spectrum is "dummy", which is perfectly according to the Nyquist theorem that says that frequencies above half the sample rate are forbidden.
I think the Nyquist theorem holds also for complex signals, which do have the "upper" half of the spectrum. In fact, the upper half are the negative frequencies. For a real signal these are the same as the positive frequencies, but with their sign flipped, which is why the real FFT doesn't store them (but conceptually they are still there).

Quote:
Originally Posted by mschnell View Post
I understand that FFT always converts a (presumed) periodic function in a periodic function. To tame that effect, we need to use the "Window" function that dims the signal at both ends.
The FFT transforms a periodic function (i.e. a single-cycle waveform) to its Fourier coefficients, or the other way around for the inverse FFT.

Note that you only need a (non-rectangular) window if you want to "glue" together a signal that is longer than your FFT size. So if your are taking the FFT of an impulse response, where you probably already know your FFT is longer than the IR, then you don't need a window. The same goes for a single-cycle waveform.

Quote:
Originally Posted by mschnell View Post
I understand that the Real FFT does not use complex numbers in time domain. So only half of the buffer is used for samples. That is why the upper half is filled with zeroes.
The real FFT simply transforms buf_size real numbers, although internally it calls the complex FFT, transforming buf_size/2 complex numbers, which is why you need fft_permute(buf_size/2). Then it reorders the bins, so it returns buf_size/2 complex numbers that represent only the positive spectrum.

Quote:
Originally Posted by mschnell View Post
Now I wonder why the fft_real function does not fill in the zeroes itself, as this action seemingly is generally necessary.
Because they are not necessary.

Quote:
Originally Posted by mschnell View Post
I understand that the spectrum created by fft_real should be the same as the lower half of that crated by a complex fft, with the imaginary part of the time domain function set to zero.
Correct.

Quote:
Originally Posted by mschnell View Post
I understand that doing a reverse real FFT from a spectrum generated by a real FFT should reproduce the real waveform.
Yes, it should.

Quote:
Originally Posted by mschnell View Post
In my example the resulting wave form is a lot smaller, nearly a sine, and the phase is not fixed.
I can't really help you there... Have you tried the two examples I posted at the beginning of this topic?
Tale is offline   Reply With Quote
Old 11-24-2016, 05:50 AM   #34
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by Tale View Post
The FFT transforms a periodic function (i.e. a single-cycle waveform) to its Fourier coefficients, or the other way around for the inverse FFT.

Note that you only need a (non-rectangular) window if you want to "glue" together a signal that is longer than your FFT size. So if your are taking the FFT of an impulse response, where you probably already know your FFT is longer than the IR, then you don't need a window. The same goes for a single-cycle waveform.
AFAIK, Theory holds that Fourier transform (both directions) of a sampled (i.e. steps with constant length) function results in a periodic Function (period = window size / sample frequency, respectively window size * sample frequency) on the other side. So digital (complex) Fourier transform assumes a periodic function at both site. So if not using a smooth Windows function in time domain, the result in frequency domain will be will
very inappropriate if the function is not periodic with a "fitting" frequency. Of course chunks of arbitrary functions can only be glued if the "forced" period is tamed by a smooth Windows function. But if not using overlapping windows, the rim-samples are weighted less than the center which seems rather inappropriate.

Quote:
Originally Posted by Tale View Post
Because they are not necessary.
That is what I initially assumed. But not setting them to zero gave a completely wrong spectrum, wile with setting them to zero the spectrum is at least similar to what I expect (the one given by the complex FFT). Am I doing something wrong with the size parameter ? (I understand that for the same window size the "size" parameter is the same for complex and real FFT, while it's different for fft_permute.) For doing a spectrum would it be unnecessary to set the upper half of the values to Zero ?

Quote:
Originally Posted by Tale View Post
I can't really help you there... Have you tried the two examples I posted at the beginning of this topic?
Of course I looked deeply at them (here I found the setting the upper part to zero:
Code:
memset(ir + size/2, 0, size/2);
fft_real(ir, size);
), but for me the first step is doing a spectrum (and this already works nicely with the complex FFT) so I try to do a "spectrum" example with the real FFT before looking at the inverse real FFT, but right now, unfortunately in my example code the spectrum is similar but not equal to the one created by the complex FFT: with a "saw" example waveform, the peaks are not sharp any more, they get lower more quickly, and they wobble slightly.

Thanks,
-Michael

Last edited by mschnell; 11-24-2016 at 06:18 AM.
mschnell is offline   Reply With Quote
Old 11-27-2016, 03:29 AM   #35
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by mschnell View Post
Of course I looked deeply at them (here I found the setting the upper part to zero:
Code:
memset(ir + size/2, 0, size/2);
fft_real(ir, size);
Ah, but that's part of the overlap-add method used in that example (I have an IR at size/2, but I use the FFT at full size, so I have to zero pad). It has nothing to do with real vs complex FFT.
Tale is offline   Reply With Quote
Old 11-27-2016, 04:13 AM   #36
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by Tale View Post
Ah, but that's part of the overlap-add method used in that example (I have an IR at size/2, but I use the FFT at full size, so I have to zero pad). It has nothing to do with real vs complex FFT.
Ah, now I understand what "overlap-add" means. I supposed the double buffer method used in the "FFT-Splitter" example was ubiquitously used.

Anyway this seems to mean that in my code (doing the overlap by dual buffers, in the same way as done in the "FFT-Splitter" example) the zero setting should not be done.

In fact when I comment out that line, the result is completely wrong.

So I'll re-check my code due to this ...

(I still suppose something is wrong with the use of size parameter. Seemingly there is no documentation on fft_real() )

With both overlapping methods, of course the FFT transform size is double to the chunks taken from the stream, as each transform uses two chunks, and each chunk is used twice. Now I want to use the same chunk size and the same FFT size for both complex and real fft functions.

-Michael

Last edited by mschnell; 11-27-2016 at 04:27 AM.
mschnell is offline   Reply With Quote
Old 11-28-2016, 11:21 PM   #37
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Hi Tale,

I can't get FFT_real working for me .

I did as pure (mono) Spectrum analyzer example that can be switched between three modes:

- Complex FFT, view full spectrum including the part higher than the Nyquist frequency
- Complex FFT, schowing only the lower half of the spectrum up to the than the Nyquist frequency
- real FFT.

I understand that the real FFT should show exactly the same as the "half" complex FFT, But I can't get this working.

UI stripped of all unnecessary code from my example project so it should be easy to debug. (Complex seems to work fine: I use a rectangle, saw ror triangle generator as an example so9urce.)

Maybe you can take a look at it.

Thanks a lot !
-Michael

Code:
desc: FFT Real Spectrum

slider1:3<0,4,1{64,128,256,1024,2048,4096,8192}>FFT Size
slider2:1<0,2,1{lin,log}>Graph
slider3:1<0,3,1{full,half,real}>Graph / Real FFT

@init 
  maxmax = 1000000000;


  fftsize = -1;
  pdc_bot_ch = 0;
  pdc_top_ch = 8;

  // array pointers
  window=0;
  buf1     = 16384;
  buf2     = buf1*2;
  result1  = buf1*3;
  result2  = buf1*4;
  bufg0    = buf1*5;
  do_display = 0;

@slider
  realfft = slider3 == 2;
  pos = 0;
  memset(buf1,    0, fftsize2);
  memset(buf2,    0, fftsize2);
  memset(result1, 0, fftsize2);
  memset(result2, 0, fftsize2);

  sliderfft = (2^(slider1+7))|0;
  fftsize != sliderfft ? (
    fftsize2 = fftsize*2;
    memset(window, 0, fftsize); 
    fftsize = sliderfft;
    w = 2.0*3.14159/fftsize;
    i = 0;
    loop(fftsize/2,
      window[i] = 0.42-0.50*cos(i*w)+0.08*cos(2.0*i*w);
      i += 1;
    ); 
    pdc_delay = fftsize;
  );  
  fftsize2 = fftsize*2; 
  fftsize_x = fftsize;
  realfft ? fftsize_x /= 2;



@sample

  pos >= fftsize_x ? (
    tmp    = buf1;
    buf1   = buf2;
    buf2   = tmp;
 
    tmp = result1;
    result1 = result2;
    result2 = tmp;    
    
    realfft ? (
//      memset(buf1+fftsize_x, 0, fftsize_x);   // ?????????????????????????????????????????
      fft_real(buf1, fftsize/2);
      fft_permute(buf1, fftsize/4);
     ) : (
      fft(buf1, fftsize);
      fft_permute(buf1, fftsize);
    );
 
//    memcpy(result1+2, buf1+2, fftsize2-2);
//    memcpy(result1, buf1, fftsize_x);
    xw = fftsize*2;
    slider3 == 1 ? (
      xw /= 2;
     ) : slider3 == 2 ? (
      xw /= 2;  
    );  
    memcpy(result1, buf1, xw);

    
    !do_display ? (
      memcpy(bufg0, result1, slider3? fftsize : fftsize2); 
      do_display = 1;
    );  

  
    pos=0;
  );


  fftreal ? (
    w1 = window[pos];
    w2 = window[fftsize_x-pos-1];
    sw = (w1+w2)*fftsize;
    buf1[pos] = w1*spl0;
    buf2[fftsize_x+pos] = w2*spl0;
    spl0r = (result1[pos]+result2[fftsize_x+pos])/sw;
    pos += 1;
   ) : (
    w1 = window[pos/2];
    w2 = window[(fftsize-pos)/2-1];
    sw = (w1+w2)*fftsize;

    buf1[pos] = w1*spl0;
    buf1[pos+1] = 0;
    buf2[fftsize+pos] = w2*spl0;
    buf2[fftsize+pos+1] = 0;
    pos += 2;
  );

@gfx 640 400


  function f0 (x) local (x2, re, im, y) (
    x2 = x * 2;
    re = bufg0[x2];
    im = bufg0[x2+1];
    y = re*re+im*im;
    y = sqrt(y);
    slider2 ? (
      y = log(y);
     ) : ( 
      y > maxmax ? y = maxmax;
    );  
  );

do_display ? (
  xgmin = 2;

  do_display = 0;
  gfx_r=gfx_g=gfx_b=0; gfx_a=1;
  gfx_x=gfx_y=0;
  gfx_rectto(gfx_w,gfx_h);

  q1 = gfx_w;
  q2 = gfx_h;


  gfx_r=gfx_g=gfx_b=0;
  gfx_g= 1;
  gfx_y = 0;
  gfx_x = 0;
  x = xgmin;
  xw = fftsize;
  slider3 == 1 ? (
    xw /= 2;
   ) : slider3 == 2 ? (
    xw /= 4;  
  );  
  yymax = 0;
  yymintest = 999999999;
  
  while (x<xw) (
    y = f0(x);
    y > yymax ? (
      yymax = y;
    );
    y < yymintest ? (
      yymintest = y;
    );
    x = x+1; 
  );
  
  
  yy = gfx_h / yymax;
  
  
  dg =  xw / gfx_w;
  dg < 1 ? dg = 1;
  x  = xgmin;
  xg = 0;
  xd = dg;
  ymax = 0;
  while (x<xw) (
    y = f0(x);
    y > ymax ? (
      ymax = y;
    ); 
    x > xd ? ( 
      xd = xd + dg;
      gfx_lineto(xg, gfx_h-ymax*yy, 1);
      ymax = 0;
      xg = xg + 1;
    );    
    x = x + 1;
  );
);
mschnell is offline   Reply With Quote
Old 11-29-2016, 03:04 AM   #38
Tale
Human being with feelings
 
Tale's Avatar
 
Join Date: Jul 2008
Location: The Netherlands
Posts: 3,645
Default

Quote:
Originally Posted by mschnell View Post
Maybe you can take a look at it.
Yeah sure.

I think you are on the right track, but there are a couple of minor issues.

Quote:
Originally Posted by mschnell View Post
Code:
      fft_real(buf1, fftsize/2);
      fft_permute(buf1, fftsize/4);
Although the real FFT does use only half the bytes in memory, it uses the same FFT size as the complex FFT, so I think this should be:

Code:
      fft_real(buf1, fftsize);
      fft_permute(buf1, fftsize/2);
Then you have a typo:

Quote:
Originally Posted by mschnell View Post
Code:
  fftreal ? (
Your code doesn't have a variable called "fftreal", so it is always zero, and thus your real FFT buffer code is never executed. Change it to this, and it will magically all start working:

Code:
  realfft ? (
And lastly:

Quote:
Originally Posted by mschnell View Post
Code:
    buf1[pos] = w1*spl0;
    buf2[fftsize_x+pos] = w2*spl0;
The real FFT uses different scaling, so you have to scale by 0.5 (even when you are just doing a forward FFT) i.e.:
Code:
    buf1[pos] = 0.5 * w1*spl0;
    buf2[fftsize_x+pos] = 0.5 * w2*spl0;
Now the real FTT will exacly match the half complex FFT.

Last edited by Tale; 11-30-2016 at 09:16 AM. Reason: Double brainfart
Tale is offline   Reply With Quote
Old 11-29-2016, 05:55 AM   #39
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Quote:
Originally Posted by Tale View Post
Yeah sure.
Your code doesn't have a variable called "fftreal",
Thanks a lot.

There needed to be a really silly but that prevented me from preceeding...

I'll take a look at your other suggestions.

-Michael
mschnell is offline   Reply With Quote
Old 11-30-2016, 07:15 AM   #40
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 14,686
Default

Success !!!!

After correcting the silly variable name error it does work fine: I get the same spectrum with real FFT as with complex FFT. It's hard to measure, but it seems like performance better with real FFT.

Just a small thing left that I don't understand: with real FFT the (absolute) values of the spectrum are exactly double than with complex FFT. Is this somehow understandable (by design) or is it a hint that my code still features a bug ?

Thanks,
-Michael
mschnell is offline   Reply With Quote
Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off

Forum Jump


All times are GMT -7. The time now is 11:03 AM.


Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2024, vBulletin Solutions Inc.