clear
close all

Im = imread('drno.jpg');
res = [480,720];
Im = imresize(Im,res);
N = 10;
fH = 15625*1000/1001;
fs = fH*(res(2)+5*N);
Tl = res(2)+5*N;

Y = double( 0.3*Im(:,:,1) + 0.59*Im(:,:,2) + 0.11*Im(:,:,3) );

Cr = double( Im(:,:,1) ) - Y ;
Cb = double( Im(:,:,3) ) - Y ;

Y =  [ zeros(size(Y,1),N), -0.3*ones(size(Y,1),N), zeros(size(Y,1),3*N), Y/255  ];
Cb = [ zeros(size(Y,1),5*N), Cb/255 ];
Cr = [ zeros(size(Y,1),5*N), Cr/255 ]; 
Y = Y';   Y = Y(:)';
Cb = Cb'; Cb = Cb(:)';
Cr = Cr'; Cr = Cr(:)';
%
%%

[bY,aY] = butter(5,4.5e6/fs);
[bC,aC] = butter(5,1e6/fs);
Y  = filter(bY,aY,Y);
Cb = filter(bC,aC,Cb);
Cr = filter(bC,aC,Cr);

fc = 455*fH/2;
N0 = length(Cb);

t = (1:N0)/fs;
chroma_mod = Cb/2.03.*cos(2*pi*t*fc) + Cr/1.14.*sin(2*pi*t*fc);

%introduce phase error
chroma_mod = real( hilbert(chroma_mod).*exp(1i*pi/2) ); 

CVBS = Y + chroma_mod;
%%
figure
plot((1:N0)/N0*fs,20*log10(abs(fft(Y))));
xlim([0,fs/2]);
hold on
plot((1:N0)/N0*fs,20*log10(abs(fft(chroma_mod))));
xlim([0,fs/2])
xlabel('f -> [Hz]')
legend('Luma(f)','Chroma(f)')
%%
figure
subplot(1,2,1)
plot(CVBS(1:4*res(2))/255)
subplot(1,2,2)
plot((1:N0)/N0*fs,20*log10(abs(fft(CVBS))));
xlim([0,fs/2])
%% Luma-chroma separation

% No separation
% Y_demod = CVBS;
% chroma_demod = -CVBS;

%Lowpass/highpass pair
% [b_lp, a_lp] = butter(3,2.5e6/fs*2);
% [b_hp, a_hp] = butter(3,2.5e6/fs*2,'high');
% Y_demod = filtfilt(b_lp, a_lp,CVBS);
% chroma_demod = -filtfilt(b_hp,a_hp,CVBS);

%Noth/bandpass pair
[b_notch, a_notch] = butter(3,(fc+[-1e6,1e6])/fs*2,'stop');
[b_bp, a_bp] = butter(3,(fc+[-1e6,1e6])/fs*2,'bandpass');
Y_demod = filtfilt(b_notch, a_notch,CVBS);
chroma_demod = -filtfilt(b_bp,a_bp,CVBS);

%Ideal comb filter
Y_demod = (CVBS + circshift(CVBS,Tl))/2;
chroma_demod = (-CVBS + circshift(CVBS,Tl))/2;

figure
plot((1:N0)/N0*fs,20*log10(abs(fft(Y_demod))));
%xlim([0,fs/2]);
hold on
plot((1:N0)/N0*fs,20*log10(abs(fft(chroma_demod))));
title('Separated components')
xlim([0,fs/2])
xlabel('f -> [Hz]')
legend('Luma(f)','Chroma(f)')
%%
% Demodulate QAM chroma
Cb_demod = -2.03*2*circshift(filter(bC,aC,chroma_demod.*cos(2*pi*t*fc)),-25);
Cr_demod = -1.14*2*circshift(filter(bC,aC,chroma_demod.*sin(2*pi*t*fc)),-25);

% Calculate RGB signals
Y_demod  = reshape(Y_demod',Tl,res(1))';
Cb_demod = reshape(Cb_demod',Tl,res(1))';
Cr_demod = reshape(Cr_demod',Tl,res(1))';

R = Cr_demod + Y_demod;
B = Cb_demod + Y_demod;
G = (Y_demod - 0.3*R - 0.11*B)/0.59;

R = uint8(double(255*R(:,5*N:end)));
G = uint8(double(255*G(:,5*N:end)));
B = uint8(double(255*B(:,5*N:end)));

Im_out = uint8(zeros(size(R,1),size(R,2),3));
Im_out(:,:,1) = R;
Im_out(:,:,2) = G;
Im_out(:,:,3) = B;

%%

figure
subplot(1,2,1)
imshow(Im)
subplot(1,2,2)
imshow(Im_out)