## digitalmars.D.learn - Autocorrelation function with ranges

• kerdemdemir (31/31) Jun 27 2015 Hi
• Rikki Cattermole (4/31) Jun 27 2015 No idea about the maths behind it are but:
• kerdemdemir (22/25) Jun 27 2015 Thanks a lot for your answer anyway. I am hoping even not related
• Timon Gehr (7/35) Jun 27 2015 One obvious problem is this:
• Timon Gehr (3/4) Jun 27 2015 (This is actually not true, your inverseFFT presumably still returns
• kerdemdemir (6/13) Jun 27 2015 You are absulately right instead
"kerdemdemir" <kerdemdemir hotmail.com> writes:
```Hi

My question is more about Maths than D lang,
I am hoping, maybe somebody worked with AutoCorrelation function
before.

auto autoCorrelation(R)(R range)
if (isRandomAccessRange!R)
{
auto residual = residualPowerOf2(range.length); // Find how
auto fftResult = range.chain(repeat(0, residual)).fft(); //
Takes FFT

//First zip each element of fft with conjagute of fft
auto autoCorrResult = fftResult.zip(fftResult.map!(a => a *
a.conj())).
map!( a=> a[0] * a[1] ). // Than multiple them
inverseFft(). // Than take inverse

auto finalResult = autoCorrResult.take(1). // Take DC element
chain(autoCorrResult[\$/2..\$]).//Take last half to
beginning
chain(autoCorrResult[1..\$/2]). // First(negative lags)
to end
map!(a => a.re); // I just need real part

return finalResult ;
}

My autocorrelation method return some crazy values(finalResult[0]
= -101652). I am mostly suspicious about calculations to set
"finalResult"  variable.

Also is there any performance issues? can I make this faster?
```
Jun 27 2015
Rikki Cattermole <alphaglosined gmail.com> writes:
```On 27/06/2015 10:29 p.m., kerdemdemir wrote:
Hi

My question is more about Maths than D lang,
I am hoping, maybe somebody worked with AutoCorrelation function before.

auto autoCorrelation(R)(R range)
if (isRandomAccessRange!R)
{
auto residual = residualPowerOf2(range.length); // Find how many
auto fftResult = range.chain(repeat(0, residual)).fft(); // Takes FFT

//First zip each element of fft with conjagute of fft
auto autoCorrResult = fftResult.zip(fftResult.map!(a => a *
a.conj())).
map!( a=> a[0] * a[1] ). // Than multiple them
inverseFft(). // Than take inverse

auto finalResult = autoCorrResult.take(1). // Take DC element
chain(autoCorrResult[\$/2..\$]).//Take last half to
beginning
chain(autoCorrResult[1..\$/2]). // First(negative lags)
to end
map!(a => a.re); // I just need real part

return finalResult ;
}

My autocorrelation method return some crazy values(finalResult[0] =
-101652). I am mostly suspicious about calculations to set
"finalResult"  variable.

Also is there any performance issues? can I make this faster?

No idea about the maths behind it are but:

chain(autoCorrResult[1..\$/2])

Should that one be a zero?
```
Jun 27 2015
"kerdemdemir" <kerdemdemir hotmail.com> writes:
```On Saturday, 27 June 2015 at 10:37:08 UTC, Rikki Cattermole wrote:
No idea about the maths behind it are but:

Thanks a lot for your answer anyway. I am hoping even not related
with D directly, this discussions may atract people from other
languages to D while looking for Domain information.

chain(autoCorrResult[1..\$/2])

Should that one be a zero?

"http://www.aip.de/groups/soe/local/numres/bookcpdf/c13-2.pdf"
Which says :
The correlation at zero lag is in r0, the first component; the
correlation at lag 1 is in r1, the second component; the
correlation at lag −1 is in rN−1, the last component; etc.

Correlation result after IFFT is like :

0 1 2 3 ..... T -1 -2 -3 .... -T

How I wanted to be :

-T -T+1 .....    -1 0 1 .... T-1 T

auto finalResult =  chain(autoCorrResult[\$/2..\$]).
chain(autoCorrResult[0..\$/2]).
map!(a => a.re);

Example above should be the way how I transform. Also now I see
inverseFft().
dropBack(residual); ==> DropBack may not be a good idea here.
```
Jun 27 2015
Timon Gehr <timon.gehr gmx.ch> writes:
```On 06/27/2015 12:29 PM, kerdemdemir wrote:
Hi

My question is more about Maths than D lang,
I am hoping, maybe somebody worked with AutoCorrelation function before.

auto autoCorrelation(R)(R range)
if (isRandomAccessRange!R)
{
auto residual = residualPowerOf2(range.length); // Find how many
auto fftResult = range.chain(repeat(0, residual)).fft(); // Takes FFT

//First zip each element of fft with conjagute of fft
auto autoCorrResult = fftResult.zip(fftResult.map!(a => a *
a.conj())).
map!( a=> a[0] * a[1] ). // Than multiple them
inverseFft(). // Than take inverse

auto finalResult = autoCorrResult.take(1). // Take DC element
chain(autoCorrResult[\$/2..\$]).//Take last half to
beginning
chain(autoCorrResult[1..\$/2]). // First(negative lags)
to end
map!(a => a.re); // I just need real part

return finalResult ;
}

My autocorrelation method return some crazy values(finalResult[0] =
-101652). I am mostly suspicious about calculations to set
"finalResult"  variable.
...

One obvious problem is this:

fftResult.zip(fftResult.map!(a => a * a.conj())).map!(a=>a[0]*a[1]).

This computes a²·a̅ instead of a·a̅.

What is the source code for residualPowerOf2?

Also is there any performance issues? can I make this faster?

Probably you should use http://dlang.org/phobos/std_complex.html#.sqAbs
instead. You then also don't need the final map to extract the real part.
```
Jun 27 2015
Timon Gehr <timon.gehr gmx.ch> writes:
```On 06/27/2015 02:17 PM, Timon Gehr wrote:
You then also don't need the final map to extract the real part.

(This is actually not true, your inverseFFT presumably still returns
complex numbers.)
```
Jun 27 2015
"kerdemdemir" <kerdemdemir hotmail.com> writes:
```On Saturday, 27 June 2015 at 12:17:31 UTC, Timon Gehr wrote:

This computes a²·a̅ instead of a·a̅.

What is the source code for residualPowerOf2?

Also is there any performance issues? can I make this faster?

Probably you should use