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
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
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
"kerdemdemir" <kerdemdemir hotmail.com> 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
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?instead. You then also don't need the final map to extract the real part.
Jun 27 2015
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
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?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 lotAlso is there any performance issues? can I make this faster?Probably you should use then also don't need the final map to extract the real part.
Jun 27 2015