CW-Robot Octave Simulation

Notes on the Project

Notes on CW Robot Matched Filter

We can have quite a bit of latency in a robotic receiver, because it isn't usually used in QSO. We could set a latency limit, and use that with a spectragram to determine the dit time and synchronization. We could use a long duration FFT to determine the frequency, and make a matched filter to the dit. Any random phase between dits won't matter, as they will just be shifted by a portion of a cycle at the audio frequency. The more latency we can have, the better the demodulation as long as the station keeps sending with the same parameters.

We would determine the synchronization and filter frequency using a lot of data, then go back and apply the matched filter to the data and decode it maybe 30 seconds after it was sent (if 30 seconds was the latency and the transmission was longer than 30 seconds).

We need to keep track of each station in a QSO separately, and when they are transmitting.

It might be advantageous to do a short term FFT to downsample.

It might be advantageous to just down convert directly to DC before doing the matched filter, because then you won't have the series of peaks in the dot product.

Octave Scripts

Here are some Octave (or MATLAB) scripts your professor put together for you to help on your project.

Morse.m

This function called morse.m will generate morse code with noise added in for specified text.

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));
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

Check the Noise Audibly

% This script allows you to use your ear to see how much noise you have.
Fs=48000;
Fnyq = Fs/2;
speed = 20;
code_f = 440;
x_length = 2^19;
x=morse(' How r u?','morsecode.wav', 4, code_f, Fs, speed, x_length);
soundsc(x,Fs);
[b,a] = cheby1(4,2,[290/(Fnyq),690/(Fnyq)]);
figure(1);
clf;
freqz(b,a,Fs);
xfiltered = filter(b,a,x);
soundsc(xfiltered,Fs);

Matched Filter Demonstration

% This script shows how a matched filter (dot product) works.
Fs=48000;
speed = 20;
code_f = 440;
dit_time = 1.2/speed;
x_length = 2^19;
x=morse(' How r u?','morsecode.wav', 4, code_f, Fs, speed, x_length);
t=0:1/Fs:1.2/speed;
burst = sin(2*pi*code_f*t);
N = length(burst);
for k=1:x_length-N
xk = x(k:1:(k+N-1));
x_filtered(k) = burst*xk;
end
t1 = 0:1/Fs:1/Fs*length(x_filtered)-1/Fs;
t2 = 0:1/Fs:1/Fs*length(x)-1/Fs;
ts = 0:dit_time:x_length/Fs-2*dit_time;
figure(2)
clf;
plot(t1,x_filtered,'b-',ts,x_filtered(1:N:end),'r*',ts,400*ones(1,length(ts)),'r-')
%figure(3)
%plot(t2,x,'r-')

Finding the Peak Frequencies Automatically

% This script demonstrates how a CW filter helps.
% It demonstrates using find.
filter_width = 100; % Hz
Fs=48000;
x=morse('How r u?','morsecode.wav', 3, 440, Fs,30,2^18); % + morse('I am just fine.','morsecode1.wav', 1, 660, Fs,35,2^18);
sound(x,Fs);
x=x';
N=length(x);
t=0:1/Fs:(N-1)/Fs;
X=fft(x(1,:));
absX=abs(X);
f=-Fs/2:Fs/N:Fs/2-Fs/N;
xlabel('Frequency (Hz)')
big_bins_only = [absX > 20 *mean(absX) ];
peaks=find( big_bins_only);
m_width = floor(400/Fs*N);
filterH = zeros(size(X));
for m = 1:length(peaks)
filterH(1,peaks(m)-m_width/2:peaks(m)+m_width/2) = hanning(m_width+1)/2;
end
filteredX = X.* filterH;
filteredx = (ifft(filteredX));
max_imag = max(imag(filteredx)) % Sanity check.  (max_imag << max_real)
max_real = max(real(filteredx))
filteredx = real(filteredx);
sound(filteredx,Fs);
figure(1)
plot(f,fftshift(absX),'b-',f,fftshift(filterH*max(absX)),'r-',f,fftshift(filteredX),'g-');
figure(2)
plot(t,x,t,filteredx)

Find the maxima of a function

function [peak_indices] = maxima(a)
% Finds the local maxima of a vector (1xn) and returns the indices
n=length(a);
r = [a(1,2:n) - a(1,1:n-1) <= 0, 1]; % derivative on the right <= 0
l = [1, a(1,2:n) - a(1,1:n-1) >= 0]; % derivative on the left >= 0
d = l & r;
peak_indices = find(d);