discrete Fourier transform

ft.cpp

// fourier.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <math.h>
#include <iostream>
using namespace std;
//#include <windows.h>

typedef float FLOAT;
const FLOAT fSampleRate = 44100;
const FLOAT fTime = 1.0f;
const FLOAT pi = 3.14159265359f;

FLOAT* GenerateSignal( FLOAT fFreq, FLOAT fAmpl )
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	FLOAT fOmega = 2.0f * pi * fFreq;
	for ( int i=0; i<nSamples; i++)
	{
		FLOAT fTime = i/fSampleRate; // cas v sekundach
		fData[i] = sin( fOmega * fTime ) * fAmpl;
	}

	return fData;
}

FLOAT* GenerateSignalCos( FLOAT fFreq)
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	FLOAT fOmega = 2.0f * pi * fFreq;
	for ( int i=0; i<nSamples; i++)
	{
		FLOAT fTime = i/fSampleRate; // cas v sekundach
		fData[i] = cos( fOmega * fTime );
	}

	return fData;
}

FLOAT* GenerateSignalSin( FLOAT fFreq)
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	FLOAT fOmega = 2.0f * pi * fFreq;
	for ( int i=0; i<nSamples; i++)
	{
		FLOAT fTime = i/fSampleRate; // cas v sekundach
		fData[i] = sin( fOmega * fTime );
	}

	return fData;
}

/* FLOAT* SignalCopy( FLOAT* arrSignal )
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	for ( int i=0; i<nSamples; i++)
		fData[i] = arrSignal[i];
	return fData;
} */

FLOAT* SignalAdd( FLOAT* arrSignal1, FLOAT* arrSignal2 )
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	for ( int i=0; i<nSamples; i++)
		fData[i] = arrSignal1[i] + arrSignal2[i];
	return fData;
}

FLOAT Sum( FLOAT* arrSignal )
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	FLOAT fSum = 0;
	for ( int i=0; i<nSamples; i++)
		fSum += arrSignal[i];
	return fSum;
}

FLOAT* SignalMul( FLOAT* arrSignal1, FLOAT* arrSignal2 )
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	for ( int i=0; i<nSamples; i++)
		fData[i] = arrSignal1[i] * arrSignal2[i];
	return fData;
}


void SignalMul( FLOAT* arrSignal, FLOAT fCoef )
{
	int nSamples = (int)(fSampleRate * fTime);
	for ( int i=0; i<nSamples; i++)
		arrSignal[i] *= fCoef;
}

FLOAT FourierCalc(FLOAT* arrData, FLOAT fFreq)
{
	FLOAT* arrHarmonicCos = GenerateSignalCos( fFreq );
	FLOAT* arrHarmonicSin = GenerateSignalSin( fFreq );

	FLOAT* arrMulCos = SignalMul( arrData, arrHarmonicCos );
	FLOAT* arrMulSin = SignalMul( arrData, arrHarmonicSin );
	
	FLOAT fSamples = fSampleRate * fTime;
	FLOAT fReal = Sum( arrMulCos ) / fSamples * 2;
	FLOAT fImag = Sum( arrMulSin ) / fSamples * 2;

	FLOAT fAmpl = sqrt(fReal*fReal + fImag*fImag);

	delete arrMulCos;
	delete arrMulSin;
	delete arrHarmonicCos;
	delete arrHarmonicSin;
	return fAmpl;
}

FLOAT* GenerateHann( FLOAT fAlpha )
{
	int nSamples = (int)(fSampleRate * fTime);
	FLOAT* fData = new FLOAT[nSamples];
	for ( int i=0; i<nSamples; i++)
	{
		FLOAT fArg = (2.0f * pi * i) / ( nSamples - 1 );
		fData[i] = 0.5f * ( 1.0f - cos( fArg ) );
	}
//		fData[i] = arrSignal1[i] * arrSignal2[i];
	return fData;
}

int _tmain(int argc, _TCHAR* argv[])
{	
	FLOAT* fSignal1 = GenerateSignal( 1000.0f, 0.5f );
	FLOAT* fSignal2 = GenerateSignal( 3500.0f, 0.3f );
	FLOAT* fSignal = SignalAdd( fSignal1, fSignal2 );  
	FLOAT* fWindow = GenerateHann(0.5f);

	cout << " vypis z fSignal1" << fSignal[0]<< endl;
	cin.get();

	FLOAT fArea = Sum( fWindow ) / ( fSampleRate * fTime );
	cout << " vypis z fSignal2" << fArea << endl;
	cin.get();
	SignalMul( fWindow, 1.0f / fArea );

	cout << "vypis z fsignal3" << fSignal[0] << endl;
	cin.get();

	printf("Obsah po prenasobeni = %f\n", Sum(fWindow) / ( fSampleRate * fTime ) ); 

	FLOAT* fWindowed = SignalMul( fSignal, fWindow);

	for (FLOAT fFreq = 100; fFreq < 5000; fFreq += 100 )
	{
		FLOAT fAmplitude = FourierCalc(fWindowed, fFreq);
		printf("freq %5.0f ampl = %2.5f\n", fFreq, fAmplitude);
	}
	printf("Test:\n");
	for (FLOAT fFreq = 1000-10; fFreq < 1000+10; fFreq += 0.7f)
	{
		FLOAT fAmplitude = FourierCalc(fWindowed, fFreq);
		printf("freq %5.2f ampl = %2.5f\n", fFreq, fAmplitude);
	}

	delete fSignal1;
	delete fSignal2;
	delete fSignal;
	delete fWindow;
	delete fWindowed;

	std::cin.get();

	return 0;
}