1

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 hereenter image description here
Improve this question

3 Answers 3

Reset to default
3

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

plot of 30 samples of a 10kHz signal sampled at 44.1hHz

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

plot of FFT amplitudes for a 10kHz signal sampled at 44.1hHz

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)
Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

>> Is it possible it contains more than one channel (left & right)? No only one channel. >>Is it realy signed 16 bit integers? Yes it is
Hi @meta4 : I think my problem is with the selection of 'n_sa' , 'T" , 'f' , 't_fft'
@Vishnu CS, Did you try plotting 30 samples of the raw (non-fft-ed) data? Can you see a "periodic" signature?
2
+25

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

Improve this answer

2 Comments

Ok .. I will check this whether it resolve my problem
Actually two issues were there. I answered my question. Thank so much
2

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 here

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

Improve this answer

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.