// FREQUENCY MODULATION
// 
// Several examples including sweeps, sine modulations, 
// oscillograms, spectra and sounds
// 
// -------------------------
// Author: Federico Miyara
// Date: 2020-03-06
//       2020-10-15

clear
close

// Sample rate
Fs = 44100;
// Durations for fast and slow sweeps
T = 2;
TT = 10;

// Time vectors 
// Fast sweep
t = 0:1/Fs:T;
// Slow sweep
tt = 0:1/Fs:TT;

// 1) Exponential sweep

// Initial and  final frequencies inicial y
f1 = 20
f2 = 20000
K = log(f2/f1)

// Frecuency vectors
f = f1 * exp(K*t/T);
ff = f1 * exp(K*t/TT);

// Phase, obtained integrating analytically the 
// frequency
// NOTE: In a difficult case to integrate analytically,
//       it could be integrated numerically, using
//       phi1 = 2*%pi*cumsum(f)/Fs; 
phi = 2*%pi*f1*T/K * (exp(K*t/T) - 1);
phii = 2*%pi*f1*TT/K * (exp(K*tt/TT) - 1);

// Sweeps
x = sin(phi);
xx = sin(phii);

// Sound playback
// Fast sewwp
playsnd(0.1*x, Fs, 1)
// Slow playback
playsnd(0.1*xx, Fs, 1)

// Sweep plots 
scf(1);
clf();

subplot(4,1,1)
plot(t,f)

subplot(4,1,2)
plot(t,x)

subplot(4,1,3)
plot(t,f)
gca().data_bounds = [0, 0; 0.4, 80]

subplot(4,1,4)
plot(t,x)
gca().data_bounds = [0, -1.2; 0.4, 1.2]

 
// 2) Sine FM

// Sample rate
Fs = 44100

// Duration
T = 20

// Carrier
fo = 400

// Modulating frequency
fm = 5

// Desviación de frecuencia
deltaf = 20

// Modulation índex
mf = deltaf/fm 

// Time vector
t = 0:1/Fs:T;

// Frequency
f = fo + deltaf*sin(2*%pi*fm*t);

// Frequencia modulated signal
x = sin(2*%pi*fo*t - mf*cos(2*%pi*fm*t));

// Playback
playsnd(0.1*x(1:2*Fs), Fs);

// Three versions of the FFT spectrum
N = 65536
fff = (0:(N/2 - 1))*Fs/N;

// a) Plain FFT ranging a single frame
X = fft(x(1:N));
X = abs([X(1)/N, X(2:N/2)]*2/N);
X = X/sqrt(2);

// b) Average over several frames, no time window
exec fftmean.sci;
Y = fftmean(x, Fs, N, "box");

// c) Average over several frames, Hann window
Z = fftmean(x, Fs, N, "hann");

// Sine modultion plots
scf(2);
clf(2);

subplot(3,1,1)
plot(t,f)
gca().data_bounds = [0, 0; 0.1, (fo + deltaf)*1.1]

subplot(3,1,2)
plot(t,x)
gca().data_bounds = [0, -1.2; 0.1, 1.2]

subplot(3,1,3)
// Single, no window
plot(fff, X*sqrt(2), 'k')
// Average, no window
plot(fff, Y*sqrt(2), 'b')
// Average, Hann window
plot(fff, Z*sqrt(2), 'r')  
legend("Single, no window", "Average, no window", "Average, Hann window")
gca().data_bounds = [0, 0; 2*fo, 0.35];
gca().grid = [0,0];

// Theoretical center band
Ao = besselj(0, mf)

// Index of carrier frequency
Nfo = round(fo/(Fs/N)) + 1

// Verification 
fff(Nfo)

// Empirical center band
// To get an approximated value of the center band partial, 
// the six closest FFT lines around the center band's peak 
// value are added energetically.
// Multiplication by sqrt(2) is because theoretical spectrum
// of FM is expressed as amplitudes and empirical spectrum is
// expressed as RMS values. 
sqrt(sum(X(Nfo-3:Nfo+3).^2))*sqrt(2)
sqrt(sum(Y(Nfo-3:Nfo+3).^2))*sqrt(2)
sqrt(sum(Z(Nfo-3:Nfo+3).^2))*sqrt(2)

besselj(1:15,mf)