Plot a fourier transform of a sin wav with matplotlib

Question

I am trying to plot a fourier transform of a sign wave based on the scipy documentation

import numpy as np
import matplotlib.pyplot as plt
import scipy.fft

def sinWav(amp, freq, time, phase=0):
    return amp * np.sin(2 * np.pi * (freq * time - phase))

def plotFFT(f, speriod, time):
    """Plots a fast fourier transform

    Args:
        f (np.arr): A signal wave
        speriod (int): Number of samples per second
        time ([type]): total seconds in wave
    """

    N = speriod * time
    # sample spacing
    T = 1.0 / 800.0
    x = np.linspace(0.0, N*T, N, endpoint=False)

    yf = scipy.fft.fft(f)
    xf = scipy.fft.fftfreq(N, T)[:N//2]
    plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
    plt.grid()
    plt.xlim([1,3])
    plt.show()


speriod = 1000
time  = {
    0: np.arange(0, 4, 1/speriod),
    1: np.arange(4, 8, 1/speriod),
    2: np.arange(8, 12, 1/speriod)
}

signal = np.concatenate([
    sinWav(amp=0.25, freq=2, time=time[0]),
    sinWav(amp=1, freq=2, time=time[1]),
    sinWav(amp=0.5, freq=2, time=time[2])
])   # generate signal

plotFFT(signal, speriod, 12)

Desired output

I want to be getting a fourier transform graph which looks like this

Current output

But instead it looks like this

Extra

This is the sin wave I am working with

Synthase · Accepted Answer · 2021-01-11T04:24:28.457

1

import numpy as np
import matplotlib.pyplot as plt
import scipy.fft

def sinWav(amp, freq, time, phase=0):
    return amp * np.sin(2 * np.pi * (freq * time - phase))

def plotFFT(f, speriod, time):
    """Plots a fast fourier transform

    Args:
        f (np.arr): A signal wave
        speriod (int): Number of samples per second
        time ([type]): total seconds in wave
    """

    N = speriod * time
    # sample spacing
    T = 1.0 / 800.0
    x = np.linspace(0.0, N*T, N, endpoint=False)

    yf = scipy.fft.fft(f)
    xf = scipy.fft.fftfreq(N, T)[:N//2]

    amplitudes = 1/speriod* np.abs(yf[:N//2])
  
    plt.plot(xf, amplitudes)
    plt.grid()
    plt.xlim([1,3])
    plt.show()


speriod = 800
time  = {
    0: np.arange(0, 4, 1/speriod),
    1: np.arange(4, 8, 1/speriod),
    2: np.arange(8, 12, 1/speriod)
}

signal = np.concatenate([
    sinWav(amp=0.25, freq=2, time=time[0]),
    sinWav(amp=1, freq=2, time=time[1]),
    sinWav(amp=0.5, freq=2, time=time[2])
])   # generate signal

plotFFT(signal, speriod, 12)

You should have what you want. Your amplitudes were not properly computed, as your resolution and speriod were inconsistent.

Longer data acquisition:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fft

def sinWav(amp, freq, time, phase=0):
    return amp * np.sin(2 * np.pi * (freq * time - phase))

def plotFFT(f, speriod, time):
    """Plots a fast fourier transform

    Args:
        f (np.arr): A signal wave
        speriod (int): Number of samples per second
        time ([type]): total seconds in wave
    """

    N = speriod * time
    # sample spacing
    T = 1.0 / 800.0
    x = np.linspace(0.0, N*T, N, endpoint=False)

    yf = scipy.fft.fft(f)
    xf = scipy.fft.fftfreq(N, T)[:N//2]

    amplitudes = 1/(speriod*4)* np.abs(yf[:N//2])
  
    plt.plot(xf, amplitudes)
    plt.grid()
    plt.xlim([1,3])
    plt.show()


speriod = 800
time  = {
    0: np.arange(0, 4*4, 1/speriod),
    1: np.arange(4*4, 8*4, 1/speriod),
    2: np.arange(8*4, 12*4, 1/speriod)
}

signal = np.concatenate([
    sinWav(amp=0.25, freq=2, time=time[0]),
    sinWav(amp=1, freq=2, time=time[1]),
    sinWav(amp=0.5, freq=2, time=time[2])
])   # generate signal

plotFFT(signal, speriod, 48)

edited Jan 11 '21 at 04:24

answered Jan 11 '21 at 03:41

Synthase

5,849
2
12
34

It should look like that graph I posted where the amplitude is 3.5, at freq 2 – Sam Jan 11 '21 at 03:43
Why doesn't this work with other sampling rates? Shouldn't the freqencies be independent of the sampling rate? I thought the sr just determined the size along the x axis that the signal was split into – Sam Jan 11 '21 at 04:07
1

Well, if you're working on continuous data you're right (so in maths). But in physics, data are not continuous. And FFT refers to Discrete FT, in that specific physical context. The reason for my changes can be found here: https://en.wikipedia.org/wiki/Nyquist_frequency It is all about this relationship. In fact, your shift is due to spectrum "folding" cause you're not enough in line with the Nyquist relation. – Synthase Jan 11 '21 at 04:09
I'm trying to make my graph look more continuous, but I keep getting the same blocky graph – Sam Jan 11 '21 at 04:11
In other words, I want the x axis split into smaller sections – Sam Jan 11 '21 at 04:15
Well, you should definitely work on your lectures :). To do that, you need more timepoints. See my edit for changes. – Synthase Jan 11 '21 at 04:16
Wouldn't I just want to increase the `speriod`? So that the time is split into smaller chunks? I've tried doing that, but receive the same results. I also noticed that the updates you just made changed the amplitude values – Sam Jan 11 '21 at 04:20
Yes sorry I corrected. In fact, amplitudes are calculated as 1/n_samples, with n_samples = speriod * nb_period. So you increase the number of samples, but not the speriod/sampling frequency. As a general rule, the longer your wave in time domain is, the better is the resolution. – Synthase Jan 11 '21 at 04:25
So this doesn't work correctly when I change the freqency of the sin waves(if say I change them to frequency 1 instead of 2) – Sam Jan 11 '21 at 04:47
I'm learning this on my own. I went through a udemy course and a textbook and it still all seems foreign – Sam Jan 11 '21 at 04:54
Well, by the way: if you set the freq to 1, the peak is at 1. which is what you expect (1Hz). So it is working :). – Synthase Jan 11 '21 at 04:59

score 0 · Answer 2 · answered Jan 15 '21 at 00:30

You can also interactively plot this. You may need to install the pip install scikit-dsp-comm

# !pip install scikit-dsp-comm
# Make an interactive version of the above
from ipywidgets import interact, interactive
import numpy as np
import matplotlib.pyplot as plt
from scipy import fftpack

plt.rcParams['figure.figsize'] = [10, 8]

font = {'weight' : 'bold',
        'size'   : 14}
plt.rc('font', **font)


def pulse_plot(fm = 1000, Fs = 2010):
    tlen = 1.0  # length in seconds
    # generate time axis
    tt = np.arange(np.round(tlen*Fs))/float(Fs)
    # generate sine
    xt = np.sin(2*np.pi*fm*tt)
    plt.subplot(211)
    plt.plot(tt[:500], xt[:500], '-b')
    plt.plot(tt[:500], xt[:500], 'or', label='xt values')
    plt.ylabel('$x(t)$')
    plt.xlabel('t [sec]')
    strt2 = 'Sinusoidal Waveform $x(t)$'
    strt2 = strt2 + ', $f_m={}$ Hz, $F_s={}$ Hz'.format(fm, Fs)
    plt.title(strt2)
    plt.legend()
    plt.grid()
    X = fftpack.fft(xt)
    freqs = fftpack.fftfreq(len(xt)) * Fs
    plt.subplot(212)
    N = xt.size
    # DFT
    X = np.fft.fft(xt)
    X_db = 20*np.log10(2*np.abs(X)/N)
    #f = np.fft.fftfreq(N, 1/Fs)
    f = np.arange(0, N)*Fs/N
    plt.plot(f, X_db, 'b')
    plt.xlabel('Frequency in Hertz [Hz]')
    plt.ylabel('Frequency Domain\n (Spectrum) Magnitude')
    plt.grid()
    plt.tight_layout()
    
    
interactive_plot = interactive(pulse_plot,fm = (1000,20000,1000), Fs = (1000,40000,10));
output = interactive_plot.children[-1]
# output.layout.height = '350px'
interactive_plot

@Sam because the y axis is just spectrum magnitude. You can make that linear if you want. — Jay Patel, Jan 15 '21 at 18:17

Plot a fourier transform of a sin wav with matplotlib

2 Answers2