%EE107: Final Project
%Elsa Ames and Bennett Brain

%all figures are commented out for runtime, see finalproj_noloop to do
%figures as it goes

%set these to pick equalizer and wave

wave = 1; %1 for SRRC, 0 for halfwave
equa = 1; %1 for MMSE, 0 for zero-forcing

close all
noisePow = .0001; %changeable parameter
%% Image Pre-Processing
image = imread('pomodoro.png');
image = rgb2gray(image);
image = im2double(image);
[m, n] = size(image);
fun = @dct2;
dct_image = blkproc(image, [8 8], fun);

%Normalize Image
zeroOffset = min(min(dct_image));
dct_image = dct_image - zeroOffset;
scaleFactor = max(max(dct_image));
dct_image = dct_image / scaleFactor;


%quantize to ints
levels = linspace(0,1,255);
B = imquantize(dct_image, levels);
B = reshape(B,[8 8 (m*n/64)]);

if (wave == 0)
hwav_B_recomp = zeros(size(B));
else
srrc_B_recomp = zeros(size(B));
end

%% make SRRC and hwav pulses to start to improve runtime
T = 1;
t = linspace(0,1,33);
A1 = 1; %amplitude



halfpulse = A1 * sin(pi*t/T);

%recalculate A1, halfpulse
A1 = 1/sqrt(sum(halfpulse.^2)); %normalize power level
halfpulse = A1*sin(pi*t/T);
halfpulse = halfpulse(1:end-1);

a = .5;
T = 1;
K = 6;

b = rcosdesign(a, K*T, 32); %already is normalized
b = b';
% figure
% plot(b)
rcosL = length(b);

%% channel Part 1

hnpad = zeros(1,31);
%initial
hn = [1,hnpad,1/2,hnpad,3/4,hnpad,-2/7,hnpad];

%indoor
%hn = [1,hnpad,.4365,hnpad,.1905,hnpad,.0832,hnpad,0,hnpad,.0158,hnpad,0,hnpad,.003,hnpad];

%outdoor
%hn = [.5,hnpad,1,hnpad,0,hnpad,.63,hnpad,zeros(1,32*4),.25,hnpad,zeros(1,32*3),.16,hnpad,zeros(1,32*12),.1,hnpad];

Hejw = fft(hn,128*64);

%% MMSE equalizer Part 1


Eb = 1; %because all signals were normalized
if (equa == 1)
Qmmse = conj(Hejw)./(abs(Hejw).^2+sigma.^2/Eb);
qnmmse = ifft(Qmmse);
else
%% Zero-Forcing equalizer part 1
Qzf = 1./Hejw;
qnzf = ifft(Qzf);
end


%%

%Regroup into N blocks and reshape
N = 6;

%% START LOOP HERE

%Note, N must divide m*n/64 without remainder for this to work well
% for our image, m*n/64 is 33750, 3 does divide it
for z = 0:N:(m*n/64)-N

currentBlocks = B(:,:,z+1:z+N); %current block will be 8 x 8 x N

dct_col = reshape(currentBlocks, 1, []);
dct_col_binary = de2bi(dct_col', 8); %this is our bitstream
bitstream = reshape(dct_col_binary', 1, []); % 2.3 conversion to a bitstream
%% Half wave

if (wave == 0)
half_wave = zeros(32*length(bitstream),1);

for k = 1:length(bitstream)
     if(bitstream(k) == 0)
        half_wave(k*32-31:k*32) = -halfpulse;
     end

     if(bitstream(k) == 1)
         half_wave(k*32-31:k*32) = halfpulse;
     end
end

%{
figure
plot(t(1:end-1), halfpulse)

figure
plot(half_wave)
xlim([1, 32*10]);
%}

else

%% SRRC

bitstream2 = bitstream*2-1;
bitstreampad = padarray(bitstream2,31,'post');
bitstreampad = reshape(bitstreampad,1,[]);

rais_cos = conv(bitstreampad,b);

%{
%testing signals
figure
plot(rais_cos)
xlim([1,600])
hold on
plot([32,32],[-.05,.2],'b')
plot([32*2,32*2],[-.05,.2], 'b')
plot([32*3,32*3],[-.05,.2], 'b')
plot([32*4,32*4],[-.05,.2], 'b')
plot([32*5,32*5],[-.05,.2], 'b')
plot([32*6,32*6],[-.05,.2], 'b')
plot([32*7,32*7],[-.05,.2], 'b')
plot([32*8,32*8],[-.05,.2], 'b')
plot([32*9,32*9],[-.05,.2], 'b')
hold off
%}

%{
eyediagram(half_wave,32,1,16)
title('eye diagram half wave mod')
eyediagram(rais_cos,32)
title('eye diagram srrc mod')
%}

end

%% channel

%hnpad = zeros(1,31);
%hn = [1,hnpad,1/2,hnpad,3/4,hnpad,-2/7,hnpad];

%{
figure
freqz(hn); %this might be wrong
%}

if(wave == 0)
halfwav_channel = conv(half_wave,hn);
else
srrc_channel = conv(rais_cos,hn);
end

%{
figure
plot(srrc_channel)
xlim([1,300])

eyediagram(halfwav_channel,32,1,16)
eyediagram(srrc_channel,32)
%}

%% adding noise

sigma = sqrt(noisePow);
if(wave == 0)
noisemat_hwav = sigma*randn(size(halfwav_channel));
hwav_noise = halfwav_channel + noisemat_hwav;
else
noisemat_srrc = sigma*randn(size(srrc_channel));
srrc_noise = srrc_channel + noisemat_srrc;
end

%{
eyediagram(hwav_noise,32,1,16)
eyediagram(srrc_noise,32)
%}

%%
%Matched filtering

if(wave == 0)
hwav_match = flip(halfpulse);
hwav_mfilt = conv(hwav_match,hwav_noise);
else
srrc_match = flip(b);
srrc_mfilt = conv(srrc_match,srrc_noise);
end

%{
figure
plot(hwav_match)
title('Impulse response half wave matched filter')
figure
freqz(hwav_match)
title('Freq response half wave matched filter')

figure
plot(srrc_match)
title('Impulse response srrc matched filter')
figure
freqz(srrc_match)
title('Freq response srrc matched filter')
%}



%{
eyediagram(hwav_mfilt,32)
eyediagram(srrc_mfilt,32)
%}

%% Zero-Forcing equalizer
% (almost) entirely commented out because MMSE equalizer is used instead
%see other code for results of ZF equalizer

%{
figure
freqz(qnzf)
%}
if(equa == 0)
    if(wave == 0)
hwav_equa = conv(qnzf,hwav_mfilt);
    else
srrc_equa = conv(qnzf,srrc_mfilt);
    end

%{
eyediagram(hwav_eq,32)
eyediagram(srrc_eq,32)
%}

else

%% MMSE equalizer

%Eb = 1; %because all signals were normalized

%Qmmse = conj(Hejw)./(abs(Hejw).^2+sigma.^2/Eb);

%qnmmse = ifft(Qmmse);

    if(wave == 0)
hwav_equa = conv(qnmmse,hwav_mfilt);
    else
srrc_equa = conv(qnmmse,srrc_mfilt);
    end


end
%{
eyediagram(hwav_mmse,32)
eyediagram(srrc_mmse,32)
%}

%% Sampling & Detection

if(wave == 0)

hwav_samplStart = 32;
hwav_samplRate = 32;
hwav_bitRecomp = zeros(size(bitstream));
for p = 0:length(bitstream)-1
    if (hwav_equa(hwav_samplStart+p*hwav_samplRate) < 0)
        hwav_bitRecomp(p+1) = 0;
    else
        hwav_bitRecomp(p+1) = 1;
    end
end

else

srrc_samplStart = 96+96;
srrc_samplRate = 32;
srrc_bitRecomp = zeros(size(bitstream));

for p = 0:length(bitstream)-1
    if (srrc_equa(srrc_samplStart+p*srrc_samplRate) < 0)
        srrc_bitRecomp(p+1) = 0;
    else
        srrc_bitRecomp(p+1) = 1;
    end

end

end

%% conversion back to image

if(wave == 0)
hwav_binary_recomp = reshape(hwav_bitRecomp,8,[]);
hwav_binary_recomp = hwav_binary_recomp';
hwav_col_recomp = bi2de(hwav_binary_recomp);
hwav_col_recomp = hwav_col_recomp';
hwav_cblock_recomp = reshape(hwav_col_recomp,8,8,[]);
hwav_B_recomp(:,:,z+1:z+N) = hwav_cblock_recomp;

else
srrc_binary_recomp = reshape(srrc_bitRecomp,8,[]);
srrc_binary_recomp = srrc_binary_recomp';
srrc_col_recomp = bi2de(srrc_binary_recomp);
srrc_col_recomp = srrc_col_recomp';
srrc_cblock_recomp = reshape(srrc_col_recomp,8,8,[]);
srrc_B_recomp(:,:,z+1:z+N) = srrc_cblock_recomp;
end

%big for loop ends here
end

fun2 = @idct2;

if(wave == 0)
hwav_B_reshape = reshape(hwav_B_recomp, m,n);
hwav_image_dct = rescale(hwav_B_reshape);
hwav_image_dct = hwav_image_dct * scaleFactor;
hwav_image_dct = hwav_image_dct + zeroOffset;
hwav_image = blkproc(hwav_image_dct, [8,8], fun2);
hwav_im = mat2gray(hwav_image);
figure
imshow(hwav_im)

else
srrc_B_reshape = reshape(srrc_B_recomp, m,n);
srrc_image_dct = rescale(srrc_B_reshape);
srrc_image_dct = srrc_image_dct * scaleFactor;
srrc_image_dct = srrc_image_dct + zeroOffset;
srrc_image = blkproc(srrc_image_dct, [8,8], fun2);
srrc_im = mat2gray(srrc_image);
figure
imshow(srrc_im)
end