I think you missed the point of the code I gave. The key line was the "x = x.astype(int8)" that converted the information to 8 bit amplitude resolution (in the C programming language, it was converted to the signed character type); this simulates the data coming from one of the instruments. The subtraction to get the derivative results in severe roundoff error. In fact, if you print out the difference array, the only numbers in the array are 0 and ±1. Scaling the difference array by any factor is useless because the information is gone.

You won't see the problem if you do it with the typical 8-byte IEEE floating point arithmetic like you did in the spreadsheet. But you can probably simulate the problem in your spreadsheet by multiplying by 128 and then taking the integer value -- then do the differencing.

But, as I said, for edge and glitch detection, that's probably not important, as it's the high frequency stuff of interest. But it is something I'll have to warn the unwary user about.

I'm struggling to get you DSP-knowledgeable folks to understand the programming constraints I'm up against here. That array x in the code I gave above is a numpy array containing 1 to 2 million floating point values. In reality, it's a chunk of memory allocated in a C library and given array semantics by numpy. But I don't have the luxury of iterating on the points in that array in C -- I can do it, but it's in a scripting language (python), which is hundreds of times slower than C. If I use a chunk of python code like

for i in xrange(len(x)):
    process(x)

to mess with individual points, I've slowed the application down a lot -- and it will be unacceptable to users. Thus, I have to use the built-in functions of numpy to do the point-by-point processing and I'm constrained by what numpy offers.

When I get back to working on this problem, I'll take a look at the suggestion by Tesla23, as it meshes with your suggestions and might be a way of getting at the interesting stuff. In fact, numpy does have primitives to do comparisons and choosing things from an array in one fell swoop, so this might be a fruitful approach.

The other significant constraint I have is that I can't throw unlimited effort at the problem, as I'm doing the work as a consultant and the customer has to be satisfied that the work is necessary. Right now, we have other issues that need to be attacked first. I recently changed to a sampled display (the waveform shown on the screen is a sample of the waveform in memory); the older version plotted every point. Thus, the old code showed every little wiggle in the waveform, but the new code doesn't. I'm anticipating that the engineers at the company I'm doing the work for will find that not displaying the glitches (i.e., interesting stuff) like before is problematic for their customers and will want some change in the code to do so. Hence this question, as I don't really have a good idea on how to attack the problem. I thought some of the AAC readers would have some good ideas, so that's why I asked.

 

I realized when I created the spread sheet that the results were precision, but the difference multiplied has to provide a derivative even if its noisy. I thought you had just looked at the difference data and concluded that all the information was hidden in the noise. This cannot be the case as the noisy sine wave still follows it recognizable path. The quantized steps being 8 bits with noise and little precision, still has to follow the path, and I believe it will still produce the desired results. The bigger question I would have is how small of a glitch should you be able to detect. You are correct in that the degree of precision will limit the useful range, but like most things there must be trade offs if the customer only supplies noisy, low precision data. I just didn't want to spend much time to program up a simulation with real random noise and lower precision, but still believe by the 2nd and 3rd order derivatives are generated, the noise will tend to equalize out.

 

Again, I thank all the folks who have provided some comments and help!

One key takeaway for me so far from this experience is that the amplitude-range-limited data provide a barrier for doing some processing. This 8-bit-limited data does have a problem in that it's essentially useless to calculate the derivative, whether or not there's noise present on the data. For example, in the python code I gave, if the line converting things to 8 bits is commented out, the derivative is calculated just fine (of course, it still needs to be scaled to see on the same graph). But if the 8-bit conversion is used, the derivative information is lost because of roundoff error. I hadn't seen this before, primarily because I hadn't worked with such amplitude limited information before.

My interest in this problem is still there and I might put some work into it -- even if the customer isn't willing to pay for it...

 

While the input precision is limited, are your calculations limited to 8 bit integer as well? I don't program in python, and I don't know what your processing time constraints are, but the input is certainly a hindrance.

I did move the data into 8 bit defined realm and it does present a problem. But if the processing of that input data is purely limited to 8 bits as well, then you are severely restricted. I still think even the difference data could identify a glitch if the definition is crafted such that a glitch is a difference of 3 or more. (In the 8 bit data for a sine wave the highest difference was 1, adding headroom for noise up to another 1, then 3 would cross the boundary certainly - but you could easily miss a "glitch" your customer is trying to identify.)

I would say it would have to be a noticeable departure on the graph to get caught. This is where the input data being pure digitized causes a problem. If you can change the sensor to a different kind of measurement, like delta-sigma, or something that identifies smaller changes with more precision than bigger changes, it would help.

 

No, all calculation can be done in 8-byte IEEE floats as is common on PC/UNIX platforms. It's just that the input is coming from an oscilloscope that has 8 bits of vertical resolution. There's no way to change that.

I'll play around with some of the ideas you folks have given here when I get another scope (the one I had quit communicating with my computer).

 

You are probably already aware of this also, but I'll throw it out just to mull over. One of the first things I had to reconcile with the digital scopes, is that the digital sampling was very fast at a point, by that I mean the aperture of acquisition was tiny (maybe 2nS) at slower speeds. This caused one time for the entire signal to become hidden as the sampling rate of the scope just happened to line up (against) with the pulse stream. This resulted in a flat scope image of a known pulse stream. I found it quickly and was able to work to get my trace, but with slower speed traces, things get lost. I say this because not only are the 8 bits going to cause a problem of glitch detection, but conceivably the acquisition window, frequency response and sampling rate can work against detection as well.
I had read an article, either in EDN or Embedded Systems where they were describing oversampling to gain higher precision from the A/D. Again you are stuck here as it appears that you just have the stream you are provided to work with. The oversampling work with an interpolation algorithm to view any change in the oversampled signal. If during that extra sample you could identify a slope change you could achieve an extra bit of resolution. There were constraints, but depending on the number of oversamples, it indicated that you could achieve an additional 2 bits. I don't know what relationship you have with the manufacturer or the customer (if they are one in the same) that might allow you to get an oversampled stream perhaps.

Just food for thought.

 

I'm not sure what "oversampling" in this context means. I skimmed a wikipedia article on it and it sounds like it's used in the case where you have a band-limited channel and it just means to increase the sampling frequency beyond the Nyquist frequency to get margin, reduce aliasing, etc. But a scope's input isn't band-limited (except for the input hardware), so I'm not making the connection. Regardless, my customer is buying their hardware from someone else and I'd doubt there would be any such features (not to mention very late in the development cycle).

Another method just occurred to me that I'm going to have to try. In numpy, taking a sample of a waveform is fast and easy. For example, to take every 10th point in a sample, the syntax is xsampled = x[::10]; that last number controls the stride of the slice through the array. In reality, it's the step size in a tight C for loop.

So, if I were to take, say, 10 samples of a waveform, each of them offset a bit from the others, then somehow OR or average them together (doing such stuff with a waveform of a few thousand points is sub-millisecond in time, so no big deal), this might stand a better chance of showing something interesting. In reality, it's similar to averaging nearby points, but it's something that can be done on a whole array at once and no 'for' loop over all the points in the waveform required.

 

Oversampling is only useful in a bandwidth limited sense. Nyquist rate would be the normal sampling rate (so the max would be at the scope bandwidth limitation), and the oversampling would be 2x or 4x to get the additional resolution. But if it is either too late in the development, or completely inaccessible, then it does not matter.

And yes you can take the jumps across the array as you indicated, but this is where I have gotten into trouble, as the things of interest [just] happen to fall in between the samples. So if you use your 10 points of jump and all the pulses or sin x/x responses fall in between, then the average is lost. Happened very frequently in image processing with gray lines or color gradiations.

I think you found one of the good problems to work on, there is certainly no easy solution.

 

Yes, I too worry about missing things because only 1/1000th to 1/100th of the waveform will be in the sample.

Again, I thank you for your help and insight. I'm hoping to get another scope in the coming weeks (my customer is a few states away, so there are delays) and play with this a bit more. When I explain the problem to them, I suspect they won't be terribly interested in investing effort into coming up with a better algorithm for displaying these big waveforms on a small screen. That's OK, but I've gotten interested in the problem. Thus, I'll probably work on it just out of curiosity.