DSPRelated.com
Forums

python psd

Started by encinoman June 8, 2009
Hi, I am trying to take a PSD of an FSK signal that I've been fooling
around with....I wrote it in python with the help of a post someone put on
this board a while back. Here is the code:

#!/usr/bin/env python
import math
import wave
import array
import random
import numpy
import pylab
from pylab import psd
from pylab import arange
from pylab import randn
from pylab import exp
from pylab import convolve
pi = pylab.pi
twopi = 2.0*pi
sin = pylab.sin
cos = pylab.cos
numBitsToTX = 1000
#data_list = [1,1,1,1,1]
#data_list = [0,1,0,1,0]
data_list=list()
def get_rand():	
	for i in xrange(numBitsToTX):
		bit =random.randrange(0,2)
		data_list.append(bit)
		print bit
	return 0
get_rand()
#init output and related vars
samplesPerSecond = 44100
Outfile = wave.open("outfile.wav", "w")
Outfile.setnchannels(1)
Outfile.setsampwidth(2) # 2 bytes for 32767 amplitude
Outfile.setframerate(samplesPerSecond)
Outfile.setcomptype("NONE", "Uncompressed")
totalArray = array.array('h')

def gen_wav(seconds, frequency, amplitudePercentOfMax):
    # calculate frequency as cycles/sample
    cyclesPerSample = float(frequency)/samplesPerSecond
    numberOfSamples = int(seconds*samplesPerSecond)
    maxAmplitude = (amplitudePercentOfMax/100.)*32767
    sampleArray = array.array('h')

    for nSamples in xrange(numberOfSamples):
        cycles = nSamples*cyclesPerSample
        cycles -= int(cycles) # makes 0.0<=cycles<1.0
        sampleValue = int(round(maxAmplitude*sin(twopi*cycles))) # round 
        sampleArray.append(sampleValue)
    return sampleArray

if __name__ == '__main__':
    
    #Outtemp = gen_wav(.01, 440, 75)
    # generate .001 sec = 44.1 samples & 2 cycles with 20 amplitude
    #sampleArray = gen_wav(.00075, 2000, (20./32767.)*100)
    for data in data_list:
		if data==1:
			sampleArray =gen_wav(1.0/2000, 2000, 75)
			totalArray=totalArray+sampleArray
	        if data==0:
			sampleArray=gen_wav(1.0/4000, 4000, 75)
			totalArray=totalArray+sampleArray

    # print first 34 samples as selfscaling ascii plot
   
    grid = list(' |%19s|%19s|' %('',''))
    scale = 1.0; maximumMagnitude = 20
    for i in range(min(1000,len(totalArray))):
        sample= totalArray[i]
        if abs(sample)>maximumMagnitude:
            maximumMagnitude=abs(sample)
            scale = 20./maximumMagnitude
        x= int(round(totalArray[i]*scale))+21
        print ''.join(grid[:x]+['*']+grid[x+1:])
    
    Outfile.writeframes(totalArray.tostring())
    Outfile.close()
print len(totalArray)
#print fft(totalArray)


#psd(totalArray, NFFT=512, Fs=44100)



I run it and get the following error:

UserWarning: psd, csd, and specgram have changed to scale their densities
by the sampling frequency for better MatLab compatibility. You can pass
scale_by_freq=False to disable this behavior.  Also, one-sided densities
are scaled by a factor of 2.
  warnings.warn("psd, csd, and specgram have changed to scale their"


I just need a second look at my code. I believe I have the right sampling
frequency and am following he new protocol for python2.6. Unfortunately, it
still wont work.

Thanks


Matt 




On Jun 8, 12:27&#4294967295;pm, "encinoman" <bajor.matt...@gmail.com> wrote:
> Hi, I am trying to take a PSD of an FSK signal that I've been fooling > around with....I wrote it in python with the help of a post someone put on > this board a while back. Here is the code: > > #!/usr/bin/env python > import math > import wave > import array > import random > import numpy > import pylab > from pylab import psd > from pylab import arange > from pylab import randn > from pylab import exp > from pylab import convolve > pi = pylab.pi > twopi = 2.0*pi > sin = pylab.sin > cos = pylab.cos > numBitsToTX = 1000 > #data_list = [1,1,1,1,1] > #data_list = [0,1,0,1,0] > data_list=list() > def get_rand(): > &#4294967295; &#4294967295; &#4294967295; &#4294967295; for i in xrange(numBitsToTX): > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; bit =random.randrange(0,2) > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; data_list.append(bit) > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; print bit > &#4294967295; &#4294967295; &#4294967295; &#4294967295; return 0 > get_rand() > #init output and related vars > samplesPerSecond = 44100 > Outfile = wave.open("outfile.wav", "w") > Outfile.setnchannels(1) > Outfile.setsampwidth(2) # 2 bytes for 32767 amplitude > Outfile.setframerate(samplesPerSecond) > Outfile.setcomptype("NONE", "Uncompressed") > totalArray = array.array('h') > > def gen_wav(seconds, frequency, amplitudePercentOfMax): > &#4294967295; &#4294967295; # calculate frequency as cycles/sample > &#4294967295; &#4294967295; cyclesPerSample = float(frequency)/samplesPerSecond > &#4294967295; &#4294967295; numberOfSamples = int(seconds*samplesPerSecond) > &#4294967295; &#4294967295; maxAmplitude = (amplitudePercentOfMax/100.)*32767 > &#4294967295; &#4294967295; sampleArray = array.array('h') > > &#4294967295; &#4294967295; for nSamples in xrange(numberOfSamples): > &#4294967295; &#4294967295; &#4294967295; &#4294967295; cycles = nSamples*cyclesPerSample > &#4294967295; &#4294967295; &#4294967295; &#4294967295; cycles -= int(cycles) # makes 0.0<=cycles<1.0 > &#4294967295; &#4294967295; &#4294967295; &#4294967295; sampleValue = int(round(maxAmplitude*sin(twopi*cycles))) # round > &#4294967295; &#4294967295; &#4294967295; &#4294967295; sampleArray.append(sampleValue) > &#4294967295; &#4294967295; return sampleArray > > if __name__ == '__main__': > > &#4294967295; &#4294967295; #Outtemp = gen_wav(.01, 440, 75) > &#4294967295; &#4294967295; # generate .001 sec = 44.1 samples & 2 cycles with 20 amplitude > &#4294967295; &#4294967295; #sampleArray = gen_wav(.00075, 2000, (20./32767.)*100) > &#4294967295; &#4294967295; for data in data_list: > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; if data==1: > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; sampleArray =gen_wav(1.0/2000, 2000, 75) > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; totalArray=totalArray+sampleArray > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; if data==0: > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; sampleArray=gen_wav(1.0/4000, 4000, 75) > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; totalArray=totalArray+sampleArray > > &#4294967295; &#4294967295; # print first 34 samples as selfscaling ascii plot > > &#4294967295; &#4294967295; grid = list(' |%19s|%19s|' %('','')) > &#4294967295; &#4294967295; scale = 1.0; maximumMagnitude = 20 > &#4294967295; &#4294967295; for i in range(min(1000,len(totalArray))): > &#4294967295; &#4294967295; &#4294967295; &#4294967295; sample= totalArray[i] > &#4294967295; &#4294967295; &#4294967295; &#4294967295; if abs(sample)>maximumMagnitude: > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; maximumMagnitude=abs(sample) > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; scale = 20./maximumMagnitude > &#4294967295; &#4294967295; &#4294967295; &#4294967295; x= int(round(totalArray[i]*scale))+21 > &#4294967295; &#4294967295; &#4294967295; &#4294967295; print ''.join(grid[:x]+['*']+grid[x+1:]) > > &#4294967295; &#4294967295; Outfile.writeframes(totalArray.tostring()) > &#4294967295; &#4294967295; Outfile.close() > print len(totalArray) > #print fft(totalArray) > > #psd(totalArray, NFFT=512, Fs=44100) > > I run it and get the following error: > > UserWarning: psd, csd, and specgram have changed to scale their densities > by the sampling frequency for better MatLab compatibility. You can pass > scale_by_freq=False to disable this behavior. &#4294967295;Also, one-sided densities > are scaled by a factor of 2. > &#4294967295; warnings.warn("psd, csd, and specgram have changed to scale their" > > I just need a second look at my code. I believe I have the right sampling > frequency and am following he new protocol for python2.6. Unfortunately, it > still wont work. > > Thanks > > Matt
Why the hell would anybody want to use such an obscure language? This is pure geeky stuff. Hardy
HardySpicer wrote:

> Why the hell would anybody want to use such an obscure language?
Python is just about mainstream. If you want an obscure langauge, try Ocaml, Haskell or Erlang. Erik -- ---------------------------------------------------------------------- Erik de Castro Lopo http://www.mega-nerd.com/
HardySpicer wrote:

[full quote and python code deleted]
> Why the hell would anybody want to use such an obscure language? This > is pure geeky stuff.
Anybody who now longer wants to take part in obfuscated C contents (meaning debugging your colleague's legacy C code). OK, the posted code lacks some empty lines to make it more readable. bye Andreas -- Andreas H&#4294967295;nnebeck | email: acmh@gmx.de ----- privat ---- | www : http://www.huennebeck-online.de Fax/Anrufbeantworter: 0721/151-284301 GPG-Key: http://www.huennebeck-online.de/public_keys/andreas.asc PGP-Key: http://www.huennebeck-online.de/public_keys/pgp_andreas.asc
> > UserWarning: psd, csd, and specgram have changed to scale their densities > by the sampling frequency for better MatLab compatibility. You can pass > scale_by_freq=False to disable this behavior. &#4294967295;Also, one-sided densities > are scaled by a factor of 2. > &#4294967295; warnings.warn("psd, csd, and specgram have changed to scale their" >
It doesn't look like an error. It looks like an informative warning. Not sure if you are referring to the printed "warning" or unintended behavior by the code. For future note, it is probably more efficient to post a simplified example that illustrates the problem than dumping a large code segment and asking others to debug it for you. Good luck.