www.digitalmars.com         C & C++   DMDScript  

digitalmars.D.learn - Autocorrelation function with ranges

reply "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 
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?
Jun 27 2015
next sibling parent reply 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
 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?
Jun 27 2015
parent "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?
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
Jun 27 2015
prev sibling parent reply 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
 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.
Jun 27 2015
next sibling parent 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
prev sibling parent "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 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
Jun 27 2015