I am using realft.f and four1.f to fourier transform a curve. I want to then use the coefficients to get back to the original curve - as a sanity check. I was trying to do this in a spreadsheet (openoffice's calc). Here I point out that I am not using NumRec to inverse transorm the curve - I am trying to do this in a spreadsheet.

To be clear, I am fourier transforming a curve and then checking that I have done things correctly by trying to use the coefficients to get back to that curve.

I have not done this successfully though. I think the problem may be with the series I am using to reform the original curve. Here it is:

curve(x) = 1/256 * [ a0 + a1*cos(2*PI*x) + a2*sin(2*PI*x) + a3*cos(2*2*PI*x) + a4*sin(2*2*PI*x) + a5*cos(3*2*PI*x) + a6*sin(3*2*PI*x) + ...+ an*cos(n*2*PI*x) + an*sin(n*2*PI*x)]

where n=255 as there are 256 datapoints.

I wasn't sure whether to use only COS, only SIN, or both as in the case shown.

When I try to plot the original curve in a spreadsheet using the fourier coefficients I get a curve around the correct magnitude but apparently not of the same shape. Admittedly I am only using the first 10 of the availavle 256 coefficients but this should be enough shouldn't it?

I would be most grateful is someone could give me some pointers, even more so if someone could implement what I am trying to do and mail it to me.

Below is the original curve data and its fourier coeffiicients that I get from realft.

Thanks,

Astro_Dave

CURVE DATA
9.45755291
9.62607002
9.41541958
9.39117432
9.37016296
9.39124203
9.34696579
9.38711452
9.36457348
9.34203339
9.36322975
9.38711262
9.38385582
9.36967945
9.35229969
9.33474636
9.41429043
9.58455181
9.36493874
9.34459877
9.34112453
9.08434486
8.77544689
9.3456583
9.33924389
9.35004616
9.33586502
9.47477341
9.3496542
9.36608887
9.38252449
9.36725998
9.35990906
9.35616302
9.37639713
9.36164188
9.38313484
9.36677933
9.37950134
9.38918877
9.38776588
9.38634396
9.41954613
9.41360092
9.41032982
9.41456127
9.41879272
9.43068314
9.43272591
9.44132328
9.45273781
9.45835781
9.47146606
9.47241497
9.47990799
9.49232864
9.49467182
9.50969124
9.5155735
9.52438927
9.53694153
9.52985477
9.56406879
9.58093166
9.58636761
9.61130333
9.62365723
9.63601017
9.64999962
9.66399193
9.67787743
9.69051933
9.70580769
9.72545052
9.74037075
9.75983334
9.77890587
9.78932285
9.81152153
9.82201862
9.83375549
9.84975529
9.86005402
9.87017727
9.87294006
9.8761692
9.87873268
9.87677765
9.87474823
9.86823368
9.86176586
9.85609913
9.84671593
9.83480835
9.81431103
9.80143452
9.78831863
9.76694489
9.75150204
9.73998928
9.72300911
9.70602798
9.68771935
9.67456532
9.6726408
9.65803337
9.63756657
9.63959789
9.61335182
9.60857677
9.58346272
9.57935429
9.57546043
9.56206989
9.54931831
9.539814
9.53035164
9.52824879
9.51671696
9.50518608
9.49293327
9.48523331
9.48058796
9.46981239
9.47700691
9.45289898
9.45382881
9.4492178
9.44504452
9.43956184
9.43259048
9.42440033
9.41771698
9.42802238
9.40560532
9.40616035
9.39897537
9.39390373
9.38931274
9.38490295
9.38225174
9.37645149
9.37421703
9.37849045
9.37527084
9.36910152
9.36735535
9.36582947
9.36591244
9.36692333
9.36865616
9.37188625
9.37511635
9.33651066
9.33363533
9.33541203
9.35033131
9.36524963
9.37273502
9.38073921
9.38606548
9.38201332
9.39547443
9.39396381
9.40096569
9.40864182
9.41763401
9.42728043
9.43514442
9.43212891
9.42689037
9.4398613
9.44766426
9.45441723
9.45598316
9.46534348
9.47952938
9.48504829
9.48683739
9.4785347
9.4829216
9.45723152
9.31463432
9.52061558
9.5274992
9.53953457
9.52950478
9.56334591
9.56492329
9.59346771
9.60110092
9.60873508
9.62848473
9.63945389
9.654562
9.66354752
9.67253304
9.69482136
9.72076797
9.71932888
9.74295712
9.76878262
9.79382324
9.80501175
9.82490635
9.88022423
9.85644817
9.87336254
9.88437843
9.8596077
9.85300922
9.91610909
9.90988159
9.91546249
9.95537663
9.91590309
9.92942047
9.9113512
9.90699196
9.90263176
9.87662983
9.8699007
9.84224033
9.83666325
9.80142212
9.80611897
9.76632404
9.76388454
9.73762703
9.72424984
9.71436882
9.66441345
9.64506054
9.6479063
9.64247513
9.65267563
9.60693455
9.56829262
9.56783104
9.52569389
9.48280907
9.4573431
9.5151968
9.4843235
9.51298141

FOURIER COEFFICIENTS
8.66784668
-10.0407591
-11.0551052
-9.02458191
-26.4664383
-10.1975155
-11.3204832
-17.3006954
-14.655899
-9.6986599
-9.99865532
-9.25029659
-7.71470308
-9.66795444
-10.0641403
-9.32540798
-10.0353327
-9.78503799
-10.4221878
-8.72947311
-10.8556919
-8.66561317
-10.5844059
-8.81226444
-9.59651661
-8.31613636
-8.35026741
-9.00021458
-8.50732803
-9.71811485
-8.99465275
-9.99221802
-9.31563473
-9.27536297
-9.00892353
-7.83052826
-8.51012897
-7.26266384
-8.12908649
-7.02258587
-6.29125071
-7.04386187
-4.83579874
-7.51653862
-4.62708998
-8.29381275
-4.81829643
-8.49766254
-4.93740368
-8.07250786
-4.9110775
-7.05817223
-4.60613585
-6.53352404
-3.839252
-6.51326895
-2.07395411
-6.16529608
-0.307168424
-6.54149008
-0.0321323164
-7.66170597
-0.111917175
-7.9978919
-0.272148043
-7.61720896
-0.667747855
-6.99189091
-0.427708298
-6.29030037
0.486132324
-5.9281683
2.06059527
-5.33631516
3.78651428
-5.17773581
4.44747448
-5.95039749
4.63988352
-6.20447636
5.04524612
-5.85038853
4.98121262
-5.16821957
4.93133354
-4.61358213
5.39141273
-3.99988627
6.47164297
-3.32578778
7.8859272
-3.06592488
8.41023159
-3.647753
8.84207726
-3.76602054
9.45459461
-3.46014237
9.55307007
-3.1047976
9.23447418
-3.0143671
9.27154636
-2.60988355
9.80605316
-2.30901837
10.377821
-2.12926435
10.5959215
-2.51077032
10.4960718
-2.62784266
10.6765432
-2.11385345
10.4980545
-1.42065346
9.93320274
-1.07053876
9.74751186
-0.828774691
9.9942894
-0.419737786
10.4992657
-0.206534937
10.2874269
-0.414682835
9.82812786
-0.760986328
9.53088379
-0.146464199
9.19252968
0.614869893
8.50602722
1.33762419
8.42090702
1.86148679
9.08857536
2.20356536
9.79204082
2.14014983
10.0342112
1.43620348
9.86219692
1.11382627
9.70605946
1.70766878
9.18973351
2.58963728
8.46724987
3.12244296
8.3094368
3.42476201
8.88469982
3.71856499
9.64365482
3.55795789
9.95253754
2.92336726
9.72990036
2.55364537
9.12665081
3.00512338
8.48521996
3.77827668
7.57858896
4.21870375
7.33317041
4.71685553
7.86352921
5.14211559
8.51159
5.3294754
8.50252724
5.03713703
8.10404968
4.92804718
7.95691872
5.21304321
7.6372304
5.71681643
7.14857292
5.88794804
6.86262941
6.33467007
7.2115922
6.76252747
7.47684336
7.10133362
7.39650869
7.02645922
7.11308289
6.69973516
6.99335241
6.8289957
6.88267851
7.14978075
6.3096652
7.29562807
6.16019058
7.30096579
6.24180174
7.68480825
6.15059757
7.85786533
5.89207172
7.7655201
5.7081337
7.54423857
5.6664362
7.75276709
5.42998028
7.96477699
4.95618153
8.01960564
4.49684286
8.0266943
4.65939188
8.55276966
4.56366682
8.69835854
4.30529594
8.55066013
3.96669149
8.39612198
3.68285322
8.68568611
3.26667452
8.98417568
2.61707497
9.09769249
2.4001193
9.38888931
2.6424458
9.71781349
3.0660038
9.69596863
2.79780149
9.13806152
2.4703033
8.88617325
1.97132027
9.1353817
1.59944129
9.78473949
0.971777081
9.94487095
0.752054751
9.89900875
1.10812867
9.92794228
1.42143166
9.57983398
1.03150284
9.01546955
0.46472615

Hey,

Im facing the exact same problem now.. If you found out the answer, please do let me know. thanks.

...facing the exact same problem now

What problem? According to the title of the thread, the Original Poster wanted, somehow, to use a Fourier Series to "reproduce the original curve." He wants to use realft.f90 to get the coefficients of the series.

Well, first of all, we don't have a "curve"; we have a set of data points.

Secondly, in general, not all "curves" have Fourier Series representations. (That is, the infinite series may not exist or may not converge.)

Thirdly, the Original Poster showed a data set of 245 points, not the 256 that he claimed, and realft only works with a number of points that is a power of 2.

Fourthly, the Original Poster wanted to know why some data manipulations with some spreadsheet formulas didn't give satisfactory results. (How the heck could we explain something like that when we don't have any access to whatever the heck he was using? I mean, guess? Or what?)

My conclusion: IF you have the exact same problem (whatever it is), then I kind of doubt that your chances of getting a helpful explanation are better than those of the Original Poster (almost three years ago).

However...

If you want to know how to use realft.f90 to get some numbers that can be used to evaluate a finite number of terms of a Fourier Series associated with a set of real data points, I'll show something to consider:

Suppose a(NP) is an array of reals that holds your NP data points.

In order to use realft, the NP must be a power of 2.

Call realft(a,+1)

When realft returns, the array holds the DFT coefficients. Now, the (two-sided) DFT actually results in (NP+1) complex coefficients, but because of symmetry conditions on the DFT of a set of real data, we will only need (NP/2+1) complex numbers, and that's what realft returns.

DFT coefficients number zero and number NP/2 are actually real numbers (imaginary parts equal to zero), so they are stored in a special way, as explained in the text.

If you want to evaluate a certain number of terms of the Fourier Series, here's a way, after calling realft:

Let OMEGA be equal to 2*pi

a(1) contains the (real-valued) coefficient for for the constant term. (The coefficient associated with zero times OMEGA.)

a(2) contains the (real-valued) coefficient associated with NP/2 * OMEGA

The other elements of a, taken in pairs, represent the real and imaginary parts of coefficients associated with the other multiples of OMEGA:

a(3) contains the real part of the (complex) coefficient associated with 1*OMEGA
a(4) contains the imaginary part of the (complex) coefficient associated with 1*OMEGA

a(5) contains the real part of the (complex) coefficient associated with 2*OMEGA
a(6) contains the imaginary part of the (complex) coefficient associated with 2*OMEGA

and so on, up to a(NP-1) and a(NP), the coefficients associated with (NP/2-1) * OMEGA

Here's a way to evaluate up to NP/2 terms of the Fourier Series for a given value of x (where x is greater than or equal to zero and less than 1). I didn't bother to use a(2) the coefficient associated with (NP/2) * OMEGA. You can add it in if you want to use all of the terms. (It gets multiplied by 2 also.)

First of all, divide all a(i) by NP. We can do it once and for all before calling the evaluation routine the first time.

Then here's the evaluation loop in Fortran:


!
! This calculates the Fourier Series approximation
! for a given value of x
!
! Data array values from realft have been divided by npoints
!
! davekw7x
!
FUNCTION fs(data,x,npoints,nterms)
USE nrtype
IMPLICIT NONE
REAL(SP),DIMENSION(:),INTENT(IN) :: data
REAL(SP),INTENT(IN) :: x
INTEGER(I4B),INTENT(IN) :: nterms,npoints
REAL(SP) :: fs
INTEGER(I4B) :: i
INTEGER(I4B) :: nt
REAL(SP),PARAMETER :: OMEGA=6.283185307_sp

! The array contains half of the number of double-sided coefficient terms.
! We count them twice. This works because of symmetry properties of
! the DFT of a real-valued function
!
if (nterms > npoints/2) then
nt = npoints/2
else
nt = nterms
endif

fs = 0

!Add the terms after the DC term
do i=3,2*nt,2
fs = fs + data(i)*cos(((i-1)/2)*OMEGA * x) + &
data(i+1)*sin(((i-1)/2)*OMEGA * x)
enddo

! The DC term appears once; the others count twice
fs = data(1) + 2.0*fs

END FUNCTION fs


Note:
As I mentioned, the Original Poster's data array didn't have 256 points as he claimed, but only 245 points. I'll use his first 128 data points so that realft can do the deed.

I have attached some plots that show the results of evaluating the Fourier Series using using the above function with 10 terms and with 64 terms.

In each plot, the blue line connects the approximate values calculated from the Fourier Series coefficients and the red line connects the original data points.

The 10-term result shows the effects of deleting frequencies above 9 times Omega.

The Original Poster wanted to see how to get original data back, and the 64-term result uses almost as much information as you can get from the data. (The ultimate would use the a(2) coefficient multiplied by 2 times the cosine of the appropriate multiple of OMEGA, but I kind of doubt that there would be any obvious difference on the plots.)

Regards,

Dave