%K Ayob, 02/Nov/2009
%binary FSK modulator/demodulator model

%summary of current set of normalised values:
%message bandwidth(bw) = Fs * .5/up
%f1 = Fs * inc1/256
%f2 = Fs * inc2/256
%modulated bw = f2-f1+2*bw

clear all;

%parameters
Fs = 1;                                %sampling frequency at upsamped stages
N = 100;                               %number of test symbols
up = 7;                                %upsample factor, symbol baseband(B) becomes Fs*.5/up 
lut = 256;                             %sine lookup table size, freq = inc/lut
inc1 = 10;                             %1st modulation frequency f1 = Fs*inc1/lut
inc2 = inc1+ inc1*200/31.25;           %2nd modulation frequency f2 = Fs*inc2/lut 
%C = 100;                              %carrier frequency = Fs*C/lut, 
%SNR = 50;                             %dB signal to noise ratio
sc = 2^15 - 1;                         %hardware scale factor

use_shaping = 0;                       %pulse shaping moves band up ??
 
%******************* Tx ********************

%generate random symbols(+/-32767), Fs=Fsym=1
symb = randsrc(1,N,[+sc -sc]);

%shaping and upsampling to Fs=Fsym*up
%message BW becomes .5/up 
 
%h = firrcos(60*up,.5/up,.3,1,'rolloff','sqrt');
%symbup = upfirdn(symb,h,up,1);
 
symbup = resample(symb,up,1);
 
%----------- binary FSK modulation ------------
%convert message amplitude to phase increments 
if use_shaping == 1
    incr = symbup + abs(min(symbup));
    incr = inc2*incr/max(incr);
    incr = incr + inc1;
else
    for i = 1:length(symbup)
       if symbup(i) > 0 incr(i) = inc2; else incr(i) = inc1; end
    end
end

%DDS structure for FSK modulation 
%based on one cycle sine/cosine data
%-----------------------------------

lut_data = exp(j*2*pi*[0:lut-1]*1/lut);

ptr = 1;   
for n = 1:length(symbup) 
    x(n) = lut_data(ptr);
    ptr = round(ptr + incr(n)); 
    if ptr > lut  ptr = ptr - lut; end;   %modulo addition
end

%note: x is normalised and independant of above symb scaling

%check baseband and modulated message spectrum
%message bw = Fs*.5/256
%modulated message width becomes [f1-f2+2*bw]
figure;hold
f = linspace(-Fs*up/2,Fs*up/2,length(x));
plot(f,20*log10(abs(fftshift(fft(symbup/32767)))),'k');
plot(f,20*log10(abs(fftshift(fft(x)))));

%-----------------------------------------
%upconvert 
%carrier = exp(j*2*pi*[0:length(x)-1]*C/lut);
%x_up = real(carrier .* x);

%check upconverted signal spectrum
%f = linspace(0,Fs,length(x_up));
%plot(f,20*log10(abs(fft(x_up))),'r');

%************************************
%channel white noise
%noise_ratio = 1/10^(SNR/20);
%x_up = x_up + noise_ratio*rand(1,length(x_up));

%**************** Rx ****************

%downconvert into I/Q
%xr = x_up .* real(carrier);
%xi = x_up .* -imag(carrier);

%downconvert into I/Q
carrier = exp(j*2*pi*[0:length(x)-1]* -(inc1+.5*(inc2-inc1))/lut);
x_rx = carrier .* x;
xr = real(x_rx);
xi = imag(x_rx);

%filter off the high sideband
%h = fir1(40,200/lut);
%xr = filtfilt(h,1,xr);
%xi = filtfilt(h,1,xi);

%check downconverted signal spectrum
plot(f,20*log10(abs(fftshift(fft(complex(xr,xi))))),'g--');
legend('baseband','FSK message','downconverted');

I = round(xr*sc);
Q = round(xi*sc);

%demodulate, get phase and convert to frequency
%method 1, using arctan
theta1 = atan2(Q,I);

%method 2
for i = 1:length(I)
    if I(i) >= 0 && Q(i) >= 0                             %quadrant I
        r = round((I(i)-Q(i)))/round((I(i)+Q(i)));
        theta2(i) = pi/4 - r*pi/4;
    elseif I(i) > 0 && Q(i) < 0                           %quadrant II
        r = round((I(i)-Q(i)))/round((I(i)+Q(i)));
        theta2(i) = -(pi/4 - r*pi/4);
    elseif I(i) < 0 && Q(i) > 0                           %quadrant III
        r = round((I(i)+Q(i)))/round((Q(i)-I(i)));
        theta2(i) = 3*pi/4 - r*pi/4;
    elseif I(i) < 0 && Q(i) < 0                           %quadrant IV
        r = round((I(i)+Q(i)))/round((Q(i)-I(i)));
        theta2(i) = -(3*pi/4 - r*pi/4);
    end;
end

theta = theta1;
%theta = theta2;

%convert phase to frequency
m = diff(unwrap(theta));
f = m * Fs/(2*pi);

%check recovered frequencies and symbols  
figure; 
plot(incr*Fs/256);hold %original frequencies
plot(f,'r--');         %recovered frequencies
legend('original frequecies','recovered frequencies');

figure;
plot(symbup);

% decimate back and slice 
m_dec = resample(m,1,up);
for i=1:length(m_dec)
    if m_dec(i) > mean(m_dec)   symb_rx(i) = sc; else symb_rx(i) = -sc; end;
end

error = symb-symb_rx;

%save data
%fid = fopen('I1.txt','w');
%for i = 1:length(I)
%    fprintf(fid,'%i\r',I(i));
%end
%fclose(fid);

%fid = fopen('Q1.txt','w');
%for i = 1:length(Q)
%    fprintf(fid,'%i\r',Q(i));
%end
%fclose(fid);

%fid = fopen('I2.txt','w');
%for i = 1:length(I)
%    fprintf(fid,'%i\r',round(real(x(i))*sc));
%end
%fclose(fid);

%fid = fopen('Q2.txt','w');
%for i = 1:length(Q)
%    fprintf(fid,'%i\r',round(imag(x(i))*sc)); 
%end
%fclose(fid);

%end of program