PSK31 Demodulation (Kurt & Michael)
Project Description
The goal of this project is to design and code a Matlab script that will encode and decode a PSK31 signal including signals with noise. The receiver should be able to read in signals (as a wav file) from other sources as well.
PSK31 is a audible text encoding that can be sent over the air by amateur radio operators. A computer's sound card can be used to send and receive the signal since the signal is audible. For more information regarding PSK31 see the Wikipedia article.[1]
Our Approach
Code Overview
Transmitter
Our code creates a PSK31 signal given an input carrier frequency and message. For testing our receiver code, the transmitter is setup to generate a random carrier frequency and phase. It also adds random noise to the signal before writing it to a wav file.
clear all
close all
clc
fc = 100+rand(1)*3900 % Carrier Frequency
%fc = 1500;
Fs = 10000; % Sampling frequency
baud = 31.25;
%% Encode the message
phase = rand(1)*2*pi; % Phase of the signal
dt = 0:1/Fs:1/baud-1/Fs; % Time between samples
one = ones(size(dt)); % A 'One' signal
zero = cos(pi*dt*baud); % A 'Zero' signal
messageString = 'Hello World! The quick brown fox jumps over the lazy dog. ABCDEFGHIJKLMNOPQRSTUVWXYZ abcdefghijklmnopqrstuvwxyz 0123456789';
[message charTable] = make_bits(messageString);
baseband = [];
m=1;
for k = 1:length(message)
if(message(k)==0)
baseband = [baseband, zero*m];
m=-m;
else
baseband = [baseband, one*m];
end
end
t = 0:1/(Fs):length(baseband)/Fs-1/Fs;
v = baseband.*cos(2*pi*fc*t+phase);
% Noise
v = v+randn(size(v))*10^(-10/20);
wavwrite(v,Fs,'signal.wav');
figure(1)
subplot(2,1,1)
plot(t,baseband)
axis([0 2 -1.2 1.2])
title('\fontsize{14}Baseband')
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,1,2)
plot(t,v)
axis([0 2 -1.2 1.2])
title('\fontsize{14}Baseband with Carrier')
xlabel('Time (s)')
ylabel('Amplitude')
movegui('northwest')
Receiver
Our receiver is split into several steps:
- We read in the signal from a wav file and generate utility variables and matrices.
- We run the signal through an FFT and take an average over the range of the FFT spike to obtain a carrier frequency guess.
- A PID controller was implemented to help offset the fact that a carrier frequency derived from the FFT produces an inaccurate carrier frequency. An offset in the carrier frequency causes the constellation diagram to rotate in a circle. The PID tries to compensate for this by adding to or subtracting from the phase of the signal while trying to drive the Q-component of the signal to zero. The feedback loop is described by the code below.
- We multiply the signal with a cosine wave at our guess frequency. This splits the signal into two parts: a high frequency component (with frequency equal to the sum of the actual carrier frequency and our guess frequency) and a low frequency component (due to the difference of the two frequencies).
- We filter the multiplied signal through a butterworth filter with a 75Hz cutoff frequency to remove the high frequency component.
- We mark the locations where the filtered signal changes sign. This identifies possible "sync" points where zeros in the signal may occur.
- We compare the length of time that occurs between each sync point with the baud rate to determine how many baud periods occur between each sync point. A single baud period indicates a zero. Each additional period over the first baud period indicates a one. (Ex. Three baud periods -> '110')
- We send the resulting bit stream through a text decoder which translates the ones and zeros into a text message.
%% Decode Message
clear all
baud = 31.25;
% [v,Fs] = wavread('PSK31_sample.wav');
% [v,Fs] = wavread('capture.wav');
[v,Fs] = wavread('signal.wav');
dt = 0:1/Fs:1/baud-1/Fs; % Time between samples
one = ones(size(dt)); % A 'One' signal
zero = cos(pi*dt*baud); % A 'Zero' signal
t = 0:1/(Fs):length(v)/Fs-1/Fs;
[sout charTable] = load_alpha();
upperNoiseCutOff = .75;
lowerNoiseCutOff = 1-upperNoiseCutOff;
% figure(1)
% plot(v)
%% Determine Carrier Frequency
vFftShift = abs(fftshift(fft(v)));
freq = linspace(0,length(v)-1,length(v))*(Fs/length(v));
freqShift = (freq - Fs/2);
% remove noise components
noise = (vFftShift < lowerNoiseCutOff*max(vFftShift));
vFftShift(noise) = 0;
% uses the center(mean) of the fft spike as the frequency
freqRange = find(vFftShift(length(vFftShift)/2:length(vFftShift)) > 0);
carrierFreq = freqShift(round(mean([max(freqRange),min(freqRange)]))+length(vFftShift)/2-1)
figure(2)
plot(freqShift,vFftShift)
movegui('center')
% PID controller constants
Kp = 25;
Ki = 50;
Kd = 1;
output = zeros(1, length(v)); % Init
setpoint = 0; % Sets the feedback point to control to.
previous_error = 0; % In this case, you want to set yt to 0.
integral = 0;
xt_filt = [];
yt_filt = [];
for k=51:length(v)
xt(k) = v(k)*cos(2*pi*u_fc*t(k)+output(k));
yt(k) = v(k)*sin(2*pi*u_fc*t(k)+output(k));
xt_filt(k:-1:k-50) = filter(b, a, xt(k:-1:k-50));
yt_filt(k:-1:k-50) = filter(b, a, yt(k:-1:k-50));
if(sign(xt(k)) == sign(yt(k)))
error = setpoint - yt_filt(k);
else
error = setpoint + yt_filt(k);
end
% PID Feedback Controller compensates for inaccurate fc guess.
% error = setpoint - yt_filt(k);
integral = integral + error*10;
derivative = (error - previous_error)/10;
output(k+1) = Kp*error + Ki*integral + Kd*derivative;
previous_error = error;
end
%% Mix and Filter
mix = v.*cos(2*pi*(carrierFreq)*t)';
% mix = v.*cos(2*pi*(fc)*t)';
figure(4)
subplot(2,2,1)
plot(t,mix)
axis([0 2 -1.2 1.2])
title('\fontsize{14}Signal Multiplied by Cosine')
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,2,3)
mixFft = abs(fftshift(fft(mix)));
plot(freqShift,mixFft)
title('\fontsize{14}FFT of Multiplied Signal')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
text(150,.8*max(mixFft),sprintf('Low Frequency\n<- Component'))
text(300,.25*max(mixFft),sprintf('High Frequency\nComponent ->'))
[b,a] = butter(5,(75)/Fs);
mixFilter = filter(b,a,mix);
mixFilterFft = abs(fftshift(fft(mixFilter)));
subplot(2,2,2)
plot(t,mixFilter)
axis([0 2 -1.2 1.2])
title('\fontsize{14}Multiplied Signal Filtered')
xlabel('Time (s)')
ylabel('Amplitude')
subplot(2,2,4)
plot(freqShift,mixFilterFft)
title('\fontsize{14}FFT of Filtered Signal')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
text(150,.8*max(mixFft),sprintf('Low Frequency\n<- Component'))
movegui('west')
%% Sign Change Method
k = 1
for n = 1:length(v)-1
if (sign(mixFilter(n+1)) ~= sign(mixFilter(n)))
zeroLoc(k) = n;
k = k+1;
end
end
period = Fs/baud;
bitstream = []
for n = 1:length(zeroLoc)-1
periodDiff = zeroLoc(n+1)-zeroLoc(n);
cross = round(periodDiff/period) - 1;
if (cross == 0)
bitstream = [bitstream 0];
else if (cross > 0)
bitstream = [bitstream ones(1,cross) 0];
end
end
end
bitstream = [bitstream ones(1,10)]; %pad with ones for end of message
sample = zeros(1,length(v));
sample(zeroLoc) = 1;
figure(5)
plot(t,mixFilter,'b')
hold on
plot(t,sample,'k')
hold off
axis([0 2 -1.2 1.2])
title('\fontsize{14}Sign Change Sync Points')
xlabel('Time (s)')
ylabel('Amplitude')
legend('Filtered Signal','Sync Points','location','best')
%% Bitstream to Characters breakLoc = []; k = 1; for n = 1:length(bitstream)-2 % Finds '00' bits indicating new character if(isequal(find(bitstream(n:n+1) == [0 0]),[1 2])) breakLoc(k) = n; k = k+1; end; end; messageText = ''; for n = 1:length(breakLoc)-1 messageChunks = bitstream(breakLoc(n):breakLoc(n+1)-1); temp = mat2str(messageChunks); temp = strrep(temp,'[',''); temp = strrep(temp,'0 0 ',''); temp = strrep(temp,']',''); if (((breakLoc(n+1) - breakLoc(n)) <= 12) && (~isequal(temp,'0')) &&(bin2dec(fliplr(temp)) <= length(charTable))) messageText = [messageText charTable(bin2dec(fliplr(temp)))]; end; end messageText
Results/Problems
Transmitter
To our knowledge, our transmitter generates a signal without any problems.
Receiver
As long as our guess frequency is close to the actual frequency, (less than 0.1Hz off) our receiver can easily handle a 1dB SNR.
Frequency Errors
In most cases, our guess frequency is within 0.2Hz of the actual frequency. As the guess frequency exceeds these bounds, the text decode becomes increasingly garbled.
PID Controller Problems
The first iteration did not take into account the I component of the signal and was constantly trying to drive the signal to the positive I coordinates. A simple if-else statement solved this problem. However there are some bugs. There might be a singularity at (0,0) which causes the controller to perform in a strange fashion. Also, the code has very bad noise rejection. It can produce baseband but not very well once noise is added. Finally, it doesn't decode signals that we generate ourselves. When we produce our signal in MATLAB, the code works fine. However, the signal from wikipedia.com and Dr. Frohne's signal don't work. One problem with the wikipedia.com signal is that it's native format is .ogg which is lossy. This could've presented problems when we convert the signal to a wave file as high frequencies get chopped in lossy formats.