Probabilistic Neural Network Classifier for Morse Code

One of the challenges in Morse Code decoder is how to build a well working adaptive classifier for incoming symbols given the variability in  "dit/dah" timing in real world CW traffic received in ham bands. I have discussed these problems in my previous blog posting.

OVERVIEW 

I have collected data from different ham QSOs in order to understand the underlying patterns and timing relationships.  One good way to represent this data is to think Morse Code as "symbols" instead of  "dits" and "dahs".   A symbol  is  a  ["tone",  " silence"]  duration pair.  For example the  letter  A   ( . - )  could be represented as two symbols   ["dit" - "ele"]  + ["dah", "chr"]  where "dit" and "dah"  represent short and long tone; respectively "ele" represent inter-element space  and "chr" represents inter-character space.  Morse code can then be represented as a sequence of symbols - letter "A" would be {S0, S4}  -  I am using the following definitions as shorthand for symbols:

• S0 = [dit, ele]  //  dit  and inter-element space 
• S1 = [dit, chr]  //  dit  and inter-character space 
• S2 = [dit, wrd]  //  dit  and inter-word space 
• S3 = [dah, ele]  //  dah  and inter-element space 
• S4 = [dah, chr]  //  dah  and inter-character space 
• S5 = [dah, wrd]  //  dah  and inter-word space 

CLASSIFICATION EXAMPLES 

If we plot  thousands of these symbols as  data vectors  in  (x, y) chart where x-axis  is  tone duration ("mark")  and y-axis is silence duration ("space") we get a picture like figure 1 below.  "Dit" duration is scaled to 0.1 in the picture and "dah" duration is  0.3 respectively.  Figure 1. was created from a recorded test audio file with  -3 dB  SNR by courtesy of Dave W1HKJ.

You can easily see different symbol clusters and there is clean separation between them.  In the "mark" dimension  (horizontal x-axis)  "dit" and "dah"  follow roughly  1:3  ratio as expected.

Dark blue and green areas represent S0 and S3  symbols - these are symbols within Morse letters as "space" value is around 0.1.

Red and brown areas represent S1 and S4 symbols  - longer inter-character "space" value ~0.3 means that these symbols are the last symbol in a letter.

Yellow and light blue areas represent S2 and S5 symbols - long "space" value around 0.6 - 0.7 shows that these are the last symbol in a word.  

Figure 1.  Morse code symbols  ( -3 dB SNR test audio file) 

While figure 1. looks nice and neat please take a look at figure 2. below.   Using the same exact classifier settings  I recorded several stations working on 7 Mhz band  and collected some 1800 symbols. This picture shows data from multiple stations - therefore the variation for each symbol is much larger.

There is much more variation in the "space" dimension.  When listening the Morse code  this is consistent with people having "thinking pauses" as they form an idea of the next words to send. Sometimes these pauses can be even longer - you can see those when a new subject is introduced in the QSO such as asking a question or responding to something that other station sent.
I have capped the duration to 1.0  ( this is 10 *  "dit" duration) for implementation reasons. In real life these pauses can be few seconds.

I created another symbol S6  that I have used for  noise spike detection. These random spikes are typically less than 1/2 "dit" duration.  Having a separate symbol for noise makes it easier to remove them when decoding characters.  

Figure 2.  Morse code symbols (7 Mhz band ham traffic)

PNN AS A CLASSIFIER 

Finding a classifier algorithm that would be able to  identify each symbol in a consistent manner was not very easy.  I used Weka software  to test different algorithms trying to find a classifier that would classify the above symbols with over 90% accuracy.  After spending many hours of running different data sets through Weka  I was quite disappointed.  None of the built-in classifiers seemed to handle the task well or I was not able to provide good enough parameters as a starting point.

I was exchanging some emails with Dave, N7AIG on how to build a classifier for Morse code symbols while he suggested that I should look at Probabilistic Neural Network algorithm.

I started reading articles on PNN and discovered that it has several great benefits.  PNN is very fast in learning new probability distributions, in fact  it is so called "one-pass" learning algorithm.  This is a great benefit as many of the more traditional neural network algorithms require  hundreds or even thousands of training cycles before they converge to a solution.  With PNN you just have to show a few examples of each symbol class and it does the classification almost immediately.

PNN networks are also relatively insensitive to outliers.  So few bad examples does not kill the classification performance. PNN networks also approach Bayes optimal classification and are guaranteed to converge to an optimal classifier as the size of the representative training set increases. Training samples can also be added or removed without extensive retraining.

PROTOTYPING WITH EXCEL 

PNN has only one parameter "sigma"  that needs tuning.  I built an Excel sheet (copy in Google Sheets) to validate the PNN algorithm  and to see how different "sigma" values  impact the classification.  By changing the "sigma" value I could easily see how the classifier handles overlapping distributions. Sigma basically determines the Gaussian shape of the PNN equation:

Figure 3.   PNN equation

Experimenting with different "sigma" values  I was able to discover that  0.03 seems to be a good compromise value.  If you increase the value much higher you tend to get more classification errors due to overlaps.  If you make it much smaller you get more classification errors due to too narrow distribution.

See figure 4.  on how the probability distribution function (PDF) looks like for "dit", "dah" and "word space" with sigma = 0.03.  The scale has been set so that "dit" duration corresponds to 0.1  like in the previous pictures.

Figure 4.  PNN tuning

IMPLEMENTATION IN C++ FOR FLDIGI

The final part of this work was to create PNN algorithm in C++  and integrate it with  FLDIGI software.  The alpha version of the classifier code is in Appendix below - I have also included the Examples[]  data I have used to classify results in Figure 1 and 2.

Figure 5. Probabilistic Neural Network 

Figure 5.  demonstrates  the role of different "layers"  in PNN algorithm.  Input layer gets the symbol  ("mark", "space") values that are scaled  between 0.0 and 1.0. The pattern layers calculates the match against each symbol class training example. In this example there are two training examples per symbol class.  The summation layer calculates the sum of probability distribution function (PDF - see Fig 3) for each class.  Finally the decision layer compares all these results and selects the largest sum value that is the best matching symbol class.

Part of my current work has been also to completely rewrite the state machine in FLDIGI  CW.CXX module that handles the conversion from KEY_DOWN and KEY_UP events to  "mark" and "space" durations.



I have also rewritten the Morse decoder part - new decoder gets stream of classified symbols from PNN() function and runs through a state machine to come up with the correct character. I am still in process of getting this state machine to behave like a Viterbi  decoder.

For this I would need the a-priori probabilities for each state and traverse through the states using the symbols and their probabilities to calculate the optimal decoding path.  I am struggling a bit with this step currently - it is not obvious to me how I should use the constraints in Morse code to derive the best path.  If you have experience in Viterbi decoders I would like to chat with you to get better insight how this algorithm should really work.

73
Mauri  AG1LE

APPENDIX  -  C++ SOURCE CODE 

#define Classes   7 // Symbols S0 ... S5  - 6 different classes + S6 for noise 
#define NEPC 2 // Number of Examples per Class
#define Dimensions 2 // dimensions - use mark & space here - could add more features like AGC 

const float Examples[Classes][NEPC][Dimensions] =  {
{{0.1, 0.1}, // S0 = dit-ele
{0.08, 0.12}}, // S0 = dit-ele

{{0.1, 0.3}, // S1 = dit-chr
{0.08, 0.23}}, // S1 = dit-chr

{{0.1, 0.7}, // S2 = dit-wrd
{0.1,  0.6}}, // S2 = dit-wrd

{{0.3, 0.1}, // S3 = dah-ele
{0.34, 0.16}}, // S3 = dah-ele

{{0.34, 0.4}, // S4 = dah-chr
{0.29, 0.22 }}, // S4 = dah-chr

{{0.23, 0.54}, // S5 = dah-wrd
{0.23,  0.9}}, // S5 = dah-wrd

{{0.01, 0.01}, // S6 = noise-noise
{0.05,  0.05}} // S6 = noise-noise

};



//=======================================================================
// (C) Mauri Niininen AG1LE 
// PNN()  is a Probabilistic Neural Network algorithm to find best matching Symbol 
// given received Mark / Space pair. Mark and Space are scaled to [0 ... 1.0] where 
// Mark: standard  'dit' length is  0.1, 'dah' is 0.3  
// Space: standard element space is 0.1, character space 0.3 and word space 0.7 
// Examples[] contains 2 examples for each symbol class
//=======================================================================
int symbol::PNN (float mark, float space) 
{
float sigma = 0.03; // SIGMA determines the Gaussian shape  z = f(x,y) = exp (-((x-x0)^2+(y-y0)^2)/2*sigma^2)
int classify = -1;
float largest = 0;
float sum [Classes];
float test_example[2];


if (abs(mark) > 1) mark = 1;
if (abs(space) > 1) space = 1; 

test_example[0] = mark;
test_example[1] = space;

// OUTPUT layer  - computer PDF for each class c  
for (int k=0; k < Classes; k++) 
{
sum[k] = 0;

// SUMMATION layer - accumulated PDF for each example from particular class k 
for (int i=0;i < NEPC; i++) 
{
float product = 0; 
// PATTERN layer - multiplies test examples by the weights 
for (int j=0; j < Dimensions; j++)
{
product += (test_example[j] - Examples[k][i][j])*(test_example[j] - Examples[k][i][j]);
}
product = -product/(2*(sigma * sigma));
product= exp(product); 
sum[k] += product;
}

sum[k] /= NEPC;
}
// sum[k] has accumulated PDF for each class k 
for (int k=0; k<Classes; k++) 
{
if ( sum[k] > largest) 
{
largest = sum[k];
classify = k;
}

}
 
learn (classify, mark, space);

return classify;
} 

Towards Bayesian Morse Decoder

I had an email exchange with Alex, VE3NEA who is the author of the popular CW Skimmer software. Since CW Skimmer is well known for its ability to decode Morse code even in tough pile-up, noise and contest situations I was asking Alex for any advice how to improve FLDIGI CW decoder functionality. In his response Alex stated that he had worked on the problem for 8 years and tried hundreds of algorithms before he got the first version of the software that worked the way he wanted.

 Alex said: "The best advice I can give you is to try to approach every problem in the Bayesian framework. Express all your prior knowledge in the form of probabilities, and use observed data to update those probabilities. For example, instead of making a hard decision at every input sample whether the signal is present or not, compute the probability that the signal is present. At a later stage you will combine that probability with other probabilities using the Bayes formula and pass the new probabilities to the subsequent stages, all the way to the word recognition unit."

Since this was written by a man who really knows how to make great Morse decoder software I started digging deeper in what this Bayesian Framework is all about. I needed to freshen up my memories about probability theory as it has been some 30+ years since I studied this at the university. Perhaps the most clear and concrete examples on conditional probabilities and how to use Bayesian Rule I found from this University of Michigan Open Textbook.

Quick Bayesian Primer 

Before applying Bayesian Rule for FLDIGI Morse decoder we need to define some terms first. In the context of Morse code we are dealing with temporal elements like 'dits' and 'dahs' and time between them - these are the short and long audio pulses that Morse code characters consists of when we listen CW on ham bands. Morse decoder would have to break signals into these basic elements to figure out what character was just received. The elements and timing rules are well explained in this Wikipedia article.

Bayes' theorem is named for Thomas Bayes (1701–1761), who first suggested using the theorem to update beliefs.  To illustrate  the Bayes Rule in ham radio context I created a simple picture below in figure 1.  I collected some 'dits' and 'dahs' as well as inter-element gaps in between  from noisy Morse code signals and used histogram to plot the distribution of the timing of these.  The horizontal axis is in milliseconds and vertical is number of occurrences captured in this session.  On the right side we see two peaks - at 60 milliseconds and at 180 milliseconds.  These corresponds to 'dit'  (60 ms) and  'dah' (180 ms). On the negative side I also plotted the inter-character gap  - peak is at 60 milliseconds as expected and number of occurrences is roughly double as expected.

Fig 1.  Bayes Rule used in Morse Code decoder

To really make use of the Bayesian Framework we need to collect these numbers to make the probability calculations on the elements. I took existing CW.CXX as basis and completely re-wrote the lower level state machine that does the heavy lifting in dealing with the above timing rules. My focus was to enable collection of these timing elements -- both signals and gaps in between. I also added some code to build statistics on these elements while the decoder is running.

I called the new state machine "edge_recorder()" and it is responsible for detecting KEY DOWN / KEY UP edges and capturing 'dit','dah' and inter-character gap durations into a timecode buffer. It also keeps statistics of these durations in a histogram for Bayesian analysis of the Morse code. GNU Scientific Library (GSL) has well documented libraries for histograms so I used GSL Histogram package.

Bayesian Classifier 

I was struggling when writing the classifier software. Baysian Rule P(dit|x) =  P(x|dit)*P(dit) / P(x) basically requires you to compute the probability of 'dit' given observed element "x" - this probability could then be passed to the next stage for character classification. The following diagram helped me to formulate how to approach this.  We get a continuous stream of numbers representing possible 'dits' and 'dahs'.  We need to calculate first a conditional probability of  P(x |dit)  - like the shaded intersection in figure 2.

Fig 2. Conditional Probability

We need to define what 'dit'  really means.   In my case I set the limits so that 'dit' duration can be anything between  0.5 *  DIT  and  2 * DIT   where DIT is the expected duration based on the speed of the received   Morse code.  For example  at 20 WPM  (words per minute)  DIT would correspond to T = 1200/WPM so 1200/20 = 60 milliseconds.   Therefore in my classifier any received element "x" between  30 msec to 120 msec would be classified as 'dit'.  To calculate the probability P(x|dit)  I need to just look up from histogram values from all the bins between 30 and 120 and divide their total by the sum of  bin values from 0 to 300.  That gives us the conditional probability of  P(x|dit).  Same logic applies for 'dahs'.  In my classifier any element "x" bigger or equal to 2 * DIT and smaller than 4 * DIT is considered a 'dah'.

Obviously if the Morse speed changes you need to collect a new histogram and re-calculate the probabilities.  We have still couple additional items to deal with before we can apply Bayes rule.  P(x) is the overall probability of  element "x"  across all occurrences.  I limited this duration between  2*DIT and 4*DIT to remove outliers and impact due to noise.  We can calculate P(x) from histogram like we did above.  The final missing piece is  P(dit).   Based on my limited observations  P(dit) and P(dah)  have roughly equal probabilies in normal (English) ham Morse traffic so I gave them value  0.5.    Now we can use the Bayes formula to calculate P(dit|x) =  P(x|dit)*P(dit) / P(x).  In the software I calculate probabilities P(dit|x) and P(dah|x)  and pass these values to the next stage which does character matching.

To illustrate the value of this approach please look at the following figure 3.  This is a plot of the demodulated noisy signal coming from FFT filtering.  It is hard even for humans to decide whether the characters below match with "I N E T"  or "I R E A"  or "Z N U" or something else.  Humans have amazing capability to build beliefs  and updating these beliefs based on observations.    We are trying to emulate this human behavior.  Bayesian rule is frequently used in machine learning systems.

Figure 3.  Noisy Morse code

Character Classifier 

Next step was to build a Morse character classifier that uses the probabilities coming from the previous stage as described above.  To simplify the problem I created a simple codebook table that encodes the 'dit' and 'dah' probabilities of each character.  Now I just need to multiply incoming number vector representing 'dit'/'dah' probabilities against each codebook entry and then find the best match.



//  26 ASCII 7bit letters
{'A', "A", {{1,0},{0,1},{0,0},{0,0},{0,0},{0,0},{0,0}}},
{'B', "B", {{0,1},{1,0},{1,0},{1,0},{0,0},{0,0},{0,0}}},
{'C', "C", {{0,1},{1,0},{0,1},{1,0},{0,0},{0,0},{0,0}}},
{'D', "D", {{0,1},{1,0},{1,0},{0,0},{0,0},{0,0},{0,0}}},
{'E', "E", {{1,0},{0,0},{0,0},{0,0},{0,0},{0,0},{0,0}}},
{'F', "F", {{1,0},{1,0},{0,1},{1,0},{0,0},{0,0},{0,0}}},
{'G', "G", {{0,1},{0,1},{1,0},{0,0},{0,0},{0,0},{0,0}}},
        .......


Above is a small section of the codebook. Each character has a set of two numbers  - (1,0) means probability of  dit is 100%  and dah is 0%.   Character "E" that has a single 'dit'  has this in the first column.  This way we can represent the characters as a sequence of probabilities  of various 'dit' / 'dah' patterns.

We could also apply the Bayes  Rule at this level if we collect additional data on character statistics. For  example if we have statistics on the usage frequency of each character we could use that as an additional parameter to multiply the calculated probability coming from the character matching algorithm. This would just require additional frequency histogram, some code to pre-populate this based on selected language and perhaps code to update histogram when we are using the FLDIGI program. I have not implemented these yet.

Word Classifier 

The last stage would be the word classifier.   As we calculate the probabilities of each character in the previous stage  we could send these probabilities upstream to the word classifier.   There are a commonly used phrases, acronyms and ham radio jargon that could be used as a corpus  to match against received characters.

For example if  we receive BALAN  where the possible 4th letter in order of probabilities are  A = 0.2345, U = 0.2340, V=0.1234,  R=0.10235   we could use Bayes Rule to figure out probability of  BALUN  vs. BALAN  and emit the most probably word.   Similarly we could use a database of contest ham stations  and try to match probabilities of received characters to most probably station call sign using Bayes Rule.   This is yet another area that would be nice to implement in the future.

Conclusions 

I want to thank  Alex VE3NEA  for opening a whole new avenue of learning about Mr. Bayes and his ideas.   Based on my initial results the Bayesian Framework provides new rich data on how to model and improve the performance of the FLDIGI Morse Decoder.

The new software I created produces a ton of useful data to analyze  what is really happening under extremely noisy  conditions and also helps to separate different problems like  signal processing & filtering, signal power detection, signal edge detection, classification of  'dits' and 'dahs'  and finally  also character and word classification.

I have a alpha quality version of  Bayesian classifier  coded for FLDIGI  v.3.21.64 and I will work with Dave W1HKJ and any volunteers  to use these ideas to improve FLDIGI Morse decoder even further.

You can leave comments here or  you can also find me from linuxham Yahoo mailing list or from http://eham.net


73
Mauri  AG1LE

Morse Decoder SNR vs CER Testing

FLDIGI Morse Decoder  CER vs SNR Testing

Summary 

I wanted to get more factual data on how different FLDIGI  CW options influence the character error rate (CER)  at different signal-to-noise ratios (SNR).  In this test I focused on comparing SOM decoder to legacy decoder option and Matched Filter to  FFT filter option. I also tested a few cases where I used ~ 2x wider FFT filter  (68 Hz  vs. 35 Hz bandwidth).  Otherwise I used the default options of FLDIGI.

Based on the test results below the following appears to be true:

• Adjusting FFT filter bandwidth has much bigger impact on CER than changing between legacy and SOM decoder. 
• Matched filter automatically sets the filter bandwidth to optimal for given CW speed. 
• SOM decoder gives about 5% better CER compared to legacy decoder @ -13dB SNR, but this advantage gets smaller as SNR decreases. 

The test results are summarized in the figure 1 below. With 3 kHz wide audio signals you can get near perfect copy (CER below 0.02) at  SNR  -10 dB  using  Matched filter or setting FFT filter at 34 Hz  corresponding to 20 WPM Morse speed.

Fig 1.  FLDIGI  Morse Code decoder CER vs. SNR

Test Cases and Setup

My goals for this this experiment were to

• Compare FLDIGI legacy  vs. SOM decoder CER at different SNR levels. 
• Compare  FLDIGI Matched Filter  vs.  FFT filter (bandwidth set to 35Hz)  CER  at different SNR levels. 
• Compare CER with  FFT filter bandwith  set to 68 Hz  vs.  35 Hz.

I used the following software versions to conduct this test:

• FLDIGI v. 3.21.64 
• RTTY Text Comparison
• WinMorse 2.0
• PathSim 1.0 
• Audacity 2.0.2
• Graph  v.4.2.2

First  I generated a test audio file with 20 WPM Morse code using WinMorse and converted it to 16 bit  8 Khz Mono Wav file using Audacity.  I played the audio file with PathSim  configured for AWGN S/N Test simulation.  Both AWGN Noise source and "Input Source" have 3 kHZ band pass filter after them as demonstrated in the figure below.

Fig 2. PathSim Configuration

Sound output from PathSim  (8 kHz 16 Bit Mono) was connected to  FLDIGI  input using Virtual Audio cable  and level was adjusted at 10 in Windows 7 volume mixer so that FLDIGI waterfall is not saturated.  Audio signal had 3 Khz bandwidth as verified on FLDIGI  waterfall display.

Audio file length was 11 min 27 seconds.  PathSim "RunTime" was adjusted to 11 minutes - this provided only 869  of 904 total character in the file.  The lack of last 27+ seconds creates a systematic error of 35 missing characters causing baseline error of CER 0.03871 -  this is already subtracted from results below.

I changed audio SNR  in PathSim in steps from -20 dB, -15 dB, -14dB, -13 dB,  -10 dB and 0 dB also did a reference test case without any injected noise. I compared FLDIGI decoded characters against sent characters using RTTY Text Comparison tool by Alex VE3NEA.  As FLDIGI sometimes emitted prosigns (like <AR> ) I edited these to a single  '*' character before entering to 'received' field in the RTTY Text Comparison tool. Note that this tool was designed for RTTY characters so BER numbers are not applicable for Morse code as RTTY has fixed bit length baudot code whereas Morse code bit length  varies by character. The character error rate (CER)  is still be applicable.

I plotted the following measured CER results against SNR values with Graph.

Test Results 

Reference case 

No noise    CER 0.0011     First character incorrectly coded 

Test cases 

SNR            CER      Test
  0 db          0.00001   SOM Decoder with FFT Filter @68 Hz

-10 dB         0.01219   Legacy decoder  with Matched Filter @35 Hz
-10 dB         0.01439   SOM decoder    with Matched Filter @35 Hz
-10 dB         0.01439   Legacy decoder with FFT Filter @35 Hz
-10 dB         0.00999   SOM Decoder with FFT Filter @35 Hz
-10 db          0.20799   SOM Decoder with FFT Filter @68 Hz

-13 dB         0.19469   Legacy decoder  with Matched Filter @35 Hz
-13 dB         0.18589   SOM decoder    with Matched Filter @35 Hz
-13 dB         0.16149   Legacy decoder with FFT Filter @35 Hz
-13 dB         0.15379   SOM Decoder with FFT Filter @35 Hz
-13 db          0.66479   SOM Decoder with FFT Filter @68 Hz


-14 dB         0.37499   Legacy decoder  with Matched Filter @35 Hz
-14 dB         0.35729   SOM decoder    with Matched Filter @35 Hz
-14 dB         0.30199   Legacy decoder with FFT Filter @35 Hz
-14 dB         0.30089   SOM Decoder with FFT Filter @35 Hz

-15 dB         0.50669   Legacy decoder  with Matched Filter @35 Hz
-15 dB         0.49779   SOM decoder    with Matched Filter @35 Hz
-15 dB         0.47899   Legacy decoder  with FFT Filter @35 Hz
-15 dB         0.46569   SOM Decoder   with FFT Filter @35 Hz
-15 db          0.71789   SOM Decoder with FFT Filter @68 Hz

-20 dB         0.74229   Legacy decoder  with Matched Filter @35 Hz
-20 dB         0.68479   SOM decoder    with Matched Filter @35 Hz
-20 dB         0.72569   Legacy decoder with FFT Filter @35 Hz
-20 dB         0.68029   SOM Decoder with FFT Filter @35 Hz
-20 db          0.75439   SOM Decoder with FFT Filter @68 Hz

Other Observations on FLDIGI Performance

At low SNR values (-20 dB) FLDIGI seemed to lose CW speed tracking, especially if Matched filter was used.  The RxWPM display showed received speed variation between 23 and 30 WPM  while sent signal was at steady 20 WPM.   This may explain some of the decoding errors.  In real life cases if you know the Morse speed it is probably better to set it manually when trying to decode very noisy signals.


73
Mauri  AG1LE

Real time Morse decoder - New Ideas

Below is a link to a video clip from CQ WW WPX ham radio contest. I am using CW Skimmer software to decode 100+ stations in this 24 KHz section of the 7 Mhz radio band.

On the left panel the waterfall display shows about the frequency spectrum on vertical axis - horizontal axis is real time.   Blue color is the background noise level,  yellow / orange shows the received Morse code signals. There is also quite a lot of  noise and interference as shown by yellow/green dots.

On the right panel the Morse decoder attempts to decode the signals in raw text mode. There is some 100+ decoder instances active on this demo video. Noise seems to generate a lot of "false" decodes - visible as lines with many "E" or "I" characters.  In Morse code  a single "dit"  is  letter "E"  so it is quite difficult to determine whether noise spike was a real signal or not.  It seems that CW Skimmer software handles this problem at higher level of the decoder chain with Bayesian algorithms.   You can see how the software corrects decoded characters once additional data becomes available.

 

As I was studying alternative and new ways to tackle the Morse decoding problem from noisy signals I stumbled across this sparse distributed models of sequence memory - see video recognition example by Dr. Rod Rinkus.

I wonder if  this technology could be used for Morse code recognition in noisy HF bands?  Sequence memory would probably help with random "dits"  problem explained above. Once the system learns typical Morse character sequences,  words and phrases that are used in ham radio communications it should be able to recognize real signals from noise.

The system Dr. Rinkus has developed has many advantages such as:

• Can learn sequences with a single trial 
• Use Sparse Distributed Representation for individual sequence items and sequences 
• Can recall individual sequences as well as recognize novel sequences 
• Constant  computational time complexity -  O(1)  sequence comparison,  O(1)  learning - may lend well in real time processing? 
• Recognition of noisy and time-warped sequences  

There is a video in  here where  Dr. Rinkus explains in detail how the system works. See also this paper.

73
Mauri  AG1LE

Simple Elecraft KX3 and PowerSDR configuration

I have played with HDSDR and Elecraft KX3 for a while.  While  HDSDR is an excellent piece of software I am more familiar with PowerSDR as I have used it over 2 years.  I wanted to see how KX3  would work with my Flex3000 / PowerSDR  setup.  This turned out to be a fairly simple configuration - in fact you can run multiple instances of PowerSDR connected to different radios simultaneously.  Note that I am using PowerSDR version v2.4.4 on Windows 7.

The first step was to connect the Elecraft KX3  RX I/Q  interface to a sound card in the computer.  A quick visit to Radio Shack was required as I did not have a 2.5 mm to 3.5 mm stereo Y-adapter (part 274-945).  Using a normal male-to-male  3.5 stereo audio cord I connected KX3 to my home brew computer  Line Interface  (Light Blue connector in ASUS P8Z68-V PRO/GEN3 motherboard as shown in the ASUS P8Z68 Manual ).  See fig 1. below  - don't use the pink microphone input jack.

Fig 1.  ASUS P8Z68-V  audio connectors

Next step was to configure PowerSDR  properly.   When you start PowerSDR  without turning Flex3000 on first it will pop up a window - note the button "Add Legacy Radios" below.

Fig 2. PowerSDR Radio Interfaces

You can add SoftRock 40 legacy radio - no need to add anything on Serial Number field.

Fig 3.  Adding Legacy Radios

You will get back to Fig 2. window  and select  "Use" button on the right. Main window of PowerSDR software should come up.  

Next step is to select  Setup menu from the top and General / Hardware Config.  You can setup the Center Frequency here - see Fig 4. below.

Fig 4.   Set Center Frequency

Next you configure the Audio / Primary settings.  In my case I am using the built-in sound card of the ASUS motherboard which is not on the supported sound cards list. There I selected "Unsupported Card" as shown in Fig 5. below.  Also,  I configured Sample Rate to 192,000 to get maximum frequency coverage from KX3  RX I/Q signals.  I also used "MME" Driver of Windows 7  and selected Input  and Mixer to "Line In High Definition Audio"   and selected Output to my video monitor "EQ276W DP-1 (NVIDIA...)".

Fig 5. Audio Hardware Configuration

Time to  press the "Start" button on  PowerSDR.  You should see some signals on panadapter if evrrything is configured correctly. 

Fig 6.  PowerSDR showing 96 Khz of signal coming from Elecraft KX3

Note that PowerSDR expects the center frequency be at 14.200 MHz   (see Fig 4. above) and
to get the frequency display correct you need to tune KX3  at that frequency.   You can also see that PowerSDR covers approximately  96 kHz  bandwidth centered around 14.2 MHz.  Of course you can change that by just going to Setup menu.  

Fig 7. Elecraft KX3  tuned at 14.200 Mhz

I was running Flex3000 and KX3   in parallel listening some stations while switching the antenna between KX3 and Flex3000.   I did not do any scientific A/B comparison study  but just by  using two instances of the same PowerSDR version connected to different radios I could not really tell much difference between these two radios.  Both picked up the faint  DX stations equally well.  The only difference was really the few kHz noise band around the KX3 center frequency  - outside of that the sensitivity and sound quality was very similar in both of these radios. 

I was also playing with the CAT interface of PowerSDR  - however, it does not recognize Elecraft as an option so I did not spend much time on that.  In comparison HDSDR uses Omnirig  CAT interface which has more choices, including Elecraft K3  that works well with KX3.    If you can figure out a way to get KX3  CAT interface working also with PowerSDR please let me know.

In conclusion I wanted to test if I can make PowerSDR software to use Elecraft KX3  as a receiver.  As you can see above this was a quick and painless configuration effort.


73  
Mauri   AG1LE

Raspberry Pi: APRS tracker mash-up

I wanted to test how Raspberry Pi would work as a small web server  so I installed NGINX  and configured it using the instructions in here. 

I also created a small Javascript page that uses  aprs.fi  JSON API  and Nokia Maps API to track  APRS enabled station on a map.  The source code is posted in github  - this is basically just one file 'index.html'  that is copied to /var/www directory on Raspberry PI.   You need to get the APRS.FI API key from  here.

Here is how the page looks like as served by Raspberry Pi  - see W1MRA-C station  near Prudential Center in Boston.  As you can see from the source code this is a very simplistic Javascript mash-up. 


However, Nokia Maps APIs  do provide a rich set of functionality that can be used to extend the capabilities.  The neat thing here  is that even a tiny web server like Raspberry Pi  can provide this rich functionality because all the heavy lifting is done by Nokia on the maps API side and on the browser doing the page rendering.  

Raspberry Pi Arrived!

I placed the Raspberry Pi order to RS Components  in early March and this credit card sized  $35  Linux computer finally arrived on Monday.  

I created a new page for Raspberry Pi related projects and experiments. Please check in here.  

Morse code detection using modified Morlet Wavelet transformation

In my previous blog post I shared some experiments with wavelets using available online tools.  In order to build better understanding on how to apply Morlet wavelets in detecting Morse code in noisy signals I wrote a little test application using Octave.   The application creates noisy Morse code and does Continuous Wavelet Transform (CWT) using modified Morlet  wavelet as well as Short Term Fourier Transform (STFT).  The Octave application allows changing various parameters and visualizing the impact on the plotted graphs and images.  See discussion on modification to Morlet wavelet below. 

The starting point was original noisy Morse code signal  that has -12.2 dB SNR as shown in the figure 1 below.  It is virtually impossible to detect the Morse signal buried in the noise. When listening the audio I can hear a faint sound in the noise but I have much difficulties recognizing any of the characters.
  
For reference, the signal-to-noise ratio calculation is done the following way: 

SNR=20*log10(norm(morsecode)/norm(morsecode-noisey_morsecode));

Figure 1.   Morse code with -12.2 dB SNR (signal-to-noise ratio)

To discover the Morse code signal I used Short Time Fourier Transformation (STFT)  with the following parameters ( x contains noisy audio,  sampling rate is Fs=4000 )

z = stft(x,1/Fs,4,128,1,512);

A very faint Morse signal is now visible at 600 Hz frequency as a horizontal pattern of "dits and "dahs" in figure 2. below.  It is very difficult to "read" from this figure 3 below what the message is.  Some "dits" and "dahs" are more visible but noise makes it difficult to detect the pattern or to decode the message. 

Figure 2.  Waterfall spectrogram of the Morse code with -12.2dB SNR.

Modified Morlet Wavelet power spectrum is shown on figure 3.  On Scale axis (vertical) you can see  at S = 8 as a horizontal pattern of lighter "dits" and "dahs".  Looking at the pattern you can almost see "dah-dit-dah-dit"  "dah dah dit dah"  "dah dit dit" "dit"  "dit dah" "dah dah dit" "dit dah dah dah dah" "dit dah dit dit"  "dit"   aka  "CQ DE AG1LE".  This represents the peak energy of the power spectrum  after Wavelet transformation. I am taking absolute value to show the envelope of the wavelet better  (see morlet.m file below, calculation is done by this line:    coeffs(k,:)=abs(fftshift(ifft(fft(w).*fft(sigin)))); )

Figure 3.  Wavelet Power Spectrum of the Morse code with -12.2 dB SNR

I plotted Wavelet coefficient C(t,8) values  (corresponding Scale = 8 above) on  x-y graph below in figure 4.  Morse code "dits" and "dahs"  are quite visible as signal peaks above threshold value of 2000.  Note that these are  abs(C(t,8))  values and low pass filtered to show envelope better. See morletdemo.m below for details.

Figure 4.  Modified Morlet coefficient C(8) values.

Original Morlet  wavelets (scales 1 to 16 shown below in figure 5.) have variable wavelet pulse length and frequency as follows: 



for k=1:scale,
    t=(-M/2:M/2-1);
   
    % Calculate Morelet Wavelet w=e-(at^2)*cos(2*pi*fo*t)
    const = 1/(sigma*K*sqrt(k)); % k impacts relative amplitude
    e = exp(-((sigma*t/k).^2));    % k impacts pulse length (t/k)
    phase = cos(2*pi*fo*t/k);    % k impacts frequency 
    w = const*e.*phase;
    plot(w) 
end

While experimenting with these wavelets it was quite difficult find the optimal wavelet to extract signal from noise. 

Figure 5.  Morlet wavelets 1...16

Looking at the impact of various parameters it became obvious that by modifying the wavelet to keep the duration constant improves the situation a lot.  I modified a single line and the corresponding wavelet graph is below in figure 6.   

  e = exp(-((sigma*t).^2));    % removed k to keep the wavelet duration constant.

Note that wavelet bandwidth sigma depends on Morse speed - I did several experiments and established the following relationship  sigma  = (1.2/speed)/w   where w is the number of the wavelet.

Figure 6.   Morlet Wavelet 

CONCLUSIONS 

With a small modification to Morlet wavelet the CWT works better than STFT in extracting the signal from the noise even at -12 dB SNR.  There are many similarities to Matched Filter method that I described in this blog post. Perhaps the main difference is the selected wavelet shape (Morlet) and the fact we use FFT to make the convolution very fast.  

Further work could include the following tasks:

• program  the modified  Morlet Wavelet algorithm in C++ 
• implement this functionality in FLDIGI  CW decoder module
• test Wavelet based decoder with real life signals from HF bands

SOFTWARE 

The Octave scripts are listed below.  The above results were created by running the following command on Octave: 

%  noise level,  signal freq, sampling rate, morsespeed
morletdemo(2,600,4000,20)

The software prints these lines

text = CQ de AG1LE
file = CQ.wav
SNR =-12.226566
Fo =  1.2000
sigma =  0.0075000

and plots the figures.


File Morlet.m


% Project 2
% Time-Frequency Representations
% Andy Doran, modified by  AG1LE Mauri Niininen 
% 
% function coeffs = cwvt(sigin,scale,quiet,sigma,fo)
%
% sigin  = sampled input signal (Should be a row vector)
% scale  = number of real, positive scale you want (1:scales)
% quiet  = plot suppression
%          1 -> suppress all plots
%          2 -> suppress wavelet plots only
%          3 -> suppress scalogram only
% sigma =0.015625;     Morlet Wavelet bandwidth
% fo = 0.25;           Center frequency of Wavelet
% coeffs = scales-by-length(sigin) matrix returning CWT of sigin at
%          each scale from 1 to scale
%
% This function takes an input signal and computes the Continuous Wavelet
% Transform at different scales using a sampled Morlet Wavelet
%
% Morelet Wavelet w(t) = (1/sigma*K)*exp-((sigma*t)^2)*cos(2*pi*fo*t)


function coeffs = morlet(sigin,scale,quiet,sigma,fo)




K = 1;           % Not sure what this is, so set to 1


M = length(sigin);
coeffs = zeros(scale,M);


for k=1:scale,
    t=(-M/2:M/2-1);


    % Calculate Morelet Wavelet w=e-(at^2)*cos(2*pi*fo*t)
    const = 1/(sigma*K*sqrt(k)); % k impacts relative amplitude
    e = exp(-((sigma*t).^2));    % removed k to keep the wavelet duration constant
    phase = cos(2*pi*fo*t/k);    % k impacts frequency 
    w = const*e.*phase;


    % Plot wavelet in time domain and frequency domain
    if ((quiet ~= 1) & (quiet ~= 2))
      figure(3)
      if (k == 1)  % Clear plot on initial run-through
        clf
      end
      subplot(scale,2,(2*k)-1)
      plot(w)
      txt = ['Modified Morlet Wavelet at scale ', num2str(k)];
      title(txt)
      %figure(4)
      %if (k == 1)  % Clear plot on initial run-through
      %  clf
      %end
      subplot(scale,2,2*k)
      plot(abs(fft(w)))
      txt = ['Frequency Spectra of Morlet Wavelet at scale ', num2str(k)];
      title(txt)
   end
    % Calculate CWT of sigin using circular convolution
%    coeffs(k,:)=ifft(fft(w).*fft(sigin));
    coeffs(k,:)=abs(fftshift(ifft(fft(w).*fft(sigin))));
end


% Coeffs should be real anyway, this just accounts for numerical error
% in circular convolution that may give small imaginary parts
coeffs = real(coeffs);


% Plot scalogram and check against MATLAB's CWT routine
if ((quiet ~= 1) & (quiet ~= 3))
  figure(1)
  %clf
  map = jet();
  colormap(map);
  imagesc(coeffs);
  axis xy;
  txt = ['abs|C(t,s)| for s = 1 to ' num2str(scale)];
  title(txt)
  ylabel('s')
  xlabel('t')
  figure(2);
  plot(sigin);  
  title('original signal');
 %figure(4);
  %clf
  %cwt(sigin,1:scale,'morl','plot');  % Call MATLAB's CWT routine
  %title('CWT Output from MATLAB')
end


File  stft.m


function y = STFT(x, sampling_rate, window, window_length, step_dist, padding)
%
%  y = STFT(x, sampling_rate, window, window_length, step_dist, padding)
%
%  STFT produces a TF image of "x".
%  The output is also stored in "y".
%
%  For "window", use one of the following inputs:
%  rectangular    = 1
%  Hamming        = 2
%  Hanning        = 3
%  Blackman-Tukey = 4
%
%  The time scale is associated with the center of the window,
%  if the window is of odd length.  Otherwise, the window_length/2
%  is used.  "Step_dist" determines the stepping distance between the number
%  of samples, and is arranged to maintain the proper time index
%  provided by "sampling_rate" in seconds.  "Padding" is the
%  total length of the windowed signal before the fft, which is
%  accomplished by zero padding.
%
%  Developed by Timothy D. Dorney
%               Rice University
%               April, 1999
%               tdorney@ieee.org
%
%  Coded using MATLAB 5.X.X.
%  See http://www.clear.rice.edu/elec631/Projects99/mit/index2.htm
%
% REVISION HISTORY
%
% VERSION 1.0.0 APR. 21, 1999 TIM DORNEY
%


if (nargin ~= 6)
        disp('STFT requires 6 input arguments!')
return;
end
if ((window < 1) | (window > 4))
window = 1;
disp('The argument "window" must be between 1-4, inclusively.  Window set to 1!');
end
if ((step_dist < 1) | (round(step_dist) ~= step_dist))
step_dist = 1;
disp('The argument "step_dist" must be an integer greater than 0.  Step_dist set to 1!');
end
if (sampling_rate <= 0)
disp('The argument "sampling_rate" must be greater than 0.');
break;
end
if (padding < window_length)
padding = window_length;
disp('The argument "padding" must be non-negative.  Padding set to "window_length"!');
end


if (window == 1)
WIN = ones(1,window_length);
elseif (window == 2)
WIN = hamming(window_length)';
elseif (window == 3)
WIN = hanning(window_length)';
elseif (window == 4)
WIN = blackman(window_length)';
end


[m,n] = size(x);
if (m ~= 1)
X = x';
else
X = x;
end
[m,n] = size(X);
if (m ~= 1)
disp('X must be a vector, not a matrix!');
break;
end


LENX = length(X);
IMGX = ceil(LENX/step_dist);
if (padding/2 == round(padding/2))
IMGY = (padding/2) + 1;
else
IMGY = ceil(padding/2);
end


y = zeros(IMGX,IMGY);


if (window_length/2 == round(window_length/2))
CENTER = window_length/2;
x_pad_st = window_length - CENTER - 1;
x_pad_fi = window_length - CENTER;
else
CENTER = (window_length+1)/2;
x_pad_st = window_length - CENTER;
x_pad_fi = window_length - CENTER;
end


X = [zeros(1,x_pad_st) X zeros(1,x_pad_fi)];


iter = 0;
for kk = 1:step_dist:LENX
iter = iter + 1;
XX = X(kk:(kk + window_length - 1));
YY = XX .* WIN;
ZZ = abs(fft(YY, padding));
y(iter,:) = ZZ(1:IMGY);
end
figure(6);
freq = (1/sampling_rate)/2;
imagesc([0:(step_dist*sampling_rate):(sampling_rate*(LENX-1))], ...
[0:(freq/(IMGY-1)):freq],y');
xlabel('Time (seconds)');
ylabel('Frequency (Hz)');
axis('xy')

File  morletdemo.m


function morletdemo(noisy,freq,Fs,speed);


x = morse('CQ de AG1LE','CQ.wav',noisy,freq,Fs,speed);


w = 8;         %  peak will be at wavelet # w
Fo = freq / (Fs/w)     % tell wavelet transform where wavelet center frequency is 
sigma  = (1.2/speed)/w  % wavelet bandwidth - impacts time resolution


c = morlet(x',16,2,sigma,Fo);         % do Morlet wavelet transform
y = filter(ones(1,20)/20,1,c(w,:));   % y = low pass filter C(t,w) wavelet 
figure(4)
plot(y);                              % plot C(t,w) envelope


z = stft(x,1/Fs,4,128,1,512); % plot spectrogram of the signal using Short Term FFT 
end;




File  morse.m

function code=morse(varargin)

% MORSE converts text to playable morse code in wav format

%

% SYNTAX

% morse(text)

% morse(text,file_name);

% morse(text,file_name,noise_multiplier);

% morse(text, file_name,noise_multiplier,code_frequency);

% morse(text, file_name,noise_multiplier,code_frequency,sample_rate);

% morse(text, file_name,noise_multiplier,code_frequency,sample_rate, code_speed_wpm, zero_fill_to_N);

% morse(text, file_name,noise_multiplier,code_frequency,sample_rate, code_speed_wpm, zero_fill_to_N, play_sound);

%

% Description:

%

%   If the wave file name is specified, then the funtion will output a wav

%   file with that file name.  If only text is specified, then the function

%   will only play the morse code wav file without saving it to a wav file.

%   If a noise multiplier is specified, zero mean addative white Gaussian

%   noise is added with 'amplitude' noise_multiplier.

%

% Examples:

%

%   morse('Hello');

%   morse('How are you doing my friend?','morsecode.wav');

%   morse('How are you doing my friend?','morsecode.wav', 0.01);

%   morse('How are you doing my friend?','morsecode.wav', 0.01, 440, ,20, Fs);

%   x = morse('How are you doing my friend?','morsecode.wav', 0.01, 440, 20, Fs, 2^20,1); %(to play the file, and make the length 2^20)

%

%   Copyright 2005 Fahad Al Mahmood

%   Version: 1.1 $  $Date: 08-Jul-2010

%   Modifications: Rob Frohne, KL7NA

%Defualt values

Fs=48000;

noise_multiplier = 0;

f_code = 375;

code_speed = 20;

text = varargin{1}

if nargin>=2

file = varargin{2}

end

if nargin>=3

noise_multiplier = varargin{3};

end

if nargin>=4

f_code = varargin{4};

end

if nargin>=5

Fs = varargin{5};

end

if nargin>=6

code_speed = varargin{6};

end

if nargin>=7

length_N = varargin{7};

end

if nargin>=8

playsound = varargin{8};

end

t=0:1/Fs:1.2/code_speed; %One dit of time at w wpm is 1.2/w.

t=t';

Dit = sin(2*pi*f_code*t);

ssp = zeros(size(Dit));

#Dah fixed by Zach Swena 

t2=0:1/Fs:3*1.2/code_speed;  # one Dah of time is 3 times  dit time

t2=t2';

Dah = sin(2*pi*f_code*t2);

lsp = zeros(size(Dah));    # changed size argument to function of Dah 

#Dah = [Dit;Dit;Dit];

#lsp = zeros(size([Dit;Dit;Dit]));

% Defining Characters & Numbers

A = [Dit;ssp;Dah];

B = [Dah;ssp;Dit;ssp;Dit;ssp;Dit];

C = [Dah;ssp;Dit;ssp;Dah;ssp;Dit];

D = [Dah;ssp;Dit;ssp;Dit];

E = [Dit];

F = [Dit;ssp;Dit;ssp;Dah;ssp;Dit];

G = [Dah;ssp;Dah;ssp;Dit];

H = [Dit;ssp;Dit;ssp;Dit;ssp;Dit];

I = [Dit;ssp;Dit];

J = [Dit;ssp;Dah;ssp;Dah;ssp;Dah];

K = [Dah;ssp;Dit;ssp;Dah];

L = [Dit;ssp;Dah;ssp;Dit;ssp;Dit];

M = [Dah;ssp;Dah];

N = [Dah;ssp;Dit];

O = [Dah;ssp;Dah;ssp;Dah];

P = [Dit;ssp;Dah;ssp;Dah;ssp;Dit];

Q = [Dah;ssp;Dah;ssp;Dit;ssp;Dah];

R = [Dit;ssp;Dah;ssp;Dit];

S = [Dit;ssp;Dit;ssp;Dit];

T = [Dah];

U = [Dit;ssp;Dit;ssp;Dah];

V = [Dit;ssp;Dit;ssp;Dit;ssp;Dah];

W = [Dit;ssp;Dah;ssp;Dah];

X = [Dah;ssp;Dit;ssp;Dit;ssp;Dah];

Y = [Dah;ssp;Dit;ssp;Dah;ssp;Dah];

Z = [Dah;ssp;Dah;ssp;Dit;ssp;Dit];

period = [Dit;ssp;Dah;ssp;Dit;ssp;Dah;ssp;Dit;ssp;Dah];

comma = [Dah;ssp;Dah;ssp;Dit;ssp;Dit;ssp;Dah;ssp;Dah];

question = [Dit;ssp;Dit;ssp;Dah;ssp;Dah;ssp;Dit;ssp;Dit];

slash_ = [Dah;ssp;Dit;ssp;Dit;ssp;Dah;ssp;Dit];

n1 = [Dit;ssp;Dah;ssp;Dah;ssp;Dah;ssp;Dah];

n2 = [Dit;ssp;Dit;ssp;Dah;ssp;Dah;ssp;Dah];

n3 = [Dit;ssp;Dit;ssp;Dit;ssp;Dah;ssp;Dah];

n4 = [Dit;ssp;Dit;ssp;Dit;ssp;Dit;ssp;Dah];

n5 = [Dit;ssp;Dit;ssp;Dit;ssp;Dit;ssp;Dit];

n6 = [Dah;ssp;Dit;ssp;Dit;ssp;Dit;ssp;Dit];

n7 = [Dah;ssp;Dah;ssp;Dit;ssp;Dit;ssp;Dit];

n8 = [Dah;ssp;Dah;ssp;Dah;ssp;Dit;ssp;Dit];

n9 = [Dah;ssp;Dah;ssp;Dah;ssp;Dah;ssp;Dit];

n0 = [Dah;ssp;Dah;ssp;Dah;ssp;Dah;ssp;Dah];

text = upper(text);

vars ={'period','comma','question','slash_'};

morsecode=[];

for i=1:length(text)

if isvarname(text(i))

morsecode = [morsecode;eval(text(i))];

elseif ismember(text(i),'.,?/')

x = findstr(text(i),'.,?/');

morsecode = [morsecode;eval(vars{x})];

elseif ~isempty(str2num(text(i)))

morsecode = [morsecode;eval(['n' text(i)])];

elseif text(i)==' '

morsecode = [morsecode;ssp;ssp;ssp;ssp];

end

morsecode = [morsecode;lsp];

end

if exist('length_N','var')

append_length = length_N - length(morsecode);

if (append_length < 0)

printf("Length %d isn't large enough for your message; it must be > %d.\n",length_N,length(morsecode));

return;

else

morsecode = [morsecode; zeros(append_length,1)];

end

end

noisey_morsecode = morsecode + noise_multiplier*randn(size(morsecode));

SNR=20*log10(norm(morsecode)/norm(morsecode-noisey_morsecode));

printf('SNR =%f\n',SNR);

if exist('file','var')

wavwrite(noisey_morsecode,Fs,16,file);

if exist('playsound')

system(['aplay ',file]);

end

else

soundsc(noisey_morsecode,Fs);

% wavplay(morsecode);

end

code = noisey_morsecode;

endfunction