PERL Modules - time-series analysis using Math::FFT

Author

time-series analysis using Math::FFT

Steve Bickerton

2006-10-07, 7:00 pm

Hi All,

I'm trying to generate a real-valued time-series with a power-spectrum
having a slope of -1 (in log-log) (called 1 over f noise). So, I start
in the frequency domain and define the real and imaginary spectral
components with the slope I want, and I use the Math::FFT invcdft()
function to invert that complex spectrum to the time domain. It works
quite well, but so far I've been generating complex time-series and
tossing away the imaginary part. What I should be doing is generating a
real-valued time-series by setting the negative frequency components to
be the complex conjugates of their positive counterparts (conjugate
symmetry). When I do this, although the resulting time-series should be
entirely real-valued, I get a time-series that still has an imaginary
part. Something I'm doing is wrong, but I don't know what. I've
appended a subroutine that I use to do this in the hope that someone
more familiar with the module than myself will spot the problem.

I'm suspicious of the way that I've filled the array containing the
spectral components. I've put the real and imaginary parts for the
negative frequencies in the upper half of the array in an alternating
sequence; which, to the best of my knowledge is the way to do it.

I apologize for posting a lengthy snippit of code for what might well be
a somewhat involved problem. Any suggestions anyone might have would be
welcome.

regards,
steve

########################################
#########################
my $number = '[+-]?\d+\.?\d*';

sub makeNoise ($$$$) {

my $usage = "Use: \$ts_ref = makeNoise(\$beta, \$dt, \$points,
\$rms);\n".
" \$ts_ref = \\\@ts with [t, real, f, power] as entries\n";

my ($beta, $dt, $Nreq, $rms) = @_;
croak $usage unless ( $beta =~ /^$number$/ &&
$dt =~ /^$number$/ &&
$Nreq =~ /^\d+$/ &&
$rms =~ /^$number$/ );

# make sure N is a power of 2
my $log2Nreq = log2($Nreq);
my $lowGuess = int($log2Nreq);
$Nreq = 2**($lowGuess+1) if ( $lowGuess != $log2Nreq );
my $N = 2*$Nreq;

# generate the fake power spectrum
my @FT = (0)x(2*$N);

my ($mean, $std) = (0.0, 1.0);
my @gaussian = random_normal($N, $mean, $std);
my @phase = random_uniform($N, 0, $TWOPI);

my $T = $N*$dt;
my $df = 1.0 / $T;
my $dw = $TWOPI*$df;
for (my $k=0; $k<$N/2; $k++) {

my $f = $k*$df;
my $mag = ($k) ? (($f)**(-$beta/2.0)*$gaussian[$k] ) : (0.0);

# assign the positive real and imaginary parts
$FT[2*$k] = $mag * cos($phase[$k]);
$FT[2*$k + 1] = $mag * sin($phase[$k]);

# assign the negative real and imaginary parts
$FT[2*$N - 2*$k - 2] = $FT[2*$k];
$FT[2*$N - 2*$k - 1] = -$FT[2*$k + 1];

}

# invert to get the time-series
my $cspec = new Math::FFT(\@FT);
my $cts = $cspec->invcdft(\@FT);

# get the stats to normalize the amplitude to an RMS of $rms
my @real;
for (my $k=0; $k<$N; $k++) { push @real, $cts->[2*$k]; }
my $stdev = rms(\@real);
my $scale = $rms/$stdev;

# return the time-series and the spectrum
my @timeseries;
for (my $k=0; $k<$N; $k++) {

my $t = $k*$dt;
my $I = 1.0 + $scale*$cts->[2*$k];
my $f = $k*$df;
my $power = $scale*($FT[2*$k]**2 + $FT[2*$k+1]**2);
push @timeseries, [ $t, $I, $cts->[2*$k], $cts->[2*$k+1],
$f, $power, $FT[2*$k], $FT[2*$k+1] ];
}

return \@timeseries;
}
########################################
##

Steve Bickerton

2006-10-12, 6:58 pm

Hi All,
I found the problem. It was an indexing error. If anyone's
interested, the following solution was required:

> #####
> for (my $k=0; $k<$N/2; $k++) {
>
> ...
> # assign the negative real and imaginary parts
> $FT[2*$N - 2*$k - 2] = $FT[2*$k];
> $FT[2*$N - 2*$k - 1] = -$FT[2*$k + 1];
>
> }
> #####

That code should have been:

#####
for (my $k=1; $k<$N/2; $k++) {

...
$FT[2*$N - 2*$k] = $FT[2*$k];
$FT[2*$N - 2*$k + 1] = -$FT[2*$k + 1];

}
@FT[0,1, $N, $N+1] = (0,0,0,0);
#####

Sorry to post prematurely.
steve

harryfmudd [AT] comcast [DOT] net

2006-10-12, 6:58 pm

Steve Bickerton wrote:
> Hi All,
> I found the problem. It was an indexing error. If anyone's
> interested, the following solution was required:

<snip content="Obi-wan error" />

>
> Sorry to post prematurely.
> steve

Thank you for posting the solution for the record.

Tom Wyant

Dr.Ruud

2006-10-12, 6:58 pm

Steve Bickerton schreef:

> for (my $k=1; $k<$N/2; $k++) {

Alternative:

for my $k (1 .. $N/2-1) {

--
Affijn, Ruud

"Gewoon is een tijger."

Sponsored Links