Old 07-17-2017, 05:43 AM   #1
weblordpepe
Human being with feelings
 
Join Date: Jul 2016
Posts: 29
Default FFT oscillating - perhaps windowing issue

Hi Guys.

I'm new to coding FFT and I have some example code I am modifying for learning.
The following code, sure enough, boosts 4KHz.

The problem is that it oscillates. If I increase the FFT size (lots of bins) to 32,768 then it oscillates very badly.

By Oscillating I mean this 4kHz will wobble around maximum amplitude, not staying still.
In addition, any frequency I attenuate down will wobble instead of staying down.

I believe its related to windowing, however I don't understand the code enough to fix this problem.

Please help!

Code:
desc:pepe FFT spectrum sandpit



@init

fft_size =2048*8;     //CHANGE THIS FFT SIZE TO LARGE TO SEE WOBBLY ACTION
delay_milliseconds = 1/srate * fft_size * 1000;

fft_buf_size = fft_size * 2;
fft_buf = 0;


out_buf = fft_buf + fft_buf_size;      //why is this not zero like the fft_buf?
out_buf_size = fft_size;

NUM_BINS=fft_size;        //number of positive bins. we have another NUM_BINS spaced after for our negative frequencies 
BIN_SIZE=srate/NUM_BINS;
NYQUIST=srate/2;

// During the first fft_size sample we don't have any output yet, so
// prefill the output buffer with silence.
memset(out_buf, 0, fft_size);

in_idx = 0;
out_idx = 0;

function freqToBin(f)(
 2*(  1 +(f/ srate * fft_size));

);
function setBin(bin,g)(
    fft_buf[bin] =g;
    fft_buf[fft_buf_size - bin] =g;
);

function norm(g)(
  g/fft_size;
);
@slider

@sample
  
  // Mono input.
  in = 0.5 * (spl0 + spl1);
  
  // Convert to complex number.
  re = in;
  im = 0.0;
  
  // Add to FFT buffer.
  fft_buf[in_idx * 2 + 0] = re;
  fft_buf[in_idx * 2 + 1] = im;
  in_idx += 1;
  
  // Do we have enough samples for FFT?
  in_idx >= fft_size ? (
  
    // FFT 
    fft(fft_buf, fft_size);
    
    //sort FFT output by frequency (permute?!)
    fft_permute(fft_buf, fft_size);
   
    
    b = freqToBin(4000);
    setBin(b,fft_size);   //'fft_size' is maximum amplitude
   
    //why does this make my 4KHz signal cry?

    
    // Prepare for inverse FFT (undo previous sort).
    fft_ipermute(fft_buf, fft_size);
    
    //inverse the FFT operation
    ifft(fft_buf, fft_size);
  
    // Scale and copy to output buffer.
    i = 0;
    loop(fft_size,
  
      // Convert to real number.
      re = fft_buf[i * 2 + 0];
      im = fft_buf[i * 2 + 1];
      
      
      // Scale and add to output buffer.  only putting real numbers into output buffer
      out_buf[i] = re / fft_size;
      
      
    
      i += 1;
    );
  
    in_idx = out_idx = 0;
  );
  
  // Output sample from output buffer.
  
  out=out_buf[out_idx];
  spl0 = out;
  spl1 = out;
  
  out_idx += 1;

Last edited by weblordpepe; 07-17-2017 at 05:44 AM. Reason: More, clearer info
weblordpepe is offline   Reply With Quote
Old 07-17-2017, 02:02 PM   #2
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 7,897
Default

You might want to take a look here: -> https://forum.cockos.com/showthread.php?t=181209

-Michael

Last edited by mschnell; 07-17-2017 at 10:09 PM.
mschnell is offline   Reply With Quote
Old 07-18-2017, 03:56 PM   #3
weblordpepe
Human being with feelings
 
Join Date: Jul 2016
Posts: 29
Default

thank you for the link.

Man I really wish people would use proper variable names instead of single letters. Thats instrumental in why I can't follow along & learn what anything does.

I understand that this is using the real_fft function and not the regular FFT function. I don't understand enough to take what was in those examples & use it to fix my code.

Are you able to look at the windowing in my code and perhaps spot what is wrong?
weblordpepe is offline   Reply With Quote
Old 07-18-2017, 04:00 PM   #4
weblordpepe
Human being with feelings
 
Join Date: Jul 2016
Posts: 29
Default

Oh wow its your thread with lots of information. I'm just seeing page 2 onwards.

I second your complaint about documentation. Its been very difficult to learn even with buffers.

I would love a whole new site with documentation on everything, and an online context help where you can F1 statements in the editor.

I'm trying to use WDL-OL and iplug but its just a nightmare to follow along and nothing works out of the box.

I've ended up going over to using JUCE because its polished & has tutorials & documentation. But I don't like its licencing and 9mb default file sizes.

I really tries to use WDL and WDL-OL but I could never get anything to work.
weblordpepe is offline   Reply With Quote
Old 07-18-2017, 10:09 PM   #5
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 7,897
Default

With FFT (and real FFT, in fact "digital Fourier Transform" in general), a problem is that it does not behave as you would suppose on the first sight, even if coded correctly, because the the mathematical construct is not that obvious.

You first need to accept that with every Window, a periodic function from t=-inf to t=+inf is transformed to a periodic function from f=-inf to f=+inf (and vice versa), and that with real (e.g. usual <time domain> samples) input data (imaginary part = zero, or using real FFT), the output function is symmetrical (i.e. showing "aliasing" above half of the sampling frequency).

And this has nothing to do with coding or commenting, but just math.

Making use of the dummy imaginary part of the input samples by a stereo channel or by "real FFT" is yet another funny math thingy.

-Michael
mschnell is offline   Reply With Quote
Old 07-19-2017, 11:13 AM   #6
geraintluff
Human being with feelings
 
geraintluff's Avatar
 
Join Date: Nov 2009
Location: mostly inside my own head
Posts: 313
Default

1) Should there be some rounding in freqToBin()? You might want something like:

Code:
function freqToBin(f)(
  2*floor(f/srate*fft_size + 0.5); // round(x) = floor(x + 0.5)
);
Otherwise, it might end up pointing to an odd number (i.e. be the imaginary part of one of the bins, instead of the real part), which might break the index calculation you're doing in setBin().

2) If you're transforming incoming audio, yeah, you should be windowing. It's pretty important for getting a clean sound with no clicks or distortion (plus you need some padding).

3) If you're modifying audio, you also shouldn't expect the input audio to ever produce a frequency in a single bin, even if it's a pure sine oscillator. It's always a smudge of frequencies near the one you expect, so trying to cut just one frequency bin will leave some residue on either side. (This is related to windowing/padding).

Last edited by geraintluff; 07-22-2017 at 04:52 AM.
geraintluff is offline   Reply With Quote
Old 07-25-2017, 08:40 PM   #7
weblordpepe
Human being with feelings
 
Join Date: Jul 2016
Posts: 29
Default

Hi guys. Thanks for ya replies

Yeah regarding binToFrequency(), I accept that it should return a range. It does in my other code - it was just a kinda baby-toys first step to getting a conversion. Secondly yeah I also accept that the frequencies should be symmetrical. However the variables are still not descriptive. They make sense if you already know what they do. I don't - because I'm learning how to do it.

Reaper has quirks regarding its buffers and the documentation is very poor. I understand that it is a translation from amplitude domain to frequency domain within a certain range - and I lose all time information inside my 2048 samples.

What I need help with is the windowing.
Can someone please look at the code provided and tell me (or show me) how to implement windowing?

I am certain the windowing is my problem.


I just want to have a working simple FFT & IFFT for my learning. I just need a working example with some variable names I can read. I can take it from there. I've been coding for 20 years and I seldom ever ask for code. Please halp. I do not understand enough of REAPER or FFT to fix this on my own. I'm learning it on my own as much as I can - but I just need this help. I understand sliding windows I just can't apply it to this. I almost can but I need help.
weblordpepe is offline   Reply With Quote
Old 07-25-2017, 10:23 PM   #8
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 7,897
Default

Quote:
Originally Posted by weblordpepe View Post
Reaper has quirks regarding its buffers
I suppose you mean "EEL array accessing notation is different from what is known from many other languages".

Here, you are right. You need to accept that there is just a single (local) "basic" array of reals in a JSFX plugin, and you need to manage the "indexes" manually, if you need multiple arrays or complex numbers.

the "convenient" term a[b] just accesses the (basic) array element with index (a+b), and so it's the same as (a+b)[0] or 0[a+b].

complex values need to be accessed as

real_part = a[2*i];
imaginary part = a[2*i+1];

-Michael
mschnell is offline   Reply With Quote
Old 07-29-2017, 10:45 PM   #9
SaulT
Human being with feelings
 
Join Date: Oct 2013
Location: Seattle, WA
Posts: 766
Default

I know I have FFT code somewhere but I think it's on another laptop. I did post this, does it help at all?


https://forum.cockos.com/showthread.php?t=151974


If you look into the JSFX file on this one I have a couple of different window methods. I can't say it's great coding, but it does work.


https://forum.cockos.com/showthread.php?t=136765
SaulT is offline   Reply With Quote
Old 08-18-2017, 06:53 AM   #10
weblordpepe
Human being with feelings
 
Join Date: Jul 2016
Posts: 29
Default reaper's memory is retarded

Quote:
Originally Posted by mschnell View Post
I suppose you mean "EEL array accessing notation is different from what is known from many other languages".

Here, you are right. You need to accept that there is just a single (local) "basic" array of reals in a JSFX plugin, and you need to manage the "indexes" manually, if you need multiple arrays or complex numbers.

the "convenient" term a[b] just accesses the (basic) array element with index (a+b), and so it's the same as (a+b)[0] or 0[a+b].

complex values need to be accessed as

real_part = a[2*i];
imaginary part = a[2*i+1];

-Michael
Oh my god this is awful. How is this way of dealing with memory a help to anybody?
Thank you for explaining it.

a[b]= basicarray[a+b] = basicarray[0+a+b]

a[b] = 0[a+b]

This means if I want a 500 length array, and a 200 length array, I literally have to the following:

firstArray=0;
secondArray=500;

firstArray[499] = someValue; //this is safe to do
firstArray[500] = someValue; //this is not as it overrides secondArray's first value.

Basically the entire thing is this:
There's only one big array.

The index number in the brackets is added to the variableName's value to create the actual offset in this grand big arse array.

It really seems like for the hand-holding JSFX does with variables, being able to handle array memory locations should be something it does.

Thanks though I finally understand this. I might actually make a video about it - its that retarded.
weblordpepe is offline   Reply With Quote
Old 08-18-2017, 08:08 AM   #11
weblordpepe
Human being with feelings
 
Join Date: Jul 2016
Posts: 29
Default

Oh man this just isn't working. I need help. Can someone please demonstrate how to do windowing?

Its apparent that yes, my code does FFT into frequency space, then iFFT back into amplitude space. It does that. However I can't manipulate anything. I need to implement windowing somehow.

Please can somebody please just show me how to implement a kind of window using the code I've got provided? The code in other places has too much going on & I can't decipher whats happening - but I understand this code here.

Code:
desc:pepe simple FFT that works


@init

fft_size =4096;

fft_buf_size = fft_size * 2;
fft_buf = 0;


out_buf = fft_buf + fft_buf_size;      //why is this not zero like the fft_buf?
out_buf_size = fft_size;



memset(out_buf, 0, fft_size);    // During the first fft_size sample we don't have any output yet, so prefill the output buffer with silence.

in_idx = 0;
out_idx = 0;

@slider

@sample
  
  // Mono input.
  in = 0.5 * (spl0 + spl1);
  
  // Convert to complex number.
  re = in;
  im = 0.0;
  
  // Add to FFT buffer.
  fft_buf[in_idx * 2 + 0] = re;
  fft_buf[in_idx * 2 + 1] = im;

  in_idx += 1;
  
  // Do we have enough samples for FFT?
  in_idx >= fft_size ? (
  
    // FFT 
    fft(fft_buf, fft_size);
    
    //sort FFT output by frequency (permute?!)
    fft_permute(fft_buf, fft_size);
    
    
    ////  --- START OF GOOFY FUN TIMES IN FREQUENCY SPACE
    
    //any manipulations here just oscillate wildly depending on the fft size.
    // all i can do is a global attenuation otherwise it just bounces around.

    
    ///----- END OF GOOFY FUN TIMES. ONTO REVERSAL TO AMPLITUDE SPACE:
    
    
    // Prepare for inverse FFT (undo previous sort).
    fft_ipermute(fft_buf, fft_size);
    
    //inverse the FFT operation
    ifft(fft_buf, fft_size);
  
    // Scale and copy to output buffer.
    i = 0;
    loop(fft_size,
  
      // Convert to real number.
      re = fft_buf[i * 2 + 0];
      im = fft_buf[i * 2 + 1];
      
      
      // Scale and add to output buffer.  only putting real numbers into output buffer
      out_buf[i] = re / fft_size;
    
      i += 1;
    );
  
    in_idx = 0;
    out_idx = 0;
  );
  
  // Output sample from output buffer.
  spl0 =out_buf[out_idx];
  spl1=spl0; //mono
  out_idx += 1;  //increment index for next time
weblordpepe is offline   Reply With Quote
Old 08-24-2017, 04:50 AM   #12
Aesis
Human being with feelings
 
Join Date: Jan 2011
Posts: 401
Default

Quote:
Originally Posted by weblordpepe View Post
Oh man this just isn't working. I need help. Can someone please demonstrate how to do windowing?

Its apparent that yes, my code does FFT into frequency space, then iFFT back into amplitude space. It does that. However I can't manipulate anything. I need to implement windowing somehow.

Please can somebody please just show me how to implement a kind of window using the code I've got provided? The code in other places has too much going on & I can't decipher whats happening - but I understand this code here.

Code:
desc:pepe simple FFT that works


@init

fft_size =4096;

fft_buf_size = fft_size * 2;
fft_buf = 0;


out_buf = fft_buf + fft_buf_size;      //why is this not zero like the fft_buf?
out_buf_size = fft_size;



memset(out_buf, 0, fft_size);    // During the first fft_size sample we don't have any output yet, so prefill the output buffer with silence.

in_idx = 0;
out_idx = 0;

@slider

@sample
  
  // Mono input.
  in = 0.5 * (spl0 + spl1);
  
  // Convert to complex number.
  re = in;
  im = 0.0;
  
  // Add to FFT buffer.
  fft_buf[in_idx * 2 + 0] = re;
  fft_buf[in_idx * 2 + 1] = im;

  in_idx += 1;
  
  // Do we have enough samples for FFT?
  in_idx >= fft_size ? (
  
    // FFT 
    fft(fft_buf, fft_size);
    
    //sort FFT output by frequency (permute?!)
    fft_permute(fft_buf, fft_size);
    
    
    ////  --- START OF GOOFY FUN TIMES IN FREQUENCY SPACE
    
    //any manipulations here just oscillate wildly depending on the fft size.
    // all i can do is a global attenuation otherwise it just bounces around.

    
    ///----- END OF GOOFY FUN TIMES. ONTO REVERSAL TO AMPLITUDE SPACE:
    
    
    // Prepare for inverse FFT (undo previous sort).
    fft_ipermute(fft_buf, fft_size);
    
    //inverse the FFT operation
    ifft(fft_buf, fft_size);
  
    // Scale and copy to output buffer.
    i = 0;
    loop(fft_size,
  
      // Convert to real number.
      re = fft_buf[i * 2 + 0];
      im = fft_buf[i * 2 + 1];
      
      
      // Scale and add to output buffer.  only putting real numbers into output buffer
      out_buf[i] = re / fft_size;
    
      i += 1;
    );
  
    in_idx = 0;
    out_idx = 0;
  );
  
  // Output sample from output buffer.
  spl0 =out_buf[out_idx];
  spl1=spl0; //mono
  out_idx += 1;  //increment index for next time

Windowing is simple, just multiply the signal, exactly like fades. Technically you are already windowing, using the "rectangular" window. Anyway the start of the buffer you would multiply by close to zero, and the end, the same. In the middle, one. Fade in, fade out...
Aesis is offline   Reply With Quote
Old 08-24-2017, 05:08 AM   #13
mschnell
Human being with feelings
 
mschnell's Avatar
 
Join Date: Jun 2013
Location: Krefeld, Germany
Posts: 7,897
Default

To get decent results you need to do "overlapping" windowing.

You collect the count of an FFT window in a buffer, then apply the windowing function (e.g. "sinc") onto same and then apply FFT.

You use two buffers with a distance of half the FFT window size for the same algorithm.

By that, the middle half (1/4 .. 3/4) of each buffer contains the (most) relevant information, while the rest is allowed to be dedicated to side-effects of the FFT (that always does virtually periodic signals).

After FFT and iFFT the buffers are windowed again and added appropriately to get a decent result.

There are examples for this among the stock JSFXes of Reapers (such as FFTSplitter).

-Michael

Last edited by mschnell; 08-24-2017 at 10:09 PM.
mschnell is offline   Reply With Quote
Old 08-24-2017, 03:53 PM   #14
jcjr
Human being with feelings
 
Join Date: Dec 2015
Location: SE TN USA
Posts: 77
Default

Quote:
Originally Posted by weblordpepe View Post
Oh my god this is awful. How is this way of dealing with memory a help to anybody?
Thank you for explaining it.

a[b]= basicarray[a+b] = basicarray[0+a+b]

a[b] = 0[a+b]

This means if I want a 500 length array, and a 200 length array, I literally have to the following:

firstArray=0;
secondArray=500;

firstArray[499] = someValue; //this is safe to do
firstArray[500] = someValue; //this is not as it overrides secondArray's first value.

Basically the entire thing is this:
There's only one big array.
Hi Pepe

Many people have made memory manager functions for js. I didn't "invent" the following, but made it stripped-down simple for my own uses. Some of the js memory manager codes I saw are more elaborate, fancier than I have needed so far. I like the KISS motto "keep it simple stupid!"

Code:
@init
  //should only call js_malloc within js @init section
  //js_malloc() preferably the first function in @init
  //example usage--
  //  Array_1 = js_malloc(500, JS_MALLOC_NO_CLEARMEM);
  //    malloc a 500 element array and don't clear the memory
  //  a = 20;
  //  b = 10;
  //  Array_2_Size = a * b; //demonstrate calculate array size
  //  Array_2 = js_malloc(Array_2_Size, JS_MALLOC_CLEARMEM);
  //    malloc a calculated-size array and clear the memory
  //  n = 41;
  //  Array_1[n] = n; //various silly trivial assignments
  //  Array_2[n + 1] = Array_1[n] * $pi;

  JS_MALLOC_NO_CLEARMEM = 0;
  JS_MALLOC_CLEARMEM = 1;
  next_js_malloc = 0;  
  function js_malloc(a_size, a_clearmem)
  local (l_result)
  (
    l_result = next_js_malloc;
    (a_clearmem != 0) ?
      memset(l_result, 0, a_size);
    next_js_malloc += a_size;
    l_result;
  );
Though computer hardware mem management is "fancy" nowadays, comceptually any computer memory is just a huge array of consecutive values. Memory management functions impose "order" on the wilderness of bytes, though one can't assume exactly how the OS and hardware are cooperating to allocate memory resources.

There may be ways to safely call a user-written js malloc function from other places than in @init, haven't studied it. My uses so far are setting things up in @init, so haven't needed to investigate the complications of calling it from other code sections.

Calling from multiple code sections could mean calling from different threads, and maybe two threads would unfortunately decide to alloc memory the same time, causing a train wreck. But it could probably be worked-around if necessary.

You could add realloc (resize block) and dealloc (free memory) functions and more. I don't need it so far. Usually @init gets called rather frequently especially if sample rate or whatever would change. If arrays need resizing because of changed samplerate then any allocation code in @init would "run from scratch" on every major change. Negating so far my needs to "on the fly" resize memory blocks, delete them, or compact the js memory heap after such memory-churning actions.

I usually remember the size of an array every time I alloc the array. But if you need to churn the memory, or want a convenient way to automatically know how big each created array happens to be--

Maybe there are smarter ways but I remember from long ago some memory managers would encode the size of each memory block in the first part of the array. If you called alloc, it would get the memory block from OS or hardware sized as RequestedSize + SizeOfInt32 (or on 64 bit system, RequestedSize + SizeOfInt64).

Because js is all float doubles, it would just be RequestedSize + 1.

Anyway, a modified js malloc remembering the size of each block would alloc RequestedSize + 1, save the RequestedSize in NewArray[0], then return NewArray + 1 as the new malloc'd array "pointer".

So if you malloc a 500 element array, it would alloc 501 elements, save the value 500 in My_Array[0], and return My_Array + 1 as the alloc result.

So you would address My_Array[0] thru MyArray[499] as expected, but if you want to know how big the array is (or you want to write realloc, dealloc and garbage collection), you can find the size of each array at (My_Array - 1)[0]. Or maybe just My_Array[-1]. I don't recall offhand if js is necessarily always friendly to negative array indices.

A possibly buggy js malloc remembering each array's size--
Code:
@init
  //should only call js_malloc within js @init section
  //js_malloc() preferably the first function in @init
  JS_MALLOC_NO_CLEARMEM = 0;
  JS_MALLOC_CLEARMEM = 1;
  next_js_malloc = 0;  
  function js_malloc(a_size, a_clearmem)
  local (l_result)
  (
    l_result = next_js_malloc;
    (a_clearmem != 0) ?
      memset(l_result, 0, a_size + 1);
    next_js_malloc += (a_size + 1);
    l_result[0] = a_size; //save the size of array in element[0]
    l_result += 1;
    l_result; //returns the "pointer" to memory right after the size field
  );
So that's enough dumbness for one message, maybe there are far superior ways to do it.
jcjr 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 06:46 AM.


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