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.

## 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:

Subscribe to:
Post Comments (Atom)

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.

ReplyDeleteTry with ff = 5.5 to see what I mean;

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

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

ReplyDeleteHello 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.

ReplyDeleteThe 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

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?

ReplyDeleteI don't know :)

ReplyDeleteHave you tried to take a look at the struct function?

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

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

ReplyDeleteOK, I solved the problem using struct pack:

ReplyDeletep = 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!

Great job Mark! you're always welcome.

ReplyDeleteabsolutly 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

ReplyDeleteHi, you can select a piece of an array in this way:

ReplyDeletemypiece = myarray[100:350]

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

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.

ReplyDeleteThank 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

ReplyDeleteplot(frq,abs(Y),'r') # plotting the spectrum by

stem(frq,abs(Y),'r') # plotting the spectrum

orcaja

and also changed

Deletet = arange(0,1,Ts) # time vector.... by

t = arange(0,1+Ts,Ts)

jut to have full cycles.

orcaja

another thing, the normalization should be twice that value:

DeleteY = 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.

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

ReplyDeleteThis 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.

ReplyDelete- 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.

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

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

ReplyDeleteHi Black, you need the inverse Fourier Transform: numpy.fft.ifft

DeleteIn 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---------------------------------------------------------------------------

ReplyDeleteTypeError 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

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

Deletefrq = 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

Hi ,

ReplyDeleteI'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

Thanks for the example!

ReplyDeleteOne 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

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

DeleteOr just...

ReplyDelete# 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

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

ReplyDeleteAmplitude is off by a factor of 2

ReplyDeleteVery good information, thanks for sharing.

ReplyDeleteWebsite Design Bangalore|Web Designing Company BangalorePython: How to plot a period of one square signal and press ENTER and second period shows on the same graph? Help please. :)

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

DeleteThanks :)

DeleteHow would this work if we had to feed bulk data in csv or txt format?

ReplyDeleteYou first need to parse the data. I recommend you to look into pandas.

DeleteThis comment has been removed by the author.

ReplyDeleteHow to plot the psd using the given example?

ReplyDeleteThis post shares excellent resources. Great!

ReplyDeleteThis 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).

ReplyDeleteAs 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?