% FM2.m
%
% A Matlab script to demonstrate the difference between using the standard
% FM algorithm to calculate an FM synthesis result versus phase increment
% approaches (including an instantaneous frequency).
%
% by Gary Scavone, McGill University, 2020.

% Signal parameters
fs = 44100;
T = 1/fs;
dur = 2.0;      % seconds
t = 0:T:dur;    % time vector

% FM parameters
H = 2;         % harmonicity ratio
fc = 200;      % center frequency
fm = H * fc;   % modulation frequency
I = 5;         % modulation index

% Standard FM algorithm
y1 = sin(2*pi*fc*t - I.*cos(2*pi*fm*t));

% Iterative approaches with phase increments (per sample period)
N = dur * fs + 1;
y2 = zeros(1, N);
y3 = zeros(1, N);
theta1 = 0;
theta2 = 0;
theta3 = 0;
theta4 = 0;
deltaTheta1 = 2 * pi * fc * T;
deltaTheta2 = 2 * pi * fm * T;
for n = 1:N
  y2(n) = cos(theta1 + theta3);   % effectively the same as above
  theta1 = theta1 + deltaTheta1;  % carrier oscillator phase increment
  theta2 = theta2 + deltaTheta2;  % modulator oscillator phase increment
  theta3 = -I*cos(theta2);        % modulated phase increment
  
  y3(n) = sin(theta4);            % calculation with instanteous frequency
  ifreq = fc + I*fm*sin(theta2);  % instantaneous frequency
  theta4 = theta4 + 2*pi*ifreq*T; % phase increment
end

Y1 = fft( y1 );
Y2 = fft( y2 );
Y3 = fft( y3 );
f = fs*(0:N-1)/N;
plot(f, 20*log10(abs([Y1; Y2; Y3].')));
xlim([0 10000])
xlabel('Frequency (Hz)');
ylabel('Log Magnitude (dB)');
sound(0.25*[y1 zeros(1, 10000) y2 zeros(1, 10000) y3], fs)