close all
clear
set(0,'DefaultFigureVisible','on'); % Turn off to suppress all plot windows.

%% Select 'role' of this script.
% 'Both' transmits and records/receives.
% TX and RX transmit only and recieve only, resp.
BOTH_MODE = 0;
TX_MODE   = 1;
RX_MODE   = 2;
MODE = BOTH_MODE;

%% Input parameters
rx_fs = 48000; % Receiver sampling frequency [Sample/s]

Rb = rx_fs/2;   % bit rate [bit/sec]
fs = rx_fs; % Transmitter sampling frequency [Sample/s]
Ts = 1/fs; % Sample time [s/Sample]
fc = 19000;  % Carrier frequency [Hz]

decimation_factor = rx_fs/fs;

span = 6;   % Pulse-width in symbol times.
rolloff = 0.4; % Roll-off factor (ie. alpha) for Raised Cosine pulse.

%% Transmitter
%% Constellation (normalized to unit energy E=d^2)
% Implemented constellations:
% 'bpsk'
% '4qam'
% '8psk'
% '32psk'
% '16qam'
% '256qam'
% '1024qam'
constellation = generateConstellation('64qam');

if length(constellation) > 2 % Voronoi need at least 3 constellation points.
    voronoi(real(constellation), imag(constellation))
    title('TX: Constellation plot')
    xlabel("Phase"); ylabel("Quadrature");
    xlim([-1.5 1.5])
    ylim([-1.5 1.5])
end

% Preamble (including initial delay where nothing is tx'd)
delay_time = randi([1,20], 1);       % Generate a random delay between 1 and 20 samples long.
delay = zeros(1, delay_time);

B2 = [1 -1]; %Barker-2
B3 = [1 1 -1]; %Barker-3
B4 = [1 1 -1 1]; %Barker-4
B5 = [1 1 1 -1 1]; %Barker-5%
B7 = [1 1 1 -1 -1 1 -1]; %Barker-7
B11 = [1 1 1 -1 -1 -1 1 -1 -1 1 -1]; %Barker 11
B13 = [1 1 1 1 1 -1 -1 1 1 -1 1 -1 1]; % Barker-13
preamble = kron(B7, B7);   % Select preamble here. if 1-dimensional barker, set first argument to 1.


%preamble_symbs = preamble.*(1+1i)./sqrt(2); % Map preamble to symbols, so that they are 180deg apart! MUST BE MODIFIED IF USING DIFFERENT CONSTELLATION.
preamble_symbs = preamble; % Do NOT need to map to valid symbols


%% Data source
% Choose data source!
% Simply uncomment the one you want.


%% Generate all unique symbols for QAM configurations
%b = genAllASCII();     %64-QAM
%b = genAllTenBit();    % 1024-QAM
%b = genAllTwelveBit(); % 4096-QAM
%N = length(b);    % Number of bits to transmit

%% Random message of any length.
%N = length(genAllTwelveBit());
%b = randsrc(1,N,[0 1]); % Random message bits.

%% Custom ASCII string
string = 'To be, or not to be: that is the question: Whether  "tis nobler in the mind to suffer The slings and arrows of outrageous fortune, Or to take arms against a sea of troubles, And by opposing end them? To die: to sleep; No more; and by a sleep to say we end The heart-ache and the thousand natural shocks That flesh is heir to,  "tis a consummation Devoutly to be wish "d. To die, to sleep; To sleep: perchance to dream: ay, there "s the rub; For in that sleep of death what dreams may come When we have shuffled off this mortal coil, Must give us pause: there "s the respect That makes calamity of so long life; For who would bear the whips and scorns of time, The oppressor "s wrong, the proud man "s contumely, The pangs of despised love, the law "s delay, The insolence of office and the spurns That patient merit of the unworthy takes, When he himself might his quietus make With a bare bodkin? who would fardels bear, To grunt and sweat under a weary life, But that the dread of something after death, The undiscover "d country from whose bourn No traveller returns, puzzles the will And makes us rather bear those ills we have Than fly to others that we know not of? Thus conscience does make cowards of us all; And thus the native hue of resolution Is sicklied o "er with the pale cast of thought, And enterprises of great pith and moment With this regard their currents turn awry, And lose the name of action.--Soft you now! The fair Ophelia! Nymph, in thy orisons Be all my sins remember "d!';
%string = 'SSY121 Group 7 says: The capital of Hawaii is Honolulu.';

testImage = imread('lego-johnny-thunder-minifigure-24_small.jpg');
string = reshape(testImage, 1, []);


b = string2bitStream(string);
%disp(['Sent message: ' string])

N = length(b);
if length(b) > N
    disp('Source too long! Truncating...')
    b = b(1:N);
end
b = [b zeros(1, N-length(b))]; % Pad string with zeros to fill length 432 packet.


%% Bits to Messages
M = length(constellation);      % M unique messages/symbols.
bpsymb = log2(M);   % Bits per symbol
Rsymb = Rb/bpsymb;  % Symbol rate [baud]
Tsymb = 1/Rsymb;    % Symbol time [s]
sps   = fs/Rsymb;   % Samples per symbol. [Sample/s / Symbol/s = Sample/Symbol] 
if floor(sps) ~= sps % Warn if sps is not integer, as it will probably cause issues.
    disp('WARNING: sps not an integer!');
    disp(sps);
    sps = fix(sps); % Convert explicitly to integer !!UGLY HACK!!
end

m = buffer(b, bpsymb)';          % Group information bits into bits per symbol
m_idx = bi2de(m, 'left-msb')'+1; % Bits to symbol index, msb: the Most Significant Bit. The +1 is unexplained...

%% Message to symbol
message_symbs = constellation(m_idx);           % Map messages to constellation symbols using indices
symbs         = [preamble_symbs message_symbs]; % Total symbol vector including preamble, message etc...

% Generate basic pulse
[pulse, pulse_t] = rtrcpuls(rolloff, Tsymb, fs, span); % Generate basis pulse.
figure
plot(pulse_t, pulse); grid on; title('Basic pulse');

% Pulse-modulation -> create baseband signal
x_upsample = upsample(symbs, sps);      % Upsample symbols from Rsymb to fs.
s          = conv(pulse, x_upsample);   % Pulse shaping the symbol-train by convolving with basis pulse. This gives baseband signal.
t_s        = (0:length(s)-1).*Ts;       % Signal time-axis ie. x[n]*n*Ts.

% Normalize baseband amplitude to 1 (because pulse-convolution attenuates signal).
s = s.*(1/max(abs(s)));

% Plot baseband signal and symbol-train.
figure
subplot(221) 
plot(real(s))
grid on; title('Real pulse-modulated baseband');

subplot(223)
stem(real(x_upsample))
grid on; title('Real symbol baseband');

subplot(222) 
plot(imag(s))
grid on; title('Imag pulse-modulated baseband');

subplot(224)
stem(imag(x_upsample))
grid on; title('Imag symbol baseband');

%% Upmix with carrier.
I = real(s);
Q = imag(s);
tx_sig = I.*cos(2*pi*fc*(t_s)) - Q.*sin(2*pi*fc*(t_s));

% Plot carrier-modulated signal.
figure
plot(t_s, tx_sig); 
title('Carrier signal')
grid on;

%% Channel
if MODE == RX_MODE || MODE == BOTH_MODE 
    %% Record!
    my_sound = audiorecorder(fs,24,1);
    record(my_sound);
    if MODE == RX_MODE
        pause(10)
    end
end

% Play carrier wave in speakers!
% Create an audioplayer object to play the noise at a given fs.
if MODE == TX_MODE || MODE == BOTH_MODE 
    audio_player = audioplayer(tx_sig, fs);  
    playblocking(audio_player);   %Play the noise. Why blocking?
    stop(my_sound);
end

if MODE == RX_MODE || MODE == BOTH_MODE 
    figure
    plot(my_sound.getaudiodata); title('Recorded Audio');

    s_passband = my_sound.getaudiodata.';

    %% Simulate oversampled signal by interpolating
    % tx_sig_interp = upsample(tx_sig, decimation_factor);
    % t_s = (0:length(s)-1).*Ts; % Signal time-axis ie. x[n]*n*Ts.
    % figure 
    % subplot(211)
    % plot(t_s, tx_sig)
    % subplot(212)
    % t_s_upsample = (0:length(tx_sig_interp)-1).*(1/rx_fs);
    % plot(t_s_upsample, tx_sig_interp)

    % %% Channel (AWGN)
    % SNRdB = 100;
    % s_passband = awgn(tx_sig_interp, SNRdB, 'measured');
    % figure; 
    % plot(s_passband)
    % title('Carrier signal with Noise')

    %% Receiver
    phzOffs = 0; % DEBUG: insert phase offset.

    % Downmix with carrier.
    t_s = (0:length(s_passband)-1).*(1/rx_fs); % Redefine time axis for our received (oversampled) signal % Signal time-axis ie. x[n]*n*Ts.
    rx_I =  s_passband.*cos(2*pi*fc*t_s + phzOffs);
    rx_Q = -s_passband.*sin(2*pi*fc*t_s + phzOffs);

    % Plot spectrum of passband signal.
    figure
    plot_fft(s_passband, rx_fs);
    title('spectrum of passband signal')
    
    % Plot downmixed signal.
    figure
    subplot(121)
    plot(t_s, rx_I')
    title('RX-I, downmixed')
    subplot(122)
    plot(t_s, rx_Q)
    title('RX-Q, downmixed')

    % Lowpass filter to give baseband signal.
    rx_I_filtered = lowpass(rx_I, Rsymb, rx_fs, 'ImpulseResponse','iir','Steepness',0.75);
    rx_Q_filtered = lowpass(rx_Q, Rsymb, rx_fs, 'ImpulseResponse','iir','Steepness',0.75);
   
    % Plot LP-filtered bandpass signal.
    figure
    subplot(121)
    plot(t_s, rx_I_filtered)
    title('RX-I, filtered')
    subplot(122)
    plot(t_s, rx_Q_filtered)
    title('RX-Q, filtered')

    rx_I_filtered_down = downsample(rx_I_filtered, decimation_factor);
    rx_Q_filtered_down = downsample(rx_Q_filtered, decimation_factor);



    % Matched filter
    MF = fliplr(conj(pulse));        %create matched filter impulse response
    MF_I_output = filter(MF, 1, rx_I_filtered_down); 
    MF_I_output = MF_I_output(length(MF):end); %remove transient
    MF_Q_output = filter(MF, 1, rx_Q_filtered_down); 
    MF_Q_output = MF_Q_output(length(MF):end); %remove transient

    % Plot matched-filtered signal.
    figure
    subplot(121)
    plot(MF_I_output)
    hold on
    grid on; title('RX-I, matched-filtered')
    subplot(122)
    plot(MF_Q_output)
    grid on; title('RX-Q, matched-filtered')

    %% Symbol Synchronizer (based on early-late gate)
    avg = zeros(sps, 1);
    MF_output_complex = MF_I_output + 1i.*MF_Q_output;
    
%         figure
%         hold off
    for delta = 1:sps
        shifted_MF_output = circshift(MF_output_complex, delta-1);          % Shift signal by one sample.
        shifted_samples = shifted_MF_output(1:sps:end);                     % Resample shifted signal at symbol-times.
        avg(delta) = sum(abs(shifted_samples))./length(shifted_samples);    % Check average energy (max at correct symbol-times).

        % Visual animation of constellation points converging (if
        % synchronisation possible)
% 
%         plot(shifted_samples, 'k*')
%         xlabel('In Phase'); ylabel('Quadrature');
%         title(['Sync offset: ',  num2str(delta-1)])
%         pause(0.1)
    end

    [max_avg, max_idx] = max(avg);  % Max average energy gives correct symbol-time (and corresponding signal time-shift).
    sync_offset = max_idx-1;        

    % Plot average energy iterated over one symbol-time.
    figure
    plot(0:sps-1,avg)
    xlabel('sync offset [samples]'); ylabel('Average energy');

    MF_output_shifted = circshift(MF_output_complex, sync_offset);
    rx_symbs = MF_output_shifted(1:sps:end);
    %% End of synchronizer


    % Plot sample points on baseband signal.
    figure
    subplot(121)
    plot(real(MF_output_shifted))
    hold on
    stem((0:length(rx_symbs)-1)*sps, real(rx_symbs))
    grid on; title('RX-I, sampled')
    subplot(122)
    plot(imag(MF_output_shifted))
    hold on
    stem((0:length(rx_symbs)-1)*sps, imag(rx_symbs))
    grid on; title('RX-Q, sampled')
    hold off

    % Unmapping to symbols
    samples = rx_symbs.'; % OBS! Transposing complex values requires dot operator!! ie. cmplx.'
    
    % Plot received samples on decision regions.
    figure
    hold on
    if length(constellation)>2
        voronoi(real(constellation), imag(constellation))
    end
    plot(samples, 'o')
    hold off
    axis square; grid on
     xlim([-1 1]);
     ylim([-1 1]);
    title('Received samples');

    % Preamble correlation
    sample_corr = conv(samples, fliplr(preamble));
    [~,peak_idx] = max(sample_corr);
    
    
    figure
    subplot(311)
    plot(abs(sample_corr))
    title(['Preamble correlation: Abs.'])
    hold on; 
    plot(peak_idx, abs(sample_corr(peak_idx)), 'or') % Plot peak value.
    grid on
    subplot(312)
    plot( real(sample_corr) )
    grid on
    title('Preamble correlation: I')
    subplot(313)
    plot(imag(sample_corr))
    grid on
    title('Preamble correlation: Q')
    
    % Get ALL received preamble symbols (based on preamble detection). 
    % eg. If B13 is used, this vector will contain all 13 symbols.
    rx_preamble_symbs =  samples(peak_idx-(length(preamble)-1):peak_idx);
    % Find average angle.
    angles = angle(rx_preamble_symbs);  % Compute each individual angle.
    neg_angles = angles < 0 ;           % Identify indicies of neg angles...
    angles(neg_angles) =angles(neg_angles)+pi;  % ...phase shift 180deg.
    rx_preamble_avg_angle = mean(angles);       % Compute average angle overall.
    
    % Phase error still ambiguous, extract coarse phase error from preamble
    % peak and adjust accurate phase error (above) to lie in correct
    % quadrant.
    % Attempt to use preamble to find phase errors/ambiguities.
    % eg. if there is a phase error of 90deg, there will be a correlation, but
    % negative in one of the I/Q streams.
    if angle(samples(peak_idx)) < 0
        rx_preamble_avg_angle = rx_preamble_avg_angle+pi;
    end
    
    % Plot phase error angle. 
%     figure
%     plot([exp(1i*rx_preamble_avg_angle), 1*exp(1i*(rx_preamble_avg_angle+pi))])
%     title('Phase offset')
    
    phzError = rx_preamble_avg_angle + 0;               % Finally, adjust phase error by nominal rotation of preamble (0 in this case, since it is real).
    samples_phzCorrected = samples.*exp(1i*-phzError);  % Compensate by counterrotating samples.

    % Normalize to average signal energy.
    samples_phzCorrected = samples_phzCorrected./( sum(abs(samples))./length(samples)); 

    % Normalize to average preamble symbol energy
    avg_preamble_energy = mean(abs(samples_phzCorrected(peak_idx-(length(preamble)-1):peak_idx)));
    samples_phzCorrected = samples_phzCorrected.*(1./avg_preamble_energy);
    
    % Extract message symbols and discard everything else.
    rx_message_symbs = samples_phzCorrected(peak_idx+1:peak_idx+ceil(N/bpsymb));
    
    % Plot phase-corrected samples.
    figure(100)
    hold on
    plot(rx_message_symbs, 'ob'); title('Phase-corrected rx. Decision regions overlaid.')
    if length(constellation)>2
        voronoi(real(constellation), imag(constellation))
    end
    grid on; axis square;
    xlim([-2 2])
    ylim([-2 2])
    hold off
    
    % Minimum euclidean distance. Unmap symbols to messages
    d = abs(repmat(rx_message_symbs, 1, length(constellation)) - repmat(constellation, length(rx_message_symbs), 1)).^2; % Compute distance of every sample symbol to every constellation point.
    [~, nearest_idx] = min(d, [], 2);  % Get index of value with minimum distance to a constellation point. Traverse along column.

    rx_mapped = constellation(nearest_idx); % Choose corresponding symbol.



    % Symbols to messages
    m_hat = rx_mapped; 
    if (length(m_hat) < length(message_symbs)) % Check if we've received at least enough symbols for 1 message.
        disp('Not enough symbols read! Check preamble detection.');
        return
    end
    m_hat = m_hat(1:round(N/bpsymb));

    [SER, SER_ratio, SER_loc] = symerr(message_symbs, m_hat);
    
    % If symbol errors detected, plot locations of offending symbols
    if SER > 0
        figure(100)
        hold on
        if length(constellation)>2
            voronoi(real(constellation), imag(constellation))
        end
        hold on
        
        rx_message_symbs_error_idx = zeros(1, length(SER));
        jj = 1;
        for idx = 1:length(rx_message_symbs)
            if SER_loc(idx) == [1]
                rx_message_symbs_error_idx(jj) = rx_message_symbs(idx);
                jj = jj + 1;
            end
        end
        
        plot(rx_message_symbs_error_idx, '*r')
        title('Incorrect Symbols')
        axis square
    end


    % Messages to bits
    b_hat_buffer = de2bi(nearest_idx-1, bpsymb, 'left-msb')';
    b_hat = b_hat_buffer(:)'; %write as a vector
    b_hat = b_hat(1:N);
    [BER, BER_ratio] = biterr(b, b_hat); %count of bit errors

    % b_hat_trunc = b_hat(1:end-rem(length(b_hat), 8)); % Truncate received bitstream so that it's divisible by 7 (ASCII length).
    disp(['Received message (quote marks added): "' bitStream2string(b_hat) '"'])

    disp('RX_symbs, RX_bits')
    disp([length(m_hat), length(b_hat)])
    
    disp('SER, SER_ratio')
    disp([SER, SER_ratio*100])
    
    disp('BER, BER_ratio')
    disp([BER, BER_ratio*100])
    
    % Plot power of passband signal
    % First peak should be at -fc: our baseband signal(?).
%     figure
%     [pxx, f] = pwelch(rx_I + rx_Q*1i ,1024,768,1024, rx_fs); % note that pwr_spect.f will be normalized frequencies
%     f = fftshift(f); %shift to be centered around fs
%     f(1:length(f)/2) = f(1:length(f)/2) - rx_fs; % center to be around zero
%     p = fftshift(10*log10(pxx/max(pxx))); % shift, normalize and convert PSD to dB
%     plot(f,p)
%     xlim([-400 400])
    

    %% Eye diagram
    % WARNING! Heavy Computation required.
%     eye_vec = MF_output_shifted+i*MF_output_shifted;
%     eye_vec = eye_vec.*exp(i*-phzError);
%     eyediagram(eye_vec, sps, Tsymb)

% eyediagram(s, sps,Tsymb)
    
    decode_bmp(b_hat, testImage)
        

    figure
    binscatter( real(rx_message_symbs), imag(rx_message_symbs))
    title('Received constellation point density')
end
