ZynAddSubFX/PADsynth/CPP Implementation

This page contains a C++ class that implements PADsynth Synthesis Algorithm. Please note that you must implement the IFFT routine.

PADsynth.h

/*
    PADsynth implementation as ready-to-use C++ class.
    By: Nasca O. Paul, Tg. Mures, Romania
    This implementation and the algorithm are released under Public Domain
    Feel free to use it into your projects or your products ;-)

    This implementation is tested under GCC/Linux, but it's 
    very easy to port to other compiler/OS.

    Please notice that you have to implement the virtual method IFFT() with your IFFT routine.
    
*/

#ifndef PADSYNTH_H
#define PADSYNTH_H

#ifndef REALTYPE
#define REALTYPE float
#endif

class PADsynth{
    public:
	/*  PADsynth:
		N                - is the samplesize (eg: 262144)
		samplerate 	 - samplerate (eg. 44100)
		number_harmonics - the number of harmonics that are computed */
	PADsynth(int N_,int samplerate_,int number_harmonics_);

	~PADsynth();

	/* set the amplitude of the n'th harmonic */
	void setharmonic(int n,REALTYPE value);

	/* get the amplitude of the n'th harmonic */
	REALTYPE getharmonic(int n);

	/*  synth() generates the wavetable
	    f		- the fundamental frequency (eg. 440 Hz)
	    bw		- bandwidth in cents of the fundamental frequency (eg. 25 cents)
	    bwscale	- how the bandwidth increase on the higher harmonics (recomanded value: 1.0)
	    *smp	- a pointer to allocated memory that can hold N samples */
	void synth(REALTYPE f,REALTYPE bw,
	    REALTYPE bwscale,REALTYPE *smp);
    protected:
	int N;			//Size of the sample

	/* IFFT() - inverse fast fourier transform
	   YOU MUST IMPLEMENT THIS METHOD!
	   *freq_real and *freq_imaginary represents the real and the imaginary part of the spectrum, 
	   The result should be in *smp array.
	   The size of the *smp array is N and the size of the freq_real and freq_imaginary is N/2 */
	virtual void IFFT(REALTYPE *freq_real,REALTYPE *freq_imaginary,REALTYPE *smp)=0;


	/* relF():
	    This method returns the N'th overtone's position relative 
	    to the fundamental frequency.
	    By default it returns N.
	    You may override it to make metallic sounds or other 
	    instruments where the overtones are not harmonic.  */
	virtual REALTYPE relF(int N);

	/* profile():
	    This is the profile of one harmonic
	    In this case is a Gaussian distribution (e^(-x^2))
            The amplitude is divided by the bandwidth to ensure that the harmonic
	    keeps the same amplitude regardless of the bandwidth */
	virtual REALTYPE profile(REALTYPE fi, REALTYPE bwi);

	/* RND() - a random number generator that 
	    returns values between 0 and 1
	*/
	virtual REALTYPE RND();

    private:
	REALTYPE *A;		//Amplitude of the harmonics
	REALTYPE *freq_amp;	//Amplitude spectrum
	int samplerate;
	int number_harmonics;
};

#endif

PADsynth.cpp

/*
    PADsynth implementation as ready-to-use C++ class.
    By: Nasca O. Paul, Tg. Mures, Romania
    This implementation and the algorithm are released under Public Domain
    Feel free to use it into your projects or your products ;-)

    This implementation is tested under GCC/Linux, but it's 
    very easy to port to other compiler/OS. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "PADsynth.h"

PADsynth::PADsynth(int N_, int samplerate_, int number_harmonics_){
    N=N_;
    samplerate=samplerate_;
    number_harmonics=number_harmonics_;
    A=new REALTYPE [number_harmonics];
    for (int i=0;i<number_harmonics;i++) A[i]=0.0;
    A[1]=1.0;//default, the first harmonic has the amplitude 1.0

    freq_amp=new REALTYPE[N/2];
};

PADsynth::~PADsynth(){
    delete[] A;
    delete[] freq_amp;
};

REALTYPE PADsynth::relF(int N){
    return N;
};

void PADsynth::setharmonic(int n,REALTYPE value){
    if ((n<1)||(n>=number_harmonics)) return;
    A[n]=value;
};

REALTYPE PADsynth::getharmonic(int n){
    if ((n<1)||(n>=number_harmonics)) return 0.0;
    return A[n];
};

REALTYPE PADsynth::profile(REALTYPE fi, REALTYPE bwi){
    REALTYPE x=fi/bwi;
    x*=x;
    if (x>14.71280603) return 0.0;//this avoids computing the e^(-x^2) where it's results are very close to zero
    return exp(-x)/bwi;
};

void PADsynth::synth(REALTYPE f,REALTYPE bw,REALTYPE bwscale,REALTYPE *smp){
    int i,nh;
    
    for (i=0;i<N/2;i++) freq_amp[i]=0.0;//default, all the frequency amplitudes are zero

    for (nh=1;nh<number_harmonics;nh++){//for each harmonic
	REALTYPE bw_Hz;//bandwidth of the current harmonic measured in Hz
        REALTYPE bwi;
	REALTYPE fi;
	REALTYPE rF=f*relF(nh);
	
        bw_Hz=(pow(2.0,bw/1200.0)-1.0)*f*pow(relF(nh),bwscale);
	
	bwi=bw_Hz/(2.0*samplerate);
	fi=rF/samplerate;
	for (i=0;i<N/2;i++){//here you can optimize, by avoiding to compute the profile for the full frequency (usually it's zero or very close to zero)
	    REALTYPE hprofile;
	    hprofile=profile((i/(REALTYPE)N)-fi,bwi);
	    freq_amp[i]+=hprofile*A[nh];
	};
    };

    REALTYPE *freq_real=new REALTYPE[N/2];
    REALTYPE *freq_imaginary=new REALTYPE[N/2];

    //Convert the freq_amp array to complex array (real/imaginary) by making the phases random
    for (i=0;i<N/2;i++){
	REALTYPE phase=RND()*2.0*3.14159265358979;
	freq_real[i]=freq_amp[i]*cos(phase);
	freq_imaginary[i]=freq_amp[i]*sin(phase);
    };
    IFFT(freq_real,freq_imaginary,smp);
    delete [] freq_real;
    delete [] freq_imaginary;

    //normalize the output
    REALTYPE max=0.0;
    for (i=0;i<N;i++) if (fabs(smp[i])>max) max=fabs(smp[i]);
    if (max<1e-5) max=1e-5;
    for (i=0;i<N;i++) smp[i]/=max*1.4142;
    
};

REALTYPE PADsynth::RND(){
    return (rand()/(RAND_MAX+1.0));
};