% dsfdj.m
% 
% Matlab implementation of Dodge and Jerse DSF formula for a
% bandlimited impulse train.
%
% By Gary P. Scavone, McGill University, 2005.

fs = 44100;   % sample rate
Ts = 1 / fs;  % sample period
N = 40000;     % signal length
F = 400.0;    % signal frequency
M = 40;       % number of harmonics
A = 1.0;      % amplitude

y = zeros(1, N);

epsilon = eps;
epsilon = 0.000001;
for n = 1:N
  den = sin( pi*F*Ts*n );
  if abs( den ) < epsilon
    y(n) = (0.5 * A / M) * ( (2*M+1) * cos( (2*M+1)*pi*F*Ts*n ) / cos( pi*F*Ts*n ) - 1 );
  else
    y(n) = (0.5 * A / M) * ( ( sin( (2*M+1)*pi*F*Ts*n ) / den ) - 1 );
  end
end
  
subplot(2, 1, 1)
plot(y)
axis([1 length(y) -0.3 1.1])
grid
xlabel('Time (samples)')
ylabel('Amplitude');

subplot(2, 1, 2)
fftplot( fft(y), fs, 1, 1 )

sound(y, fs)
%print -deps 'dsfdj.eps'
