Understanding Digital Signal Processing Chapter 3

Chapter 3

The continuous Fourier transform is given by.

$X(f) = \int^{\infty}_{-\infty} x(t) e^{-j \, 2 \pi \, f \, t}\mathrm{dt}$

DFT: In the discrete case, this becomes a sum, where $t \rightarrow n$, and $f \rightarrow m/N$.

$X(m) = \sum^{N-1}_{n = 0} x(n) e^{-j \, 2 \pi \, n \, m/N}$

Another way to get to the discrete Fourier transform: Pick a test-frequency $m/N$. Multiply a cosine and a sine (or complex exponential, your choice) that has that test-frequency with your signal data, and sum all $N$ terms up. If the original signal had a component at that frequency, the net-sum will be non-zero (basically a DC-term when multiplying signals that have the same frequency). Or, in the words of the original author:

Any individual X(m) output value is nothing more than the sum of the term-by-term products, a correlation, of an input signal sample sequence with a cosine and a sinewave whose frequencies are m complete cycles in the total sample interval of N samples.

The (discrete) Fourier Transform of a real-valued signal is symmetric, as can be seen by sending $m \rightarrow N-m$ in the equation above, and then dropping the $2 \pi n$ term in the exponential, which leads to

$X(N-m) = X^*(m)\quad .$

The DFT of a symmetric signal, where $x(n) = x(N-n)$, is real. You can either see this by reversing the index of summation in the DFT's definition, or intuitively: A symmetric signal can be described with cosines only, so the amplitudes of the sines are all 0.

In our definition of the DFT the calculated phases of the frequency components are always relative to cosines.

If you have a real-valued signal and want to know the amplitude of one of its frequency components $f_n < f_{\mathrm{s}}/2$, then you have to divide the calculated (DFT) amplitude by $N/2$. You can see why this factor is $N/2$ when you consider what happens for some test signal. Each point in time contributes with 1/2. (Again, this is only true for real-valued signals, where the missing 1/2 is at the other end of the spectrum!) Side note: the author is back at all-positive frequencies here... There are other definitions of DFT's where the scale factor $N$ is distributed somewhere else, but that's not important to us.

The "fundamental frequency" or frequency resolution when computing a DFT of a signal that has length $N$ and sample rate $f_{\mathrm{s}}$ is given by

$f_{\mathrm{res}} = \frac{f_{\mathrm{s}}}{N}$

The inverse DFT follows from the view that "the signal is just a bunch of harmonic signals added together" and that we already know the magnitude and phases of those signals. It differs from the DFT only in the sign of the exponential, and in the scaling factor $1/n$.

$x(n) = \frac{1}{N}\sum^{N-1}_{m = 0} X(m) e^{j \, 2 \pi \, n \, m/N}$

A time-shift of the input signal merely results in a phase-shift of the DFT.

Leakage: If the signal-frequency is not one of the DFT-frequency "bins" $f_{\mathrm{bin}}(m) = m f_{\mathrm{s}} /N$, signal energy "leaks" into all bins. Example: A real cosine signal with $k \in \Re^+$ cycles in an $N$-point sample and amplitude $A_0$ has DFT-amplitudes given by the sinc-function

$X(m) = \frac{A_0 N}{2} \cdot \frac{\sin[\pi(k-m)]}{\pi(k-m)} \quad .$

This makes sense, because due to only looking at N points, we are looking at a "burst" that has a FT given by a sinc.

Leakage can be surprisingly asymmetric, especially if $f_{\mathrm{signal}}$ is near 0 or $f_{\mathrm{s}}/2$, because the leaking spectrum wraps around (more like - is reflected? No!) at these corner frequencies.

Windowing is a technique to mitigate leakage. The input signal is multiplied (scalar product) with an input window, for instance a triangle function or a "raised cosine" (Hamming window, better than triangle) function. Tradeoff: Peaks become broader, and overall you lose signal (because start and end of signal tend towards 0). However, it's now easier to spot small signals in the vicinity of big ones. Other window functions: Chebyshev, Kaiser.

Scalloping is somewhat related to leakage: For signal frequencies directly between two DFT bins, the max(magnitude) drops by more than 30% (-4 dB), because the signal leaks into other bins. Windowing improves this effect as well (example: -2 dB for Hamming). Real-world signals are often wider than 1 bin, making this not very relevant.

Zero-Padding is used to increase the frequency resolution by appending 0's to the sample. One case where this makes sense if the signal really is 0 outside of the sampled interval. But you can also think about it in a different way: When zero-padding, you are basically "probing" your signal with non-integer frequencies, which is of course absolutely ok. In the extreme limit, this becomes the discrete time Fourier transform (DTFT), which has a continuous output. Careful: When combining zero-padding with windowing, only window the original $L$ non-zero samples. Note that the calculated DFT amplitudes will go down with the total, padded length $N$.

DTFT stands for "discrete time Fourier transform" and is basically zero-padding taken to the extreme, such that the transform's output is continuous. This does not improve resolution, which due to the width of the (rectangular or funky looking) "input window" surrounded by 0's is limited by sidelobes to roughly $f_{\mathrm{res}} = 1/T_\mathrm{sample}$.

SNR is output signal power / average noise power.

You can think of a DFT's frequency bins as narrow bandpass filters. When you increase the sample length $N$, the average noise amplitude increases by $\sqrt{N}$, but the measured amplitude of a (coherent!) signal increases by $N$. Thus, SNR increases by $\sqrt{N}^2 = N$ (assuming that the power in the signal bin is dominated by the signal, not the noise). The book does not give a concise definition of processing gain - only that due to the effects described above it increases by 3dB per doubling of $N$.

The DFT of a rectangular pulse is always a (sampled) sinc-function. The region around the central peak is called the Dirichlet kernel. Assuming a pulse with $K$ non-zero values of amplitude $A_0$, the peak in the DFT has a height of $A_0 K$ an a full width (to 0 amplitude) of $2 N/K$. In the book's example, the author unexptectedly uses negative indexes and frequencies. Trying the example with Matlab, I get the same amplitudes, but the phases are quite different, which is due to his sneaky redefinition of the DFT in equation 3-36 (the n in the exponent should be shifted, if he decides to shift the summation index as he did!).

Further, the DFT of a rectangular pulse has a phase-factor and an amplitude factor. The amplitude factor only depends on the amplitude, sample-length, and number of (consecutive) non-zero samples. In other words, if you shift the onset of the pulse, you only get phase shifts in the frequency domain.

The frequency axis can be expressed in many different ways: In frequency $f$ with a resolution $f_{\mathrm{s}}/m$, in normalized frequency or "cycles per sample" $f/f_{\mathrm{s}}$ with resolution $1/N$, or in radians per sample $\omega$ with resolution $2 \pi/N$. Note that in the latter case $\omega$ has the non-standard dimension radians/sample!!

Understanding Digital Signal Processing Chapter 1 - 2

Chapter 1

Analog usually means "continuous". The word probably originates from analog computers, where each signal was "analog" to some real-world value.

Discrete signals are usually discrete in both time and value.

Nonlinear Systems lead to intermodulation. (Always?)

LTI systems are commutative. (Always!)

Important LTI system: The time delay.

If the linear system's parameters do not change over time, it does not mean it's LTI. What is important is that it is shift-invariant.

Magnitude = Abs(Amplitude)

Chapter 2

Aliasing: The spectrum of a discrete signal repeats in intervals of the sampling frequency $f_{\mathrm{s}}/2$, from "DC to Daylight".

Always use a low-pass filter with $f_{\mathrm{corner}} < f_{\mathrm{s}}/2$ to prevent aliasing.

Bandpass-sampling = IF-sampling = under-sampling = deliberate abuse of aliasing to look at some signal with $f_{\mathrm{signal}} > f_{\mathrm{s}}/2$. This degrades SNR, as all noise is folded into the $[-f_{\mathrm{s}}/2; f_{\mathrm{s}}/2]$ band.

Apparently, negative frequencies are all the rage. Does he not want to use negative amplitudes? Explanation is promised in chapters 3 and 8.

Bash scripts for Rohde&Schwarz ZVL network/spectrum analyzer

These are part of my efforts on remote controlling LXI enabled instruments using bash scripts. I started by writing the following wrapper for netcat - that way I can easily pipe stuff in and out. I named it "zvlcat" and placed it in /usr/local/bin/

#!/bin/bash
# Script for writing commands to and reading responses from
# a ZVL network analyzer. 
# Reading responses might fail if the response time is 
# greater than 1s.
#
# Usage: 
# zvlcat "some_command"  
# Sends "some_command" to ZVL. If cmd_str contains a "?",
# then it must be a query and the script waits for up to 
# 1 second for the reponse.
# The argument is optional and defaults to "*IDN?".

cmd_str="*IDN?"
opt_str=

# If argument was given, use it as command
if [ $# -gt 0 ]; then
    cmd_str=$1
fi

# Count number of question marks in command string
#  grep -o: put each find in new line
#  wc -l: count number of lines
num_quest=$(echo $cmd_str | grep -o "?" | wc -l)

# If there is a ? in the command, it's a query
# -> We should wait for 1 second for the reply.
if [ $num_quest -gt 0 ]; then
    opt_str="-q 1"
fi

echo $cmd_str | netcat $opt_str zvl.dhcp.mpi-hd.mpg.de 5025

I can then control the instrument using scripts such as this one:

#!/bin/bash

# This script performs a measurement of the ion amplitude
# using the ZVL network analyzer.

# Parameters
center="4.0565MHz"
rbw=1
aver_coun=50


# Reset instrument, activate spectrum analyser mode
zvlcat "*RST; INST:SEL SAN"

# single sweep mode, only sweep current channel
zvlcat "INIT:CONT OFF; SCOP SING"

# center frequency, span
zvlcat "FREQ:CENT $center; SPAN 100 Hz"

# filter type fft
zvlcat "BAND:TYPE FFT"

# bandwidth (seperate line because fft must be active at this point)
zvlcat "BAND $rbw"

# Set trace mode to average (Trace->Trace Mode->Average)
zvlcat "AVER ON; AVER:COUN $aver_coun"

# Set Range and Reference level
zvlcat "DISP:TRAC:Y 20; Y:RLEV -60"

# Calculate measurement time
swe_time=$(zvlcat "SWE:TIME?")
# This will be a float, but that's ok
meas_time=$(echo "($aver_coun * $swe_time)" | bc)

# Start single sweep, wait until it's finished
# (Note: This is the only "overlapping" command; needs *WAI)
zvlcat "INIT:IMM; *WAI"
sleep $meas_time

# Print to "Printer 1"
# zvlcat "HCOP 1"

# Read out the data (takes 1 second), write to file
# zvlcat "TRAC? TRACE1" > trace.dat

# I want to measure signal to noise. This is complicated, because
# I am not allowed to average over logarithms! Therefore I am using
# ZVL functions to do that

# Switch on marker
zvlcat "CALC:MARK:STAT ON"

# Put marker on peak
zvlcat "CALC:MARK:MAX"

# Read out marker
max_pow=$(./zvlcat "CALC:MARK:Y?")

# Activate noise measurement
zvlcat "CALC:MARK:FUNC:NOIS ON"

# Set marker lower by three divisions
zvlcat "CALC:MARK:X DOWN"
zvlcat "CALC:MARK:X DOWN"
zvlcat "CALC:MARK:X DOWN"

# Read out noise
nois_pow=$(zvlcat "CALC:MARK:FUNC:NOIS:RES?")

# Append values to file, including a timestamp
echo "RBW was $rbw"
echo "signal: $max_pow dBm"
echo "noise: $nois_pow dBm/Hz"
echo "$(date); $max_pow; $nois_pow" >> log_ion_vs_noise.txt

# Back this up.
# TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO TO DO

# Turn on continous measurement again
zvlcat "INIT:CONT ON"

# Return to local control
zvlcat "@LOC"


###################
## DATA ANALYSIS ##
###################

# sed s/,/\\n/g trace.dat > rows.dat

PTS Frequency Synthesizers

Amongst the dinosaurs in our labe are "Programmed Test Sources" (PTS) frequency synthesizers from the late 1970's - early 1980's. Earlier this week I observed that for certain frequency settings (26 612 123.818 000 000 Hz +- a few mHz, for example), their phase noise becomes pretty big. So I dug into the (excellent!) documentation to find out how this can be. Did we miss this behaviour entirely for the last 4 years?

No. As it turns out, the problem was intermittent and went away after two days... But I did learn how they operate:

A 10 MHz reference signal (either internal or external source) is transformed (analog) 14 MHz, 16 MHz, 18 MHz, 20 MHz, 22 MHz, 112 MHz, and 113 MHz. These are then used to form a 14 MHz carrier signal that carries all the fine digits up to the 100kHz dial using several identical sections. Each section takes the 14.X MHz carrier of the previous section (first one takes pristine 14 MHz), and then adds a combination of the 14 MHz, ... 113 MHz signals such that the total signal is between 140.X MHz and 149.X MHz. The specific combination depends on the dial setting for the particular digit, so that you end up with 143.X MHz if the current dial was set to 3. Finally, the section divides the carrier by 10 so that the output is between 14 MHz and 15 MHz again (in this case, 14.3X MHz). They daisy-chain 7 of these modules, with an option to connect an external source to the 14MHz input of the first section, so that you can daisy-chain two or more synthesizers and thereby get arbitrary precision.

PTS200 Overview

PTS200 Fine Section

As of 2014, the company still exists and sells frequency synthesizers with impressive specs (just look at that phase noise!), but by now they mostly switched to direct digital synthesis for the lower digits (although their FAQ mentions the analog mixing as an option).

As a side note: The other, popular way to synthesize arbitrary frequencies is to use a phase locked loop (PLL): You run a voltage controlled oscillator (VCO), divide it's output by a user-set number so that the divider-output coincides with the frequency of the reference oscillator. You compare the phase of the ouput with the reference and servo the VCO to keep the phase difference constant.

Agilent 33600A: Output off means output off

Lower-end agilent synthesizers such as the 33250A have the quirk that they don't actually switch the output off - they just attenuate it by 50dB. I got to play around with the new 33600A series (up to 120 MHz, 2 channel), and this one is different: Off means attenuate by 134dBm, which is "off" in all practical cases I can think off.

Instrument control using SCPI and Bash scripts

You don't need a bloated LabView installation for controllilng your lab instruments. If they have a LAN connector, you can probably talk to them using simple bash (or cmd or pyton or ...) scripts.

Instrument Requirements

• LAN-port
• SCPI support (see here for the shortest introduction possible)
• Working network connection
  • Try to find out the instrument's IP-address and on what ports it listens on. Standard port is 5025.
  • Careful: Often there is no username/password, so anyone in your network can control your instrument! I suggest to set up a dedicated lab subnet.

High level overview

• Hook up instrument to LAN
• Connect to port 5025
• Send SCPI commands
• Receive data and redirect it to a file

Hooking up instrument to LAN

Plug a LAN connector into your lab instrument and power it on. Many lab instruments support DHCP and will automatically get an IP address. Find the IP-adress by going through all instrument menus that are related to "SYSTEM" or "REMOTE".

Connecting to port 5025: Ubuntu

You can use telnet for a test connection (see Windows section below for an example). It is available on most Windows and Linux installations. But on Linux (including Ubuntu), the prettier alternative is netcat.

Open a terminal and install netcat (should be installed already):

sudo apt-get install netcat

Open a netcat session to the lab instrument:

netcat 149.100.100.5 5025

Here, 149.100.100.5 is the instrument's IP adress, and 5025 ist the port that we want to connect to on the instrument. Port 5025 is the standard port for SCPI commands.

Query the instrument for its ID number by typing

*IDN?

Press enter and the device will return its ID. Press CTRL+C to exit the netcat session.

You can also run netcat non-interactively by sending input commands through a "pipe" |:

echo "*IDN?" | netcat -q 1 149.100.100.5 502

The option "-q 1" makes netcat wait for 1 second for the answer from the instrument. You only need this option for queries. (You don't strictly need it - but I had a few queries fail before using that wait option.)

Redirect netcat's output to a file:

echo "*IDN?" | netcat -q 1 149.100.100.5 5025 > ans.txt

Then type

cat ans.txt

and you will see the contents of the file.

To set something, you can use

echo "FREQUENCY:CENTER 100MHz" | netcat 149.100.100.5 5025

if your device understands that command, its center frequency will now be set to 100MHz. This command has no output, so there is no point in redirecting it to a file.

Connecting to port 5025: Windows

On Windows, this is a bit harder. PuTTY comes with a command line utility called plink - but I can't get it to open raw connections. Installing netcat is tricky because antivirus software classifies it as risky (it's often included in malware). Telnet always works (also in corporate environments with strict antivirus settings), but telnet does not support scripting. We'll script it anyways.

Note: On Windows 7, you have to activate telnet first.

In a cmd prompt, type

telnet

It will greet you by saying something like "Welcome - the escape combination is CTRL++. (This is different for each Locale.) You will need this escape sequence later.

Type

o 149.100.100.5 5025

The bottom of the screen will now keep saying "Connecting". Don't be fooled by this - most likely you are already connected. Type

*IDN?

and press enter. The device should answer with its ID number, but it may take two or three tries, because backspace doesn't work in Windows telnet.

Quitting: Press your escape key combination (in my case pressing CTRL and + simultaneously) to get to the telnet prompt. Type q, then enter, to quit.

Let's script this. Telnet, by default, does not accept input pipes and cannot be scripted. We'll be using a visual basic script that simulates key presses. It's easier than it sounds.

Create a file called "telnet_commands.vbs". Paste the following lines into it:

set OBJECT=WScript.CreateObject("WScript.Shell")
WScript.sleep 50 
OBJECT.SendKeys "*IDN?{ENTER}" 
WScript.sleep 500 
OBJECT.SendKeys "^({+})" 
WScript.sleep 50 
OBJECT.SendKeys "q{ENTER}" 

• set OJBECT=... Create a "shell" object that can do system stuff. This can be used to simulate key presses.
• WScript.sleep 50 Wait 50ms.
• OBJECT.SendKeys ... Simulate button presses. These go to whatever program that has "focus" (is active) right now.
• OBJECT.SendKeys "^({+})" Press control key ^, keep it pressed () while pressing plus key {+}. The plus sign can have special meaning, which is why it's escaped by curly braces.

Create another file called query_idn.bat. Paste the following lines into it:

@echo OFF
start telnet 149.100.100.5 5025 -f output.txt
cscript telnet_commands.vbs

• @echo OFF Don't echo (repeat) this command in the output (the @ at the beginning), and switch echoing off for all subsequent commands.
• start telnet... Start a telnet session in its own window, write the output to a file
• cscript ... Execute the vbs file

Double-click the bat file and rejoice. The first line of the output file output.txt is bogus junk, but the second line has the interesting stuff.

Escaping WordPress Shortcodes

Shortcodes, which allow you to easily insert galleries, videos, ... into a blog post, are written in square brackets:

[gallery]

When writing about shortcodes, they are unfortunately parsed even when inside a <code> </code> environment. You can escape this by enclosing a valid shortcode in double brackets:

[[gallery]]

MathJax and Wordpress

Install the mathjax-latex plugin.

In visual editor, put

[mathjax]

anywhere on the page. This will enable LaTeX parsing so that you can use LaTeX style formula formatting.

Formulas that stand on their own lines can be writte LaTeX style \[ c = \frac{a}{b} \] or TeX style $$c = \frac{a}{b}$$. Both render like this:

$c = \frac{a}{b}$

Due to MathJax limitations, inline formulas can olny be written in LaTeX style \( a^2 + b^2 = c^2 \), which renders $a^2 + b^2 = c^2$ in the same line.

Note: The author of the mathjax-latex plugin prefers to have all formulas enclosed by [latex] [/latex] shortcodes, but I find the LaTeX style easier to write.