FLDIGI: more test results on experimental features

In my previous blog post I shared two experimental features that I had developed on FLDIGI.
Matched Filter was to improve signal-to-noise ratio and an algorithm borrowed from Self Organizing Map (SOM) neural network to improve morse decoding.


With help of Rene I posted experimental FLDIGI software in here
FLDIGI Experimental features

Dave W1HKJ was kind enough to provide a set of test files.  I copied the wav files at http://www.w1hkj.com/developers/
to my local hard disk.

I opened AlsaPlayer and added the wav files in the playlist as shown in Figure 1.


On my system the software version crashes after 2-3 memory replay events when experimental features are enabled. Therefore I had to take the tests with only two files at the time when experimental features were enabled.  I have not yet found the reason for these Sig11 events. I used memory replay feature in the tests below. 

Here are the tests I performed: 

1)  Legacy FLDIGI  detection and decoder - all files

I disabled  matched filter  (check box "off") and disabled SOM detector ("Use Farnsworth timing"  "off").  I started AlsaPlayer and let it cycle through all files in order shown below.  FLDIGI morse detection result is shown on Figure 1.  Errors are reduced from 1db onwards. Signals between -3dB to 0dB produce decoding errors.

Figure 1.  Legacy FLDIGI detection and decoder - test signals -3dB....6dB

2)  Using Matched Filter and SOM detection for -3dB and -2DB signals

For  files cw.-3db.wav and cw.-2db.wav  I enabled Matched Filter and SOM detection. I let the whole 36 seconds run and clicked waterfall to cause memory replay event for  -3dB signal. I did the same for -2dB signal.  The result is shown on Figure 2. below. 

Figure 2.  -3dB and -2dB detection with Matched Filter and SOM

3)  Using Matched Filter and SOM detection for -1dB and 0DB signals

 For  files cw.-1db.wav and cw.0db.wav  I enabled Matched Filter and SOM detection. I let the whole 36 seconds run and clicked waterfall to cause memory replay event for  -1dB signal. I did the same for 0dB signal.  The result is shown on Figure 3. below.

Figure 3. -1dB and 0dB detection with Matched Filter and SOM

4)  Using Matched Filter and SOM detection for 0dB and 1DB signals

For  files cw.0db.wav and cw.1db.wav  I enabled Matched Filter and SOM detection. I let the whole 36 seconds run and clicked waterfall to cause memory replay event for  0 dB signal. I did the same for 1dB signal.  The result is shown on Figure 4. below.

Figure 4. 0dB and 1 dB detection using Matched Filter and SOM

5)  Using Matched Filter and SOM detection for 1dB and 2DB signals

For  files cw.1db.wav and cw.2db.wav  I enabled Matched Filter and SOM detection. I let the whole 36 seconds run and clicked waterfall to cause memory replay event for  1 dB signal. I did the same for 2dB signal.  The result is shown on Figure 5. below.

Figure 5.  1dB and 2dB detection with Matched Filter and SOM

6)  Using Matched Filter and SOM detection for 3dB and 6DB signals

For  files cw.3db.wav and cw.6db.wav  I enabled Matched Filter and SOM detection. I let the whole 36 seconds run and clicked waterfall to cause memory replay event for  3 dB signal. I did the same for 6dB signal.  The result is shown on Figure 6. below.

Figure 6.  3dB and 6 dB detection with Matched Filter and SOM

7)  Legacy FLDIGI detection  for  -3db and -2dB 

In this test I focused on -3dB and -2dB signals.  Figure 7 below was taken using legacy FLDIGI decoding.  One CQ and  parts of call sign  HKJ  was detected. 

Figure 7. FLDIGI legacy decoder for -3dB and -2dB signals

8)  FLDIGI SOM detection  for  -3db and -2dB 

In this test I focused on -3dB and -2dB signals.  Figure 7 below was taken using SOM  decoding (Matched Filter is off).  One CQ and  parts of call sign  HKJ  was detected.   SOM algorithm has returned best matching characters - but most of these are garbage due to noise.

Figure 8.  FLDIGI SOM decoding for -3dB and -2dB signals

9)  FLDIGI SOM & Matched Filter and FLDIGI legacy detection comparison for  -3db  signal

In this test I focused on -3dB signal.  Figure 9 below was taken using SOM  decoding (Matched Filter is on) first 2 rounds.  Memory replay event was triggered by clicking waterfall. Then SOM and Matched Filter were turned off and 2 rounds of legacy decoder was run.   SOM decoder shows in green text, legacy decoder in black text below.

Figure 9.  FLDIGI  SOM and legacy decoder for -3db signal

10) CONCLUSIONS

• Matched filter works best when a long buffer is processed. Using memory replay events (clicking Waterfall) you can get the whole 36 second long signal run through the Matched Filter, like in Figure 9.  
• Matched filter is sensitive to correct frequency. Using 4x zoom feature on waterfall allows you to center the bandwith on exactly right frequency.
• Matched Filter is sensitive to correct Morse speed.  With manual adjustment to 20 WPM you can reduce errors compared to default 18 WPM setting. Automatic speed adjustment may not work correctly during memory replay events - how to verify this? 
•  SOM decoder calculates Euclidian minimum distance of the detected lengths of "dits" and "dahs" durations compared to codebook values. 
• SOM decoder tries to find best matching characters from the codebook. Noise can cause incorrect coding  - a "dit" is interpreted as "dah"  or vice versa, or some "dits" or "dah" is missing.  See Figure 9.  for examples:
  W1DKR is miscoded  -  D  (-..)  instead if H(....)  and R (.-.) instead of J (.---).   WLUAJ   is miscoded -  L(.-..) instead of 1 (.----)  and U(..-)  instead of H(....) and A(.-) instead of  K(.-.)  
  RQ (.-. --.-) instead of CQ (-.-.  --.-)

FLDIGI: Adding matched filter feature to CW mode

In my previous blog post I was experimenting with some algorithms to detect and decode morse code from noisy RF band.   Since the results looked promising I have spent some time implementing matched filter algorithm in  FLDIGI.  I also implemented a neural network decoder based on Self Organized Maps (SOM) algorithms.

I started by downloading fldigi-3.21.41 source code from W1HKJ download page.  Once I got the software to compile properly on my Linux system I studied the internals how CW detection & decoding really works in FLDIGI.

1) FLDIGI  CW MODULE INTRODUCTION

Main modules are in CW.CXX  and MORSE.CXX.  After initialization the  cw:rx_process() function receives new audio data 512 samples on continuous basis (this block size is actually configured in sound.h  - see #define SCBLOCKSIZE  ).

Once a block of data is received, there are some checks done if user has changed the low pass filter bandwidth. Then the  algorithm creates a baseband signal by mixing with audio carrier frequency of selected signal. The next step is to run a low pass filter,  demodulate by taking the magnitude of the signal and finally running the envelope through a moving average filter.  There is also AGC (automatic gain control) function with fast attack and slow decay. A check is done if squelch is on and if signal exceeds squelch value (this is to reduce noise created errors).  If the signal has upward trend the control gets passed to handle_event() function and same is done with downward trend.  Handle_event() keeps track of  "dits" and "dahs"  and eventually returns decoded character using morse::rx_lookup() function.

There is actually much more going on than the above simple explanation. The software keeps track of morse speed, automatic gain control,  dit/dah ratio and various other parameters.  Also, the user interface is updated - if you have signal scope on the signal amplitude value is updated between characters etc.

2) EXPERIMENTAL FEATURES ADDED 

Matched Filter 

I implemented Matched Filter algorithm using a simple convolution method. The implementation is in MFILT.CXX and it has few key functions.

Mfilt::create_filter() calculates "dit" time based on current morse speed. It also allocates memory for 3 buffers for the convolution. Finally, it creates a convolution kernel - this is a "dit" long sine wave burst at selected CW frequency.

Mfilt::convolve() performs linear convolution with very straight forward calculation. It has not been speed optimized but seems to be fast enough on my desktop computer.

Mfilt::run()  collects input samples until selected input buffer has been filled. Once enough samples are collected it will call the convolve() function.

Compared to other FLDIGI filters this is a very straight forward filter and no speed optimization has been done. I am sure that the gurus who wrote FLDIGI filters could improve this a lot.

SOM (Self Organizing Map) decoding

Last year I spent some time trying to understand how to use Self Organizing Map algorithm to decode noisy signals.  I collected quite a lot of data from my Flexradio 3000 trying to cover multiple signals from different bands (40m, 20m, 17m) and from different stations with varying signal strengths and keying styles.  The idea behind this work was to figure out if SOM algorithm would learn unsupervised how to decode morse code when presented a variety of signals as learning data.

Self Organizing Maps are different from other artificial neural networks as they use a neighborhood function to preserve topological properties of the input space. Like most artificial neural networks, SOMs operate in two modes: training and mapping. Training builds the map using input examples. It is a competitive process, also called vector quantization. Mapping automatically classifies a new input vector.

I used the SOM Toolbox 2.0 from my alma mater.  This excellent software package written for Matlab 5 did also work with Octave after some tweaking.
I basically created one large morse learning dataset and then run SOM training algorithm to create 20 x 20, 10 x 10 and 7 x 7 rectangular neuron grids with morse characters clustered based on their topological properties.  It was very interesting to see how morse code characters with similar properties ("E" = "dit" and "I" = "dit dit") automatically converged to nearby neuron cells.

I was first planning to implement this type of self learning algorithm as part of FLDIGI but realized soon that coding the entire SOM toolbox functions with C++ would be quite a large project.

After pondering this problem for a while I realized that the essence of the problem was perhaps not the learning part but using the "Best Matching Unit" algorithm against a known codebook. After all, morse code is well defined and every character should match to one codebook entry.  The BMU algorithm  basically iterates through all the nodes, calculates Euclidian distance between each node weigth vector and current input vector. The node with weight vector closest to input vector is tagged as the "best matching unit".

The Euclidian distance is given as:

where V is the current input vector and W is the weight vector.

I created the following  "find_winner()" function. It takes a buffer with detected "dit" and "dah" lengths and the duration of two dots as input parameters. It then computes the distance between codebook and input vector. The winner is the closest match with smallest distance. Since taking square root is computationally expensive function we simplify by calculating Manhattan distance that works as well here:


const char *find_winner (float *inbuf, int twodots)
{
    int i;
    SOM_TABLE *som, *winner;
    float diffsf = 999999999999.0;
    float difference = 0.0;
    float diff;

    if ( normalize (inbuf, 7, twodots) == 0) return NULL;

    winner = NULL;
    for (som = som_table; som->chr != 0; som++) {
         /* Compute the distance between codebook and input entry */
        difference = 0.0;
           for (i = 0; i < 7; i++) {
            diff = (inbuf[i] - som->wgt[i]);
                    difference += diff * diff;
                    if (difference > diffsf) break;
              }

     /* If distance is smaller than previous distances */
            if (difference < diffsf) {
                  winner = som;
                  diffsf = difference;
            }
    }

    if (winner != NULL)
        return winner->prt;
    else
        return NULL;

}


In order to normalize the durations of "dits" and "dahs"  I used the following function:

int normalize(float *v, int n, int twodots)
{
    float max = -99999999999999.9;
    float min = 999999999999999.9;
    int j;
    /* find max and min values */   
    for (j=0; j<n; j++) {
            if (v[j] > max)    max = v[j];
            if (v[j] < min)    min = v[j];
    }
    /* all values 0 - no need to normalize or decode */
    if (max == 0.0) return 0;

    /* scale values between  [0,1] -- if Max longer than 2 dots it was "dah" and should be 1.0, otherwise it was "dit" and should be 0.33 */
    if (max > twodots)
        for (j=0; j<n; j++) v[j] = v[j] / max;
    else
        for (j=0; j<n; j++) v[j] = 0.33 * v[j]/ max;

    return (1);
}

Finally,  I created a codebook that has normalized durations - "dah" is 1.0  and "dit is  0.33  following the standard morse code conventions.

struct SOM_TABLE {
    char chr;    /* The character(s) represented */
    const char *prt;    /* The printable representation of the character */
    float wgt[7];    /* Dot-dash weight vector */
};


static SOM_TABLE som_table[] = {
    /* Prosigns */
    {'=',    "<BT>",   {1.0,  0.33,  0.33,  0.33, 1.0,   0, 0}    }, // 0
    {'~',    "<AA>",   { 0.33, 1.0,  0.33, 1.0,   0,   0, 0}    }, // 1
    {'%',    "<AS>",   { 0.33, 1.0,  0.33,  0.33,  0.33,   0, 0}     }, // 2
    {'+',    "<AR>",   { 0.33, 1.0,  0.33, 1.0,  0.33,   0, 0}     }, // 3
    {'>',    "<SK>",   { 0.33,  0.33,  0.33, 1.0,  0.33, 1.0, 0}    }, // 4
    {'<',    "<KN>",   {1.0,  0.33, 1.0, 1.0,  0.33,   0, 0}     }, // 5
    {'&',    "<INT>",  { 0.33,  0.33, 1.0,  0.33, 1.0,   0, 0}    }, // 6
    {'}',    "<HM>",   { 0.33,  0.33,  0.33,  0.33, 1.0, 1.0, 0}    }, // 7
    {'{',    "<VE>",   { 0.33,  0.33,  0.33, 1.0,  0.33,   0, 0}    }, // 8
    /* ASCII 7bit letters */
    {'A',    "A",    { 0.33, 1.0,   0,   0,   0,   0, 0}    },
    {'B',    "B",    {1.0,  0.33,  0.33,  0.33,   0,   0, 0}    },
    {'C',    "C",    {1.0,  0.33, 1.0,  0.33,   0,   0, 0}    },
    {'D',    "D",    {1.0,  0.33,  0.33,   0,   0,   0, 0}     },
    {'E',    "E",    { 0.33,   0,   0,   0,   0,   0, 0}    },
    {'F',    "F",    { 0.33,  0.33, 1.0,  0.33,   0,   0, 0}    },
    {'G',    "G",    {1.0, 1.0,  0.33,   0,   0,   0, 0}    },
    {'H',    "H",    { 0.33,  0.33,  0.33,  0.33,   0,   0, 0}    },
    {'I',    "I",    { 0.33,  0.33,   0,   0,   0,   0, 0}    },
    {'J',    "J",    { 0.33, 1.0, 1.0, 1.0,   0,   0, 0}    },
    {'K',    "K",    {1.0,  0.33, 1.0,   0,   0,   0, 0}    },
    {'L',    "L",    { 0.33, 1.0,  0.33,  0.33,   0,   0, 0}    },
    {'M',    "M",    {1.0, 1.0,   0,   0,   0,   0, 0}    },
    {'N',    "N",    {1.0,  0.33,   0,   0,   0,   0, 0}    },
    {'O',    "O",    {1.0, 1.0, 1.0,   0,   0,   0, 0}    },
    {'P',    "P",    { 0.33, 1.0, 1.0,  0.33,   0,   0, 0}    },
    {'Q',    "Q",    {1.0, 1.0,  0.33, 1.0,   0,   0, 0}    },
    {'R',    "R",    { 0.33, 1.0,  0.33,   0,   0,   0, 0}    },
    {'S',    "S",    { 0.33,  0.33,  0.33,   0,   0,   0, 0}    },
    {'T',    "T",    {1.0,   0,   0,   0,   0,   0, 0}    },
    {'U',    "U",    { 0.33,  0.33, 1.0,   0,   0,   0, 0}    },
    {'V',    "V",    { 0.33,  0.33,  0.33, 1.0,   0,   0, 0}    },
    {'W',    "W",    { 0.33, 1.0, 1.0,   0,   0,   0, 0}    },
    {'X',    "X",    {1.0,  0.33,  0.33, 1.0,   0,   0, 0}    },
    {'Y',    "Y",    {1.0,  0.33, 1.0, 1.0,   0,   0, 0}    },
    {'Z',    "Z",    {1.0, 1.0,  0.33,  0.33,   0,   0, 0}    },
    /* Numerals */
    {'0',    "0",    {1.0, 1.0, 1.0, 1.0, 1.0,   0, 0}    },
    {'1',    "1",    { 0.33, 1.0, 1.0, 1.0, 1.0,   0, 0}    },
    {'2',    "2",    { 0.33,  0.33, 1.0, 1.0, 1.0,   0, 0}    },
    {'3',    "3",    { 0.33,  0.33,  0.33, 1.0, 1.0,   0, 0}    },
    {'4',    "4",    { 0.33,  0.33,  0.33,  0.33, 1.0,   0, 0}    },
    {'5',    "5",    { 0.33,  0.33,  0.33,  0.33,  0.33,   0, 0}    },
    {'6',    "6",    {1.0,  0.33,  0.33,  0.33,  0.33,   0, 0}    },
    {'7',    "7",    {1.0, 1.0,  0.33,  0.33,  0.33,   0, 0}    },
    {'8',    "8",    {1.0, 1.0, 1.0,  0.33,  0.33,   0, 0}    },
    {'9',    "9",    {1.0, 1.0, 1.0, 1.0,  0.33,   0, 0}    },
    /* Punctuation */
    {'\\',    "\\",    { 0.33, 1.0,  0.33,  0.33, 1.0,  0.33, 0}    },
    {'\'',    "'",    { 0.33, 1.0, 1.0, 1.0, 1.0,  0.33, 0}    },
    {'$',    "$",    { 0.33,  0.33,  0.33, 1.0,  0.33,  0.33,1.0}    },
    {'(',    "(",    {1.0,  0.33, 1.0, 1.0,  0.33,   0, 0}    },
    {')',    ")",    {1.0,  0.33, 1.0, 1.0,  0.33, 1.0, 0}    },
    {',',    ",",    {1.0, 1.0,  0.33,  0.33, 1.0, 1.0, 0}    },
    {'-',    "-",    {1.0,  0.33,  0.33,  0.33,  0.33, 1.0, 0}    },
    {'.',    ".",    { 0.33, 1.0,  0.33, 1.0,  0.33, 1.0, 0}    },
    {'/',    "/",    {1.0,  0.33,  0.33, 1.0,  0.33,   0, 0}    },
    {':',    ":",    {1.0, 1.0, 1.0,  0.33,  0.33,  0.33, 0}    },
    {';',    ";",    {1.0,  0.33, 1.0,  0.33, 1.0,  0.33, 0}    },
    {'?',    "?",    { 0.33,  0.33, 1.0, 1.0,  0.33,  0.33, 0}    },
    {'_',    "_",    { 0.33,  0.33, 1.0, 1.0,  0.33, 1.0, 0}    },
    {'@',    "@",    { 0.33, 1.0, 1.0,  0.33, 1.0,  0.33, 0}    },
    {'!',    "!",    {1.0,  0.33, 1.0,  0.33, 1.0, 1.0, 0}    },
    {0, NULL, {NULL}}
};

I implemented above SOM feature in FLDIGI to see if the morse decoding works any better than the token lookup based algorithm in MORSE.CXX.

It is hard to quantify the improvement dealing with real life noisy morse signals as there are so many variables.  In order to compare the algorithms I decided to create a test bed where I can run a known audio file with known morse code content and use FLDIGI  to decode the signals.

3) TEST RESULTS  

To test these new features  I created a very noisy morse signal using morse.m  by Rob, KL7NA.  The text was "THE QUICK BROWN FOX JUMPS OVER THE LAZY DOG - 1234567890". The noise level is 2  making it quite hard to hear the actual signal  from the noise. As expected the original FLDIGI  CW module has difficulties detecting the morse code. See figure 1 below.

Figure 1. FLDIGI - errors detecting morse code from noise

Using the new features I implemented (Matched Filter and SOM) I run the same noisy signal. See figure 2. below - now the text is almost readable and the amount of errors has significantly reduced.

Figure 2. FLDIGI - with Matched Filter and SOM features enabled

Looking at the FLDIGI signal scope below you can see the difference on noise level. Note that these were captured at different times so they represent different morse letters.  Figure 3 is without matched filter and figure 4 has matched filter enabled.  Reduction of noise is clearly visible.

Figure 3. Without Matched Filter

Figure 4. With Matched Filter

To enable testing Matched Filter I modified the FLDIGI Modems / CW / General configuration screen.  I added a checkbox and a "Matched Filter length" slider  that represents buffer length in 10ms units.  So  "227" in figure 5 below represents 2.27 second input buffer for the matched filter. The results improve the longer the buffer is set  but this also causes longer latency between the received audio and decoded characters appearing on the screen.  At 18 WPM morse speed  1...2 second delay seemed like a reasonable compromise. I set the range between  100ms and 3 seconds in the user interface. For longer messages such as bulletins even longer buffer could be applicable.

I am also using "Use Farnsworth timing" tick box to enable/disable SOM detection feature (I was lazy and did not want to create another tick box user interface). This enables to test SOM feature in real time and see if it makes any difference.

Figure 5.  Modified FLDIGI configuration screen

4) CONCLUSIONS

I implemented two new experimental features to FLDIGI software package, namely  Matched Filter  and Self Organizing Maps decoding.  My objective was to improve FLDIGI morse code detection and decoding capabilities.

Based on above testing  it seems that both detection and decoding capabilities improved  compared to baseline algorithms in FLDIGI.  However,  the testing above was done with artificially generated noisy morse signal, not with real life audio from noisy RF bands. I am planning to run more tests to verify how well these features work with real life cases.

These two experimental features could potentially be beneficial in the following areas: 

• Decoding morse signals with low signal-to-noise ratio such as in VHF / UHF / 50 Mhz bands  - combined with PSK reporter type functionality we could get automatic alerts on band openings,  such as Es or tropo, from FLDIGI 
• Monitoring weak CW beacons   
• Working with QRP stations 
• Perhaps even  working EME  QSOs  with smaller antennas

The software is alpha quality and not ready for broader distribution. Since I have only limited amount of time available I would like to collaborate with other hams who could help me to improve this software and hopefully make it ready for mainline FLDIGI distribution in future.

If you are able to compile FLDIGI from source tar files on Linux  and you don't need any support to resolve problems I can send the tar file via email.

5) ACKNOWLEGEMENT

I would like to express my gratitude to David H Freese, W1HKJ for making FLDIGI available as open source software.  This great software package is not only very useful piece of ham radio software but also a fantastic treasure trove of digital signal processing gems. 

I would also like to thank  Rob Frohne, KL7NA for making matched filter and many other algorithms publicly available for his students and providing great insights over multiple emails we have exchanged. 
 

73  de Mauri,  AG1LE

Experiment: Decoding multiple Morse code signals automatically on a noisy RF band

Library of Congress,
Prints & Photographs Division,
FSA-OWI Collection

While experienced CW operators can easily copy morse code from multiple stations in a pile-up even on a noisy RF band,  writing software to enable computers to accomplish similar task is quite difficult. CW Skimmer is an example of such software and it has gained popularity and some good reviews. However, CW Skimmer is not open source software and it has some limitations.

Rob Frohne, KL7NA has provided both excellent papers on the topic (1,2) as well as  working code examples to create, detect and decode morse code from noisy audio signals. I wanted to experiment and expand Rob's software written for Octave, open source software similar to Matlab.  I have been using Octave v3.4.2 on Ubuntu V11.1 running on my Thinkpad T43 laptop.

In order to do some experiments  and learn more how this problem could be solved I decided to break the problem into 4 different tasks, namely the following:

1)  Create a test case with 9 simultaneous morse code signals in different frequencies with noise.
2)  Create an algorithm to find all morse signals by frequency.
3)  Utilize  Rob's  matched filter code to improve signal to noise on selected frequency.
4)  Apply morse code decoding algorithm to filtered signals.

1)  CREATING A TEST CASE 

I used Rob's morse.m  version to create a test audio file simulating a snapshot from a contest.  The code below simulates 9 stations with different call signs in various phases sending "CQ TEST DE  XXXXX"  with different speeds on different audio frequencies. These signals are then scaled and added up to create roughly S1...S9 signals between 400Hz and 1200Hz. Resulting 10 second long audio signal is then saved on 'cwcombo.wav' file.


function morse_file()
% create a test audio file with multiple morse stations in a pile-up


Fs = 48000; % Fs is sampling frequency - 48 Khz
Ts = 10*Fs; % Total sample time is 10 seconds


% create 9 different parallel morse sessions - 10 seconds each at 20-35 WPM speed
%         TEXT         audio file  noiselevel Hz    speed WPM
x1=morse('CQ TEST DE AG1LE','cw1.wav', 0.03125,1200,Fs,20, Ts);
x2=morse('TEST DE SP3RQ CQ','cw2.wav', 0.0625, 1100,Fs,35, Ts);
x3=morse('DE W3RQS CQ TEST','cw3.wav', 0.125, 1000,Fs,30, Ts);
x4=morse('SM0LXW CQ TEST DE','cw4.wav',0.25, 900,Fs, 25, Ts);
x5=morse('CQ TEST DE HS1DX','cw5.wav', 0.5, 800,Fs, 20, Ts);
x6=morse('TEST DE JA1DX CQ','cw6.wav', 1, 700,Fs, 20, Ts);
x7=morse('DE JA2ATA CQ TEST','cw7.wav',2, 600,Fs, 20, Ts);
x8=morse('UA2HH CQ TEST DE','cw8.wav', 4, 500,Fs, 20, Ts);
x9=morse('CQ TEST DE CT1CX','cw9.wav', 8, 400,Fs, 20, Ts);


% weighted sum - merge all the audio streams together 
% 2x signal strength corresponds to 6 dB (one S-unit)
% 9 signals arranged S9 to S1 in frequency order 1200Hz ... 400Hz
y = 256*x1 + 128*x2 + 64*x3 + 32*x4 + 16*x5 + 8*x6 + 4*x7 + 2*x8 + x9;


% write to cwcombo.wav file 
wavwrite(y,Fs,'cwcombo.wav');


% look at the spectrum of combined signals - display 100Hz...2500Hz 
spektri(y,Fs,100,2500);


endfunction;

2) ALGORITHM TO FIND MORSE SIGNALS BY FREQUENCY

I created a small function "spektri"  to calculate FFT (fast fourier transform) and plot the audio spectrum in a given frequency range. The input is the above 10 second combined audio clip. The resulting spectrum is shown in Figure 1. below. The 9 simulated "stations" are 100 Hz apart from each other. Noise level is at around 90 dB and the  1200Hz  peak signal is at about 144 dB.  The range is about 54 dB corresponding to scale of roughly 9  S units if this signal would be coming from a real ham radio.  Signal x1 would peak at S9 and signal x9 would be S1 at almost noise level.

Figure 1.  Morse code  test case - spectrum 100 - 2500 Hz.

Visually it is relatively easy to find the 7 to 8 strongest peaks from the picture above. However, finding those frequency peaks automatically is more challenging due to noise and some splatter components.

After testing various peak detection algorithms I found the following that seems to work quite well:



function [maxtab, mintab]=peakdet(v, delta, x)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, and column 2 the found values.
%      
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA.


% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.


maxtab = [];
mintab = [];


v = v(:); % Just in case this wasn't a proper vector


if nargin < 3
  x = (1:length(v))';
else 
  x = x(:);
  if length(v)~= length(x)
    error('Input vectors v and x must have same length');
  end
end
  
if (length(delta(:)))>1
  error('Input argument DELTA must be a scalar');
end


if delta <= 0
  error('Input argument DELTA must be positive');
end


mn = Inf; mx = -Inf;
mnpos = NaN; mxpos = NaN;


lookformax = 1;


for i=1:length(v)
  this = v(i);
  if this > mx, mx = this; mxpos = x(i); end
  if this < mn, mn = this; mnpos = x(i); end
  
  if lookformax
    if this < mx-delta
      maxtab = [maxtab ; mxpos mx];
      mn = this; mnpos = x(i);
      lookformax = 0;
    end  
  else
    if this > mn+delta
      mintab = [mintab ; mnpos mn];
      mx = this; mxpos = x(i);
      lookformax = 1;
    end
  end
end



I created another function that plots the spectrum and detected peaks on the same graph using the peak detection function above.

function  fpeaks = spektri_p(data, Fs, str_f, stp_f, delta);
% plot spectrum of data 
% Fs = sampling frequency (48000) 
% str_f = start frequency range to plot
% stp_f = stop frequency range to plot
% delta  =  A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by delta




N = length(data);
spec = fft(data); % do fourier transform
df = Fs/N; % frequency bin size
minf = -Fs/2;
maxf = Fs/2 -df;
i = round(N/2+(str_f*N/2)/(Fs/2)); % start index of freq range
j = round(N/2+(stp_f*N/2)/(Fs/2)); % stop index of freq range


f = [minf:df:maxf]; % frequency axis i.e. [ -24kHz...+24kHz ] 
y = 20*log10(abs(fftshift(spec))); % dB magnitude
[maxtab,mintab] = peakdet(y(i:j,1),delta); % detect max peaks in range w/ delta dB difference
nr = length(maxtab(:,1))   % nr of peaks found



figure(2);
plot(f(1,i:j),y(i:j,1)); % plot audio spectrum within range
hold on; 
plot(f(1,maxtab(:,1)+i),maxtab(:,2),'r*');  % plot found peaks within range
hold off;
title('Spectrum');
xlabel('Frequency(Hz)');
ylabel('Volume dB');
fpeaks = f(1,maxtab(:,1)+i)'; % return found frequency peaks within range


Using above plotting function pks = spektri_p(y,Fs,300,1300,45);
creates the following figure 2. Adjusting the last parameter (delta) will increase/decrease the number of peaks detected. Peaks are marked with red stars in the graph.

Figure 2.  Frequency peaks marked with red star

 Looking at the result  with delta = 45  this algorithm detects all peaks correctly within few Hz accuracy.
The pks variable contains the detected values, see below:


pks =


    308.30
    400.00
    500.00
    600.00
    700.00
    799.90
    895.70
    999.80
   1103.50
   1199.90

If the program would be more interactive this parameter "delta" could be adjustable variable, similar to squelch function in some radios. Lowering this parameter value would work like squelch - more noise peaks would be detected. Increasing the value would  allow to detect only the strongest stations.

3) MATCHED FILTER TO IMPROVE SIGNAL-TO-NOISE RATIO

Now that we have automatically found the peak audio frequencies of the morse signals we can apply a matched filter to improve the signal-to-noise ratio.

Below is the matched filter algorithm. The parameters include the audio signal, estimated morse speed in WPM,  sampling frequency Fs, and audio frequency of the morse signal.



% This script shows how a matched filter (dot product) works.
% See http://en.wikipedia.org/wiki/Matched_filter for theory


function x_f = mfilter(x,speed,Fs,codef);
%x  is audio signal (read from wavfile) 
%Fs=48000; sampling frequency
%speed = 20; morse code speed in WPM
%code_f = 440; morse code audio frequency


dit_time = 1.2/speed;
x_length = length(x);


t=0:1/Fs:dit_time;


burst = sin(2*pi*codef*t);   % template of a "dit"
N = length(burst);


for k=1:x_length-N
xk = x(k:1:(k+N-1)); % time reversed signal
x_f1(k) = burst*xk;  % dot (inner) product with template
end
x_f = x_f1; %return filtered signal




Matched filter produces pretty dramatic effect - you can see it from the "before and after"  figure 3  below.
All the graphs show the actual audio signals - 480,000 samples over 10 second period.

Red graph shows the original 10 seconds audio clip produced by this line of code:
x6=morse('TEST DE JA1DX CQ','cw6.wav', 1, 700,Fs, 20, Ts);
Morse signal is barely visible from the noise. 


Blue  graph shows the combined 10 second audio clip produced by this line of code:   
y=256*x1+128*x2+64*x3+32*x4+16*x5+8*x6+4*x7+2*x8+x9;
The  strong S9 signal x1 dominates and the morse signal is clearly visible. Can you read the morse code produced by this line of code below by looking at the blue signal alone? x1=morse('CQ TEST DE AG1LE','cw1.wav', 0.03125,1200,Fs,20, Ts);
As there are 8 other signals embedded in the blue graph it is not as clean as in figure 4 where only x1 signal is visible.




Green  graph shows the output of  the matched filter:  y6 = mfilter(y,20,Fs,700). 
The matched filter extracts from the combined audio the original noisy S4 level x6 signal  at 700Hz and using the morse code template "dit"  makes the original morse code visible.   Can you read 'TEST DE JA1DX CQ' by looking at green graph alone? There are some extra noise peaks making it more difficult.

Figure 3. Matched Filter  demonstration

By using the results from section 2 above we can now run the matched filter on all peak frequencies where we found morse signals. This way we can extract automatically each "station" for the final morse code decoding step.

One area that is still open here is the morse speed. The speed is an essential parameter for the matched filter to create a "dit" template with proper duration  for the dot product.  In this example above we used a known speed. In real life morse speed has to be estimated. Once a few proper  "dits" and "dahs"  are received the speed estimation is pretty simple. Speed is  1.2 / dit_time. I have not implemented automatic speed tracking yet.

Matched filter has also some other interesting properties - more details in here.

4) MORSE CODE DECODING 

Decoding morse code seems intuitively a pretty simple task, but when the signal has significant amount of noise present it becomes more difficult to produce reliable results. 

Rob Frohne, KL7NA provides morse decoding software written for Octave in here.
The software reads the audio wavfile,  creates rectified and filtered morse signal and then detects short / long pulses, uses tokenized table lookup to decode morse code to text.

However,  Rob's version has some built-in assumptions about the audio signals.  He has a fixed threshold (0.05) and in the presence of noise this assumption can produce a lot of errors. I experimented with the threshold with many different kind of signals  and the following works a little bit better.


agc = max(y);
threshold = agc/2; 

Proper AGC  (automatic gain control) would be needed if the signal strength varies a lot over time. Otherwise decoder will use incorrect threshold to determine when the signal starts and ends. This leads to incorrect timing of "dits" and "dahs".

Rob's  slow wave filter  (y = filter(ones(1,20)/20,1, x2);)  has also a fixed assumption of  the window size.   If the morse speed is very high the slow wave filter can smoothen fast "dits" too much.  On the other hand noise spikes need to be filtered, otherwise decoder misinterprets those as "dits".  There should be some mechanism to adjust filter window size based on morse speed.  I tried the following code that made some improvement. It requires  "speed" parameter though. This is a "chicken and egg" problem of the morse code  - you need to receive at least one "dit" and one "dah" to determine the actual speed. Starting from some reasonable default value like 20 WPM is also possible.  



Fs =48000;
dit = 1.2 / speed;
dit_samples = Fs*dit;

flt = dit_samples;

y = filter(ones(1,flt)/flt,1, x2);

With relatively noise free signals this algorithm seems to work OK.
See figure 4. below for  S9  signal - original x1 version, produced by this line of code
x1=morse('CQ TEST DE AG1LE','cw1.wav', 0.03125,1200,Fs,20, Ts);
The upper red graph is the original signal - 480,000 samples over 10 seconds.  The bottom blue graph shows rectified & filtered envelope of the signal, and red overlapping line with threshold at 0.32 level determining "dits" and "dahs". 

Figure 4. Morse decoding process - original S9 signal

With somewhat noisy signals the matched filter helps - this S7 signal below in figure 5 is extracted using
 y3 = mfilter(y,30,Fs,1000)  - text "DE W3RQS CQ TEST" is still readable from the bottom red graph below. 

Figure 5.  Morse decoding process - S7 signal y3 extracted with matched filter.

With  noisy signals the matched filter helps a lot - this S5 signal in figure 6 below is extracted using
y5 = mfilter(y,20,Fs,800)  - text "CQ TEST DE HS1DX" is still feasible to decode despite noise spikes.

Figure 6. Morse decoding process - S5 signal y5 extracted with matched filter

With  noisy signals the matched filter helps but is not enough  - this S3 signal below is extracted using y7 = mfilter(y,20,Fs,600)  - text "DE JA2ATA CQ TEST" is missing some "dits". If threshold would be set lower then extra "dits" would be created from noise spikes. This will show up as errors in the decoding process.

Figure 7. Morse decoding process - S3 signal y7 extracted with matched filter

Finally with signals buried under noise the matched filter helps but is not enough especially in pile-up or contest situation.  This S1 signal below is extracted using y9 = mfilter(y,20,Fs,400)  - text "CQ TEST DE CT1CX" is no longer decodeable. Human eye is good for pattern matching but noise distorts the signal below to almost unreadable.  

Figure 8. Morse decoding process - S1 signal  y9 extracted with matched filter

CONCLUSIONS

Detecting and decoding multiple morse signals simultaneously in a ham radio contest or pile-up situation on a noisy RF band  is not an easy task.  Even for experienced human operators it is hard to copy weak S1..S3 signals in the presence of strong S7...S9 signals.

With some advanced digital signal processing algorithms  it is possible to automatically detect multiple stations with S1...S9 signal strengths close to each others on a noisy RF band.  The algorithm presented in section 2 above will correctly determine peak signal frequencies for further processing.

Matched Filter has some unique properties making it suitable for extracting maximum signal-to-noise ratio for the signal of interest.  Using the matched filter for peak frequencies found with the algorithm above enables  automatic decoding of S5...S9 level signals virtually error free.  However, for the weaker S1...S3 level signals the presented method will produce errors.

On the decoding software some minor improvements were implemented, such as threshold determination based on overall signal level, as well as  slow wave filtering using morse code speed dependent filter windows. These minor improvements made the decoding software a bit more robust.

Experimenting with Gnu Octave software is an excellent way to gain insight on signal processing concepts and also to understand better the limitations of morse code in the presence of noise.


 73 de AG1LE

PS. See also the latest experiment where above functionality was added on FLDIGI

DX pedition to KP2 / St Croix USVI - IOTA NA-106 in April 2012

I am active from  St Croix, USVI  island  (IOTA NA-106)
between April 4 - 11, 2012 using call sign  KP2/AG1LE.

Location is Radio Reef  (also known as KP2M contest station).

Radio Reef

Apr 8, 2012 update:  weather has been perfect, lot of sunshine and temperature near 30 C.   I have 1005 contacts in the logbook so far, propagation has been good on 20m, 17m, 15m, 12m and 10m bands. I will try to be on air from 6 AM local time as well as in the evenings. We have had also great time scuba diving & snorkeling so I won't be all the time working on radios ;-).

Apr 12, 2012 update:  All 1845 QSO are now uploaded to LotW and eQSL. 


73

Mauri  AG1LE