Re: date-compensated fourier transform

[Rob de Graaf <rdegraaf@rdegraaf.nl>]

Hi Justin,

you can use the following script using octave. 
The data input is a file and does a fourier transform and pipes it out.

Regards,

Rob de Graaf

[ code ]

#!/usr/bin/octave -qf
# dftdata var from to

# read data
fi = popen(sprintf("cat %s | grep -v '#' ", argv(1,:) ),"r");  # open pipe with a cat and filter lines on '#'
d = fscanf(fi,"%f",[2 Inf])';                                  # read the content into matrix,
                                                               #   notice the transpose!!!
pclose(fi);                                                    # close pipe
n = size(d,1);                                                 # determine size
if n < 3                                                       # check size
  if size(d,2) == 0
    n = 0;
  endif
  fprintf(stderr,"%d data points - not enough!\n",n);
  exit(1);
endif

# transform parameters
dt = d(2:n,1)-d(1:n-1,1);
mdt = mean(dt);
maxfreq = 1/(2*mdt);
nf = floor(n/2);
df = maxfreq/nf;
T = 1/df; # period for the transformed function

# remove points beyond T
# we rely that time for the 1st point is 0: 'get_data' property
j = find(d(:,1) < T);
tim = d(j,1); # time coordinate
val = d(j,2); # value
val = val - mean(val); # remove constant
# make source function periodic by repeating 1st point at T
tim = [tim; T];
val = [val; val(1)];

# prepare weights for trapezoidal integration
n = size(tim,1);
dt = tim(2:n)-tim(1:n-1);
w = zeros(n,1);
w(2:n-1) = (dt(1:n-2)+dt(2:n-1))/2;
w(1) = dt(1)/2;
w(n) = dt(n-1)/2;
wval = w.*val;

# prepare storage for result
# we omit f = 0;
freq = (1:nf)'*df;
ftres = zeros(nf,1);

# do actual integration
for j = 1:nf
  x = 2*pi*freq(j)*tim;
  cosum = sum(wval.*cos(x));
  sisum = sum(wval.*sin(x));
  ftres(j) = sqrt(cosum*cosum+sisum*sisum);
endfor

# normalize
ftres = ftres/T;

# output
printf("%g %g\n",[freq ftres]');

[ code ]





On Thu, 17 Nov 2005 01:27:48 -0500
Justin Pryzby <justinpryzby@users.sourceforge.net> wrote:

> Does anyone have experience with date-compensated fourier transforms?
> This is a transform where the input data is spaced unevenly in time.
> I've glanced at some papers, but haven't seen any library code that
> implements this.  I'm interested in finding such code, but it would
> also be nice to understand the algorithm in intuitive terms.
> 
> -- 
> Clear skies,
> Justin
> 
> 
> -- 
> To UNSUBSCRIBE, email to debian-science-request@lists.debian.org
> with a subject of "unsubscribe". Trouble? Contact listmaster@lists.debian.org
>