FFT spectrogram in python

Question

I have python 3.4.

I transmitted a 2MHz (for example) frequency and received the cavitation over the time (until I stopped the measurement). I want to get a spectrogram (cavitation vs frequency) and more interesting is a spectrogram of cavitation over the time of the sub-harmonic (1MHz) frequency.

The data is saved in sdataA (=cavitation), and t (=measurement time)

I tried to save fft in FFTA

FFTA = np.array([])
FFTA = np.fft.fft(dataA)
FFTA = np.append(FFTA, dataA)

I got real and complex numbers Then I took only half (from 0 to 1MHz) and save the real and complex data.

nA = int(len(FFTA)/2)
yAre = FFTA[range(nA)].real
yAim = FFTA[range(nA)].imag

I tried to get the frequencies by:

FFTAfreqs = np.fft.fftfreq(len(yAre))

But it is totally wrong (I printed the data by print (FFTAfreqs))

I also plotted the data and again it's wrong:

plt.plot(t, FFTA[range(n)].real, 'b-', t, FFTA[range(n)].imag, 'r--')
plt.legend(('real', 'imaginary'))
plt.show()

How can I output a spectrogram of cavitation over the time of the sub-harmonic (1MHz) frequency?

EDIT:

Data example:

see a sample of 'dataA' and 'time':

dataA = [6.08E-04,2.78E-04,3.64E-04,3.64E-04,4.37E-04,4.09E-04,4.49E-04,4.09E-04,3.52E-04,3.24E-04,3.92E-04,3.24E-04,2.67E-04,3.24E-04,2.95E-04,2.95E-04,4.94E-04,4.09E-04,3.64E-04,3.07E-04]
time = [0.00E+00,4.96E-07,9.92E-07,1.49E-06,1.98E-06,2.48E-06,2.98E-06,3.47E-06,3.97E-06,4.46E-06,4.96E-06,5.46E-06,5.95E-06,6.45E-06,6.94E-06,7.44E-06,7.94E-06,8.43E-06,8.93E-06,9.42E-06]

EDIT II: From @Martin example I tried the following code, please let me know if I did it right.

In the case that dataA and Time are saved as h5 files (or the data that I posted already)

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

dfdata = pd.read_hdf("C:\\data_python\\DataA.h5")
dft = pd.read_hdf("C:\\data_python\\time.h5")
dft_cor = int((len(dft)-2)*4.96E-6)   # calculating the measured time

fs = 2000000 #sampling frequency 2MHz
CHUNK = 10000
signal_time = dft_cor # seconds

def sine(freq,fs,secs):
    data=dfdata
    wave = np.sin(freq*2*np.pi*data)
    return wave

a1 = sine(fs,fs,120)
a2 = sine(fs/2,fs,120)

signal = a1+a2
afft = np.abs(np.fft.fft(signal[0:CHUNK]))
freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
spectrogram_chunk = freqs/np.amax(freqs*1.0)

# Plot spectral analysis
plt.plot(freqs[0:1000000],afft[0:1000000]) # 0-1MHz
plt.show()

number_of_chunks = 1000

# Empty spectrogram
Spectrogram = np.zeros(shape = [CHUNK,number_of_chunks])

for i in range(number_of_chunks):
    afft = np.abs(np.fft.fft(signal[i*CHUNK:(1+i)*CHUNK]))
    freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
    spectrogram_chunk = afft/np.amax(afft*1.0)
    try:
        Spectrogram[:,i]=spectrogram_chunk
    except:
        break


import cv2
Spectrogram = Spectrogram[0:1000000,:]
cv2.imshow('spectrogram',np.uint8(255*Spectrogram/np.amax(Spectrogram)))
cv2.waitKey()
cv2.destroyAllWindows()

Answer 1

It seems your problem is not in Python but in understanding what is Spectrogram.

---

Spectrogram is sequences of spectral analysis of a signal.

1) You need to cut your signal in CHUNKS.

2) Do spectral analysis of these CHUNKS and stick it together.

Example:

You have 1 second of audio recoding (44100 HZ sampling). That means the recording will have 1s * 44100 -> 44100 samples. You define CHUNK size = 1024 (for example).

For each chunk you will do FFT, and stick it together into 2D matrix (X axis - FFT of the CHUNK, Y axis - CHUNK number,). 44100 samples / CHUNK ~ 44 FFTs, each of the FFT covers 1024/44100~0.023 seconds of the signal

The bigger the CHUNK, the more accurate Spectrogram is, but less 'realtime'.

The smaller the CHUNK is, the less acurate is the Spectrogram, but you have more measurements as you measure frequencies 'more often'.

If you need 1MHZ - actually you cannot use anything higher than 1MHZ, you just take half of the resulting FFT array - and it doesnt matter which half, because 1MHZ is just the half of your sampling frequency, and the FFT is mirroring anything that is higher than 1/2 of sampling frequency.

---

About FFT, you dont want complex numbers. You want to do

FFT = np.abs(FFT) # Edit - I just noticed you use '.real', but I will keep it here

because you want real numbers.

Preparation for Spectrogram - example of Spectrogram

Audio Signal with 150HZ wave and 300HZ Wave

import numpy as np
import matplotlib.pyplot as plt


fs = 44100#sampling frequency
CHUNK = 10000
signal_time = 20 # seconds

def sine(freq,fs,secs):
    data=np.arange(fs*secs)/(fs*1.0)
    wave = np.sin(freq*2*np.pi*data)
    return wave

a1 = sine(150,fs,120)
a2 = sine(300,fs,120)

signal = a1+a2
afft = np.abs(np.fft.fft(signal[0:CHUNK]))
freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
spectrogram_chunk = freqs/np.amax(freqs*1.0)

# Plot spectral analysis
plt.plot(freqs[0:250],afft[0:250])
plt.show()

number_of_chunks = 1000

# Empty spectrogram
Spectrogram = np.zeros(shape = [CHUNK,number_of_chunks])

for i in range(number_of_chunks):
    afft = np.abs(np.fft.fft(signal[i*CHUNK:(1+i)*CHUNK]))
    freqs = np.linspace(0,fs,CHUNK)[0:int(fs/2)]
    #plt.plot(spectrogram_chunk[0:250],afft[0:250])
    #plt.show()
    spectrogram_chunk = afft/np.amax(afft*1.0)
    #print(signal[i*CHUNK:(1+i)*CHUNK].shape)
    try:
        Spectrogram[:,i]=spectrogram_chunk
    except:
        break


import cv2
Spectrogram = Spectrogram[0:250,:]
cv2.imshow('spectrogram',np.uint8(255*Spectrogram/np.amax(Spectrogram)))
cv2.waitKey()
cv2.destroyAllWindows()

Spectral analysis of single CHUNK

Spectrogram