The personal website of Scott W Harden
July 22nd, 2010

Idea: vdFSK modulation

My goal is to create a QRPP (extremely low power) transmitter and modulation method to send QRSS (extremely slow, frequency shifting data) efficiently, able to be decoded visually or with automated image analysis software. This evolving post will document the thought process and development behind AJ4VD's Frequency Shift Keying method, vdFSK.

Briefly, this is what my idea is. Rather than standard 2-frequencies (low for space, high for tone) QRSS3 (3 seconds per dot), I eliminate the need for pauses between dots by using 3 frequencies (low for a space between letters, medium for dot, high for dash). The following images compare my call sign (AJ4VD) being sent with the old method, and the vdFSK method.

Again, both of these images say the same thing: AJ4VD, .- .--- ....- ...- -.. However, note that the above image has greater than a 3 second dot, so it's unfairly long if you look at the time scale. Until I get a more fairly representative image, just appreciate it graphically. It's obviously faster to send 3 frequencies rather than two. In my case, it's over 200% faster.

This is the code to generate audio files converting a string of text into vdFSK audio, saving the output as a WAV file. Spectrographs can be created from these WAV files.

generate_audio.py

# converts a string into vdFSK audio saved as a WAV file

import numpy
import wave
from morse import *


def makeTone(freq, duration=1, samplerate=5000, shape=True):
    signal = numpy.arange(duration*samplerate) / \
        float(samplerate)*float(freq)*3.14*2
    signal = numpy.sin(signal)*16384
    if shape == True:  # soften edges
        for i in range(100):
            signal[i] = signal[i]*(i/100.0)
            signal[-i] = signal[-i]*(i/100.0)
    ssignal = ''
    for i in range(len(signal)):  # make it binary
        ssignal += wave.struct.pack('h', signal[i])
    return ssignal


def text2tone(msg, base=800, sep=5):
    audio = ''
    mult = 3  # secs per beep
    msg = " "+msg+" "
    for char in msg.lower():
        morse = lookup[char]
        print char, morse
        audio += makeTone(base, mult)
        for step in lookup[char]:
            if step[0] == ".":
                audio += makeTone(base+sep, int(step[1])*mult)
            if step[0] == "-":
                audio += makeTone(base+sep*2, int(step[1])*mult)
            if step[0] == "|":
                audio += makeTone(base, 3*mult)
    return audio


msg = "aj4vd"
file = wave.open('test.wav', 'wb')
file.setparams((1, 2, 5000, 5000*4, 'NONE', 'noncompressed'))
file.writeframes(text2tone(msg))
file.close()

print 'file written'

morse.py

# library for converting between text and Morse code
raw_lookup="""
a.- b-... c-.-. d-.. e. f..-. g--. h.... i.. j.--- k-- l.-.. m--
n-. o--- p.--. q--.- r.-. s... t- u.- v...- w.-- x-..- y-.-- z--..
0----- 1.---- 2..--- 3...-- 4....- 5..... 6-.... 7--... 8---.. 9----.
..-.-.- =-...- :---... ,--..-- /-..-. --....-
""".replace("n","").split(" ")

lookup={}
lookup[" "]=["|1"]
for char in raw_lookup:
    """This is a silly way to do it, but it works."""
    char,code=char[0],char[1:]
    code=code.replace("-----","x15 ")
    code=code.replace("----","x14 ")
    code=code.replace("---","x13 ")
    code=code.replace("--","x12 ")
    code=code.replace("-","x11 ")
    code=code.replace(".....","x05 ")
    code=code.replace("....","x04 ")
    code=code.replace("...","x03 ")
    code=code.replace("..","x02 ")
    code=code.replace(".","x01 ")
    code=code.replace("x0",'.')
    code=code.replace("x1",'-')
    code=code.split(" ")[:-1]
    #print char,code
    lookup[char]=code

Automated decoding is trivial. The image above was analyzed, turned into the image below, and the string (AJ4VD) was extracted:

decode.py

# given an image, it finds peaks and pulls data out
from PIL import Image
from PIL import ImageDraw
import pylab
import numpy

pixelSeek = 10
pixelShift = 15


def findPeak(data):
    maxVal = 0
    maxX = 0
    for x in range(len(data)):
        if data[x] > maxVal:
            maxVal, maxX = data[x], x
    return maxX


def peaks2morse(peaks):
    baseFreq = peaks[0]
    lastSignal = peaks[0]
    lastChange = 0
    directions = []
    for i in range(len(peaks)):
        if abs(peaks[i]-baseFreq) < pixelSeek:
            baseFreq = peaks[i]
        if abs(peaks[i]-lastSignal) < pixelSeek and i < len(peaks)-1:
            lastChange += 1
        else:
            if abs(baseFreq-lastSignal) < pixelSeek:
                c = " "
            if abs(baseFreq-lastSignal) < pixelSeek:
                c = " "
            if abs(baseFreq-lastSignal) < pixelSeek:
                c = " "
            directions.append(
                [lastSignal, lastChange, baseFreq, baseFreq-lastSignal])
            lastChange = 0
        lastSignal = peaks[i]
    return directions


def morse2image(directions):
    im = Image.new("L", (300, 100), 0)
    draw = ImageDraw.Draw(im)
    lastx = 0
    for d in directions:
        print d
        draw.line((lastx, d[0], lastx+d[1], d[0]), width=5, fill=255)
        lastx = lastx+d[1]
    im.show()


im = Image.open('raw.png')
pix = im.load()
data = numpy.zeros(im.size)
for x in range(im.size[0]):
    for y in range(im.size[1]):
        data[x][y] = pix[x, y]

peaks = []
for i in range(im.size[0]):
    peaks.append(findPeak(data[i]))

morse = peaks2morse(peaks)
morse2image(morse)
print morse
Markdown source code last modified on January 18th, 2021
---
title: Idea: vdFSK modulation
date: 2010-07-22 12:39:54
tags: python, qrss, old
---

# Idea: vdFSK modulation

<blockquote class="wp-block-quote"><p>My goal is to create a QRPP (extremely low power) transmitter and modulation method to send QRSS (extremely slow, frequency shifting data) efficiently, able to be decoded visually or with automated image analysis software. This evolving post will document the thought process and development behind AJ4VD's Frequency Shift Keying method, <b>vdFSK</b>.</p></blockquote>

__Briefly, this is what my idea is.__ Rather than standard 2-frequencies (low for space, high for tone) QRSS3 (3 seconds per dot), I eliminate the need for pauses between dots by using 3 frequencies (low for a space between letters, medium for dot, high for dash). The following images compare my call sign (AJ4VD) being sent with the old method, and the vdFSK method.

<div class="text-center img-border">

[![](traditional_thumb.jpg)](traditional.png)

</div>

__Again,__ both of these images say the same thing: AJ4VD, `.- .--- ....- ...- -..` However, note that the above image has greater than a 3 second dot, so it's unfairly long if you look at the time scale. Until I get a more fairly representative image, just appreciate it graphically. It's obviously faster to send 3 frequencies rather than two. In my case, it's over 200% faster.

<div class="text-center img-border">

[![](modulation_thumb.jpg)](modulation.png)

</div>

__This is the code to generate audio files__ converting a string of text into vdFSK audio, saving the output as a WAV file. Spectrographs can be created from these WAV files.

### generate_audio.py

```python
# converts a string into vdFSK audio saved as a WAV file

import numpy
import wave
from morse import *


def makeTone(freq, duration=1, samplerate=5000, shape=True):
    signal = numpy.arange(duration*samplerate) / \
        float(samplerate)*float(freq)*3.14*2
    signal = numpy.sin(signal)*16384
    if shape == True:  # soften edges
        for i in range(100):
            signal[i] = signal[i]*(i/100.0)
            signal[-i] = signal[-i]*(i/100.0)
    ssignal = ''
    for i in range(len(signal)):  # make it binary
        ssignal += wave.struct.pack('h', signal[i])
    return ssignal


def text2tone(msg, base=800, sep=5):
    audio = ''
    mult = 3  # secs per beep
    msg = " "+msg+" "
    for char in msg.lower():
        morse = lookup[char]
        print char, morse
        audio += makeTone(base, mult)
        for step in lookup[char]:
            if step[0] == ".":
                audio += makeTone(base+sep, int(step[1])*mult)
            if step[0] == "-":
                audio += makeTone(base+sep*2, int(step[1])*mult)
            if step[0] == "|":
                audio += makeTone(base, 3*mult)
    return audio


msg = "aj4vd"
file = wave.open('test.wav', 'wb')
file.setparams((1, 2, 5000, 5000*4, 'NONE', 'noncompressed'))
file.writeframes(text2tone(msg))
file.close()

print 'file written'
```

### morse.py

```python
# library for converting between text and Morse code
raw_lookup="""
a.- b-... c-.-. d-.. e. f..-. g--. h.... i.. j.--- k-- l.-.. m--
n-. o--- p.--. q--.- r.-. s... t- u.- v...- w.-- x-..- y-.-- z--..
0----- 1.---- 2..--- 3...-- 4....- 5..... 6-.... 7--... 8---.. 9----.
..-.-.- =-...- :---... ,--..-- /-..-. --....-
""".replace("n","").split(" ")

lookup={}
lookup[" "]=["|1"]
for char in raw_lookup:
    """This is a silly way to do it, but it works."""
    char,code=char[0],char[1:]
    code=code.replace("-----","x15 ")
    code=code.replace("----","x14 ")
    code=code.replace("---","x13 ")
    code=code.replace("--","x12 ")
    code=code.replace("-","x11 ")
    code=code.replace(".....","x05 ")
    code=code.replace("....","x04 ")
    code=code.replace("...","x03 ")
    code=code.replace("..","x02 ")
    code=code.replace(".","x01 ")
    code=code.replace("x0",'.')
    code=code.replace("x1",'-')
    code=code.split(" ")[:-1]
    #print char,code
    lookup[char]=code

```

<div class="text-center img-border">

[![](produced_thumb.jpg)](produced.png)

</div>

__Automated decoding__ is trivial. The image above was analyzed, turned into the image below, and the string (AJ4VD) was extracted:

### decode.py

```python
# given an image, it finds peaks and pulls data out
from PIL import Image
from PIL import ImageDraw
import pylab
import numpy

pixelSeek = 10
pixelShift = 15


def findPeak(data):
    maxVal = 0
    maxX = 0
    for x in range(len(data)):
        if data[x] > maxVal:
            maxVal, maxX = data[x], x
    return maxX


def peaks2morse(peaks):
    baseFreq = peaks[0]
    lastSignal = peaks[0]
    lastChange = 0
    directions = []
    for i in range(len(peaks)):
        if abs(peaks[i]-baseFreq) < pixelSeek:
            baseFreq = peaks[i]
        if abs(peaks[i]-lastSignal) < pixelSeek and i < len(peaks)-1:
            lastChange += 1
        else:
            if abs(baseFreq-lastSignal) < pixelSeek:
                c = " "
            if abs(baseFreq-lastSignal) < pixelSeek:
                c = " "
            if abs(baseFreq-lastSignal) < pixelSeek:
                c = " "
            directions.append(
                [lastSignal, lastChange, baseFreq, baseFreq-lastSignal])
            lastChange = 0
        lastSignal = peaks[i]
    return directions


def morse2image(directions):
    im = Image.new("L", (300, 100), 0)
    draw = ImageDraw.Draw(im)
    lastx = 0
    for d in directions:
        print d
        draw.line((lastx, d[0], lastx+d[1], d[0]), width=5, fill=255)
        lastx = lastx+d[1]
    im.show()


im = Image.open('raw.png')
pix = im.load()
data = numpy.zeros(im.size)
for x in range(im.size[0]):
    for y in range(im.size[1]):
        data[x][y] = pix[x, y]

peaks = []
for i in range(im.size[0]):
    peaks.append(findPeak(data[i]))

morse = peaks2morse(peaks)
morse2image(morse)
print morse
```

June 28th, 2010

Getting GTK, Glade, and Python to Work in Windows

These screenshots show me running the Py2EXE-compiled script I wrote last weekend on a Windows 7 machine. Additionally there is a screenshot of the "Add/Remove Programs" window demonstrating which versions of which libraries were required.

Markdown source code last modified on January 18th, 2021
---
title: Getting GTK, Glade, and Python to Work in Windows
date: 2010-06-28 09:46:26
tags: python, old
---

# Getting GTK, Glade, and Python to Work in Windows

__These screenshots__ show me running the Py2EXE-compiled script I wrote last weekend on a Windows 7 machine. Additionally there is a screenshot of the "Add/Remove Programs" window demonstrating which versions of which libraries were required.

<div class="text-center img-border">

[![](glade_exe_thumb.jpg)](glade_exe.png)

[![](needToInstall_thumb.jpg)](needToInstall.png)

</div>

June 24th, 2010

Fast TK Pixelmap generation from 2D Numpy Arrays in Python

I had TKinter all wrong! While my initial tests with PyGame's rapid ability to render Numpy arrays in the form of pixel maps proved impressive, it was only because I was comparing it to poor TK code. I don't know what I was doing wrong, but when I decided to give TKinter one more shot I was blown away -- it's as smooth or smoother as PyGame. Forget PyGame! I'm rendering everything in raw TK from now on. This utilizes the Python Imaging Library (PIL) so it's EXTREMELY flexible (supports fancy operations, alpha channels, etc).

The screenshot shows me running the script (below) generating random noise and "scrolling" it horizontally (like my spectrograph software does) at a fast rate smoothly (almost 90 FPS!). Basically, it launches a window, creates a canvas widget (which I'm told is faster to update than a label and reduces flickering that's often associated with rapid redraws because it uses double-buffering). Also, it uses threading to handle the calculations/redraws without lagging the GUI. The code speaks for itself.

import Tkinter
from PIL import Image, ImageTk
import numpy
import time


class mainWindow():
    times = 1
    timestart = time.clock()
    data = numpy.array(numpy.random.random((400, 500))*100, dtype=int)

    def __init__(self):
        self.root = Tkinter.Tk()
        self.frame = Tkinter.Frame(self.root, width=500, height=400)
        self.frame.pack()
        self.canvas = Tkinter.Canvas(self.frame, width=500, height=400)
        self.canvas.place(x=-2, y=-2)
        self.root.after(0, self.start)  # INCREASE THE 0 TO SLOW IT DOWN
        self.root.mainloop()

    def start(self):
        global data
        self.im = Image.fromstring('L', (self.data.shape[1],
                                         self.data.shape[0]), self.data.astype('b').tostring())
        self.photo = ImageTk.PhotoImage(image=self.im)
        self.canvas.create_image(0, 0, image=self.photo, anchor=Tkinter.NW)
        self.root.update()
        self.times += 1
        if self.times % 33 == 0:
            print "%.02f FPS" % (self.times/(time.clock()-self.timestart))
        self.root.after(10, self.start)
        self.data = numpy.roll(self.data, -1, 1)


if __name__ == '__main__':
    x = mainWindow()
Markdown source code last modified on January 18th, 2021
---
title: Fast TK Pixelmap generation from 2D Numpy Arrays in Python
date: 2010-06-24 21:29:29
tags: python, old
---

# Fast TK Pixelmap generation from 2D Numpy Arrays in Python

__I had TKinter all wrong!__ While my initial tests with PyGame's rapid ability to render Numpy arrays in the form of pixel maps proved impressive, it was only because I was comparing it to poor TK code. I don't know what I was doing wrong, but when I decided to give TKinter one more shot I was blown away -- it's as smooth or smoother as PyGame. Forget PyGame! I'm rendering everything in raw TK from now on. This utilizes the Python Imaging Library (PIL) so it's EXTREMELY flexible (supports fancy operations, alpha channels, etc).

<div class="text-center img-border">

[![](glade_python_improving_thumb.jpg)](glade_python_improving.png)

</div>

__The screenshot shows__ me running the script (below) generating random noise and "scrolling" it horizontally (like my spectrograph software does) at a fast rate smoothly (almost 90 FPS!). Basically, it launches a window, creates a canvas widget (which I'm told is faster to update than a label and reduces flickering that's often associated with rapid redraws because it uses double-buffering). Also, it uses threading to handle the calculations/redraws without lagging the GUI. The code speaks for itself.

```python
import Tkinter
from PIL import Image, ImageTk
import numpy
import time


class mainWindow():
    times = 1
    timestart = time.clock()
    data = numpy.array(numpy.random.random((400, 500))*100, dtype=int)

    def __init__(self):
        self.root = Tkinter.Tk()
        self.frame = Tkinter.Frame(self.root, width=500, height=400)
        self.frame.pack()
        self.canvas = Tkinter.Canvas(self.frame, width=500, height=400)
        self.canvas.place(x=-2, y=-2)
        self.root.after(0, self.start)  # INCREASE THE 0 TO SLOW IT DOWN
        self.root.mainloop()

    def start(self):
        global data
        self.im = Image.fromstring('L', (self.data.shape[1],
                                         self.data.shape[0]), self.data.astype('b').tostring())
        self.photo = ImageTk.PhotoImage(image=self.im)
        self.canvas.create_image(0, 0, image=self.photo, anchor=Tkinter.NW)
        self.root.update()
        self.times += 1
        if self.times % 33 == 0:
            print "%.02f FPS" % (self.times/(time.clock()-self.timestart))
        self.root.after(10, self.start)
        self.data = numpy.roll(self.data, -1, 1)


if __name__ == '__main__':
    x = mainWindow()
```

June 24th, 2010

Detrending Data in Python with Numpy

⚠️ SEE UPDATED POST: Signal Filtering in Python

While continuing my quest into the world of linear data analysis and signal processing, I came to a point where I wanted to emphasize variations in FFT traces. While I am keeping my original data for scientific reference, visually I want to represent it emphasizing variations rather than concentrating on trends. I wrote a detrending function which I'm sure will be useful for many applications:

def detrend(data,degree=10):
    detrended=[None]*degree
    for i in range(degree,len(data)-degree):
        chunk=data[i-degree:i+degree]
        chunk=sum(chunk)/len(chunk)
        detrended.append(data[i]-chunk)
    return detrended+[None]*degree

However, this method is extremely slow. I need to think of a way to accomplish this same thing much faster. [ponders]

UPDATE: It looks like I've once again re-invented the wheel. All of this has been done already, and FAR more efficiently I might add. For more see scipy.signal.detrend.html

import scipy.signal
ffty=scipy.signal.detrend(ffty)
Markdown source code last modified on January 18th, 2021
---
title: Detrending Data in Python with Numpy
date: 2010-06-24 08:38:52
tags: python, old
---

# Detrending Data in Python with Numpy

> **⚠️ SEE UPDATED POST:** [**Signal Filtering in Python**](https://swharden.com/blog/2020-09-23-signal-filtering-in-python/)

__While continuing my quest__ into the world of linear data analysis and signal processing, I came to a point where I wanted to emphasize variations in FFT traces. While I am keeping my original data for scientific reference, visually I want to represent it emphasizing variations rather than concentrating on trends. I wrote a detrending function which I'm sure will be useful for many applications:

```python
def detrend(data,degree=10):
	detrended=[None]*degree
	for i in range(degree,len(data)-degree):
		chunk=data[i-degree:i+degree]
		chunk=sum(chunk)/len(chunk)
		detrended.append(data[i]-chunk)
	return detrended+[None]*degree
```

<div class="text-center">

[![](detrend_fft_thumb.jpg)](detrend_fft.png)

</div>

However, this method is extremely slow. I need to think of a way to accomplish this same thing much faster. \[ponders\]

__UPDATE:__ It looks like I've once again re-invented the wheel. All of this has been done already, and FAR more efficiently I might add. For more see [scipy.signal.detrend.html](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.detrend.html)

```python
import scipy.signal
ffty=scipy.signal.detrend(ffty)
```

June 23rd, 2010

Insights Into FFTs, Imaginary Numbers, and Accurate Spectrographs

I'm attempting to thoroughly re-write the data assessment portions of my QRSS VD software, and rather than rushing to code it (like I did last time) I'm working hard on every step trying to optimize the code. I came across some notes I made about Fast Fourier Transformations from the first time I coded the software, and though I'd post some code I found helpful. Of particular satisfaction is an email I received from Alberto, I2PHD, the creator of Argo (the "gold standard" QRSS spectrograph software for Windows). In it he notes:

I think that [it is a mistake to] throw away the imaginary part of the FFT. What I do in Argo, in Spectran, in Winrad, in SDRadio and in all of my other programs is compute the magnitude of the [FFT] signal, then compute the logarithm of it, and only then I do a mapping of the colors on the screen with the result of this last computation.

Alberto, I2PHD (the creator of Argo)

UPDATE IN SEPTEMBER, 2020 (10 years later): I now understand that magnitude = sqrt(real^2 + imag^2) and this post is a bit embarrassing to read! Check out my .NET FFT library FftSharp for a more advanced discussion on this topic.

These concepts are simple to visualize when graphed. Here I've written a short Python script to listen to the microphone (which is being fed a 2kHz sine wave), perform the FFT, and graph the real FFT component, imaginary FFT component, and their sum. The output is:

Of particular interest to me is the beautiful complementary of the two curves. It makes me wonder what types of data can be extracted by the individual curves (or perhaps their difference?) down the road. I wonder if phase measurements would be useful in extracting weak carries from beneath the noise floor?

Here's the code I used to generate the image above. Note that my microphone device was set to listen to my stereo output, and I generated a 2kHz sine wave using the command speaker-test -t sine -f 2000 on a PC running Linux. I hope you find it useful!

import numpy
import pyaudio
import pylab
import numpy

### RECORD AUDIO FROM MICROPHONE ###
rate = 44100
soundcard = 1  # CUSTOMIZE THIS!!!
p = pyaudio.PyAudio()
strm = p.open(format=pyaudio.paInt16, channels=1, rate=rate,
              input_device_index=soundcard, input=True)
strm.read(1024)  # prime the sound card this way
pcm = numpy.fromstring(strm.read(1024), dtype=numpy.int16)

### DO THE FFT ANALYSIS ###
fft = numpy.fft.fft(pcm)
fftr = 10*numpy.log10(abs(fft.real))[:len(pcm)/2]
ffti = 10*numpy.log10(abs(fft.imag))[:len(pcm)/2]
fftb = 10*numpy.log10(numpy.sqrt(fft.imag**2+fft.real**2))[:len(pcm)/2]
freq = numpy.fft.fftfreq(numpy.arange(len(pcm)).shape[-1])[:len(pcm)/2]
freq = freq*rate/1000  # make the frequency scale

### GRAPH THIS STUFF ###
pylab.subplot(411)
pylab.title("Original Data")
pylab.grid()
pylab.plot(numpy.arange(len(pcm))/float(rate)*1000, pcm, 'r-', alpha=1)
pylab.xlabel("Time (milliseconds)")
pylab.ylabel("Amplitude")
pylab.subplot(412)
pylab.title("Real FFT")
pylab.xlabel("Frequency (kHz)")
pylab.ylabel("Power")
pylab.grid()
pylab.plot(freq, fftr, 'b-', alpha=1)
pylab.subplot(413)
pylab.title("Imaginary FFT")
pylab.xlabel("Frequency (kHz)")
pylab.ylabel("Power")
pylab.grid()
pylab.plot(freq, ffti, 'g-', alpha=1)
pylab.subplot(414)
pylab.title("Real+Imaginary FFT")
pylab.xlabel("Frequency (kHz)")
pylab.ylabel("Power")
pylab.grid()
pylab.plot(freq, fftb, 'k-', alpha=1)
pylab.show()

After fighting for a while long with a "shifty baseline" of the FFT, I came to another understanding. Let me first address the problem. Taking the FFT of different regions of the 2kHz wave I got traces with the peak in the identical location, but the "baselines" completely different.

Like many things, I re-invented the wheel. Since I knew the PCM values weren't changing, the only variable was the starting/stopping point of the linear sample. "Hard edges", I imagined, must be the problem. I then wrote the following function to shape the PCM audio like a triangle, silencing the edges and sweeping the volume up toward the middle of the sample:

def shapeTriangle(data):
    triangle=numpy.array(range(len(data)/2)+range(len(data)/2)[::-1])+1
    return data*triangle

After shaping the data BEFORE I applied the FFT, I made the subsequent traces MUCH more acceptable. Observe:

Now that I've done all this experimentation/thinking, I remembered that this is nothing new! Everyone talks about shaping the wave to minimize hard edges before taking the FFT. They call it windowing. Another case of me re-inventing the wheel because I'm too lazy to read others' work. However, in my defense, I learned a lot by trying all this stuff -- far more than I would have learned simply by copying someone else's code into my script. Experimentation is the key to discovery!

Markdown source code last modified on January 18th, 2021
---
title: Insights Into FFTs, Imaginary Numbers, and Accurate Spectrographs
date: 2010-06-23 22:21:00
tags: qrss, python, old
---

# Insights Into FFTs, Imaginary Numbers, and Accurate Spectrographs

__I'm attempting to thoroughly re-write the data assessment__ portions of my QRSS VD software, and rather than rushing to code it (like I did last time) I'm working hard on every step trying to optimize the code. I came across some notes I made about Fast Fourier Transformations from the first time I coded the software, and though I'd post some code I found helpful. Of particular satisfaction is an email I received from Alberto, I2PHD, the creator of Argo (the "gold standard" QRSS spectrograph software for Windows). In it he notes:

<blockquote class="wp-block-quote"><p>I think that [it is a mistake to] throw away the imaginary part of the FFT. What I do in Argo, in Spectran, in Winrad, in SDRadio and in all of my other programs is compute the magnitude of the [FFT] signal, then compute the logarithm of it, and only then I do a mapping of the colors on the screen with the result of this last computation.</p><cite> Alberto, I2PHD (the creator of Argo)</cite></blockquote>

> __UPDATE IN SEPTEMBER, 2020 (10 years later):__ I now understand that `magnitude = sqrt(real^2 + imag^2)` and this post is a bit embarrassing to read! Check out my .NET FFT library [FftSharp](https://github.com/swharden/FftSharp) for a more advanced discussion on this topic.

__These concepts are simple__ to visualize when graphed. Here I've written a short Python script to listen to the microphone (which is being fed a 2kHz sine wave), perform the FFT, and graph the real FFT component, imaginary FFT component, and their sum. The output is:

<div class="text-center">

[![](real_imaginary_fft_pcm_thumb.jpg)](real_imaginary_fft_pcm.png)

</div>

__Of particular interest__ to me is the beautiful complementary of the two curves. It makes me wonder what types of data can be extracted by the individual curves (or perhaps their difference?) down the road. I wonder if phase measurements would be useful in extracting weak carries from beneath the noise floor?

<div class="text-center">

[![](fft_base2_thumb.jpg)](fft_base2.png)

</div>

__Here's the code I used to generate the image above.__ Note that my microphone device was set to listen to my stereo output, and I generated a 2kHz sine wave using the command `` speaker-test -t sine -f 2000 `` on a PC running Linux. I hope you find it useful!

```python
import numpy
import pyaudio
import pylab
import numpy

### RECORD AUDIO FROM MICROPHONE ###
rate = 44100
soundcard = 1  # CUSTOMIZE THIS!!!
p = pyaudio.PyAudio()
strm = p.open(format=pyaudio.paInt16, channels=1, rate=rate,
              input_device_index=soundcard, input=True)
strm.read(1024)  # prime the sound card this way
pcm = numpy.fromstring(strm.read(1024), dtype=numpy.int16)

### DO THE FFT ANALYSIS ###
fft = numpy.fft.fft(pcm)
fftr = 10*numpy.log10(abs(fft.real))[:len(pcm)/2]
ffti = 10*numpy.log10(abs(fft.imag))[:len(pcm)/2]
fftb = 10*numpy.log10(numpy.sqrt(fft.imag**2+fft.real**2))[:len(pcm)/2]
freq = numpy.fft.fftfreq(numpy.arange(len(pcm)).shape[-1])[:len(pcm)/2]
freq = freq*rate/1000  # make the frequency scale

### GRAPH THIS STUFF ###
pylab.subplot(411)
pylab.title("Original Data")
pylab.grid()
pylab.plot(numpy.arange(len(pcm))/float(rate)*1000, pcm, 'r-', alpha=1)
pylab.xlabel("Time (milliseconds)")
pylab.ylabel("Amplitude")
pylab.subplot(412)
pylab.title("Real FFT")
pylab.xlabel("Frequency (kHz)")
pylab.ylabel("Power")
pylab.grid()
pylab.plot(freq, fftr, 'b-', alpha=1)
pylab.subplot(413)
pylab.title("Imaginary FFT")
pylab.xlabel("Frequency (kHz)")
pylab.ylabel("Power")
pylab.grid()
pylab.plot(freq, ffti, 'g-', alpha=1)
pylab.subplot(414)
pylab.title("Real+Imaginary FFT")
pylab.xlabel("Frequency (kHz)")
pylab.ylabel("Power")
pylab.grid()
pylab.plot(freq, fftb, 'k-', alpha=1)
pylab.show()
```

__After fighting for a while long with__ a "shifty baseline" of the FFT, I came to another understanding. Let me first address the problem. Taking the FFT of different regions of the 2kHz wave I got traces with the peak in the identical location, but the "baselines" completely different.

<div class="text-center">

[![](fft_base3_thumb.jpg)](fft_base3.png)

</div>

__Like many things, I re-invented the wheel.__ Since I knew the PCM values weren't changing, the only variable was the starting/stopping point of the linear sample. "Hard edges", I imagined, must be the problem. I then wrote the following function to shape the PCM audio like a triangle, silencing the edges and sweeping the volume up toward the middle of the sample:

```python
def shapeTriangle(data):
    triangle=numpy.array(range(len(data)/2)+range(len(data)/2)[::-1])+1
    return data*triangle
```

__After shaping the data BEFORE I applied the FFT,__ I made the subsequent traces MUCH more acceptable. Observe:

__Now that I've done all this experimentation/thinking,__ I remembered that this is nothing new! Everyone talks about shaping the wave to minimize hard edges before taking the FFT. They call it _windowing._ Another case of me re-inventing the wheel because I'm too lazy to read others' work. However, in my defense, I learned a lot by trying all this stuff -- far more than I would have learned simply by copying someone else's code into my script. Experimentation is the key to discovery!
Pages