%Script for simulating wavetable synthesis and seeing aliasing.

%The script allows you to specify a sample rate, a playback frequency of

% a wavetable, and to specify an amount of aliasing to see. Increasing

% amounts of aliasing can be seen by increasing the aliasing factor.

SAMPLE_RATE = 32768;%Set the sampling rate

NYQUIST_RATE = SAMPLE_RATE/2;%Figure out the Nyquist Rate

NUM_OF_SAMPLES = 256;%This will be the number of samples in the wavetable.

PLAYBACK_FREQUENCY = 2000; %The frequency of the note we are simulating

NUMBER_OF_PLAYBACK_CYCLES = 25;%The number of cycles of playback of the wavetable

ALIASING_FACTOR = 1;

NUM_OF_ALIASES = ALIASING_FACTOR*NYQUIST_RATE/PLAYBACK_FREQUENCY;%The wavetable is built from summing together sine wave aliases

%Make large wavetable containing the 128 samples of each of the first

%2000 harmonics of a 1Hz sine wave

sin_table = zeros(2000,NUM_OF_SAMPLES);

p = 1:NUM_OF_SAMPLES;

for z = 1:2000

for y = 1:NUM_OF_SAMPLES

sin_table(z,y) = sin(2.*pi.*z.*y./NUM_OF_SAMPLES);

end

end

%plot(p,sin_table);

wavetable = zeros(1, NUM_OF_SAMPLES);

%This may seem unnecessary. We are creating a sawtooth wavetable to make

%a square wave. But, this is how it's done in my synth and that's

%particularly what I wanted to simulate. We can make square

%waves of various pulse widths by adding together square waves with

%different phase angles between them.

for k = 1:NUM_OF_SAMPLES

for n = 1:NUM_OF_ALIASES

wavetable(1,k) = wavetable(1,k) + (1/n).*(sin_table(n,k).*(cos((n-1).*pi./(NUM_OF_ALIASES*2))).^2);

end

end

%This process below simulates a wavetable lookup oscillator for generating

%a square wave with 50% pulse width, or duty cycle.

phase_accumulator = 0;

wavetable_lookup_index = 0;

wavetable_reverse_lookup_index = 0;

wavetable_divisor = SAMPLE_RATE/NUM_OF_SAMPLES;

%Make a matrix to hold the samples of the wavetable playback

NUMBER_OF_PLAYBACK_SAMPLES = ceil(SAMPLE_RATE*NUMBER_OF_PLAYBACK_CYCLES/PLAYBACK_FREQUENCY);

playback = zeros(1, NUMBER_OF_PLAYBACK_SAMPLES);

index = 1;

for f = 1:PLAYBACK_FREQUENCY:(SAMPLE_RATE*NUMBER_OF_PLAYBACK_CYCLES)

phase_accumulator = phase_accumulator + PLAYBACK_FREQUENCY;

if (phase_accumulator > SAMPLE_RATE)

phase_accumulator = phase_accumulator - SAMPLE_RATE;

end

wavetable_lookup_index = floor(phase_accumulator/wavetable_divisor);

if(wavetable_lookup_index < 1)

wavetable_lookup_index = 1;

end

if(wavetable_lookup_index > NUM_OF_SAMPLES/2)

wavetable_reverse_lookup_index = wavetable_lookup_index - NUM_OF_SAMPLES/2;

else

wavetable_reverse_lookup_index = wavetable_lookup_index + NUM_OF_SAMPLES/2;

end

if(wavetable_reverse_lookup_index < 1)

wavetable_reverse_lookup_index = 1;

end

playback(1, index) = wavetable(1, wavetable_lookup_index) - wavetable(1, wavetable_reverse_lookup_index);

index = index + 1;

end

%Plot the wavetable and its FFT.

k = 1:NUMBER_OF_PLAYBACK_SAMPLES;

T = 1/SAMPLE_RATE;

subplot(2,1,1);

plot(T*k,playback);

NFFT = 2^nextpow2(NUMBER_OF_PLAYBACK_SAMPLES*2); % Next power of 2 from length of y

Y = fft(playback, NFFT)/NUMBER_OF_PLAYBACK_SAMPLES;

f = SAMPLE_RATE/2*linspace(0,1,NFFT/2+1);

subplot(2,1,2);

plot(f,2*abs(Y(1:NFFT/2+1)))

print('-dpng', 'test.png');

%plot(T*k,fft(playback,k));