%% EE 351 - Week 7 (optional): television video chain in plain MATLAB
% From an image row to a sampled video waveform, its spectrum, and DSB
% versus VSB amplitude modulation. No toolboxes, no Simulink.
%
%   1. Turn one image row into a video line (sync, back porch, active video,
%      front porch) sampled at fs.
%   2. Check f = N / (2 T_active) against the FFT of a stripe line.
%   3. Amplitude-modulate the video (DSB) and shape it with a vestigial
%      filter (VSB); compare the recovered video.
%
% The RF carrier is scaled down to 8 MHz so that everything fits inside the
% sampling rate. Real NTSC carriers are in the VHF/UHF bands (54-806 MHz in
% the textbook's chart), but the baseband mathematics is identical.

clear; close all; clc;

%% Parameters (NTSC System M, colour)
fs   = 40e6;                       % sample rate, Hz
fH   = 4.5e6/286;                  % line frequency, Hz (15 734.27)
TH   = 1/fH;                       % line period, s (63.556 us)
Tsync = 4.7e-6;  Tbp = 4.7e-6;  Tfp = 1.5e-6;
Tact = TH - Tsync - Tbp - Tfp;     % active video, s (52.66 us)
B    = 4.2e6;                      % nominal video bandwidth, Hz
fc   = 8e6;                        % scaled-down "picture carrier", Hz
fv   = 0.75e6;                     % vestige width, Hz

nLine = round(TH*fs);              % samples per line
t = (0:nLine-1)/fs;
iSync = t < Tsync;
iAct  = t >= Tsync+Tbp & t < Tsync+Tbp+Tact;
xAct  = (t(iAct) - (Tsync+Tbp))/Tact;          % 0..1 across the active line

%% 1. Image row -> video line (levels in IRE: sync -40, blank 0, black 7.5, white 100)
Hh = 48; Ww = 64;
img = ones(Hh, Ww);                                % white background
img(8:40, 20:27)  = 0.04;                          % letter F: stem
img(8:15, 20:46)  = 0.04;                          % top bar
img(22:29, 20:40) = 0.04;                          % middle bar
row = 10;                                          % row crossing the top bar

makeLine = @(rowvals) buildLine(rowvals, nLine, iSync, iAct, xAct);
vF = makeLine(img(row, :));

figure('Name','Image row to video line');
subplot(1,2,1); imagesc(img); colormap(gray); axis image; hold on;
plot([0.5 Ww+0.5], [row row], 'r-', 'LineWidth', 1.5); title(sprintf('Scanned row %d', row));
subplot(1,2,2); plot(t*1e6, vF, 'LineWidth', 1.4); grid on;
xlabel('time within one line (\mus)'); ylabel('video (IRE)');
title(sprintf('One line, T_H = %.2f \\mus (active %.2f \\mus)', TH*1e6, Tact*1e6));

%% 2. Stripes: frequency of the fundamental, f = N/(2 T_active)
N = 100;                                           % alternating elements per line
% Build the stripe pattern directly from position along the active line (white = 1):
vS = zeros(1, nLine); vS(iSync) = -40;
vS(iAct) = 7.5 + 92.5*double(mod(floor(xAct*N), 2) == 1);

seg = vS(iAct) - mean(vS(iAct));
w = 0.5 - 0.5*cos(2*pi*(0:numel(seg)-1)/(numel(seg)-1));        % Hann window
Nfft = 2^18;
Sg = abs(fft(seg.*w, Nfft));
fax = (0:Nfft-1)*fs/Nfft;
half = fax <= B;
[~, ipk] = max(Sg(half));
f0_measured = fax(ipk);
f0_formula  = N/(2*Tact);
fprintf('Stripe fundamental: FFT peak = %.4f MHz, N/(2T_active) = %.4f MHz\n', f0_measured/1e6, f0_formula/1e6);
assert(abs(f0_measured - f0_formula) < 0.02*f0_formula, 'FFT peak does not match N/(2T)');

fprintf('Largest N with f <= %.1f MHz: N_max = 2*T_active*B = %.1f elements\n', B/1e6, 2*Tact*B);

figure('Name','Stripe line spectrum');
subplot(2,1,1); plot(t*1e6, vS, 'LineWidth', 1.2); grid on; xlabel('time (\mus)'); ylabel('video (IRE)');
title(sprintf('N = %d stripes', N));
subplot(2,1,2); plot(fax(half)/1e6, 20*log10(Sg(half)/max(Sg(half)) + eps)); grid on; hold on;
xline(f0_formula/1e6, 'r--', 'N/(2T_{active})'); ylim([-80 5]);
xlabel('frequency (MHz)'); ylabel('dB'); title('Spectrum of one active line (Hann window)');

%% 3. AM of a multi-line video signal: DSB versus VSB
nLines = 8;
vMulti = repmat(vS, 1, nLines);                    % eight identical lines
tt = (0:numel(vMulti)-1)/fs;
Nf = numel(vMulti);
f  = (0:Nf-1)*fs/Nf;  f(f > fs/2) = f(f > fs/2) - fs;          % two-sided axis
fa = abs(f);
df = fa - fc;                                      % offset from the carrier
lpf = double(fa <= B);                             % ideal low-pass at B

m  = real(ifft(fft(vMulti - mean(vMulti)) .* lpf));% band-limit the video to B first (as a transmitter would)
m  = m/max(abs(m));                                % zero-mean message, |m| <= 1
sDSB = (1 + 0.8*m) .* cos(2*pi*fc*tt);

nyq = zeros(size(f));                              % odd-symmetric (Nyquist) vestigial slope
nyq(df >= fv) = 1;
mid = abs(df) < fv;
nyq(mid) = 0.5 + 0.5*sin(pi*df(mid)/(2*fv));
nyq(df > B) = 0;                                   % nothing above the video band
nyq(fa < fc - 2*B | fa > fc + 2*B) = 0;

abrupt = double(df >= -fv);                        % vestige kept flat, then cut
abrupt(df > B) = 0;
abrupt(fa < fc - 2*B | fa > fc + 2*B) = 0;

demod = @(sRF) real(ifft(fft(sRF .* (2*cos(2*pi*fc*tt))) .* lpf));   % coherent detector
ref   = 1 + 0.8*m;                                 % what a perfect receiver returns

% DSB passes both sidebands (gain 2 after detection, absorbed in the factor 2 above).
% A vestigial filter with H = 1 in the passband halves the baseband gain, so the
% receiver applies a gain of 2 to compare like with like.
yDSB = demod(sDSB);
yNyq = 2*demod(real(ifft(fft(sDSB) .* nyq)));
yAbr = 2*demod(real(ifft(fft(sDSB) .* abrupt)));

trim = round(fs*TH) : numel(tt) - round(fs*TH);    % ignore the first/last line (filter edges)
rms  = @(y) sqrt(mean((y(trim) - ref(trim)).^2));
fprintf('RMS error vs ideal video band: DSB %.4f | VSB Nyquist slope %.4f | VSB abrupt vestige %.4f\n', ...
        rms(yDSB), rms(yNyq), rms(yAbr));
assert(rms(yDSB) < 0.02, 'DSB recovery should be nearly exact');
assert(rms(yNyq) < 0.05, 'VSB with a Nyquist slope should recover the video faithfully');
assert(rms(yAbr) > 3*rms(yNyq), 'Abrupt vestige should distort the low frequencies');

% Nyquist condition H(fc+f) + H(fc-f) on a grid of offsets
off = linspace(0, B, 400);
Hn = @(d) (d >= fv) + (abs(d) < fv).*(0.5 + 0.5*sin(pi*d/(2*fv)));
Ha = @(d) double(d >= -fv);
sumN = Hn(off) + Hn(-off);   sumA = Ha(off) + Ha(-off);
fprintf('H(fc+f)+H(fc-f): Nyquist slope min %.3f max %.3f | abrupt vestige min %.3f max %.3f\n', ...
        min(sumN), max(sumN), min(sumA), max(sumA));
assert(max(abs(sumN - 1)) < 1e-9, 'Nyquist slope must sum to a constant');

% Bandwidth the two schemes are designed for (definitions, not measured from this test pattern)
fprintf('Designed RF bandwidth for B = %.1f MHz: DSB 2B = %.1f MHz, SSB B = %.1f MHz, VSB B + vestige = %.2f MHz\n', ...
        B/1e6, 2*B/1e6, B/1e6, (B+fv)/1e6);

figure('Name','DSB versus VSB');
S1 = fftshift(abs(fft(sDSB))/Nf);  S2 = fftshift(abs(fft(real(ifft(fft(sDSB).*nyq))))/Nf);  fs_ = fftshift(f);
subplot(3,1,1); plot(fs_/1e6, S1); xlim([0 15]); grid on; title('DSB-AM spectrum'); ylabel('|S|');
subplot(3,1,2); plot(fs_/1e6, S2); xlim([0 15]); grid on; title('VSB-AM spectrum (Nyquist vestigial slope)'); ylabel('|S|');
subplot(3,1,3); plot(off/1e6, sumN, off/1e6, sumA, '--', 'LineWidth', 1.4); grid on; ylim([0 2.2]);
legend('Nyquist slope','abrupt vestige'); xlabel('video frequency f (MHz)'); ylabel('H(f_c+f)+H(f_c-f)');

%% Local functions ----------------------------------------------------------
function v = buildLine(rowvals, nLine, iSync, iAct, xAct)
    Ww = numel(rowvals);
    v = zeros(1, nLine);                           % blanking level
    v(iSync) = -40;                                % sync tip
    col = min(Ww, floor(xAct*Ww) + 1);
    v(iAct) = 7.5 + 92.5*rowvals(col);             % black = 7.5, white = 100 IRE
end