FFT plot of raw PCM comes wrong for higher frequency in python

Question

Here I am using fft function of numpy to plot the fft of PCM wave generated from a 10000Hz sine wave. But the amplitude of the plot I am getting is wrong.

The frequency is coming correct using fftfreq function which I am printing in the console itself. My python code is here.

import numpy as np
import matplotlib.pyplot as plt 


frate = 44100
filename = 'Sine_10000Hz.bin' #signed16 bit PCM of a 10000Hz sine wave 


f = open('Sine_10000Hz.bin','rb') 
y = np.fromfile(f,dtype='int16') #Extract the signed 16 bit PCM value of 10000Hz Sine wave
f.close()

####### Spectral Analysis #########

fft_value = np.fft.fft(y)
freqs = np.fft.fftfreq(len(fft_value))  # frequencies associated with the coefficients:
print("freqs.min(), freqs.max()")

idx = np.argmax(np.abs(fft_value))  # Find the peak in the coefficients
freq = freqs[idx]
freq_in_hertz = abs(freq * frate)
print("\n\n\n\n\n\nfreq_in_hertz")
print(freq_in_hertz)


for i in range(2):
    print("Value at index {}:\t{}".format(i, fft_value[i + 1]), "\nValue at index {}:\t{}".format(fft_value.size -1 - i, fft_value[-1 - i]))

#####


n_sa = 8 * int(freq_in_hertz)
t_fft = np.linspace(0, 1, n_sa)
T = t_fft[1] - t_fft[0]  # sampling interval 
N = n_sa   #Here it is n_sample
print("\nN value=")
print(N)

# 1/T = frequency
f = np.linspace(0, 1 / T, N)

plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")
plt.xlim(0,15000)
# 2 / N is a normalization factor  Here second half of the sequence gives us no new information that the half of the FFT sequence is the output we need.
plt.bar(f[:N // 2], np.abs(fft_value)[:N // 2] * 2 / N, width=15,color="red") 

Output comes in the console (Only minimal prints I am pasting here)

freqs.min(), freqs.max()
-0.5 0.49997732426303854






freq_in_hertz
10000.0
Value at index 0:       (19.949569768991054-17.456031216294324j)
Value at index 44099:   (19.949569768991157+17.45603121629439j)
Value at index 1:       (9.216783424692835-13.477631008179145j)
Value at index 44098:   (9.216783424692792+13.477631008179262j)

N value=
80000

The frequency extraction is coming correctly but in the plot something I am doing is incorrect which I don't know.

Updating the work:

1. When I am change the multiplication factor 10 in the line n_sa = 10 * int(freq_in_hertz) to 5 gives me correct plot. Whether its correct or not I am not able to understand
2. In the line plt.xlim(0,15000) if I increase max value to 20000 again is not plotting. Till 15000 it is plotting correctly.
3. I generated this Sine_10000Hz.bin using Audacity tool where I generate a sine wave of freq 10000Hz of 1sec duration and a sampling rate of 44100. Then I exported this audio to signed 16bit with headerless (means raw PCM). I could able to regenerate this sine wave using this script. Also I want to calculate the FFT of this. So I expect a peak at 10000Hz with amplitude 32767. You can see i changed the multiplication factor 8 instead of 10 in the line n_sa = 8 * int(freq_in_hertz). Hence it worked. But the amplitude is showing incorrect. I will attach my new figure here

Answer 1

I'm not sure exactly what you are trying to do, but my suspicion is that the Sine_10000Hz.bin file isn't what you think it is.

Is it possible it contains more than one channel (left & right)?

Is it realy signed 16 bit integers?

It's not hard to create a 10kHz sine wave in 16 bit integers in numpy.

import numpy as np
import matplotlib.pyplot as plt

n_samples = 2000
f_signal = 10000 # (Hz) Signal Frequency
f_sample = 44100 # (Hz) Sample Rate
amplitude = 2**3 # Arbitrary. Must be > 1. Should be > 2. Larger makes FFT results better
time = np.arange(n_samples) / f_sample # sample times
# The signal
y = (np.sin(time * f_signal * 2 * np.pi) * amplitude).astype('int16')

If you plot 30 points of the signal you can see there are about 5 points per cycle.

plt.plot(time[:30], y[:30], marker='o')
plt.xlabel('Time (s)')
plt.yticks([]); # Amplitude value is artificial. hide it

If you plot 30 samples of the data from Sine_10000Hz.bin does it have about 5 points per cycle?

This is my attempt to recreate the FFT work as I understand it.

fft_value = np.fft.fft(y)                          # compute the FFT
freqs = np.fft.fftfreq(len(fft_value)) * f_sample  # frequencies for each FFT bin

N = len(y)
plt.plot(freqs[:N//2], np.abs(fft_value[:N//2]))
plt.yscale('log')
plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")

I get the following plot

The y-axis of this plot is on a log scale. Notice that the amplitude of the peak is in the thousands. The amplitude of most of the rest of the data points are around 100.

idx_max = np.argmax(np.abs(fft_value))  # Find the peak in the coefficients
idx_min = np.argmin(np.abs(fft_value))  # Find the peak in the coefficients
print(f'idx_max = {idx_max}, idx_min = {idx_min}')
print(f'f_max = {freqs[idx_max]}, f_min = {freqs[idx_min]}')
print(f'fft_value[idx_max] {fft_value[idx_max]}')
print(f'fft_value[idx_min] {fft_value[idx_min]}')

produces:

idx_max = 1546, idx_min = 1738
f_max = -10010.7, f_min = -5777.1
fft_value[idx_max] (-4733.232076236707+219.11718299533203j)
fft_value[idx_min] (-0.17017443966211232+0.9557200531465061j)

Answer 2

I'm adding a link to a script I've build that outputs the FFT with ACTUAL amplitude (for real signals - e.g. your signal). Have a go and see if it works:

dt=1/frate in your constellation....

https://stackoverflow.com/a/53925342/4879610

Answer 3

After a long home work I could able to find my issue. As I mentioned in the Updating the work: the reason was with the number of samples which I took was wrong.

I changed the two lines in the code

n_sa = 8 * int(freq_in_hertz)
t_fft = np.linspace(0, 1, n_sa)

to

n_sa = y.size //number of samples directly taken from the raw 16bits
t_fft = np.arange(n_sa)/frate //Here we need to divide each samples by the sampling rate

This solved my issue.

My spectral output is

Special thanks to @meta4 and @YoniChechik for giving me some suggestions.