time-series analysis using Math::FFT
am 07.10.2006 22:41:33 von Steve BickertonHi 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;
}
##########################################