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 zeros to add
    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
				dropBack(residual);//Drop the additional zeros
		
    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?

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
> zeros to add
>      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
>                  dropBack(residual);//Drop the additional zeros
>
>      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?

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?

I found this link below.
"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

After reading the link I think you are right,

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. I will think about it

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
> zeros to add
>      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
>                  dropBack(residual);//Drop the additional zeros
>
>      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.

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.)

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 http://dlang.org/phobos/std_complex.html#.sqAbs instead. You then also don't need the final map to extract the real part.

You are absulately right instead
fftResult.zip(fftResult.map!(a => a * a.conj()))
I corrected it to
fftResult.zip(fftResult.map!(a => a.conj()))

Thanks a lot