%{

PADsynth MATLAB implementation
Matt Montag - October 23, 2010

%}
clear all; clc;

% Input:
N = 131072;
	% wavetable size. Recommended to be a power of 2.
	% usually a big number, like 262144.
samplerate = 22050; 
	% the samplerate (eg. 44100)
f = [50 250 300 375 450]; % an arbitrary set of tones
shift = 0; f = f .* 2^(shift/12); % transpose if you want
	% frequency of the the fundamental note(s) (eg. 440) 
bw = 100;
	% bandwidth of first harmonic in cents (eg. 50 cents)
number_harmonics = 4; 
	% the number of harmonics. Of course, number_harmonics<(samplerate/f)
A = [1 rand(1, number_harmonics-1)]; 
	% amplitude of each overtone - randomize or invent your own

% Output:
smp = zeros(1, N); % [0..N-1] stores the generated wavetable

% Internal variables:
freq_amp = zeros(1, N/2); % [0..N/2-1] = {0,0,0,0,...,0}
freq_complex = zeros(1, N/2);
freq_phase = rand(1,N/2) * 2*pi; % randomize phases

% Harmonic generation:
for j=1:length(f) % For each frequency...
	for nh = 1:number_harmonics % For each harmonic...   
		bw_Hz=((2^(bw/1200))-1.0)*f(j)*nh;
		bwi=bw_Hz/(2.0*samplerate);
		fi=f(j)*nh/samplerate;
		for i=1:N/2
			% Add the gaussian harmonic profile centered on fi:
			% (profile can be a different function if we want)
			hprofile=profile((i/N)-fi,bwi);
			freq_amp(i) = freq_amp(i)+hprofile*A(nh);
			freq_complex(i) = pol2cart(freq_phase(i), freq_amp(i));
		end
	end
end

% Frequency to time domain:
smp = real(ifft(freq_complex,N));
%smp = [smp smp smp smp];

% Apply amplitude envelope:
attack = linspace(0,1,length(smp)/4);
smp(1:length(smp)/4) = smp(1:length(smp)/4) .* attack;
release = linspace(1,0,length(smp)/16);
smp(length(smp)*15/16+1:length(smp)) = ...
	smp(length(smp)*15/16+1:length(smp)) .* release;

% Apply convolution reverb:
reverb = wavread('church_ir.wav');
wav(:,1) = fconv(smp, reverb(:,1)');
wav(:,2) = fconv(smp, reverb(:,2)');

% 50% wet/dry mix:
wav(1:length(smp),1) = wav(1:length(smp),1) + smp';
wav(1:length(smp),2) = wav(1:length(smp),2) + smp';

% Normalize and output:
wav = .9*wav./max(max(abs(wav)));
wavwrite(wav, samplerate, 16, 'padsynth.wav');
soundsc(wav, samplerate);