?? dft.c
字號(hào):
/* 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]);
}
}
?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -