SSB_modulation_FFT_IFFT.m

% Input signal: s(t) = 0.5 * cos(2*pi*800*t) - 0.5 * cos(2*pi*1200*t);
% mixer wave: c(t) = cos(2*pi*5000*t)
f1 = 800
f2 = 1200
fc = 5000

% FFT size
N = 1024


% input signal sampling
% assume at time t0, input signal angular value = w0
% assume cos(2*pi*800*t0) = cos(w_basis), we got cos(2*pi*1200*t0) = cos(1.5*w_basis)
% so, s(w0) = 0.5 * cos(w_basis) - 0.5 cos(1.5 * w_basis)
% note that, s(0) = s(2*pi), we w_basis has period = 4*pi
% Sampling interval in angular domain
ws = 4*pi/N;
% sample sequence in angular domain
xw = 0:ws:4*pi*(N-1)/N;
% s(t) = 0.5 * cos(2*pi*800*t) - 0.5 * cos(2*pi*1200*t);
s = 0.5 * cos(xw) - 0.5 * cos(f2/f1 * xw);

% carrier for input signal
c = cos(fc/f1 * xw);

figure(1);
subplot(231);
plot(s);
title("s(t) = 0.5 * cos(2*pi*800*t) - 0.5 * cos(2*pi*1200*t)");
hold;
%plot(c);
hold off;

% sample frequency
% as we sample cos(2*pi*800*t) for two periods, we can get the sampling frequency
% N*ws = 4*pi
%   => N*ws/2 = 2*pi, because T = 1/800
%   => Ts * N/2 = 1/800
%   => Fs = 400 * N
fs = 400*N;
S = fft(s, N);
xf = [-ceil((N-1)/2) : floor((N-1)/2)] * (fs/N);

subplot(234);
range = 20
sel = (N/2 + 1 - range) : (N/2 + 1 + range);

stem(xf(sel), abs(fftshift(S))(sel));
hold;
scale = 8;
xf2 = [-ceil((scale*N - 1)/2) : floor((scale*N - 1)/2)] * (fs/(N*scale));
sel2 = (N/2 + 1 - range) * scale : (N/2 + 1 + range) * scale;
plot(xf2(sel2), abs(fftshift(fft(s, N * scale))(sel2)));
title("FFT of s(t)");

hold off;

% modulated signal s(t) * c(t)
subplot(232);

x = s .* c;
plot(x);
title("x(t) = s(t)*cos(2*pi*5000*t)");


% spectrum of modulated signal s(t) * c(t)
subplot(235);

stem(xf(sel), abs(fftshift(fft(x, N))(sel)));
hold;

plot(xf2(sel2), abs(fftshift(fft(x, N * scale))(sel2)));
title("FFT of x(t)");

hold off;

% modulated signal s(t) * c(t) - s_h(t) * c_h(t)
pkg load signal;
% hilbert transform of input signal
s_h = imag(hilbert(s));
% carrier for hilbert transformed signal
c_h = sin(fc/f1 * xw);
y = x - s_h .* c_h;
subplot(233);
plot(y);
title("SSB(t) = s(t)*cos(2*pi*5000*t) - s^(t)*sin(2*pi*5000*t)");

% spectrum of modulated signal s(t) * c(t) - s_h(t) * c_h(t)
subplot(236);
stem(xf(sel), abs(fftshift(fft(y, N))(sel)));
hold;
plot(xf2(sel2), abs(fftshift(fft(y, N * scale))(sel2)));
title("FFT of SSB(t)");
hold off;

% demodulation
figure(2);
subplot(131);
% receiver demod: z(t) = ssb(t) * cos(wt)
z = y .* c
plot(z);
title("z(t) = SSB(t)*cos(2*pi*5000*t)");

subplot(132);
M = N * scale;
intval = fs/M;
plot(0:intval:12000-intval , abs(fft(z, M))(1:12000/intval));
title("z(t) spectrum");
hold;
% show low-pass window
plot(0:intval:5000-intval, [ones(1, 99),0]*180);
hold off;

subplot(133);

% low-pass filtering to keep only single side band
window_size = ceil(fc/ (fs/N));
Z = fft(z, N)(1:window_size);
% recover signal by IFFT on the single side band
plot(0:N-1, real(ifft(Z, N)));
title("recovered signal s'(t)");