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

Ultimate Morse Code Decoder?

Over the last few weeks I started reading articles and papers about wavelet transformation to figure out how to build an ultimate Morse code decoder. Wavelets have been used since 1980's in digital signal processing and wavelet transforms are now being adopted for a vast number of applications, often replacing the conventional  Fourier transformation.

One particular application is for smoothing/denoising data based on wavelet coefficient thresholding, also called wavelet shrinkage. By adaptively thresholding the wavelet coefficients that correspond to undesired frequency components smoothing and/or denoising operations can be performed.

WAVELETS FOR DETECTING NOISY SIGNALS

I stumbled on this paper from Aly, Omar and Eldherbeni few days ago. The paper describes an algorithm for extracting and localizing an RF radar pulse from a noisy background. The algorithm combines two powerful tools: the wavelet packet analysis and higher-order-statistics (HOS). The use of the proposed technique makes detection and localization of RF radar pulses possible in very low signal-to-noise ratio conditions. 

The proposed algorithm is able to detect and well localize RF radar pulses without a prior knowledge of the pulse parameters (e.g., its frequency and duration). The proposed  algorithm has been tested for SNR down to −24 dB and proved to work successfully. See figures 9 and 7 below on the results they achieved.

In my previous blog post  I used matched filter to extract Morse code pulses from noisy signals. Matched filter can be implemented  either in time or frequency domain. However,  one problem with narrow matched filter is "ringing" artifacts  - see figure 8  in the blog post. This creates uncertainty and jitter in pulse width.  According to above paper  one of the advantages of  using wavelets is possibility to obtain adapted tiling in the time-frequency plane, which is automatically generated based on the signal observation. This leads to improved time & frequency resolution compared to traditional short-time windowed Fast Fourier Transformation (SFFT) based methods.  

DIFFERENT WAVELET TRANSFORMATIONS

I found this website that provides online tools to visualize various wavelet types. Since the discovery of wavelets many years ago there has been active development and many different wavelets have been created for various purposes. This site provide very nice online visualization tool.

I created  a small noisy Morse code file using  morse.m  Octave function by Rob Frohne, KL7NA.  The above website allows only 2000 samples so I copied numbers to the above site (select "Your data" section).
Figures 1..3  below show a small noisy audio section with a "dit" tone  and start of "dah" tone and corresponding wavelet transformation (Morlet, Gaussian & Paul wavelet types) underneath.  There is also a global wavelet showing the variance by period on the right side.  The wavelet transformation values are color coded - red color is showing the highest energy at Period (scale) 8 clearly identifying where the tone starts and where it ends.  Different wavelets have slightly different properties as is quite visible looking at the figures below.

Figure 1.  Morlet Wavelet - noisy Morse code

Figure 2.  Gaussian Wavelet - noisy Morse code

Figure 3.  Paul Wavelet - noisy Morse code

MORLET WAVELET TRANSFORMATIONS 

As the results from previous section demonstrate Morlet type wavelets seem to provide particularly good delineation of Morse code tones in noisy signals.   I did three more visualizations with  various level of noise included in the signals.  The first one (MORSE_SNR0) below has no noise and the signal is visible as clear red color  line on wavelet power spectrum image (period = 8). It is very easy to see where the signal ends and next one starts (see the gap between samples 1450 ...1900 ). 

The following visualization (MORSE_SNR1) has approximately - 9 dB  SNR  and quite a lot of noise. The red line is still visible in the wavelet power spectrum image (period = 8). However,  noise spikes make it a bit more difficult to determine signal end and start timing. There is also other red noise signals at periods 32, 64 and 256. As these are on different wavelet periods they can be filtered out. 

The last visualization (Morse_SNR2)  has  approximately  - 12 dB SNR  and it is quite difficult to see where the tone ends and next one starts. Looking from original signal it is almost impossible to tell where the tone is.  Global wavelet shows a peak at period 8  where most of the red dots are also aligned. There are many more red areas showing noise energy peaks.

WAVELETS FOR PULSE TRAIN

I  did another experiment with Morse code pulse train.  This is the envelope of the noisy audio signal after detection and filtering. Pulse train has some sharp edges that should correspond to high frequency components as well as longer stable plateaus corresponding low frequency components. 

I used  Haar and Gaussian wavelets in this experiment. Figure 4 and 5 below show the signal and corresponding wavelet transformation.  As expected the high frequency components are visible where the signal edges are.

Figure 4. Haar Wavelet - noisy Morse code pulse train

Notice the red high energy components between 600 - 750 ms in scale 128..512 range. This represents noise that is visible also on the time domain signal in figure 5.  In a wavelet filter implementation these  coefficients could be set to zero  to denoise the signal.

Figure 5. Gaussian Wavelet - noisy Morse code pulse train

CONCLUSIONS

Wavelet transformation  is a powerful signal processing tool to manipulate signals. Based on literature the wavelet transform gives better localization in the time-frequency domain than the discrete windowed Fourier transform.  Wavelets enable also flexible manipulation of the signals to remove noise, find signals buried under noise  etc. 

For a real time CW decoding software like in FLDIGI  wavelet transformation could open a whole new performance level dealing with noisy signals.  More experimentation is definitely needed to unlock this potential.


73
Mauri AG1LE

FLDIGI - Analysing SOM decoder errors

In my previous blog post I shared some test results on new alpha version of Fldigi.

The new version of the Fldigi software works pretty well but occasionally it still generates some decoding errors.  I spent some time today in instrumenting the CW.CXX module to collect some measurements during the normal CW decoding operation.

COLLECTING DATA 

The first place to focus was on how the software detects  "dits" and "dahs" from the incoming signal.
The cw::handle_event() function keeps track of CW_KEYDOWN_EVENTs and CW_KEYUP_EVENTs based on signal thresholds in the new  cw::decode_stream() function.
Therefore this was an obvious place to put a hook.  I created a new function cw::histogram() that collects the current duration value after CW_KEYUP_EVENT is detected. I placed the collection after the noise spike code as these values are not collected for decoding anyways.

I recorded about 3:17 into an audio file from 14 Mhz band with multiple CW stations having contacts. The band was quite noisy and I could see some decoder errors on the Fldigi instance I had connected to my Flex3000.  I replayed the processed audio file on Linux environment with instrumented version of Fldigi.  Figure 1 below shows the probability distribution of "dits" and "dahs".  There are multiple peaks visible and I also noticed that automatic Morse speed tracking changed from 16 WPM to  13 WPM, corresponding to "dit" values of 75ms and 92 ms respectively.  There are also a small number of outliers between 100 and 250 ms range.

Figure 1.  Dit & Dah timing distribution.

Since I was using the SOM decoder  for this experiment I decided to utilize a very nice feature built-in to the Best Matching Unit algorithm.  Since we are calculating Euclidian distance between incoming data vectors and codebook weight vectors and selecting the codebook vector with smallest distance we can use this distance as an error metric.  In other words if the best matching unit did not really provide a good match the distance (diffsf variable in the find_winner() function) should be higher than normally.
I  plotted figure 2  in order to demonstrate where SOM decoder indicated it had trouble matching the right codebook entry.  Not surprisingly when looking at each decoded character with high (over 0.05) error metric you can see some problem with the input vector.  The mean error value was 0.0125 and standard deviation was 0.047556. As you can see the the figure 2 the maximum was 0.50672  and there are several other peaks over 0.1 below.   Each entry on x-axis corresponds to one detected character and y axis is the SOM decoder error.

Figure 2.  SOM decoder - errors over time

Since there are so many error peaks in the figure 2  I started wondering if there is something else than just noise peaks that could be causing problems for the SOM decoder.  In most cases where SOM decoder indicates a bad match there was either additional "dit" or "dah" elements concatenated to the input vector  or some values that were not "dit" or "dah"  according to timing distribution shown in figure 1.

I let the voice file play again and plotted two other variables on the same timescale.  Agc_peak is used to determine with the threshold when signal is considered going up or down and it is a variable itself with fast attack and slow decay. If agc_peak falls down that means that signal level is also going down giving an indication on potential signal-to-noise problem.

Cw_adaptive_receive_threshold is a variable tracking the duration of 2 * "dit" length with a time constant defined by trackingfilter. This is a key variable determining whether received key down event was a "dit" or "dah".  If this variable is not able to follow the changing Morse speed then SOM decoder could get incorrectly assigned "dit" and "dah" values in the  input vector.

I normalized these 3 variables to  [0.0 1.0] range and plotted them on the same graph to see if there are any dependencies.  Figure 3  shows  SOM error rate,  agc_peak and  cw_adaptive_receive_threshold variables over time.  X axis represents time of each decoded character (sampling done when SOM decoder emits a character) and Y axis is normalized value of these variables.

Figure 3.   Error rate, Automatic Gain Control peak and CW speed over time 

By looking the figure 3  it is not obvious that agc_peak or cw_adaptive_receive_threshold values would correlate with higher SOM error rates.  Interestingly, even when the  agc_peak value goes down significantly showing that signal has some fading down to S2..S3 level, the SOM error  rate is not increasing during these fading events.

I decided to have another look at the audio file.  The original audio file recorded by PowerSDR/Flexradio 3000 sounded OK when played with Windows media player.  I usually import audio files to Linux environment where I have development environment.  Alsaplayer is also playing the audio files OK  but for some reason Fldigi  does not show waterfall when playing files originated from PowerSDR.  I have used Octave to do a format conversion:

PowerSDR format: Uncompressed 32-bit IEEE float audio, stereo, 12000Hz sample rate
Octave format:        Uncompressed 16-bit PCM audio, stereo, 8000Hz sample rate.

I copied the Octave formatted audio file back to Windows environment and heard what sounded like clipping  (there was several CW stations sending at the same time on this sample).

I decimated the audio file to the same length as measurements in the above figure (235 samples),  took absolute value and plotted the signal on the graph with other variables.  Now the culprit for these odd errors  is more visible - the signal is  clipping severly at x = 20...40,  x = 125..130, x = 190..210.  In the Fldigi waterfall display this clipping caused some visible interference spikes.

Figure 4.  Signal clipping is causing SOM decoder error rate peaks  

CONCLUSIONS 

SOM decoder has an error metric  feature that is very useful in debugging problems.  In this particular case I was able to track the problem down to incorrectly converted audio file that caused clipping and artefacts on the frequency that Fldigi was decoding.

This decoding error metric could be used for other purposes as well.  Instead of printing incorrectly decoded characters on Fldigi display we could establish an  error limit  perhaps based on mean & standard deviation of normal  SOM decoding.  If error metric goes above this limit we can either stop printing characters or show  "*"  like  the original decoder is doing when it cannot decode.

For hams who want to learn manually to send high quality Morse code these features could provide some numerical metric on progress.  Non-gaussian "dit" / "dah" distribution indicates problems in rhythm   and SOM error metrics indicates timing problems in patterns.


73
Mauri AG1LE

Morse Code decoding - machine learning

In the previous blog entry I covered an experiment using Self Organizing Maps with SOM toolbox to automatically learn Morse code from the input vectors that contained duration of the "dit" and "dah" tones.

While I have been reading academic research papers  Dave W1HKJ and  alpha testing team has been working on a new version of FLdigi.

NEW ALPHA VERSION OF FLDIGI

I was testing the brand new Fldigi alpha release during the CQ WW WPX contest  as the band was literally full of CW stations around the world. It was very difficult to find an empty frequency on the CW bands. Figure 1 below shows a screenshot from my Flex3000  PowerSDR pan adaptor / waterfall display. I recorded  a  230 MB size IQ audio Wav file for further testing purposes.

Figure 1.  CQ WW WPX  contest at 7 Mhz CW band

One of the practical challenges is how to detect the start and end point of the tone correctly in the presence of noise or interference. When the detected signal amplitude varies greatly due to fading, noise or interference having a simple threshold does not work well.  For example using the upper detector threshold in figure 2 would only produce two "dits" (aka "E E").  Having two detector thresholds and some hysteresis will produce better results.

Figure 2. Upper and lower threshold.

In the  FLdigi 3.22.0CB alpha version Dave W1HKJ has added Lower and Upper Detector thresholds parameters that are also visible in the Scope display.  See figure 3 below.  

Figure 3.  Detector Threshold settings

This simple enhancement works surprisingly well.  In figure 4 below I have used SOM detection and alternatively used FFT filtering and Matched Filtering.  As you can see from the figure there are some stations with very high signal strength (you can see the Morse code spectrum spreading over 2 kHz below from 7072.7 kHZ ) and some with much smaller signal strength. I am getting some decoding errors (mostly extra "E", "I", "T" characters)  but the new decoder works quite well compared to previous decoder. Here are some ideas how to improve tone start and end point detecting even more, such as using energy and Zero Crossing Ratio (ZCR) values as thresholds.

Figure 4.  FLDigi 3.22.0CB version in action

HOW TO ENGAGE WITH THE MACHINE LEARNING COMMUNITY

Over the last 2 weeks I have read many research papers on various topics such as Hidden Markov Models, Dynamic Time Warping, matching patterns in time with Bayesian classifiers,  Restricted Boltzman Machines (RBM), Deep Belief Networks (DBN) etc.

Some of these machine learning algorithms tend to focus on creating a set of features that one (or multiple) classifier(s) then use trying to find the best match to pre-labeled training data. This PhD theses is somewhat old (1992) but well written and it covers most commonly used techniques and algorithms.

Other alternative is "unsupervised learning" approach where a lot of unlabeled training data is used to find patterns & clusters that can then be classified and labeled.  Especially RBM and DBN camp seems to believe this is better and more universal approach. This presentation in Youtube (see demo 21:38 forward)  by Geoffrey E. Hinton  explains well this latter approach.

Just looking at the amount of research done on machine learning it is quite a challenge trying to find the best way to move forward. One idea that I was playing with was to create a database of Morse code audio wav files as training & testing material  - it looks like the machine learning community has focused so far only on broad speech, music and  environment audio categories. These type of training sample  audio files can be found in many websites of machine learning teams.  Having good quality labeled material seems to be a big problem for the machine learning community.

I think it would be a worthwhile technical challenge  to advance the state of the art in machine learning by automatically decoding all the hundreds of Morse code conversations taking place during these ham radio contests. It would push the machine learning algorithm development further and hopefully bring some bright minds into our ham community as well.

I wonder if ARRL or some other organization could setup this kind of public challenge for schools and universities engaged with machine learning ?  As we can see from the above Morse code is very much alive and being used every day.  It would be relatively easy for hams to produce such labeled audio content. We could use FLDIGI to decode the audio file content as a training reference and have different kinds of CW files recorded from real ham RF bands, like the  above CQ WW WPX contest example on 7 Mhz


73
Mauri  AG1LE

















EPD algorithms

Morse Code decoding with Self Organizing Maps

Continuing the series of articles on the Morse code decoding using Self Organizing Maps (SOM)  I am covering some new experiments with SOM based learning.

Using the experimental version of FLDIGI setup explained in my previous blog entry  I collected a dataset of  dit/dah durations from the following decoded text:

"<AS>CWDEWFHKJW1HKJW1HKJPSEMCQCQCQDEW1HKJW1HKJW1HKJCQCQCQDEW1HKJW1HKJW1HKJPSEKCQCQCQDEWFHKJW1HNJW1HKJCQCQCQDEW1HKJW1HKJW1HKJPSEKCQCQCQDEW1HKJW1HKJW1HKJCQCQCQDEW1HKJW1HKJW1HNJPSEKCQCQCQDEW1HKJW1HKJW1HKJCQCQCQDEW1HKJW1HKJW1HKJPSEKCQCQCQDEW1HKJW1HKJW1HKJCQCQCQDEW1HKJW1HKJW"  

I used the SOM toolbox to create a 7 x 7 SOM  and learn the 14 different morse characters in the text above based  purely on the dit/dah timing information.  Each morse character is vector with 7 numbers  (6 for dit/dah durations in milliseconds and 7th number as terminating zero).  Here is an example of the data:



56  160 156   0   0 0 0 W
68  162 176 152 154 0 0 1
58   60  66  46   0 0 0 H
174  64 156   0   0 0 0 K
44  154 160 146   0 0 0 J

Once the 7 x 7  rectangular SOM was created and automatic learning process was completed I created a few visualizations to look at the data.  

VISUALIZATIONS 

 Figure 1 below shows a hit histogram in colors and numbers.  For example top left corner node has 26 hits from input vectors matching letter "Q".  SOM was also able to cluster the  learned characters automatically.  Dark lines between nodes depict the clusters that were created by Ward's linkage feature of the SOM toolbox.   

Figure 1.  SOM - hit histogram and clusters

The next visualization in Figure 2 below  is a SOM nearest neighbourhood graph. The method defines graphs resulting from nearest neighbor- and radius-based distance calculations in the data space, and shows projections of these graph structures on the map. You can observe how relations between the data are preserved by the projection, yielding interesting insights into the topology of the mapping, and helping to identify outliers. 


For example you can see that letters  Q (dah-dah-dit-dah),  M  (dah-dah) and D (dah-dit-dit)  are neighbours and  connected via a graph projected with red line on the SOM map below.  The line segment D <-> M  is quite long indicating a topology violation. You can also see that different variations of H, K, J, W, C are nearest neighbours. 


The original signal contained a lot of noise which was causing jitter on the dit/dah duration, and sometimes even errors.  However,  SOM was able to learn from this dataset the Morse characters having similar patterns. 

Figure 2. SOM Neighbourhood graph k-nn

Another visualization technique is U-matrix in figure 3 below. U-matrix colors the nodes with the average distance of that unit to other adjacent units. The blueish colors represent close distance whereas red colors represent long distance. The SOM network appears to have clustered the morse characters in distinctive groups. You can compare how U-matrix and  clusters (depicted by black lines) are aligned.  

Figure 3. SOM U-matrix visualization

CONCLUSIONS 

Self Organizing Maps  is a powerful  method to visualize multidimensional data as it preserves topological properties of the input data.  In this example we used a set of Morse code patterns as training data. Each Morse code consists of  a "dit/dah"  pattern representing tone duration in milliseconds.  The original data was collected from noisy audio files using  FLDIGI  that had experimental matched filter and SOM decoder features enabled.


The author has done similar testing last year with a much larger dataset collected from actual CW contacts recorded in multiple RF bands (7 Mhz, 14 Mhz and 18 Mhz bands)  from multiple stations with different CW styles.  The dataset contained almost  40,000 characters  covering QSOs,  rag chewing, bulletins, and other types of CW communication.   Testing was performed with 20 x 20, 10x 10 and 7 x7 sized SOMs.


SOM was able to cluster similar Morse code patterns together despite noise and jitter on the original data.  We also demonstrated that automatic clustering and character classification would be feasible using the features available in the SOM toolbox.

Looking at the weight vectors created by this automatic learning process they resemble very much the SOM codebook that the author created for the SOM decoder function for FLDIGI software.  Using the  codebook entry for each cluster (see the letters Q,M, H, E, S, D, J, P, W, K, C, 1 in figure 3) and Best Matching Unit (BMU) algorithm would classify the incoming data vector correctly despite some noise and jitter.

There are still multiple unresolved problems in error free automatic detection and decoding of Morse code.  As the signal-to-noise ratio decreases the noise and jitter  makes it more difficult to produce accurate timing information.  Also, some CW stations don't comply with  dit / dah timing standards which creates more variability in the data.  Self Organizing Maps is one approach to tackle these problems  and based on preliminary testing it shows a lot of potential for further improvements.

73
Mauri AG1LE

FLDIGI: Matched Filter and SOM decoder - new ideas

In my previous blog entry I described experiments adding Self Organized Maps (SOM)  based Morse code decoder into FLDIGI.  

I have some new test results on the experimental FLDIGI code  as well as some ideas how to improve SOM decoder.

NEW TEST RESULTS 

I did a new comparison test with a FLDIGI patch on 3.22RC that Dave W1HKJ provided.

Dave is saying:
" I completed the fast-convolution-fft implementation of the matched filter for
the CW decoder. Also changed the decimate in time processing in the main Rx
process loop. No special attention is needed for audio stream replay. This
version of your mod works very well and is very cpu efficient. The combination
of matched filter and SOM decoding give excellent results!"

The  9 test wav files with s/n ranging from -3 dB to +6 dB used are in here  and Dave's test results are in here .  I added the test audio files in a playlist ordered from -3dB to +3dB signal to noise ratio (SNR)  and played them in that order with AlsaPlayer.

To demonstrate the improvements I took two screen captures of the experimental software.

Figure 1 shows the new experimental matched filter and SOM decoder features enabled.   The new features enable CW decoding from -3dB SNR upward relatively error free.

Figure 1. Matched Filter and SOM decoder enabled

Figure 2 shows the new experimental matched filter and SOM decoder features disabled. The legacy CW decoder is relatively error free  from +1dB  SNR upwards.
 

Figure 2. Matched Filter and SOM decoder disabled

IMPROVEMENT IDEAS FOR  SOM DECODER

There has been also discussion in Linuxham mailing list on the SOM decoder.
Several people have shared ideas on how we could improve CW recognition / decoding using different approaches, like  Hidden Markov Models   and Scale Invariant Transforms.

My focus in this work has been to improve CW decoding in the presence of noise. However, there are multiple other problems to solve in the CW decoding space. Here are some examples: 

• Decoding CW from those irregular fists and all the bug users that have very short dits and long dahs, way different than the standard 1:3  dit:dah timing ratio
• Decoding variable speed  CW  - I have heard stations who give other station's call sign manually at 15 WPM  and  give rest of the message from elbug memory at 25 WPM. For human operators this may be OK if the message is standard exhange like RST and serial number. However, for decoding software this rapidly varying speed may be a problem.
• Decoding multiple CW stations simultaneously - all with different SNR, speed and timing.

The timing problem is obviously also quite important as there are still many hams working CW manually and timing errors can also be generated by noise, operating practices etc.

I added some debugging code to FLDIGI to be able to extract the timing of decoded characters.  Figure 3 below shows the distribution of "dits" and "dahs" taken from the above 9 computer generated noisy wav files.  The speed in the above example is 24 WPM  so  "dit" length should be 1.2/24 = 50 ms.  Accordingly  "dah" length should be 150 ms.

As you can see from figure 3   noise creates jitter in timing.  The sample size was only 267 decoded characters but you can already observe a bell shaped normal distribution being formed below.

Figure 3.  Dit/Dah timing histogram (horizontal axis in milliseconds). 

You can also see some outliers around 250 ms and 350 ms area. We have  4/267 outliers in 250 ms range and 2/267 in 350 ms range.  We are talking about  6/267  cases so approximately   2.2%  of incorrectly decoded characters due to noise in this sample of 267 characters.

 When looking at these outliers I found  letter Q misinterpreted as W  in the first CQ  sequence:

Case 1:

== duration (ms)==== CHAR

164 62  158 60 0 0 0 C     //  dah dit dah dit = "C"
156 172 258 0  0 0 0 W     //  dit dah dah   = "W" (should be "Q")

SOM algorithm looked at this pattern and tried to find best match to what looks like   dit - dah - dah  sequence (W) after normalization.  However,   what apparently happened here was that  letter  Q  ( dah - dah - dit - dah)  got scrambled by noise and the last dit - dah was merged to one extra long dah with 258 ms duration. 

Below is another example.  W1HKJ  was misdecoded as  WFHKJ. When looking at the details it looks like a noise spike merged two 'dahs' into one extra long  'dah' with 360 ms duration.

Case 2:
== dit/dah duration (ms) = CHAR
70  164 154 0   0 0 0 W
58  152 360 158 0 0 0 F       //  dit-dit-dah-dit   = "F" (should be "1")
62  56  56  50  0 0 0 H
162 76  160 0   0 0 0 K
70  158 158 170 0 0 0 J

These examples above give some hints how to improve SOM decoder even further.  Right now  SOM decoder treats every character separately.  After normalization it just tries to find "Best Matching Unit"  (BMU) based on input vector  (dit/dah durations) by calculating Euclidian distance to codebook entries.  The codebook entry with minimum distance is the "winner". In the current version the SOM codebook is fixed and does not use SOM learning algoritms.

OUTLIER PROCESSING IDEA

We could maintain a table that learns the  frequency (or probability) of dit and dah durations.  If the incoming pattern of dits and dahs  fits under the normal distribution  we can proceed to BMU algorithm and find the winner.  However,  if the incoming pattern has an outlier  we need to figure out most likely alternative before we pass the pattern to BMU.

For example in the above case 1  - 258 ms is an outlier.  Since dit time mean is 50 ms  ( ~ space = 50ms, ~ dah = 150 ms) we have the following potential alternatives to replace the outlier:

Alt1: dit - dah         =  50ms + 50ms + 150ms  = 250 ms 

Alt2: dah - dit         = 150ms + 50ms + 50 ms  = 250 ms 

Alt3: dit - dit - dit   =  5 x 50 ms = 250 ms
etc.

Input: 156 172 258 0  0 0 0   // 258 is outlier, try alternatives 

Alt1:  156 172 50 150 0 0     //  = "Q"

Alt2:  156 172 150 50 0 0     //  = not in codebook

Alt3:  156 172 50 50 50 0     // = "7"

In the above case 2  -  360 ms is an outlier.  We have the following potential alternatives

Alt1:  dah - dah        =  150ms + 50ms + 150ms  = 350 ms

Alt2:  dit - dit - dah  =  50 + 50 + 50 + 50 + 150 = 350 ms 

Alt3:  dah - dit - dit  =  150 + 50 + 50 + 50 + 50 = 350 ms 

etc.

Input: 58  152 360 158 0   0   0   //  360 is outlier, try alternatives

Alt1:  58  152 150 150 158 0   0   // = "1"    
Alt1:  58  152 50  50  150 158 0   // = not in codebook
Alt2:  58  152 150 50  50  158 0   // = not in codebook  

 These alternatives could be arranged in some sort of probability order. I am not sure how to measure likelyhood of one noise spike  vs. multiple noise spikes impacting the same  character. Intuitively one noise spike would be more probably but this depends on morse speed, noise type and other factors I assume.

Alternatively, if we know typical outlier alternatives we can just add the most likely cases in the SOM codebook and let the BMU algorithm to match incoming vector (with outlier)  with the SOM codebook.  This makes the processing really simple and focus shifts to maintaining an up-to-date  SOM  codebook that includes error handling cases.  We just need to collect enough data to learn about most typical  error cases.  This could be done also via some kind of web service where FLDIGI clients send error cases with enough context  and in return will get updated SOM codebook.  This way  FLDIGI community would improve the decoder & codebook by exposing the SOM algorithm to a large variety of real world cases.


Update May 27, 2012:  
 I added a new histogram collecting function in the case CW_KEYUP_EVENT:  section of  cw::handle_event() function.  I also added a new histogram display to FLDIGI  to visualize cumulative distribution of detected signals. As the program runs and detects CW signals it builds a table of possible "dit" and "dah" durations.  Figure 4 below shows an example distribution collected from over 20 minutes period. Notice that current experimental prototype does not include logic  to accomodate Morse speed changes - it uses the same cumulative distribution for all stations. We should zero the distribution table every time we change the station we are listening, if the speed varies a lot.

Figure 4.  Cumulative distribution of "dit" and "dah" durations. 

73
Mauri  AG1LE