Wednesday, August 3, 2011

How to plot the frequency spectrum with scipy

Spectrum analysis is the process of determining the frequency domain representation of a time domain signal and most commonly employs the Fourier transform. The Discrete Fourier Transform (DFT) is used to determine the frequency content of signals and the Fast Fourier Transform (FFT) is an efficient method for calculating the DFT. Scipy implements FFT and in this post we will see a simple example of spectrum analysis:
from numpy import sin, linspace, pi
from pylab import plot, show, title, xlabel, ylabel, subplot
from scipy import fft, arange

def plotSpectrum(y,Fs):
 """
 Plots a Single-Sided Amplitude Spectrum of y(t)
 """
 n = len(y) # length of the signal
 k = arange(n)
 T = n/Fs
 frq = k/T # two sides frequency range
 frq = frq[range(n/2)] # one side frequency range

 Y = fft(y)/n # fft computing and normalization
 Y = Y[range(n/2)]
 
 plot(frq,abs(Y),'r') # plotting the spectrum
 xlabel('Freq (Hz)')
 ylabel('|Y(freq)|')

Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = arange(0,1,Ts) # time vector

ff = 5;   # frequency of the signal
y = sin(2*pi*ff*t)

subplot(2,1,1)
plot(t,y)
xlabel('Time')
ylabel('Amplitude')
subplot(2,1,2)
plotSpectrum(y,Fs)
show()
The program shows the following figure, on top we have a plot of the signal and on the bottom the frequency spectrum.

39 comments:

  1. you might want to add some zero filling to sample the frequencey axis more densely. the spectrum looks artificially good now because the peak maximum is exactly on a datapoint of the frequency axis.

    Try with ff = 5.5 to see what I mean;

    ReplyDelete
  2. Hi, thank for your comment. I know what you mean, I chose ff = 5 on purpose.

    ReplyDelete
  3. How should the function be called? Does Fs represent the frequency spectrum? Should Fs be a list? Isn't Fs the thing to be calculated?

    ReplyDelete
  4. Hello Mark, Fs is the number of samples per second, hence it's an integer. For example, if you have an audio signal sampled with 44100 samples per second you have to set Fs = 44100.

    The aim of this snippet is to compute the frequency spectrum, not the sampling rate. Usually the sampling rate is known.

    You can find out more about signal processing in python on this post:
    http://glowingpython.blogspot.com/2011/09/sound-synthesis.html

    ReplyDelete
  5. I'm experimenting to see how fast Python and SciPy can calculate sound. In the sound synthesis post, you output to a wave file of 16 bit signed integers. However, I'm using PyAudio.write to output directly to the Windows audio and it expects data frames of 2 byte strings in little-endian format. Is there a simple way to convert the data?

    ReplyDelete
  6. I don't know :)
    Have you tried to take a look at the struct function?

    http://docs.python.org/library/struct.html

    ReplyDelete
  7. I'll look at the struct function in the next few days. I don't use Python enough to know all the library functions.

    ReplyDelete
  8. OK, I solved the problem using struct pack:
    p = pyaudio.PyAudio()
    stream = p.open(format = 8, channels = 1, rate = 44100, output = True)
    fmt = struct.Struct('h')

    for i in range(len(sound)):
    stream.write(fmt.pack(sound[i]))

    stream.close()
    p.terminate()

    Thanks for the help!

    ReplyDelete
  9. Great job Mark! you're always welcome.

    ReplyDelete
  10. absolutly helpfull, thx... But does anybody an idea, how to "window" the Signal? What I mean is, I've a Signal with 4000000 Datasamples. How can I create a spectrum without a so much datapoints? (I ve a knot in my brain, sry^^)... I need some kind of a meanvalue without changing the Information too much :) thx

    ReplyDelete
  11. Hi, you can select a piece of an array in this way:

    mypiece = myarray[100:350]

    where 100 and 350 the indexes that extract the elements 100 through 344.

    ReplyDelete
  12. thanks :) ... by the way, I´ve take a look to thw pylab.specgram function. It does what I need... not really clean and beauty, because of all the created figures, but anyway. Never the less, I am thankfull for your help to get a better understanding. I also read some stuff about "windowing" with a "hann-window". This is a possibility with overlap add, to fouriertransform the signaal without getting to much trouble couse of the "windotransformation". At least a rectangel becames a sinc-window through FFT.

    ReplyDelete
  13. Thank you, I like it! I used it to starting my own program, however, since the spectrum is not continuous I imported stem from pylab, and changed
    plot(frq,abs(Y),'r') # plotting the spectrum by
    stem(frq,abs(Y),'r') # plotting the spectrum

    orcaja

    ReplyDelete
    Replies
    1. and also changed
      t = arange(0,1,Ts) # time vector.... by
      t = arange(0,1+Ts,Ts)
      jut to have full cycles.

      orcaja

      Delete
    2. another thing, the normalization should be twice that value:
      Y = 2*fft(y)/n # fft computing and normalization
      In this way the spectrum will have the amplitud of 1, the same as in the input.

      Delete
  14. i tried this on my xubuntu and the bottom graph (the red one) don't come up at all, it's just plane white. not really getting why as i'm a newbie you might say at this high level math (i mean FFT). can someone help, i'm really interested in finding peak frequency at any (real) time of sound file

    ReplyDelete
  15. This may be off the subject. I am trying to make a simple hearing aid. An audiologist maps the decibels at each frequency a person hears at. The decibel threshold for normal is less than 25 decibels.
    - The simplest hearing aid would just amplify the sound to the equivalent 25 decibels.
    - A slightly more complex system would amplify individual frequencies according to the hearing pattern of the individual. - Expensive ($3,000+) hearing aids can filter out background noise.
    Hearing aids are very easily $2,000 per ear for the low end. Some go down to $350, but it is uncertain what they do. The elderly need hearing aids the most, but they are also the ones who can least afford it. They are rarely covered by insurance.

    I would like to break the sound into frequencies and amplitudes for the frequencies. The system would then amplify each frequency in a sound sample (10ms?) separately and combine the resulting frequencies outputting a 10ms processed interval. For speed, this would be run parallel on a multicore system or on a GPU, using something like PyCuda. Ideally, this would be done in real time. However, it could give a delayed result if speed is not possible.

    Am I on the right track? This seems like a simple problem using fft generating a frequency spectrum. Thanks.

    ReplyDelete
    Replies
    1. I believe you want to filter your signal then amplify it. You can do it in real time without using any GPU.

      Delete
  16. How would you reverse process. Start with frequency and amplitude information and then display the resulting pattern?

    ReplyDelete
    Replies
    1. Hi Black, you need the inverse Fourier Transform: numpy.fft.ifft

      Delete
  17. In the line """Y = fft(y)/n # fft computing and normalization""", I don't understand why you divide by `n`. I would interpret that to mean that longer the duration of the signal, the smaller the amplitude values in `Y`. I would think you would want to divide by the max possible value for the amplitude data to nomalize it so it falls between 0 and 1.

    ReplyDelete
  18. ---------------------------------------------------------------------------
    TypeError Traceback (most recent call last)
    in ()
    32 ylabel('Amplitude')
    33 subplot(2,1,2)
    ---> 34 plotSpectrum(y,Fs)
    35 show()

    in plotSpectrum(y, Fs)
    11 T = n/Fs
    12 frq = k/T # two sides frequency range
    ---> 13 frq = frq[range(n/2)] # one side frequency range
    14
    15 Y = fft(y)/n # fft computing and normalization

    TypeError: 'float' object cannot be interpreted as an integer

    ReplyDelete
    Replies
    1. I believe they are looking to slice the frequency spectrum there. The correct way to do this now is:

      frq = frq[1:n/2] # one side frequency range (sliced the frq in half)

      you will also have this problem with the Y values. It should be

      Y = Y[1:n/2]

      please correct me if i am wrong

      Delete
  19. Hi ,

    I'd like to get quite the same analysis here to reproduce the spectrum while playing a wav file. Something like here : https://www.youtube.com/watch?v=Tmpl5KA02S4.
    I just need to get the 8x8 matrix for each time (let's say each 1s) and I can manage the leds after that.

    If it doesn't bother you could you point out some code please ?
    it about many days I couldn't end with something successful.

    Thanks

    ReplyDelete
  20. Thanks for the example!
    One question though,

    how can n=len(y) actually be calculated since y is a function? Also if I try to print n it says it's undefinied. But how is it possible to (finally) use it as the x axis?

    Thanks in advance!
    Jan

    ReplyDelete
    Replies
    1. hi, y must be an array. If your y is a function you need to sample it. n is the length of this array.

      Delete
  21. Or just...

    # Express FFT in the frequency domain.
    def spectrum(signal, Time):

    frq = fftfreq(signal.size, d = Time[1] - Time[0] )
    Y = rfft(signal)

    return frq, Y

    ReplyDelete
  22. Why does the fft distribution change with the length of sound sample in SciPy?

    ReplyDelete
  23. Amplitude is off by a factor of 2

    ReplyDelete
  24. Python: How to plot a period of one square signal and press ENTER and second period shows on the same graph? Help please. :)

    ReplyDelete
    Replies
    1. Hi Ivana, you have to use the interactive mode, have a look here: http://matplotlib.org/users/shell.html

      Delete
  25. How would this work if we had to feed bulk data in csv or txt format?

    ReplyDelete
    Replies
    1. You first need to parse the data. I recommend you to look into pandas.

      Delete
  26. This comment has been removed by the author.

    ReplyDelete
  27. How to plot the psd using the given example?

    ReplyDelete
  28. This post shares excellent resources. Great!

    ReplyDelete
  29. This is not a one sided spectrum, but rather just half of a 2 sided spectrum. You should multiply the remaining one-sided values by 2 to normalise to a one-sided spectrum (to preserve the energy of the spectrum).

    As f. sheng noted, its clear to see the resulting amplitudes in this plot are off by a factor of 2.

    Or maybe im missing something here?

    ReplyDelete

Note: Only a member of this blog may post a comment.