how to plot fft of signal with correct frequencies on x-axis?

Create the signal (a sine wave) using numpy. Compute the one-dimensional discrete Fourier Transform.,Set the figure size and adjust the padding between and around the subplots.,To plot the FFT (Fast Fourier Transform) of a signal with correct frequencies on the X-axis in matplotlib, we can take the following steps −,Return the Discrete Fourier Transform sample frequencies.

Example

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["figure.figsize"] = [7.00, 3.50]
plt.rcParams["figure.autolayout"] = True

N = 256
t = np.arange(N)
m = 4
nu = float(m) / N

signal = np.sin(2 * np.pi * nu * t)

ft = np.fft.fft(signal)

freq = np.fft.fftfreq(N)

plt.plot(freq, ft.real ** 2 + ft.imag ** 2)

plt.show()

Suggestion : 2

When you use welch, the returned frequency and power vectors are not sorted in ascending frequency order. You must fftshift the output before you plot. Also, the sample frequency you pass welch must be a float. Make sure to use the scaling="spectrum" option to get the power instead of the power density. To get the correct power value, you also need to scale the power to account for the windowing effect. For a hann window, you need to divide by 1.5. Here is a working example:

from scipy.signal
import welch
from numpy.fft
import fftshift, fft
from numpy
import arange, exp, pi, mean
import matplotlib.pyplot as plt

#generate a 1000 Hz complex tone singnal
sample_rate = 48000. #sample rate must be a float
t = arange(1024 * 1024) / sample_rate
signal = exp(1 j * 2000 * pi * t)
#amplitude correction factor
corr = 1.5

#calculate the psd with welch
sample_freq, power = welch(signal, fs = sample_rate, window = "hann", nperseg = 256, noverlap = 128, scaling = 'spectrum')

#fftshift the output
sample_freq = fftshift(sample_freq)
power = fftshift(power) / corr

#check that the power sum is right
print sum(power)

plt.figure(figsize = (9.84, 3.94))
plt.plot(sample_freq, power)
plt.xlabel("Frequency (MHz)")
plt.ylabel("Relative power (dB)")
plt.show()

Suggestion : 3

The fftfreq function tells us which frequencies we are looking at specifically:,Compute the DFT of a real sequence, exploiting the symmetry of the resulting spectrum for increased performance. Also used by fft internally when applicable.,Let’s examine the frequency components in a noisy image (Figure 4-7). Note that, while a static image has no time-varying component, its values do vary across space. The DFT applies equally to either case.,There is, and it is called the discrete Fourier transform, or DFT, where discrete refers to the recording consisting of time-spaced sound measurements, in contrast to a continual recording as, for example, on magnetic tape (remember cassettes?). The DFT is often computed using the FFT algorithm, a name informally used to refer to the DFT itself. The DFT tells us which frequencies or “notes” to expect in our signal.

We’ll start by setting up some plotting styles and importing the usual suspects:

# Make plots appear inline, set custom plotting style
   %
   matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

# Make plots appear inline, set custom plotting style
%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('style/elegant.mplstyle')

import numpy as np

The discrete¹ Fourier transform (DFT) is a mathematical technique used to convert temporal or spatial data into frequency domain data. Frequency is a familiar concept, due to its colloquial occurrence in the English language: the lowest notes your headphones can rumble out are around 20 Hz, whereas middle C on a piano lies around 261.6 Hz; Hertz, or oscillations per second, in this case literally refers to the number of times per second at which the membrane inside the headphone moves to-and-fro. That, in turn, creates compressed pulses of air which, upon arrival at your eardrum, induces a vibration at the same frequency. So, if you take a simple periodic function, sin(10 × 2πt), you can view it as a wave:

f = 10 # Frequency, in cycles per second, or Hertz
f_s = 100 # Sampling rate, or number of measurements per second

t = np.linspace(0, 2, 2 * f_s, endpoint = False)
x = np.sin(f * 2 * np.pi * t)

fig, ax = plt.subplots()
ax.plot(t, x)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Signal amplitude');

from scipy import fftpack

X = fftpack.fft(x)
freqs = fftpack.fftfreq(len(x)) * f_s

fig, ax = plt.subplots()

ax.stem(freqs, np.abs(X))
ax.set_xlabel('Frequency in Hertz [Hz]')
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlim(-f_s / 2, f_s / 2)
ax.set_ylim(-5, 110)

(-5, 110)

Listen to this snippet of the nightingale birdsong (released under CC BY 4.0):

from IPython.display
import Audio
Audio('data/nightingale.wav')

We load the audio file, which gives us the sampling rate (number of measurements per second) as well as audio data as an (N, 2) array—two columns because this is a stereo recording.

from scipy.io
import wavfile

rate, audio = wavfile.read('data/nightingale.wav')

We convert to mono by averaging the left and right channels.

audio = np.mean(audio, axis = 1)

N = audio.shape[0]
L = N / rate

print(f'Audio length: {L:.2f} seconds')

f, ax = plt.subplots()
ax.plot(np.arange(N) / rate, audio)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Amplitude [unknown]');

Audio length: 7.67 seconds

Start by chopping up the signal into slices of 1024 samples, each slice overlapping the previous by 100 samples. The resulting slices object contains one slice per row.

from skimage
import util

M = 1024

slices = util.view_as_windows(audio, window_shape = (M, ), step = 100)
print(f 'Audio shape: {audio.shape}, Sliced audio shape: {slices.shape}')

Let us illustrate:

import time

from scipy
import fftpack
from sympy
import factorint

K = 1000
lengths = range(250, 260)

# Calculate the smoothness
for all input lengths
smoothness = [max(factorint(i).keys()) for i in lengths]

exec_times = []
for i in lengths:
   z = np.random.random(i)

# For each input length i, execute the FFT K times
# and store the execution time

times = []
for k in range(K):
   tic = time.monotonic()
fftpack.fft(z)
toc = time.monotonic()
times.append(toc - tic)

# For each input length, remember the * minimum * execution time
exec_times.append(min(times))

f, (ax0, ax1) = plt.subplots(2, 1, sharex = True)
ax0.stem(lengths, np.array(exec_times) * 10 ** 6)
ax0.set_ylabel('Execution time (µs)')

ax1.stem(lengths, smoothness)
ax1.set_ylabel('Smoothness of input length\n(lower is better)')
ax1.set_xlabel('Length of input');

For historical reasons, most implementations return an array where frequencies vary from low to high to low (see “Discrete Fourier Transforms” for further explanation of frequencies). For example, when we do the real Fourier transform of a signal of all ones, an input that has no variation and therefore only has the slowest, constant Fourier component (also known as the “DC,” or Direct Current, component—just electronics jargon for “mean of the signal”), appearing as the first entry:

from scipy
import fftpack
N = 10

fftpack.fft(np.ones(N)) # The first component is np.mean(x) * N

from scipy import fftpack
N = 10

fftpack.fft(np.ones(N))  # The first component is np.mean(x) * N

array([10. + 0. j, 0. + 0. j, 0. + 0. j, 0. + 0. j, 0. + 0. j, 0. + 0. j,
   0. - 0. j, 0. - 0. j, 0. - 0. j, 0. - 0. j
])

When we try the FFT on a rapidly changing signal, we see a high-frequency component appear:

z = np.ones(10)
z[::2] = -1

print(f 'Applying FFT to {z}')
fftpack.fft(z)

Applying FFT to [-1.  1. -1.  1. -1.  1. -1.  1. -1.  1.]

array([0. + 0. j, 0. + 0. j, 0. + 0. j, 0. + 0. j, 0. + 0. j, -10. + 0. j,
   0. - 0. j, 0. - 0. j, 0. - 0. j, 0. - 0. j
])

Note that the FFT returns a complex spectrum that, in the case of real inputs, is conjugate symmetrical (i.e., symmetric in the real part and antisymmetric in the imaginary part):

x = np.array([1, 5, 12, 7, 3, 0, 4, 3, 2, 8])
X = fftpack.fft(x)

np.set_printoptions(precision = 2)

print("Real part:     ", X.real)
print("Imaginary part:", X.imag)

np.set_printoptions()

If we examine the Fourier transform of a rectangular pulse, we see significant side lobes in the spectrum:

x = np.zeros(500)
x[100: 150] = 1

X = fftpack.fft(x)

f, (ax0, ax1) = plt.subplots(2, 1, sharex = True)

ax0.plot(x)
ax0.set_ylim(-0.1, 1.1)

ax1.plot(fftpack.fftshift(np.abs(X)))
ax1.set_ylim(-5, 55);

We measure the signal for only a short time, labeled T_eff. The Fourier transform assumes that x(8) = x(0), and that the signal is continued as the dashed, rather than the solid, line. This introduces a big jump at the edge, with the expected oscillation in the spectrum:

t = np.linspace(0, 1, 500)
x = np.sin(49 * np.pi * t)

X = fftpack.fft(x)

f, (ax0, ax1) = plt.subplots(2, 1)

ax0.plot(x)
ax0.set_ylim(-1.1, 1.1)

ax1.plot(fftpack.fftfreq(len(t)), np.abs(X))
ax1.set_ylim(0, 190);

We can counter this effect by a process called windowing. The original function is multiplied with a window function such as the Kaiser window K(N, β). Here we visualize it for β ranging from 0 to 100:

f, ax = plt.subplots()

N = 10
beta_max = 5
colormap = plt.cm.plasma

norm = plt.Normalize(vmin = 0, vmax = beta_max)

lines = [
   ax.plot(np.kaiser(100, beta), color = colormap(norm(beta)))
   for beta in np.linspace(0, beta_max, N)
]

sm = plt.cm.ScalarMappable(cmap = colormap, norm = norm)

sm._A = []

plt.colorbar(sm).set_label(r 'Kaiser $\beta$');

First, we take the FFTs of our three signals (synthetic single target, synthetic multi-target, and real) and then display the positive frequency components (i.e., components 0 to N/2; see Figure 4-12). These are called the range traces in radar terminology.

fig, axes = plt.subplots(3, 1, sharex = True, figsize = (4.8, 2.4))

# Take FFTs of our signals.Note the convention to name FFTs with a
# capital letter.

V_single = np.fft.fft(v_single)
V_sim = np.fft.fft(v_sim)
V_actual = np.fft.fft(v_actual)

N = len(V_single)

with plt.style.context('style/thinner.mplstyle'):
   axes[0].plot(np.abs(V_single[: N // 2]))
            axes[0].set_ylabel("$|V_\mathrm{single}|$") axes[0].set_xlim(0, N // 2)
               axes[0].set_ylim(0, 1100)

               axes[1].plot(np.abs(V_sim[: N // 2]))
                        axes[1].set_ylabel("$|V_\mathrm{sim} |$") axes[1].set_ylim(0, 1000)

                        axes[2].plot(np.abs(V_actual[: N // 2]))
                                 axes[2].set_ylim(0, 750) axes[2].set_ylabel("$|V_\mathrm{actual}|$")

                                 axes[2].set_xlabel("FFT component $n$")

                                 for ax in axes:
                                 ax.grid()

This result is quite satisfying—but the dynamic range is so large that we could very easily miss some peaks. Let’s take the log as before with the spectrogram:

c = 3e8 # Approximately the speed of light and of
   # electromagnetic waves in air

fig, (ax0, ax1, ax2) = plt.subplots(3, 1)

def dB(y):
   "Calculate the log ratio of y / max(y) in decibel."

y = np.abs(y)
y /= y.max()

return 20 * np.log10(y)

def log_plot_normalized(x, y, ylabel, ax):
   ax.plot(x, dB(y))
ax.set_ylabel(ylabel)
ax.grid()

rng = np.arange(N // 2) * c / 2 / Beff

      with plt.style.context('style/thinner.mplstyle'):
      log_plot_normalized(rng, V_single[: N // 2], "$|V_0|$ [dB]", ax0)
            log_plot_normalized(rng, V_sim[: N // 2], "$|V_5|$ [dB]", ax1)
                  log_plot_normalized(rng, V_actual[: N // 2], "$|V_{\mathrm{actual}}|$ [dB]"
                        , ax2)

                     ax0.set_xlim(0, 300) # Change x limits
                     for these plots so that ax1.set_xlim(0, 300) # we are better able to see the shape of the peaks.ax2.set_xlim(0, len(V_actual) // 2)
                        ax2.set_xlabel('range')