A tiny problem ...

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

A tiny problem ...

dieter.jurzitza@t-online.de
Dear listmembers,
I am using a peace of code as shown below:

"signalfilt" is a low pass filtered vector of 44100 complex elements, "noise"
a similar but uncorrelated one of the same size.

"rmsval" computes the rms value of an arbitrary (complex) input vector

****************** SNIPP **********************
NUMEL=5000;

atten=zeros(1,NUMEL);
b=signalfilt+noise;
b=b/rmsval(b);

for i=1:NUMEL
  b=shift(b,1,2);
  atten(i)=rmsval(signalfilt-b);
endfor
****************** SNAPP **********************

What I actually do here is some kind of correlation measurement I do not plan
to dive too deep into right now. The loop takes about 6s on my computer given
the dimensions of the vectors and numbers as listed above.

The shift operation is the most costly one in the loop, say about 10% longer
than the rmsval of the vector difference below.

If I had to implement this in "C" I would put b two times one behind the other
into memory, rather than actually shifting I would only increment the start
address pointing to that new "b" as used within the loop by one for each loop,
thereby avoiding any real shift at all. This way I could replace a shift of
44100 elements by one single pointer increment.

However, I cannot see a way to do this in octave, as there are no pointers
AFAIK but maybe one much more experienced user in the list has a good idea for
me.

By the way, I am using octave 4.2, but this should not have a significant
impact IMHO.

Thank you very much for your efforts,
regards




Dieter Jurzitza
--
-----------------------------------------------------------
Dr.-Ing. Dieter Jurzitza                    76131 Karlsruhe


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: A tiny problem ...

Doug Stewart-4


On Mon, May 1, 2017 at 3:57 PM, Dr.-Ing. Dieter Jurzitza <[hidden email]> wrote:
Dear listmembers,
I am using a peace of code as shown below:

"signalfilt" is a low pass filtered vector of 44100 complex elements, "noise"
a similar but uncorrelated one of the same size.

"rmsval" computes the rms value of an arbitrary (complex) input vector

****************** SNIPP **********************
NUMEL=5000;

atten=zeros(1,NUMEL);
b=signalfilt+noise;
b=b/rmsval(b);

for i=1:NUMEL
  b=shift(b,1,2);
  atten(i)=rmsval(signalfilt-b);
endfor
****************** SNAPP **********************

What I actually do here is some kind of correlation measurement I do not plan
to dive too deep into right now. The loop takes about 6s on my computer given
the dimensions of the vectors and numbers as listed above.

The shift operation is the most costly one in the loop, say about 10% longer
than the rmsval of the vector difference below.

If I had to implement this in "C" I would put b two times one behind the other
into memory, rather than actually shifting I would only increment the start
address pointing to that new "b" as used within the loop by one for each loop,
thereby avoiding any real shift at all. This way I could replace a shift of
44100 elements by one single pointer increment.

However, I cannot see a way to do this in octave, as there are no pointers
AFAIK but maybe one much more experienced user in the list has a good idea for
me.

By the way, I am using octave 4.2, but this should not have a significant
impact IMHO.

Thank you very much for your efforts,
regards




Dieter Jurzitza
--
-----------------------------------------------------------
Dr.-Ing. Dieter Jurzitza                    76131 Karlsruhe


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave

try this
a=1:7
b=[0 0 a]
is that what you want?


--
DASCertificate for 206392


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: A tiny problem ...

Kire Pudsje


On Mon, May 1, 2017 at 10:03 PM, Doug Stewart <[hidden email]> wrote:


On Mon, May 1, 2017 at 3:57 PM, Dr.-Ing. Dieter Jurzitza <[hidden email]> wrote:
Dear listmembers,
I am using a peace of code as shown below:

"signalfilt" is a low pass filtered vector of 44100 complex elements, "noise"
a similar but uncorrelated one of the same size.

"rmsval" computes the rms value of an arbitrary (complex) input vector

****************** SNIPP **********************
NUMEL=5000;

atten=zeros(1,NUMEL);
b=signalfilt+noise;
b=b/rmsval(b);

for i=1:NUMEL
  b=shift(b,1,2);
  atten(i)=rmsval(signalfilt-b);
endfor
****************** SNAPP **********************

What I actually do here is some kind of correlation measurement I do not plan
to dive too deep into right now. The loop takes about 6s on my computer given
the dimensions of the vectors and numbers as listed above.

The shift operation is the most costly one in the loop, say about 10% longer
than the rmsval of the vector difference below.

You seem to care about performance. Like you say, you are doing a kind of correlation. Then why not actually do some math and rewrite your equation to perform a cross covariance using FFT's. You will get the modulo-N for free and perform all covariances in N*log(N) time. When performing the math, you will only be left with the square root of the cross covariance plus the auto correlation of the noise.
I am not sure as to which of the covariance functions in octave use the fft (ie have modulo-N behavior), but it is easy enough to write by hand using the fft functions yourself.



_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: A tiny problem ...

dieter.jurzitza@t-online.de
Hello Kire,
hello Doug,
thank you for the feedback and for taking the time to reply.

It does not look like I have overlooked a vectorizing option in the code I
wrote, though - according to your input Kire - there might be a way to do this
in the frequency domain what could be faster because the number of iterations
could be reduced.

Honestly speaking I did not understand what Doug was suggesting ...

But anyway, the code runs - though a little slower than it could probably be
running, and the results I get are the way I expected them to be.
Regards




Dieter Jurzitza
Am Montag, 1. Mai 2017, 23:34:33 schrieb Kire Pudsje:
> On Mon, May 1, 2017 at 10:03 PM, Doug Stewart <[hidden email]> wrote:
> > On Mon, May 1, 2017 at 3:57 PM, Dr.-Ing. Dieter Jurzitza <
> >
> > [hidden email]> wrote:
***************

--
-----------------------------------------------------------
Dr.-Ing. Dieter Jurzitza                    76131 Karlsruhe


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: A tiny problem ...

Kire Pudsje


On Wed, May 3, 2017 at 7:58 PM, Dr.-Ing. Dieter Jurzitza <[hidden email]> wrote:
Hello Kire,
hello Doug,
thank you for the feedback and for taking the time to reply.

It does not look like I have overlooked a vectorizing option in the code I
wrote, though - according to your input Kire - there might be a way to do this
in the frequency domain what could be faster because the number of iterations
could be reduced.

Honestly speaking I did not understand what Doug was suggesting ...

But anyway, the code runs - though a little slower than it could probably be
running, and the results I get are the way I expected them to be.
Regards
which shows it for autocorrelation. From there it is a small step to covariance.
There are two tricks involved:
A convolution in time domain is just a multiplication in the frequency domain.
A time shift will become an exponential in the frequency domain.


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave