Smoothing Window Data Averaging in Python – Moving Triangle Tecnique

Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

While I wrote a pervious post on linear data smoothing with python, those scripts were never fully polished. Fred (KJ4LFJ) asked me about this today and I felt bad I had nothing to send him. While I might add that the script below isn’t polished, at least it’s clean. I’ve been using this method for all of my smoothing recently. Funny enough, none of my code was clean enough to copy and paste, so I wrote this from scratch tonight. It’s a function to take a list in (any size) and smooth it with a triangle window (of any size, given by “degree”) and return the smoothed data with or without flanking copies of data to make it the identical length as before. The script also graphs the original data vs. smoothed traces of varying degrees. The output is below. I hope it helps whoever wants it!

import numpy
import pylab

def smoothTriangle(data,degree,dropVals=False):
"""performs moving triangle smoothing with a variable degree."""
"""note that if dropVals is False, output length will be identical
to input length, but with copies of data at the flanking regions"""
triangle=numpy.array(range(degree)+[degree]+range(degree)[::-1])+1
smoothed=[]
for i in range(degree,len(data)-degree*2):
point=data[i:i+len(triangle)]*triangle
smoothed.append(sum(point)/sum(triangle))
if dropVals: return smoothed
smoothed=[smoothed]*(degree+degree/2)+smoothed
while len(smoothed)<len(data):smoothed.append(smoothed[-1])
return smoothed

### CREATE SOME DATA ###
data=numpy.random.random(100) #make 100 random numbers from 0-1
data=numpy.array(data*100,dtype=int) #make them integers from 1 to 100
for i in range(100):
data[i]=data[i]+i**((150-i)/80.0) #give it a funny trend

### GRAPH ORIGINAL/SMOOTHED DATA ###
pylab.plot(data,"k.-",label="original data",alpha=.3)
pylab.plot(smoothTriangle(data,3),"-",label="smoothed d=3")
pylab.plot(smoothTriangle(data,5),"-",label="smoothed d=5")
pylab.plot(smoothTriangle(data,10),"-",label="smoothed d=10")
pylab.title("Moving Triangle Smoothing")
pylab.grid(alpha=.3)
pylab.axis([20,80,50,300])
pylab.legend()
pylab.show()

Simple Python Spectrograph with PyGame

Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

While thinking of ways to improve my QRSS VD high-definitions spectrograph software, I often wish I had a better way to display large spectrographs. Currently I’m using PIL (the Python Imaging Library) with TK and it’s slow as heck. I looked into the PyGame project, and it seems to be designed with speed in mind. I whipped-up this quick demo, and it’s a simple case audio spectrograph which takes in audio from your sound card and graphs it time vs. frequency. This method is far superior to the method I was using previously to display the data, because while QRSS VD can only update the entire GUI (500px by 8,000 px) every 3 seconds, early tests with PyGame suggests it can do it about 20 times a second (wow!). With less time/CPU going into the GUI, the program can be more responsivle and my software can be less of a drain.

import pygame
import numpy
import pyaudio
import scipy
import scipy.fftpack
import scipy.io.wavfile
import wave
rate=12000 #try 5000 for HD data, 48000 for realtime
soundcard=2
windowWidth=500
fftsize=512
currentCol=0
scooter=[]
overlap=5 #1 for raw, realtime - 8 or 16 for high-definition
def graphFFT(pcm):
global currentCol, data
ffty=scipy.fftpack.fft(pcm) #convert WAV to FFT
ffty=abs(ffty[0:len(ffty)/2])/500 #FFT is mirror-imaged
#ffty=(scipy.log(ffty))*30-50 # if you want uniform data
print "MIN:t%stMAX:t%s"%(min(ffty),max(ffty))
for i in range(len(ffty)):
if ffty[i]<0: ffty[i]=0
if ffty[i]>255: ffty[i]=255
scooter.append(ffty)
if len(scooter)<6:return
scooter.pop(0)
ffty=(scooter+scooter*2+scooter*3+scooter*2+scooter)/9
data=numpy.roll(data,-1,0)
data[-1]=ffty[::-1]
currentCol+=1
if currentCol==windowWidth: currentCol=0

def record():
p = pyaudio.PyAudio()
inStream = p.open(format=pyaudio.paInt16,channels=1,rate=rate,
input_device_index=soundcard,input=True)
linear=*fftsize
while True:
linear=linear[fftsize/overlap:]
linear=numpy.append(linear,pcm)
graphFFT(linear)

pal = [(max((x-128)*2,0),x,min(x*2,255)) for x in xrange(256)]
print max(pal),min(pal)
data=numpy.array(numpy.zeros((windowWidth,fftsize/2)),dtype=int)
#data=Numeric.array(data) # for older PyGame that requires Numeric
pygame.init() #crank up PyGame
pygame.display.set_caption("Simple Spectrograph")
screen=pygame.display.set_mode((windowWidth,fftsize/2))
world=pygame.Surface((windowWidth,fftsize/2),depth=8) # MAIN SURFACE
world.set_palette(pal)
t_rec.daemon=True # daemon mode forces thread to quit with program
clk=pygame.time.Clock()
while 1:
for event in pygame.event.get(): #check if we need to exit
if event.type == pygame.QUIT:pygame.quit();sys.exit()
pygame.surfarray.blit_array(world,data) #place data in window
screen.blit(world, (0,0))
pygame.display.flip() #RENDER WINDOW
clk.tick(30) #limit to 30FPS

Simple-Case PyGame Example

Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

I’m starting to investigate PyGame as an alternative to PIL and K for my QRSS VD spectrograph project. This sample code makes a box bounce around a window.

import pygame, sys
pygame.init() #load pygame modules
size = width, height = 320, 240 #size of window
speed = [2, 2] #speed and direction
screen = pygame.display.set_mode(size) #make window
s=pygame.Surface((100,50)) #create surface 100px by 50px
s.fill((33,66,99)) #color the surface blue
r=s.get_rect() #get the rectangle bounds for the surface
clock=pygame.time.Clock() #make a clock
while 1: #infinite loop
clock.tick(30) #limit framerate to 30 FPS
for event in pygame.event.get(): #if something clicked
if event.type == pygame.QUIT: #if EXIT clicked
sys.exit() #close cleanly
r=r.move(speed) #move the box by the "speed" coordinates
#if we hit a  wall, change direction
if r.left < 0 or r.right > width: speed = -speed
if r.top < 0 or r.bottom > height: speed = -speed
screen.fill((0,0,0)) #make redraw background black
screen.blit(s,r) #render the surface into the rectangle
pygame.display.flip() #update the screen

QRSS VD Image Assembler

Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

This minimal Python script will convert a directory filled with tiny image captures such as this into gorgeous montages as seen below! I whipped-up this script tonight because I wanted to assess the regularity of my transmitter’s embarrassing drift. I hope you find it useful.

import os
from PIL import Image

x1,y1,x2,y2=[0,0,800,534] #crop from (x,y) 0,0 to 800x534
squish=10 #how much to squish it horizontally

### LOAD LIST OF FILES ###
workwith=[]
for fname in os.listdir('./'):
if ".jpg" in fname and not "assembled" in fname:
workwith.append(fname)
workwith.sort()

### MAKE NEW IMAGE ###
im=Image.new("RGB",(x2*len(workwith),y2))
for i in range(len(workwith)):
im2=Image.open(workwith[i])
im2=im2.crop((x1,y1,x2,y2))
im.paste(im2,(i*x2,0))
print "saving BIG image"
im.save("assembled.jpg")
print "saving SQUISHED image"
im=im.resize((im.size/10,im.size),Image.ANTIALIAS)
im.save("assembled-squished.jpg")
print "DONE"

Script to download every image linked to from a webpage:

import urllib2
import os

suckFrom="http://w1bw.org/grabber/archive/2010-06-08/"

f=urllib2.urlopen(suckFrom)
f.close()

for line in s:
if ".jpg" in line and not line in download and not "thumb" in line:

fname = url.split("/")[-1].replace(":","-")
if fname in os.listdir('./'):
else:
output=open(fname,'wb')
output.close()

Generate LUTs with Python

Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

I wrote a script to generate and display LUTs with Python. There has been a lot of heated discussion in the QRSS Knights mailing list as to the use of color maps when representing QRSS data. I’ll make a separate post (perhaps later?) documenting why it’s so critical to use particular mathematically-generated color maps rather than empirical “looks good to me” color selections. Anyway, this is what I came up with:

For my QRSS needs, I desire a colormap which is aesthetically pleasing but can also be quickly reverted to its original (gray-scale) data. I accomplished this by choosing a channel (green in this case) and applying its intensity linearly with respect to the value it represents. Thus, any “final” image can be imported into an editor, split by RGB, and the green channel represents the original data. This allows adjustment of contrast/brightness and even the reassignment of a different colormap, all without losing any data!

ORIGINAL DATA:
(that’s the “flying W” and the FSK signal below it is WA5DJJ)

Note that it looks nice, shows weak signals, doesn’t get blown-out by strong signals, and it fully includes the noise floor (utilizing all available data). BLUE CHANNEL: weak signals / noise floor RED CHANNEL: strong signals / no noise They look okay, but not as great as linear

The following links are downloadable LUTs which can be applied to 8-bit grayscale images using most editors (i.e., MBF ImageJ) generated by the python script below.
Linear LUT
Nonlinear LUT

This is the Python script I wrote to generate the downloadable LUTs, graphs, and scale bars / keys / legends which are not posted. It requires python, matplotlib, and PIL.

import math
import pylab
from PIL import Image

####################### GENERATE RGB VALUES #######################

r,g,b=[],[],[]
name="Blin_Glin_Rlin"
for i in range(256):
if i>128: #LOW HALF
j=128
k=i
else: #HIGH HALF
k=128
j=i
#b.append((math.sin(3.1415926535*j/128.0/2))*256)
#r.append((1+math.sin(3.1415926535*(k-128*2)/128.0/2))*256)
r.append(k*2-255)
g.append(i)
b.append(j*2-1)

if r[-1]<0:r[-1]=0
if g[-1]<0:g[-1]=0
if b[-1]<0:b[-1]=0

if r[-1]>255:b[-1]=255
if g[-1]>255:g[-1]=255
if b[-1]>255:b[-1]=255

####################### SAVE LUT FILE #######################
im = Image.new("RGB",(256*2,10*4))
for x in range(256):
for y in range(10):
pix[x,y] = (r[x],g[x],b[x])
pix[x,y+10] = (r[x],0,0)
pix[x,y+20] = (0,g[x],0)
pix[x,y+30] = (0,0,b[x])
a=(g[x]+g[x]+g[x])/3
pix[256+x,y] = (a,a,a)
pix[256+x,y+10] = (r[x],r[x],r[x])
pix[256+x,y+20] = (g[x],g[x],g[x])
pix[256+x,y+30] = (b[x],b[x],b[x])
#im=im.resize((256/2,40),Image.ANTIALIAS)
im.save(name+"_scale.png")

####################### PLOT IT #######################
pylab.figure(figsize=(8,4))
pylab.grid(alpha=.3)
pylab.title(name)
pylab.xlabel("Data Value")
pylab.ylabel("Color Intensity")
pylab.plot(g,'g-')
pylab.plot(r,'r-')
pylab.plot(b,'b-')
pylab.axis([-10,266,-10,266])
pylab.subplots_adjust(top=0.90, bottom=0.14, left=0.1, right=0.97)
pylab.savefig(name+"_graph.png",dpi=60)
#pylab.show()

####################### SAVE LUT FILE #######################
f=open(name+".lut",'w')
out="IndextRedtGreentBluen"
for i in range(256):
out+=("t%dt%dt%dt%dn"%(i,r[i],g[i],b[i]))
f.write(out)
f.close()