dft.c

Discrete Fourier Transform and Power Spectrum Calculates Power Spectrum from a Time Series

/*	Discrete Fourier Transform and Power Spectrum
	Calculates Power Spectrum from a Time Series
	Copyright 1985 Nicholas B. Tufillaro
*/

#include 
#include 

#define PI (3.1415926536)
#define SIZE 512

double ts[SIZE], A[SIZE], B[SIZE], P[SIZE];

main()
{
	int i, k, p, N, L;
	double avg, y, sum, psmax;

	/* read in and scale data points */
	i = 0;
	while(scanf("%lf", &y) != EOF) {
		ts[i] = y/1000.0;
		i += 1;
	}
	/* get rid of last point and make sure # 
		of data points is even */
	if((i%2) == 0)
		i -= 2;
	else
		i -= 1;
	L = i; N = L/2;
	/* subtract out dc component from time series */
	for(i = 0, avg = 0; i < L; ++i) {
		avg += ts[i];
	}
	avg = avg/L;
	/* now subtract out the mean value from the time series */
	for(i = 0; i < L; ++i) {
		ts[i] = ts[i] - avg;
	}
	/* o.k. guys, ready to do Fourier transform */
	/* first do cosine series */
	for(k = 0; k <= N; ++k) {
		for(p = 0, sum = 0; p < 2*N; ++p) {
			sum += ts[p]*cos(PI*k*p/N);
		}
		A[k] = sum/N;
	}
	/* now do sine series */
	for(k = 0; k < N; ++k) {
		for(p = 0, sum = 0; p < 2*N; ++p) {
			sum += ts[p]*sin(PI*k*p/N);
		}
		B[k] = sum/N;
	}
	/* lastly, calculate the power spectrum */
	for(i = 0; i <= N; ++i) {
		P[i] = sqrt(A[i]*A[i]+B[i]*B[i]);
	}
	/* find the maximum of the power spectrum to normalize */
	for(i = 0, psmax = 0; i <= N; ++i) {
		if(P[i] > psmax)
			psmax = P[i];
	}
	for(i = 0; i <= N; ++i) {
		P[i] = P[i]/psmax;
	}
	/* o.k., print out the results: k, P(k) */
	for(k = 0; k <= N; ++k) {
		printf("%d %g\n", k, P[k]);
	}
}