PPM MarkIV project progress report

[ApBerg]

Hi Ap,

A long time since we talked together.
If you are still ON to try and build a PPM yourself, you could get my support.
Tell me why you did not yet achieve to make your PPM to work and I'll help you.
Cheers,
Willy

Hello Willy,

Thanks for your replay !

Yes, now and then I am trying to get some proton noise.. , it is very intriguing, the mystery from the earth magnetism’s and of course the possibility to detect magnetic variations !

Think that the problems I have is mostly the lack to have a place where is no 50Hz noise.
Have been training with a tone generator, to make 2Khz signal, and so simulating the proton signal. This has some small success. I wonder how the signal must look when it leaves the last opamp, think it has to be a squire wave ?

Made a pcb that has audio out….for the fun to hear the signal ..( not fully working )

Best regards.

Ap

[Willy Bayot]

Hello Konstantin.

Is it the main.c file that is with the Koehler project in CircuitCeller May 2007 ?
Have you made that project ?
Last month I bought the CC backup CD from 2007 and have read about that project..
Nice info.. as al the info from Jim Koehler is, also the info (and project) from Willy Bayot is great.

But with all that info I cannot yet get a proton mag running...to bay a mag from Willy ( I sure like that ) is for me not an option... no budget (wife kills me ....yes she told me..)

I sure would like a diy project that is on a lower budget...with a Microchip that knows what to do..

Best regards

Ap

Hi Ap,

A long time since we talked together.
If you are still ON to try and build a PPM yourself, you could get my support.
Tell me why you did not yet achieve to make your PPM to work and I'll help you.
Cheers,
Willy

[Willy Bayot]

I tested it with two compilers (GCC and Keil), on three different platforms (ARM9, CortexM3, CortexM4). Also, I tested it in Matlabe. Tested with test and real signals. For example, you can get a signal from the computer sound card. There are many programs that give a signal with the desired SNR. You could check the accuracy of the algorithm in this way?

We have done all these things a long time ago (about 4 years ago).
Jim is a Doctor in Physics and owns all the necessary equipment to make all the simulations and lab tests. There must be a mistake somewhere in our understanding of the algorithm that we called 'PHASE SLIP' algorithm.
If you wish, I am ready to talk privately in order to remove this misunderstanding.

Willy

[KonstantinT]

I do not know whether this deeply modified code of the same algorithm works but I know very well that the original algorithm worked.
It is true that this code here is looking more optimized in terms of C writing but we happily count of the good optimizations applied automatically by the Keil compiler.
Changing the original code as heavily will always give some risks of bad conversion.
Willy

I tested it with two compilers (GCC and Keil), on three different platforms (ARM9, CortexM3, CortexM4). Also, I tested it in Matlabe. Tested with test and real signals. For example, you can get a signal from the computer sound card. There are many programs that give a signal with the desired SNR. You could check the accuracy of the algorithm in this way?

This algorithm works, but gives the worst results in terms of accuracy compared with other algorithms. Achieved accuracy of 1,3-1,5 nT at SNR 30-40 dB.

You can't show your sources. I just wanted to know - my code is like you code or not? I have too much time spent on the investigation of this algorithm and want to do for himself the final conclusions.

BR
Konstantin

[Willy Bayot]

This is our piece of code of the Koehler's algorithm.

for(j=0; j<calcLength; j++) {

fS = (float)(*(uk+0))-(float)(*(uk+;
fC = (float)(*(uk+4))-(float)(*(uk+12));

if(fC == 0) {

if(fS>=0)
*(arctangents+j) = (float)PI_OVER_TWO;
else
*(arctangents+j) = (float)-PI_OVER_TWO;

} else
*(arctangents+j) = (float)atan2((float)fS,(float)fC);


fS = (float)(*(uk+1)+50)-(float)(*(uk+9));
fC = (float)(*(uk+5)+50)-(float)(*(uk+13));


if(fC == 0) {

if(fS>=0)
*(arctangents+1*calcLength+j) = (float)PI_OVER_TWO;
else
*(arctangents+1*calcLength+j) = (float)-PI_OVER_TWO;

} else
*(arctangents+1*calcLength+j) = (float)atan2((float)fS,(float)fC);

fS = (float)(*(uk+2))-(float)(*(uk+10));
fC = (float)(*(uk+6))-(float)(*(uk+14));

if(fC == 0) {
if(fS>=0)
*(arctangents+2*calcLength+j) = (float)PI_OVER_TWO;
else
*(arctangents+2*calcLength+j) = (float)-PI_OVER_TWO;
} else
*(arctangents+2*calcLength+j) = (float)atan2((float)fS,(float)fC);

fS = (float)(*(uk+3))-(float)(*(uk+11));
fC = (float)(*(uk+7))-(float)(*(uk+15));

if(fC == 0) {
if(fS>=0)
*(arctangents+3*calcLength+j) = (float)PI_OVER_TWO;
else
*(arctangents+3*calcLength+j) = (float)-PI_OVER_TWO;
} else
*(arctangents+3*calcLength+j) = (float)atan2((float)fS,(float)fC);


uk = uk + 16;

}

for (j = 0; j < 4; j++) {
offset = 0.0;

for (i = 1; i < calcLength; i++) {
diff = *(arctangents+j*calcLength+i)- *(arctangents+j*calcLength+i-1);

if(diff > PI)
for (k = i; k < calcLength; k++)
*(arctangents+j*calcLength+k) -=TWO_PI;


if(diff < - PI)
for (k = i; k < calcLength; k++)
*(arctangents+j*calcLength+k) +=TWO_PI;
}

}




for (j = 0; j < 4; j++) {
sumx = 0.0;

sumy = 0.0;
sumxy = 0.0;
sumxsq = 0.0;

for (i = 0; i < calcLength; i++) {
sumx += (float) i;
sumy += *(arctangents+j*calcLength+i);
sumxy += (*(arctangents+j*calcLength+i))* (float) i;
sumxsq += (float) i * (float)i;
}

det = (calcLength * sumxsq) - (sumx * sumx);
slope[j] = ((calcLength * sumxy) - (sumy * sumx)) / det;
intercept[j] = ((sumy * sumxsq) - (sumxy * sumx)) / det;
}

sigma = 0.0;
for (j = 0; j < 4; j++) {
for (i = 0; i < calcLength; i++) {
sumx= ((intercept[j] + (slope[j] * (float)i)) - *(arctangents+j*calcLength+i));
sigma += sumx * sumx;
}
}


sigma /= (4.0 * calcLength) - 1.0;
sigma = sqrtf((float)sigma);
sumx = 0.0;

for (j = 0; j < 4; j++) sumx += slope[j];
sumx /= 4.0;
delta_f = sumx * sample_freq / TWO_PI;
precession_freq = sample_freq + delta_f;

It is no different from the code, published in the Philips contest. I tested it on 12, 16 and 24-bit ADC. It is a mistake or not? If yes, please and I'll check its accuracy.

BR
Konstantin

I do not know whether this deeply modified code of the same algorithm works but I know very well that the original algorithm worked.
It is true that this code here is looking more optimized in terms of C writing but we happily count of the good optimizations applied automatically by the Keil compiler.
Changing the original code as heavily will always give some risks of bad conversion.

However, I have still not yet understood from your various posts whether the original code was working at all (with or without the two 'corrections' you published)or if it was working but giving less accurate results than OTHER algorithms.
Could you please clarify this?

Willy

[KonstantinT]

That could very well be.
Did you REALLY try it yourself?
Willy

This is our piece of code of the Koehler's algorithm.

for(j=0; j<calcLength; j++) {

fS = (float)(*(uk+0))-(float)(*(uk+;
fC = (float)(*(uk+4))-(float)(*(uk+12));

if(fC == 0) {

if(fS>=0)
*(arctangents+j) = (float)PI_OVER_TWO;
else
*(arctangents+j) = (float)-PI_OVER_TWO;

} else
*(arctangents+j) = (float)atan2((float)fS,(float)fC);


fS = (float)(*(uk+1)+50)-(float)(*(uk+9));
fC = (float)(*(uk+5)+50)-(float)(*(uk+13));


if(fC == 0) {

if(fS>=0)
*(arctangents+1*calcLength+j) = (float)PI_OVER_TWO;
else
*(arctangents+1*calcLength+j) = (float)-PI_OVER_TWO;

} else
*(arctangents+1*calcLength+j) = (float)atan2((float)fS,(float)fC);

fS = (float)(*(uk+2))-(float)(*(uk+10));
fC = (float)(*(uk+6))-(float)(*(uk+14));

if(fC == 0) {
if(fS>=0)
*(arctangents+2*calcLength+j) = (float)PI_OVER_TWO;
else
*(arctangents+2*calcLength+j) = (float)-PI_OVER_TWO;
} else
*(arctangents+2*calcLength+j) = (float)atan2((float)fS,(float)fC);

fS = (float)(*(uk+3))-(float)(*(uk+11));
fC = (float)(*(uk+7))-(float)(*(uk+15));

if(fC == 0) {
if(fS>=0)
*(arctangents+3*calcLength+j) = (float)PI_OVER_TWO;
else
*(arctangents+3*calcLength+j) = (float)-PI_OVER_TWO;
} else
*(arctangents+3*calcLength+j) = (float)atan2((float)fS,(float)fC);


uk = uk + 16;

}

for (j = 0; j < 4; j++) {
offset = 0.0;

for (i = 1; i < calcLength; i++) {
diff = *(arctangents+j*calcLength+i)- *(arctangents+j*calcLength+i-1);

if(diff > PI)
for (k = i; k < calcLength; k++)
*(arctangents+j*calcLength+k) -=TWO_PI;


if(diff < - PI)
for (k = i; k < calcLength; k++)
*(arctangents+j*calcLength+k) +=TWO_PI;
}

}




for (j = 0; j < 4; j++) {
sumx = 0.0;

sumy = 0.0;
sumxy = 0.0;
sumxsq = 0.0;

for (i = 0; i < calcLength; i++) {
sumx += (float) i;
sumy += *(arctangents+j*calcLength+i);
sumxy += (*(arctangents+j*calcLength+i))* (float) i;
sumxsq += (float) i * (float)i;
}

det = (calcLength * sumxsq) - (sumx * sumx);
slope[j] = ((calcLength * sumxy) - (sumy * sumx)) / det;
intercept[j] = ((sumy * sumxsq) - (sumxy * sumx)) / det;
}

sigma = 0.0;
for (j = 0; j < 4; j++) {
for (i = 0; i < calcLength; i++) {
sumx= ((intercept[j] + (slope[j] * (float)i)) - *(arctangents+j*calcLength+i));
sigma += sumx * sumx;
}
}


sigma /= (4.0 * calcLength) - 1.0;
sigma = sqrtf((float)sigma);
sumx = 0.0;

for (j = 0; j < 4; j++) sumx += slope[j];
sumx /= 4.0;
delta_f = sumx * sample_freq / TWO_PI;
precession_freq = sample_freq + delta_f;

It is no different from the code, published in the Philips contest. I tested it on 12, 16 and 24-bit ADC. It is a mistake or not? If yes, please and I'll check its accuracy.

BR
Konstantin

[Willy Bayot]

In this forum, I found threads where people have complained that this code does not work.
Konstantin

That could very well be.
Did you REALLY try it yourself?
Willy

[KonstantinT]

I can assure that the published code of Jim works (it was working and tested with the hardware platform described in the contest) and it is the essentially same as our current algorithm. By the way, Jim and I are partners on that project from the beginning and we together designed this algorithm a few years ago. He chose to publish it in the context of that contest.
We also are using an FFT algorithm but only to evaluate the precession frequency of a local region just before tuning the sensor for that region. FFT is not accurate enough for a short data capture of 1/4 second.
We also store in memory a segment of the signal and process it.
The algorithm is indeed using a sampling frequency around 16 times the expected precession frequency. See '#define SAMPLES_PER_CYCLE 16'
This means that we are making 4 frequency measurements per signal period and we, at the time, were capturing about 1000 signal periods during the 500 msec period.
Willy

In this forum, I found threads where people have complained that this code does not work.
Konstantin

[Willy Bayot]

In their instruments, we used the signal sampling frequency is 16 times greater than the expected frequency of precession. Without our changes, the code did not work.

Our device produces a Fourier transform processing of the input signal and calculates the signal to noise ratio. I can publish the source code of Koehler's algorithm with which we worked. In turn, you too can publish your code for this algorithm and then we can compare them. Unfortunately, I never saw of your code.

Impulse noise simply removed with simple digital filters. In our implementation of the algorithm, we first write the entire signal into memory and then process it. This is the only difference from the published algorithm. If you do not want to publish in open access, you can send me a compiled library and I can check it on our hardware. Or give a source under DNA.

On the Internet you can find a lot of scientific work to estimate the frequency of a sinusoidal signal in white noise. For example, you can look at the site of IEEE.

In an attachment - one of the first papers on this subject, which can begin to explore.

I can assure that the published code of Jim works (it was working and tested with the hardware platform described in the contest) and it is the essentially same as our current algorithm. By the way, Jim and I are partners on that project from the beginning and we together designed this algorithm a few years ago. He chose to publish it in the context of that contest.
We also are using an FFT algorithm but only to evaluate the precession frequency of a local region just before tuning the sensor for that region. FFT is not accurate enough for a short data capture of 1/4 second.
We also store in memory a segment of the signal and process it.
The algorithm is indeed using a sampling frequency around 16 times the expected precession frequency. See '#define SAMPLES_PER_CYCLE 16'
This means that we are making 4 frequency measurements per signal period and we, at the time, were capturing about 1000 signal periods during the 500 msec period.
Willy

[KonstantinT]

Development the overhausers sensor is even more difficult task than creating a proton magnetometer. If I give for you detailed description, you are unlikely it to repeat. It took us more than two years before we got a working sample. For a start, we tried to copy the finished commercial sensors.

For example, for tuning the RF resonator are necessary skills in radio engineering. To select the desired radical - PhD in chemistry. For the manufacture of glass ampules - a good glass-blowing workshop. To remove oxygen from the solution - a chemical laboratory. Going this way you will get an increase in SNR by only 10-15 dB.

In an attachment - a photo's of the sensor prototype.

[KonstantinT]

These corrections are indeed useful but that does not change the accuracy of the algorithm.
This code segment is just required when the reference sampling frequency is too far from the frequency to be measured, thus, in the case of very large field gradients. In those exceptional cases, the slope has to be corrected because of quadrant overlapping. The original code is still working well even if the threshold of PI is more precise to detect those overlaps than the approximate value of 5.
This algorithm is the best we could design to accurately measure a noisy signal with spikes as it usually is in the real world. It uses all the possible data redundancy available in the signal over a period of 250msec, i.e. around 1000 periods and uses the average of 1000 individual measurements to calculate the final precession frequency.
It also calculates the quality factor of the measurement (standard deviation between the 1000 measurements) and produces it together with the net nT value.
Note that the published algorithm was using a sample period of 500msec and a polarization time of 5sec. We are now doing a net reading rate of one reading per second
Our system resists to a hot 220V 50Hz cable running at 50cm from the sensor.
The sensitivity of our whole system is proven to be sub-0.1nT in real survey practice.
The algorithm does not live by itself in the PPM. The performance of the whole system also depends very much upon its sensor technology, its circuit switching and analog stages.

However, we are always ready to learn and acquire more experience in order to improve our PPM systems.
Since we have made the effort to publish our algorithm as originally implemented in code (even if it has been much improved since then), we would be very glad to get in return more than just the sentence 'His algorithm is much inferior to other algorithms'.
If Konstantin knows a better algorithm, we would be listening to his suggestions and possibly, read and study the real code of a better algorithm.
We also are interested in the practical building of an Overhauser system. Any info available?

Thanks,

Willy

In their instruments, we used the signal sampling frequency is 16 times greater than the expected frequency of precession. Without our changes, the code did not work.

Our device produces a Fourier transform processing of the input signal and calculates the signal to noise ratio. I can publish the source code of Koehler's algorithm with which we worked. In turn, you too can publish your code for this algorithm and then we can compare them. Unfortunately, I never saw of your code.

Impulse noise simply removed with simple digital filters. In our implementation of the algorithm, we first write the entire signal into memory and then process it. This is the only difference from the published algorithm. If you do not want to publish in open access, you can send me a compiled library and I can check it on our hardware. Or give a source under DNA.

On the Internet you can find a lot of scientific work to estimate the frequency of a sinusoidal signal in white noise. For example, you can look at the site of IEEE.

In an attachment - one of the first papers on this subject, which can begin to explore.

[KonstantinT]

What do the corrections you make in the code do? Is it more accurate? Did the code not work without those changes?

With the SNR of 40db the code was not accurate enough? IF you feed a strong signal into the code, does it perform better? Is the problem of the ability of the algorithm to not be able to generate accurate results? Or is the code not immune enough to the noise of the proton signal to produce a better accuracy?

Looking at the video of the bayot mag, it seems to be more accurate than you are claiming it is...using a .5 second sample time.

I've built a pretty nice proton mag. Not as many bells and whistles as Willy's but performs as well. I'd be interested in building an Overhauser mag. Do you sell sensors? a schematic? a kit?
Thanks,
SM

Without these changes the code does not work.

When SNR 30-40 dB is the ultimate precision. When a signal is applied to the SNR 50-60 db accuracy increased to 0,3 nT. Under such conditions, other algorithms work much better. On the other hand, it is quite resistant to induced noise. He also works with rectangular signals.

Maybe Bayout code changed? Maybe he uses averaging?

We compared our hardware and commercial magnetometers in the field. Koehler's algorithm showed significantly worse results. With our algorithm, we achieved much better accuracy. Therefore, the problem is not hardware.

It makes no sense to build overhauser's sensor until you have achieved in the proton sensor accuracy 0,1 nT.

Separately overhauser's sensors, we do not sell. Only by working together with our instruments. I can advise you if you have any questions.

BR
Konstantin

[Willy Bayot]

These corrections are indeed useful but that does not change the accuracy of the algorithm.
This code segment is just required when the reference sampling frequency is too far from the frequency to be measured, thus, in the case of very large field gradients. In those exceptional cases, the slope has to be corrected because of quadrant overlapping. The original code is still working well even if the threshold of PI is more precise to detect those overlaps than the approximate value of 5.
This algorithm is the best we could design to accurately measure a noisy signal with spikes as it usually is in the real world. It uses all the possible data redundancy available in the signal over a period of 250msec, i.e. around 1000 periods and uses the average of 1000 individual measurements to calculate the final precession frequency.
It also calculates the quality factor of the measurement (standard deviation between the 1000 measurements) and produces it together with the net nT value.
Note that the published algorithm was using a sample period of 500msec and a polarization time of 5sec. We are now doing a net reading rate of one reading per second
Our system resists to a hot 220V 50Hz cable running at 50cm from the sensor.
The sensitivity of our whole system is proven to be sub-0.1nT in real survey practice.
The algorithm does not live by itself in the PPM. The performance of the whole system also depends very much upon its sensor technology, its circuit switching and analog stages.

However, we are always ready to learn and acquire more experience in order to improve our PPM systems.
Since we have made the effort to publish our algorithm as originally implemented in code (even if it has been much improved since then), we would be very glad to get in return more than just the sentence 'His algorithm is much inferior to other algorithms'.
If Konstantin knows a better algorithm, we would be listening to his suggestions and possibly, read and study the real code of a better algorithm.
We also are interested in the practical building of an Overhauser system. Any info available?

Thanks,

Willy

[signman]

Questions about corrections

In an attachment to my previous post are the source of the Jim Koehler's project from Philips contest. I'm not made this project. We produce our proton magnetometers. Also, we have launched the production of our overhauser sensors. The price range for our devices - from $ 2200 to $ 8000.

We had high hopes for Jim Koehler's algorithm, because we wanted to get the accuracy of 0,1 nT with proton sensors. But unfortunately, it does not get the required accuracy with the standard signal-to-noise ratio of 30-40 dB . So we had to develop their own algorithms for computing the precession frequency.

BR
Konstantin

What do the corrections you make in the code do? Is it more accurate? Did the code not work without those changes?

With the SNR of 40db the code was not accurate enough? IF you feed a strong signal into the code, does it perform better? Is the problem of the ability of the algorithm to not be able to generate accurate results? Or is the code not immune enough to the noise of the proton signal to produce a better accuracy?

Looking at the video of the bayot mag, it seems to be more accurate than you are claiming it is...using a .5 second sample time.

I've built a pretty nice proton mag. Not as many bells and whistles as Willy's but performs as well. I'd be interested in building an Overhauser mag. Do you sell sensors? a schematic? a kit?
Thanks,
SM

[KonstantinT]

Hello Konstantin.

Is it the main.c file that is with the Koehler project in CircuitCeller May 2007 ?
Have you made that project ?
Last month I bought the CC backup CD from 2007 and have read about that project..
Nice info.. as al the info from Jim Koehler is, also the info (and project) from Willy Bayot is great.

But with all that info I cannot yet get a proton mag running...to bay a mag from Willy ( I sure like that ) is for me not an option... no budget (wife kills me ....yes she told me..)

I sure would like a diy project that is on a lower budget...with a Microchip that knows what to do..

Best regards

Ap

In an attachment to my previous post are the source of the Jim Koehler's project from Philips contest. I'm not made this project. We produce our proton magnetometers. Also, we have launched the production of our overhauser sensors. The price range for our devices - from $ 2200 to $ 8000.

We had high hopes for Jim Koehler's algorithm, because we wanted to get the accuracy of 0,1 nT with proton sensors. But unfortunately, it does not get the required accuracy with the standard signal-to-noise ratio of 30-40 dB . So we had to develop their own algorithms for computing the precession frequency.

BR
Konstantin