Monday, December 31, 2012

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








Tuesday, December 25, 2012

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









Monday, July 30, 2012

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.  




Wednesday, July 4, 2012

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.  







Friday, June 29, 2012

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

Sunday, June 24, 2012

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



Popular Posts