How to split range in frequency domain and inverse FT?

[larsien]

Hi! First of all, thank you for your fast library for FT.

I tried to split the range of frequencies and perform an inverse FT. In detail, my EEG data is sampled at 200Hz per second. Here are the steps I followed:

1. Convert from time domain to frequency domain using STFT
2. split the frequencies into three ranges: 1~4Khz, 4~16Khz, 16~64Khz
3. Perform an Inverse FT

As far as I know, a spectrogram is Frequency + Amplitude(or can be Magnitude) + Time. When I ran the Spectrogram example, I couldn't understand the result:

        final chunkSize = 300;
        final stft = STFT(chunkSize, Window.hanning(chunkSize));
        final spectrogram = <Float64List>[];
        stft.run(d, (Float64x2List freq) {
          spectrogram.add(freq.discardConjugates().magnitudes());
        });

The result of the spectrogram is a list of 151 double array and etc

stft = {STFT} 
 _fft = {CompositeFFT} CompositeFFT(300)
 _win = {_Float64List} size = 300
 _chunk = {_Float64x2List} size = 300

Where can I find frequency information from this spectrogram? and how to split?

[liamappelbe]

The frequency is implicit in the indices of the array. In other words, in stft.run(d, (Float64x2List freq) {, the freq is a an array of amplitude+phase info indexed by the frequency (tho the index is not in Hz). To convert a freq list index to a frequency in Hz, you can use the stft.frequency() method. So the frequency of freq[i] is stft.frequency(i, sampleRate) (and in your case sampleRate would be 200000).

To split freq at the bands you mentioned, you can just split the freq array at the indices corresponding to those frequencies. Try something like this (I haven't run this code, so it may have bugs :P):

const int sampleRate = 200000;
final int indexOf1kHz = stft.indexOfFrequency(1000, sampleRate).round();
final int indexOf4kHz = stft.indexOfFrequency(4000, sampleRate).round();
final int indexOf16kHz = stft.indexOfFrequency(16000, sampleRate).round();
final int indexOf64kHz = stft.indexOfFrequency(64000, sampleRate).round();

stft.run(d, (Float64x2List freq) {
  final realFreq = freq.discardConjugates();
  spectrogram1to4.add(Float64x2List.sublistView(realFreq, indexOf1kHz, indexOf4kHz).magnitudes());
  spectrogram4to16.add(Float64x2List.sublistView(realFreq, indexOf4kHz, indexOf16kHz).magnitudes());
  spectrogram16to64.add(Float64x2List.sublistView(realFreq, indexOf16kHz, indexOf64kHz).magnitudes());
});

I'm assuming here that you want 3 separate spectrograms, spectrogram1to4, spectrogram4to16, and spectrogram16to64. As I was writing that code, I realized I should have a function that's the inverse of stft.frequency. So I've added stft.indexOfFrequency. You can update to 1.3.1 to get it.

However, if you want to do an inverse FFT, STFT might not be the best choice. Inverse FFT can perfectly reverse an FFT, but not an STFT. See how you go, but if your results aren't what you expect, try doing an FFT instead of an STFT. Also, if you want to do an inverse FFT, you'll need to retain the phase info, so don't use magnitudes() to discard the phase info. If you give me some more info about what you're trying to do (what you're using the spectrogram for, and what you're using the inverse of the frequency bands for), maybe I can give you some more help.

[larsien]

Thank you for your answer above. Actually I am trying to draw a graph of Brain EEG signal and want to analyze features of them in realtime or near-real time. For example, in a scenario where the user blinks their eye, then graph will check the timing in real time domain.

Each frequency band may have a feature about 'eye blinking, close, open', + 'various noise - muscular moving, electronic power etc' So I want to split it into frequency domain and restore it to time domain again to identify the moment of eye blinking.

And I want to fix my wrong information before. the frequency ranges are 0 ~ 4Hz, 4 ~ 8Hz, 8 ~ 12Hz, 15 ~ 20Hz. Not kHz.

[liamappelbe]

Ah ok. Yeah, for real time you do need to use STFT, because it breaks the signal up into chunks. FFT works on the entire signal. At the moment my API doesn't support streaming, but I can add that easily.

For frequencies that low, it's important to keep in mind that your window size limits the frequencies you can analyse. For example, if your window size is 1 second, the lowest frequency you can measure is 1Hz (otherwise the wave won't be able to complete a cycle within the window). Of course, there's a trade off with the time resolution of your windows.

You mentioned that the sampling rate is 200kHz. Is that also actually 200Hz, like the others?

As for the inverse FFT, I can see 3 fiddly bits to think about. I'll outline them now, but we can talk about them in more detail once you get the spectrogram working.

1. Instead of "splitting" the freq array, you probably need to copy it for each frequency band, and zero out the elements you don't care about. Otherwise the inverse FFT will be wrong. You also shouldn't take the magnitude or discard the redundant conjugates.
2. The windowing function means there's a lot of overlap between adjacent chunks, so you have to overlap the time domain output of the inverse FFT by the same amount.
3. I haven't tried inverse FFTing an STFT before. When I say it probably won't perfectly invert, the way that would manifest is that the time domain of each chunk might not join up properly. But see how you go. It might be fine ;)

[larsien]

Yes, Sampling rate is 200 Hz. 200 sampling in 1 second.

I tried to implement your comments for 0~4Hz IFFT. It seems working. But I don't 100% make sure. Could you check this?

    double frequencySampleRate = 200;
    final stft = STFT(chunkSize, Window.hanning(chunkSize));
    List<Float64x2> spectrogram0to4 ;
    Float64x2List spectrogram4to8 ;
    Float64x2List spectrogram8to12 ;

    final int indexOf4Hz = stft.indexOfFrequency(4, frequencySampleRate).round(); //4 in this case
    final int indexOf8Hz = stft.indexOfFrequency(8, frequencySampleRate).round(); //8 in this case
    final int indexOf12Hz = stft.indexOfFrequency(12, frequencySampleRate).round();//12 in this case

    int timeIndex = 0;

    final fft = FFT(frequencySampleRate.round());
    final ifft0To4HzList = <Float64List>[];
    // final ifft4To8HzList = <Float64List>[];
    // final ifft8To12HzList = <Float64List>[];

    stft.run(d, (Float64x2List freq){
      // final realFreq = freq.discardConjugates();//shouldn't discard conjugates
      spectrogram0to4 = Float64x2List.sublistView(freq, 0, indexOf4Hz).toList();
      spectrogram0to4.addAll(Float64x2List(freq.length - indexOf4Hz).toList());//add 0 as padding

      ifft0To4HzList.add(fft.realInverseFft(Float64x2List.fromList(spectrogram0to4)));
      spectrogram0to4.clear();
   }

[liamappelbe]

Hmm. Close. There's a small problem tho. I'm not sure what your knowledge level about fourier transforms is, so I'll have to explain some basics. Sorry if you already know this stuff.

FFT is a fast algorithm for doing discreet fourier transforms (DFT). A DFT takes a complex time domain signal and returns the complex frequency domain representation of it. If you run a DFT on a real time domain signal, you still get a complex frequency domain output, but the result must necessarily be around 50% redundant data, because the real input has half as many degrees of freedom as the complex output. It turns out that the redundancy is that the upper half of the frequency domain output is a reflected and conjugated copy of the lower half: [sum term, ...terms..., nyquist term, ...reflected and conjugated terms...]. The sum term (aka 0Hz) is real, and the nyquist term (half your sample rate) is also real, and neither of these are duplicated. Further reading.

It's perfectly reasonable to mess with the frequency representation and then inverse FFT, like you are here. But those redundant conjugates are important for the inverse FFT. If you invert a frequency domain signal that is not conjugate symmetric you will get a complex time domain signal, which is not what you want in your case. Your code has 0s instead of the conjugates, so the output of the inverse FFT will be complex.

I added the discardConjugates() function to get rid of that redundant data. I've been meaning to add another function that fabricates redundant conjugates, so that you can discard the conjugates, mess with the frequency representation, then reconstruct the conjugates, and inverse FFT. But actually I think in your case you can do something simpler.

Your frequency domain signal looks like this: [sum (0Hz) term, ..., 4Hz term, ..., 8Hz term, ..., nyquist (100Hz) term, ..., conjugate of 8Hz term ..., conjugate of 4Hz term, ...] So for the 0 to 4Hz spectrogram, you can just make a copy of freq and set all the terms to 0 from indexOf4Hz + 1 to freq.length - indexOf4Hz - 1, inclusive. For the others you can do something similar, but any time you 0 out term i, you also need to 0 out term freq.length - i (unless i is 0).

IMPORTANT: You will know if you've done this correctly because the time domain output of the inverse FFT will be real valued (or close to it). If you write some code to find the largest imaginary part of the output, it should be 0 or very close to 0.

[larsien]

Let me study one more time about FFT, and your explanation! Really appreciate