|
Posted by Steve Bickerton on October 7, 2006, 4:41 pm
Please log in for more thread options
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;
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;
}
##########################################
|
|
Posted by Steve Bickerton on October 9, 2006, 2:25 pm
Please log in for more thread options
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
|
|
Posted by harryfmudd [AT] comcast [DOT] on October 9, 2006, 7:28 pm
Please log in for more thread options
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
|
|
Posted by Dr.Ruud on October 9, 2006, 9:20 pm
Please log in for more thread options
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."
|
| Similar Threads | Posted | | Net::Analysis Parse tcpdump for HTTP Request/Response Headers | July 29, 2007, 6:41 am |
| GMP::Math ...PLEASE HELP!! | February 7, 2005, 11:02 am |
| Math::TrulyRandom | May 22, 2006, 9:56 pm |
| Help with Math::GMP compile on AIX 5.3 | September 11, 2007, 4:47 pm |
| Math-Random-MT-Auto-4.08.00.tar.gz | September 18, 2005, 12:16 am |
| Math::Macopt version 0.01 available | January 9, 2006, 9:21 pm |
| A problem with Math::Prime::XS | February 14, 2006, 12:35 pm |
| Math::Pari headaches | August 15, 2006, 1:50 pm |
| Math::Pari in turmoil | September 21, 2008, 11:42 pm |
| Memory Leak in Math::Pari | December 16, 2004, 6:21 am |
|
XML Sitemap