Time domain interpolation using FFT with zero padding on the end I've got a situation where I'd like to use an FFT to do interpolation in time on some complex data (I need to go to the frequency domain anyways to window my data).

The notional way of doing this is to take an FFT, which, in the case of something like FFTW will put the frequency components "in-order" which can be viewed as having positive frequencies first, followed by negative.

$\mathcal{F}([\begin{array}{ccccccc} a & b & c & d & e & f & g & h \end{array}]) = \begin{array}{cccccccc} [A & B & C & D & E & F & G & H] \end{array}$

Where ABCD are positive frequency components, and EFGH are negative frequency components.

You can split the frequency response in two and zero-pad between D and E, but you have to be careful if your transform is an odd length, as well as possibly splitting the nyquist bin, and it's more cumbersome from a buffer-juggling perspective too.  It'd be better if I could pre-multiply by Fs/2, which will rotate the spectrum:

$\mathcal{F}(exp(j2\pi fs/2t) * [\begin{array}{ccccccc} a & b & c & d & e & f & g & h \end{array}]) = \begin{array}{cccccccc} [E & F & G & H & A & B & C & D] \end{array}$

Then it's just a matter of zero padding on the end, and getting back to the time domain.  This might be as simple as just transforming into a sufficiently large zeroed buffer.

But, it's the getting back to the time domain that's got me wrapped around the axle.  You can certainly inverse-FFT, and then multiply by -Fs/4 to shift things back to the way they should be, but I'd prefer to avoid the post-multiplication if I can (performance reasons).

My gut says that the shifted spectrum is now in the proper order (again in FFTW land) to just do a forward-FFT on it, but that might time-reverse the data, so perhaps a conjugate before transforming would be in order, but I'm not sure.  Can anyone help me puzzle this out?

It's important to remember that, for an even N, you have both F(0) and F(N/2) points, neither of which has a counterpart.  So, for instance, with N = 8, you have:

F(0),  F(1),  F(2),  F(3),  F(4),  F(-3),  F(-2),  F(-1)

(some people consider the N/2 point as a negative frequency ( ie: F(-4) ), but we'll consider it to be positive)

So to interpolate an 8 point FFT by a factor of 2, you'd do a 16 point inverse FFT using:

F(0),  F(1),  F(2),  F(3),  F(4)/2,  0,  0,  0,  0,  0,  0,  0,  F(4)/2,  F(-3),  F(-2),  F(-1)

So you split the N/2 point of your forward FFT, and add the appropriate number of zeroes in the center, and do the inverse.

Now, for an odd number N, you DON'T have an N/2 point (but you still have a unique 0 point), and the splitting business doesn't happen (but you still have to add the appropriate zeroes in the center).

I suggest you experiment with the FFT/IFFT method by using: a real input time domain waveform; a complex input time domain waveform; even and odd N; an integer number of cycles in the input; a non-integer number of cycles in the input.  Then compare/contrast your results from the FFT/IFFT method with to those generated by a strictly time domain based interpolation (eg: use A*cos(omega*t) for the FFT/IFFT method and compare to A*cos(omega*t/2) generated in the time domain).  Experience can be enlightening.

As for your multiplying method – I don't think it will work for an even N.  And the FFT/IFFT timing, unless you're using very small FFT's, will probably dominate the total run time (of, course, you should always profile your programs to make sure).


Is the system $y[n]=x[n]+2=T\{x[n]\}$ an LTI-System? Is the system $y[n]=x[n]+2=T\{x[n]\}$ an LTI-System?

Linearity:
$ \alpha \cdot T\{x_1[n] \} + \beta \cdot T\{x_2[n] \} = T\{\alpha \cdot x_1[n]\ + \beta \cdot x_2[n] \} \\\alpha \cdot (x_1[n]+2) + \beta  \cdot(x_2[n]+2)= \alpha \cdot (x_1[n]+2) + \beta  \cdot(x_2[n]+2)$

Time-Invariance:
$y_1[n] = y[n-n_0] = x[n-n_0]+2 \\ x_2[n]= x_1[n-n_0]\\ y_2 [n]=T\{x_2[n ]\}= x[n-n_0]+2 \\y_1[n]=y_2[n]$

So I would say it is an LTI-System, is that right?

A system that adds a constant is rarely linear. If you have two inputs $x_1$ and $x_2$, each one gets the constant, so their combination get the constant two times. While the input $x_1+x_2$ gets the constant only once.


Is a DC signal 0 frequency or extremely high? If I pass a signal with DC component through a Low Pass Filter will the output be without the DC component?

In most mathematical treatments of DSP, a "DC" signal acts like a frequency component with a frequency and phase of zero (0 Hz, cos(0)).


Cross-Correlation of two signals each with independent complex Gaussian noise The nearest question I can find is here, but it does not have an answer. Say I want to correlate the following two noisy channels: $x_1 = s_1 + n_1$ and $x_2 = s_2 + n_2$, where $s_1$ and $s_2$ are arbitrary signals and $n_1$ and $n_2$ are independent complex Gaussian noise (identically distributed for starters). Further, let either $M \lt N$ or $M \ll N$, with $x_1 \in \mathbb{C}^M$ and $x_2 \in \mathbb{C}^N$. The cross-correlation of the two channels yields

$$
x_1 \star x_2 = s_1 \star s_2 + s_1 \star n_2 + n_1 \star s_2 + n_1 \star n_2.
$$

Perhaps I am searching with the wrong keywords, but I am having trouble finding information about $n_1 \star n_2$. What kind of distribution does it have? How much power does it contribute to the output?

I don't understand why your complex vectors have different dimensions but there is an answer for a scaler case,

from:

Leo A. Goodman. On the exact variance of products. Journal of the American Statistical Association, 55(292):708–713, 1960.

Some definitions,

$$\begin{array}{c}
m_{x1}= E\{ x_1 \}\\
var(x_1) = E\{ |x_1|^2\} - |m_{x1}|^2 \\
m_{x2}= E\{ x_2 \} \\
var(x_2) = E\{ |x_2|^2\} - |m_{x2}|^2
\end{array} $$
for $x_1$ and $x_2$ independent,

$$ var(x_1 x_2^{*}) = var(x_1)|m_{x2}|^2 +var(x_2^*)|m_{x1}|^2 + var(x_2^*)var(x_1) $$

You don't specify the distribution of your signals.
assuming deterministic signals

N. O’Donoughue and J. M. F. Moura. On the product of independent complex gaussians. IEEE Transactions on Signal Processing, 60(3):1050– 1063, March 2012.


Crop all written letters from single image Above image has handwritten hindi / deavanagari letters all on a single

image.
Task is to Prepare a dataset from the given assignment.jpg image using image

processing tools
Idon't have clarity whether i can use photoshop editing tools like adobe photoshop etc.
If image processing tools mean opencv , is it really possible to create like a emnist like dataset from single image.
I am not able to come up with any way to seperate the letters from eachother

Presumably this is an academic assignment of some sort and you are expected to use a specific tool, OpenCV I am guessing. You should review the documentation for that tool and try to find a feature which matches the requirements. I suspect you would need to use the Image Segmentation feature: https://docs.opencv.org/3.4.3/d3/db4/tutorial_py_watershed.html


How to understand the energy centralizes at low frequency after making DCT on an image?What does energy mean? I've already know the low frequency represents a smooth area,and most part of an image doesn't have abrupt change which can explain natural image has more low frequency. But I just can't understand what does the term 'energy' mean.I would really appreciate if some can help me!

It's just a fancy way of saying "amplitude" or "signal strength".  The DCT splits the original signal into a bunch of outputs, each of which gives the amplitude for a particular frequency.  Add them all back together again, through an inverse DCT, and you get the original signal back.


Large Laplacian Spatial Filter on EEG? I want to know that why large laplacian spatial filtering is done on EEG signals? I tried a lot but could not find material on LLSF. And what is meant by spatial filter ( in context of Digital Signal Processing)? Why is it called "Large" and "laplacian"?

Laplacian

It is a 2nd derivative of the recorded voltage for each electrode. Suppose the recordings at position $(i,j)$ is $V(i,j)$, then Laplacian operation  at $(i,j)$ is

$(V(i,j)-2V(i-1,j)+V(i-2,j))^2 + (V(i,j)-2V(i,j-1)+V(i,j-2))^2$

And it serves as a high-pass filter that enhances localized activity while suppresses the diffusion activity. The implementation of the Laplacian in EEG filtering on the voltage at each electrode is to subtract the weighted voltages from the surrounding electrodes from the voltage recording at current electrode, where the weight is electrode distance dependent.

Large

The difference between large and small is determined by how you select the electrodes. For the small Laplacian, you select the nearest-neighbor electrodes, while for the large Laplacian, they are the next-nearest-neighbor electrodes. See the figure below. If the signal you are interested in is less highly localized and/or varies in exact location over
time, the large Laplacian is expected.

Read the paper Spatial filter selection for EEG-based communication by Dennis J. McFarland et all for more details.


To remove doppler for equalizer doppler does it need decision directed like carrier recovery? Normally equalizer is used at the start of the receive chain .So do we need decision directed to estimate the doppler here and also later for carrier recovery for carrier phase recovery?

Doppler offsets are typically small fractions on the symbol rate and as such do not significantly effect the task of the equalizer which can perform under such carrier offsets. The equalized waveform is then passed through timing recovery and carrier recovery synchronization loops.


Advantage of using FIR Polyphase Decimator over just computing FIR output every D samples Let's say you want to implement an FIR filter $h[k]$ with $L$ taps and downsample the output:

$v[k] = x[k] * h[k]$

$y[m] = v[kD]$

The "naive" way would be to compute all the samples of $v[k]$, then compute $y[m]$.  This should have complexity of $O(LD)$ to process a block of $D$ samples.

The "less naive" way would be to just compute every $D^{th}$ sample of $y[m]$ directly:

$y[m] = v[mD] = \sum_{l = 0}^{L-1} x[mD - l] h[l]$

This should have a complexity of $O(LD/D) = O(L)$ to process a block of $D$ samples.

The "smart" way, according to everything I read (e.g. Proakis & Manolakis Digital Signal Processing: Principles, Algorithms, and Applications, 4th ed. p.771, also here), is to use a polyphase filter/decimator.

This seems to have the same computational burden: every $D^{th}$ timestep, you are computing the output of $D$ filters, each of which has $L/D$ taps, for a total complexity of roughly $O(DL/D) = O(L)$.

Is my analysis correct? If so, what is the advantage of using a polyphase decimator if the "skip $D$" method has the same complexity and is much easier to implement?

One advantage of the polyphase filter is when hardware resources are limited. The polyphase approach replaces a $L$ tap filter with $N$ filter sets of $L/N$ taps. Once you've done that you can just use a single $L/N$ filter set and swap coefficients.


how can i center an audio signal? I'm trying to reconstruct an audio signal, and below is what i got after rescaling.

How can i obtain centered audio signal such as below?

SOLVED:

ruoho ruotsi's answer addresses my problem :)

i found someone asked the similar question before

How to remove or filter the drift problem in measured Strain signal?

"Center" isn't the word that I'd use to describe what you want to do. Your signal does not have a zero-mean, so like @Dole suggested you should use a DC-Blocker (DC-Offset removing) filter to get the job done:

Here are your references:

https://ccrma.stanford.edu/~jos/filters/DC_Blocker.html
https://ccrma.stanford.edu/~jos/filters/DC_Blocker_Frequency_Response.html
https://ccrma.stanford.edu/~jos/filters/DC_Blocker_Software_Implementations.html

the code (from the above links) can be as simple as this, where x is the current sample and y is the filtered sample, while xm1 and ym1 are the x,y from the previous step (i.e. x[i-1], y[i-1]):

  y = x - xm1 + 0.995 * ym1;
  xm1 = x;
  ym1 = y;


Quantify the difference between two audio spectrograms I have two speech signals of the same length. They contain the exact same speech content but are recorded from two different microphones. I want to compare the differences between them.

1) Assuming that the two signals are perfectly aligned, I can compute the spectrograms for both of them and calculate a distance/coherence metric (L2 distance or magnitude-squared coherence) between the two spectrograms. Are there any other ways of quantifying spectrogram differences?

2) If the signals are not perfectly aligned - even a misalignment of a few milliseconds - can mess up the L2 distance between spectrograms. What could be some strategies in this scenario?

Thanks for your inputs.

You could do a cross correlation between signals to even out the time offset.
But why go for spectrogram in the first place? Do you have acoustic conditions varying over time? If not, go for the spectrum, should do the same job but without the time alignement problem.


Is there any algorithm to decompose a wave into a set of given frequencies? I heard about Fast Fourier Transform can decompose any waves.

But it seems like it was hitting for "all frequencies", namely, n HAVE TO be from 0 to N-1!

What if I want to decompose the wave into the frequencies I want!?

Such as, for example, 5Hz, 14Hz, 37Hz, 42Hz, 59Hz, etc.!?

Is there any algorithm that could do the job!?

PS. My understanding of FFT is from Google, and the degree of understanding is VERY limited! So PLEASE forgive me for my lack of knowledge and any stupid comments!

You could use the Goertzel algorithm to compute the discrete Fourier transformation for just the frequencies you want to detect.

The Goertzel algorithm is used to detect single bins from the DFT.  You implement a bunch of them, each set to one of your desired frequencies, and calculate them for each audio sample.

Anything that doesn't fit your selected bins is ignored.  If you implemented it for the frequencies given in your example (5Hz, 14Hz, 37Hz, 42Hz, and 59Hz) but your audio was a sine wave at 100Hz, then your software could not detect that 100Hz signal.

The Goertzel algirithm can be "tuned" to make the bins narrower.  You would probably need to do that to separate the frequencies given in your example.  They are very close together, so you would need to use very large $N$ to calculate the coefficients for each bin.  This has the side effect of making the algorithm respond slowly - quickly changing frequencies would not be detected properly.


what is happening when i run this code line? I have the following MATLAB code:

x=[1 2 3; 4 5 6; 7 8 9;];
u=(0:255)/(255-1).^1.05;
x2=u(x+1)

when I run this code I get the following as a result in x2:

x2 =

0.0030    0.0060    0.0090
0.0119    0.0149    0.0179
0.0209    0.0239    0.0269

what is happening when I call this line x2=u(x+1)? and how to write it mathematically?

Thank you very much in advance

You are putting the values of

u(2),...,u(10)

in a $3\times 3$ matrix with the shape of the index matrix $x$. This could be written simply as $$x2_{(i,j)}= u (x_{(i,j)}+1),\;(i,j)\in[1,2,3]\times[1,2,3]\,. $$

For further matrix shape conversions, you can have a look at vectorization or matricization, and you can find notations in chapter 3.4 Matricization of a tensor of Multilinear operators for higher-order decompositions, T. Kolda, 2006.


What is the set of most predictable or in other words most self-correlated signals? Which signals or set of signals is the most predictable or in other words most self-correlated (autocorrelation)?

I think one of them is $x(t)=c$, for all $t$, where $c$ is a scalar constant. The other one can be $x(t)=c$, for $t=t_0$ and $x(t)=0$, otherwise (a discrete signal with only one value). $t_0$ is a constant.

Here the book says: "the most predictable signal" is the sum of discrete cosine functions with random phases. Why? My two examples above seem to me more predictable!

Sums of discrete cosines are good candidates because they're bound to repeat after a finite number of samples. In fact,  anything that repeats after some finite number of samples is the most predictable thing in the world because you can predict the entire thing from just a handful of samples. In fact, the converse of this statement is true, i.e. anything that repeats after a finite number of samples can be written down as a sum of discrete cosines with different phases. In fact, so can your own candidate:
$$c = \sum_{k=0}^{\infty}c_k\cos(\omega_k + \phi_k)$$
where $c_k=c\delta(k)$, $\omega_0=0$, $\phi_0=0$. In other words a constant signal is just a cosine of zero frequency. So there's no tautology here.

You may say of course, that your example is still better, but if we're trying to establish classes of predictability, then both your example and the generic form fall into the same class. I don't see a contradiction here.


FIR Filter - Low Pass Filter coefficients to High Pass filter coefficients I need to convert this 12 Tap Low Pass Filter to a High Pass Filter.

The coefficients of the Low Pass Filter are as follows:

h[k] = [0.022 -0.083 -0.083 0.083 0.311 0.417 0.311 0.083 -0.083 -0.083 0.022 0.083]

My understanding is that if its an even-tap - We Multiply the coeffs by -1.

However the answer is:

hT[k] = [-0.022 -0.083 0.083 -0.083 -0.083 0.311 -0.417 0.311 -0.083 0.083 -0.083 -0.022 0.083]

I do not know how this came about. I looked for answers online alot but did not get the answer. Any help is appreciated

HINT:

It looks like they shifted the frequency response by half the sampling frequency, i.e.,

$$H_{HP}\left(e^{j\omega}\right)=H_{LP}\big(e^{j(\omega-\pi)}\big)\tag{1}$$

Frequency shifting corresponds to modulation (multiplication) in the time domain. Now you just have to figure out the modulation sequence that achieves the correct frequency shift.


Red Blood Cell Color Image Segmentation I want to segment this Red Blood Cell Image on HSI color space.I want to try the Otsu Method but from what I read that method work only on Grayscale image.Do you have any suggestion for suitable image segmentation on those red blood cell image ?
[EDIT]
This is what i want to achieve(Segmented Blood Cell)

Using Matlab:

You might try converting your image to greyscale first with

greycells = rgb2gray(rgbcells);

where rgbcells contains the color image of your cells (see here).

Then use Otsu's threshold on the greyscale image like so

threshcells = graythresh(greycells);

like so (see here).

EDIT
After OP explained what image should be the resulting image.

Using: FIJI (formerly ImageJ)

I am still not 100% sure as to which images you finally want to obtain. It seems like you only want to see the cells with the many purple dots (sorry, I am not a biologist myself) inside them and completely remove the others. For this attempt at a solution I use FIJI, but the principle can be easily applied in Python / Matlab, you name it...

Your images are somewhat different, because in the first image the other cells pretty much have the same color as the one in the middle (that you want to segment probably). In the second image the other cells are much lighter.

Let us look at the second image first:

Converting to HSI

Use the Color Transformer 2 Plugin to convert the images to HSI using the options "From Colorspace: sRGB" and "To Colorspace: HSI". Looking at the channels we can see that S and I carry pretty much the same information. So let's stick with the I channel. It looks like this: .

Selecting a Threshold

Now we use the useful "Auto Threshold" (Image->Adjust->Auto Threshold) feature of Fiji on this channel. Clip the image from the HSI stack and convert it to 16 bit first. Use the option "Try All", also check "Ignore White" and uncheck "white objects on black background". Now we can see a whole panel of different thresholding algorithms applied to this image. It looks like this:

You can see that many thresholding methods do what you want. Among them Otsus method (as you suggested) and MaxEnt. Now that we know MaxEnt or Otsus thresholding works for this image lets do the same thing again (meaning AutoThreshold), but instead of "Try All" we select "MaxEnt" as our threshold of choice (you'll see below why MaxEnt). Then we invert the threshold and overlay it on the image and we end up with:

Now for the first image

Here we can do exactly the same. Repeating the steps, we end up with an I channel like this:

The other cells are very well visible here. This is also true for the other channels. I have also plotted the histogram. The left hump corresponds to the wihte background, the right hump to pretty much all cells. This will be important in the next step.

Now try all thresholds again like above and we will see this:

You see that Otsus threshold separates the cells from the background, but does not give you the cell in the middle that you want. Look at the Histogram: due to the way Otsus method is designed (see here) it will seperate these two clusters. Which is not what you want (as I understand) in this case. But MaxEnt(Maximum Entropy) does. So use the MaxEnt threshold as above.

You will end up with this:

SUMMARY
This is my processing chain:

Do HSI transform
Apply MaxEnt threshold to the I channel
Use the binary mask on the RGB image again to clip out the relevant parts

OUTLOOK
If you are not satisfied with the results, then play around with different thresholds. Ultimately changing the colorspance and thresholding are powerful but very low tech ways of image segmentation. You might be interested in more sophisticated ways. For a software suite catered especially to cell biologists check out Cell Profiler. Former colleagoues have used it and were very happy with it.

LITERATURE
Maximum Entropy Threshold: Wong et al: "A gray-level threshold selection method based on maximum entropy principle",  https://ieeexplore.ieee.org/document/35351


Finding minimum-phase/allpass of transfer function I am trying to find the minimum-phase system of the transfer function $H(z) = \frac{2 + 3.125z^{-2}}{1-0.9z^{-1}+0.81z^{-2}}$.

I know I need to find and remove the allpass (basically reflect the two zeroes on the imaginary axis inside the unit circle?), and I've gotten so far as this: $H(z) = 2\frac{1 + \frac{3.125}{2}z^{-2}}{1-0.9z^{-1}+0.81z^{-2}} = 2\frac{1 + \frac{3.125}{2}z^{-2}}{1-0.9z^{-1}+0.81z^{-2}} \frac{\frac{3.125}{2} + z^{-2}}{1 + \frac{3.125}{2}z^{-2}}=,$

and then I'm not really sure where to go next. Is my thinking correct thus far? What step should I do next?

Thank you.

Edit:
This is the kind of thing I'm trying to do, as per the example in my book:

The allpass part must have its zeros outside the unit circle, so it gets the zeros of your system, with the corresponding poles inside the unit circle. You get the correct pole locations by simply swapping the coefficients of the numerator polynomial:

$$H_a(z)=\frac{2+3.125z^{-2}}{3.125+2z^{-2}}\tag{1}$$

The minimum phase part is obtained by taking the system's poles and cancelling the poles of the allpass part:

$$H_m(z)=\frac{3.125+2z^{-2}}{1-0.9z^{-1}+0.81z^{-2}}\tag{2}$$

The transfer function of the original system is given by

$$H(z)=H_a(z)\cdot H_m(z)\tag{3}$$


Z-Transform of $x(n) = 3^n$ First of all, thank you all for your answers. I know the z transform for

$$ x(n)=3^n \space ; \space  n\geqslant  3 $$

or rather

$$ x(n)= 3^n u(n-3) $$$$\begin{align}X(z)&=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\\&=\sum_{n=-\infty}^\infty3^nu(n-3)z^{-n}\\& \text{change of variables: } n'=n-3 \\& =\sum_{n=-\infty}^\infty 3^{(n'+3)} u(n') z^{(-n'-3)}\\& \text{recall:  } n'-> n \\&=3^3z^{-3}\sum_{n=0}^{\infty}(3 z^{-1})^{n}\\& X(z) =3^3 z^{-3}\frac{1}{1-3{z^{-1}}},\quad |z|>3  \tag{1}\end{align} $$\begin{align} \\&  \text{instead, for the sequence: } x(n)=3^{n} \ ; n\leqslant -3  \\&=\sum_{n=-\infty}^\infty3^nu(-n-3)z^{-n}\\& \text{change of variables: } n'=-n-3 \\& =\sum_{n=-\infty}^\infty 3^{(-n'-3)} u(n') z^{-(-n'-3)}\\& \text{recall:  } n'-> n \\&=3^{-3} z^{3}\sum_{n=0}^{\infty}(3^{-1} z)^{n}\\&  X(Z)= 3^{-3}z^{3}\frac{1}{1- \frac{z}{3}} \quad |z|<3\tag{2}\end{align}
Is it correct?
And for @Matt L. I don't understand your mathematical step initial when you write:
$$\begin {align} \sum_{n=-\infty}^3 3^nz^{-n}  =\sum_{n=-\infty}^0 3^{(n+3)} u(n') z^{-(n+3)} \end {align} $$

Can you write me the generic formula?

The signal $x[n]=3^n$ for $n\le -3$ (and zero otherwise, I assume) can be written as

$$x[n]=3^nu[-n-3]\tag{1}$$

To compute its $\mathcal{Z}$-transform you simply need to use the formula:

$$\begin{align}X(z)&=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\\&=\sum_{n=-\infty}^{-3}3^nz^{-n}\\&=\sum_{n=-\infty}^03^{(n-3)}z^{-(n-3)}\\&=3^{-3}z^{3}\sum_{n=0}^{\infty}\left(\frac{z}{3}\right)^n\\&=3^{-3}z^{3}\frac{1}{1-\frac13 z},\quad |z|<3\tag{2}\end{align}$$

It's important to specify the region of convergence ($|z|<3$), otherwise the $\mathcal{Z}$-transform does not uniquely identify the corresponding sequence $x[n]$.


Fastest Algorithm for Distance Transform I am looking for the fastest available algorithm for distance transform.

According to this site http://homepages.inf.ed.ac.uk/rbf/HIPR2/distance.htm, it describes:

  The distance transform can be calculated much more efficiently using
  clever algorithms in only two passes (e.g. Rosenfeld and Pfaltz
  1968).

Searching around, I found: "Rosenfeld, A and Pfaltz, J L. 1968. Distance Functions on Digital Pictures. Pattern Recognition, 1, 33-61."

But I believe we should have a better and faster algorithm than the one in 1968 already? In fact, I could not find the source from 1968, so any help is highly appreciated.

Pedro F. Felzenszwalb and Daniel P. Huttenlocher have published their implementation for the distance transform. You cannot use it for volumetric images, but maybe you can extend it to support 3d data. I have only used it as a black box.


Guitar String frequency I'm making a guitar tuner for iOS.Due to the fact I'm Beginner I'm struggling a bit to gather all the resources and information about it. I know the theory like ( correct me if I'm wrong ) :-

First I need to get the input from microphone.
Then need to apply apply FFT algorithm to get the frequency
Then compare the frequency with the fundamental frequency of notes.

according to my research the fundamental frequency of the guitar are measured as:-

String  Frequency
1 (E)   329.63 Hz
2 (B)   246.94 Hz
3 (G)   196.00 Hz
4 (D)   146.83 Hz
5 (A)   110.00 Hz
6 (E)   82.41  Hz

and from an test app example I'm getting the frequency:-

String  Frequency
1 (E)   333.02  Hz
2 (B)   247.60  Hz
3 (G)   398-193 Hz    (398 when start and 193 when end)
4 (D)   290-150 Hz    (290 when start and 150 when end)
5 (A)   333-215 Hz    (333 when start and 215 when end)
6 (E)   247-161 Hz    (247 when start and 161 when end)

  One thing to notice the example was showing the Max Frequency.

So what I'm asking here is there anybody who have implemented this before
and give me the direction and some detail information about the topic so In future nobody have to research more and more to gather all the resources.I need help with the following topics:-

What is the best possible way to get the accurate frequency of guitar string.
Which type of frequency need to tune a guitar like maxFrequency, is there a role of altitude, magnitude (i have less information about this topic).
What to do after getting the right frequency.

Any help would be truly appreciated and my apologies in advance if i'm asking something stupid .

The fft provides the spectrum of the audio signal, which tells you the frequencies your signal is made of. In the case of harmonic signals like guitar, violin and other musical instruments, the spectrum exhibits peaks that are (almost) equally separated. For example, 110, 220, 330, 440... Hz. These peaks are called harmonics and they have different amplitudes, which contribute to the perception of timbre.

The perceived note (aka the pitch) depends on the frequency of the first peak, that is the fundamental frequency. It could happen that the fundamental frequency, which is what you want to retrieve with your tuner, is not the strongest one in terms of amplitude. For this reason, you need more complex techniques than just retrieving the maximum of the spectrum. I suggest you to google pitch detection, which is the problem you are working on. There is a lot of material and possible solutions that you can compare and implement.


Changepoint Detection of a 2D signal Consider the 2D signal below

I would like to find the changepoints in the signal (using a hierarchical approach). Based on our visual observation, we can see there are 2 main changepoints where the signal 'starts' and 'stops' and we could argue there is potentially another one where the red starts to come in. We could continue to look in more and more granularity if we so wished.

My question is, how could we get a computer to detect these automatically. My intuition from the way I did it was to do as follows:

Select the first point as our starting point.
Move along the time axis with another point, finding the Fourier transform of the data between the first point and this point.
As long as the transformed signal is not become too much more complex (in the information sense), we know that we have not hit any big change.
Once the signal starts to become more complex (above some threshold level), we select this as our new starting point and repeat.

After this process is finished we now have selected an array of points where our signal changes drastically w.r.t. its spectral make up.

My question then is:

Does this make sense?
What metric should I use for the complexity of the Fourier transform?
Has this already been tackled successfully and if so any references?

One approach would be to use the spectral flux of a sequence of STFT frames $X(n,k)$, namely

$$
\mathit{SF}(n) = \sum_{k = -N/2}^{N/2 - 1}{H(\left|X(n,k)\right| - \left|X(n-1,k)\right|)},
$$

where $n$ is the frame number, $N$ is the size (in samples) of each frame, and $H(x) = \frac{x + \left|x\right|}{2}$ is the half-wave rectifier function. Intuitively, it attempts to measure how much the spectrum changes from frame to frame.

The overall area of research is called onset detection. A good overview of techniques can be found in the following paper (which also describes spectral flux):

  S. Dixon. Onset detection revisited. In Proceedings of the 9th
  International Conference on Digital Audio Effects, volume 120, pages
  133-137, 2006.

An example implementation of spectral flux can be found in the NSGToolbox.


FFT radix-2 in audio - resolution issues - logarithmic scale - sample rate VS buffer size I would like to ask how to deal with audio DSP. Is FFT good at all to deal with audio? I am working on an app, that generates a single impulse, rout it to an EQ plugin, and the output of the plugin is run through the FFT to see what exact changes were made by the plugin.
And I have found some issue with FFT, that I'm not sure how to solve:
FFT gives you the same resolution for the whole frequency range. If I have a sampling rate of 48 kHz, and a buffer (and DTFT) size of 1024 samples, then the resolution is 48000/1024 = 46.875 Hz. And the same resolution is for lower, medium and highest freq, in one word for the whole range. But in audio, we often work on logarithmic scales. So we need better resolution on low freq, than high freq. How do you deal with that?

Of course I can store data from buffer size until I have 48000 samples, and then run it through the FFT, then I would have resolution 1 Hz. Great, but at first: then I can refresh at least only one time per second, at second: I still think I don't need such a high resolution for high frequencies, and at third it didn't solve my problem if we have other sample rate like 44100 which is not 2^L.

So, is FFT not a good choice for that purpose, please tell me in which situations it is good to use radix-2 FFT in audio industry? As I know for filtering there are also better choices like FIR filters. So what is about all that FFT? :)

If you overlap, you can refresh an N long FFT every sample, or every M samples, no need to wait an entire 1 second, or for N samples.

No need to find a sample rate that is 2^n, as most modern FFT libraries support any reasonable product of several small primes.

For analysis, no need to use the same size FFT for both low and high frequencies (or mid, or upper mid, or an entire octave, etc.)

Why FFT?  O(NLogN) provided a speedup that in ancient pre-history allowed DFTs to be run on chips with far far less than a few million transistors, and today still allows real-time performance on power sources that fit in a watch or hearing aid or implant.  Could even allow some Fourier analysis of data when computers were women who used paper and pencil.  The square matrix derivation of a DFT requires less chalk on the classroom chalkboard (today: less iPad scrolling) than more modern audio tricks.

But if you have a larger power source with modern hardware/chips/iPhones, you can use piles of wavelets, large filter banks, psycho-acoustic models, Goertzels, and several other algorithms that are far less efficient than an O(NLogN) sped-up DFT, but together can still easily run faster that real-time requires.


Visualising a Z-transformed Transfer Function? For designing any analog filter and various other outputs of filter we use  laplace transform,I can visualise a laplace transform like for ex.
s[X(s)] can be implemented as differentiator fetched with signal x(t)while implementing  differentiators we generally use
capacitors or by using OPAMP similarly,
[X(s)]/s can be implemented as integrators using OPAMP.

But I cannot visualize the z transform in designing of digital filters,
while the implementing the digital filters we often map the s plane to z plane i cannot understand the significance

There are various method of implementing but in end we compare the s plane with z-plane ,
and deduce a equation where the value of z in found in place of s,

Why do we need to map the S plane to z plane?

in continuous-time functions, like $x(t)$, the operation of the derivative makes sense, it is well defined.  the three components to a continuous-time LTI filter are adders (two or more signals being added), scalers (a signal is simply multiplied by a constant), and integrators (the $s^{-1}$ operators).  the first two components do not discriminate with regard to frequency, so with just those two components, you cannot make a "filter" that filters out some frequencies more so than others.  but the integrator does act differently on sinusoids of different frequencies.  higher frequencies come out of the integrator reduced in amplitude more than lower frequencies.

in discrete-time functions (or "sequences"), like $x[n]$, the derivative operator does not makes sense, it is not defined.  the three components to a discrete-time LTI filter are adders (two or more signals being added), scalers (a signal is simply multiplied by a constant), and delay elements (the $z^{-1}$ operators).  the first two components do not discriminate with regard to frequency, so with just those two components, you cannot make a "filter" that filters out some frequencies more so than others.  but the delay element does act differently on sinusoids of different frequencies.  higher frequencies come out of the delay element shifted more in phase than lower frequencies.

analog filters act on physical quantities in time.  digital filters act on numbers.  these numbers are samples of a continuous-time function ($x[n]=x(nT)$) and make up a sequence.  there is no way to do derivatives or integrals directly.  but we can delay any of these sequences by use of computer memory.  that's why we talk of $H(z)$ more so than $H(s)$ when designing and implementing digital filters.


how to find sample rate of sum of 2 cos signals I have an exercise to find a sample rate of that signal:

while the filter is:

while A << 1.
I need to set the sampling rate so that at the output of the filter I'll get s[n] and reduce v[n] according to A.
can someone give me a guidance?
thanks alot

Your signal has maximum frequency component of 10Hz which means a sampling rate of 20Hz according to Nyquist theorem should suffice. Then the cut-off frequency of 5Hz in the ideal filter should be enough to retain s(t) and eliminate v(t). However in the digital domain, typically a much higher sample rate is needed to process and display continuous time signals. I used a sampling rate of 100Hz in a simple python script below. Cutoff frequency of the selected FIR filter is still 5Hz to match s(t). It is customary to normalize sampling frequency to \pi in which case cut-off frequency becomes \pi/20 instead of \pi/4 as you showed for the theoretical case.
Below python script:
import numpy as np
from scipy.signal import firwin, lfilter
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt

# sampling frequency in Hz
Fs = 100

# center frequencies of a linearly combined signal
fc1 = 5
fc2 = 10

# Observe duration in seconds, forumulate a time-axis
T = 5
tvals = np.arange(0, T, 1/Fs)

fig=plt.figure(figsize=(10,10))
s_t = np.cos(2*np.pi*fc1*tvals)
v_t = np.cos(2*np.pi*fc2*tvals)

plt.subplot(3,1,1)
plt.plot(tvals,s_t)
plt.ylabel('low-pass signal')
plt.grid()
x_t = s_t + v_t

plt.subplot(3,1,2)

plt.plot(tvals, x_t)
plt.ylabel('combined signal')
plt.grid()

# Filter Length
N = 31
taps = firwin(N, fc1/Fs)

# Simple FIR filtering
filtered_signal = lfilter(taps, 1.0, x_t)
ll = len(filtered_signal)
plt.subplot(3,1,3)

# Plot final result
plt.plot(tvals[int((N -1)/2):ll], filtered_signal[int((N - 1)/2):ll])
plt.grid()
plt.xlabel('Time')
plt.ylabel('filtered signal')
plt.tight_layout()
plt.show()

will generate a figure like:
Observe the recovered signal after filtering is the same (minus the filter delay effect) as the s(t)


Is a notch filter for 50 Hz without delay possible? I am processing an EEG signal (2 kHz sampling frequency) from hardware and the built-in filter doesn't work as intended. FIR gives a delay, and IIR is unstable.

The signal of interest is from 4 Hz to 360 Hz, with the goal to clear the utility signal (AC 50 Hz) from the signal. The main constraint is to get as many untouched (no noise, perturbance, and leaks) samples as possible. The tradeoff here is that gamma brain waves starts from about 40 Hz and goes up through 50 Hz and further.

Are there any options available to get no delay or at least a very small one (a 357 tap FIR gives more than 130 dB attenuation, but the delay with 2 kHz is about 180 ms. I aim at less than a 30 ms delay with at least 60 dB)?

The start of the pass band is 49 Hz and the end is at 51 Hz (in fact, it looks like this: https://www.swissgrid.ch/en/home/operation/regulation/grid-stability.html), so it varies in time and is non-stationary.

You could use a 2nd order IIR notch filter as I describe in this post Transfer function of second order notch filter - That post demonstrates a 50 Hz IIR notch with 1 KHz sampling.

[Update: As @user47050 astutely points out in the comments, the IIR notch would also have minimal delay regardless of notch bandwidth, since the dominat delay in the IIR notch filter is specific to the signal components in the rejection bandwidth of the filter, not the signals that you want to pass through. The advantage of what follows is not necessarily then to eliminate a delay issue, but to provide tracking rejection of just the interference that is approximately 50Hz (but slowly varying in phase and amplitude) while minimizing impact to adjacent frequencies of interest.]

However to provide cancellation of the interference with zero delay, consider a feedback signal cancellation servo such as I depict in the diagram below. This would have no delay to your signal but will have a delayed response to changes in your interfering signal. Results will depend on the relative stability of the interference. If the interference is a relatively slowly changing sinusoidal tone, this should work very well.

This is an all digital nulling loop that can null down to the rms power level of your signal (at which point it can no longer estimate the cancellation coeeficient). The bold lines shown in green represent complex signals:The 50 Hz PLL is a NCO with sine/cosine outputs, the multiplier following it is a full complex multiplier (essentially a vector modulator providing complete control of amplitude or phase), the multiplier after the second BPF is essentially a quadrature down-converter providing a measure of the amplitude and phase of the residual.

The way this works can be described as follows: The input signal is assumed to have a strong 50 Hz signal, strong enough such that combined with the filtering of the first 50 Hz BPF the PLL can lock in phase to the tone (the tone needs to be 6 dB stronger than the noise at the PLL input).  This will provide a phase locked copy of the 50 Hz signal as a quadrature NCO output. The quadrature PLL is basically a very tight tracking filter, providing a clean replica for purposes of cancellation.

Next we have a cancellation loop provided with a post cancellation band pass filter to detect the residual and adjust a complex coefficient to adjust the gain and phase of the 50 Hz PLL output prior to subtraction. I show the approach of quadrature downconverting the 50 Hz BPF output to complex baseband (by multiplying with the complex PLL output), and accumulating the result which will find the null with the convergence gain coefficient k. This is a cancellation loop using the method of stochastic gradient descent. The accumulator accumulates residual error as measured on the I and Q output of the quadature downcoverter resulting in the correct setting to the "Vector Modulator" multiplier resulting in a real sine wave output that will be equal in magnitude and opposite in phase to the interfering signal at the input.

The diagram can be simplified in that the first 50 Hz PLL BPF can be eliminated and the PLL replaced with a constant 50 Hz sine/cosine generator, since the vector modulator provides complete control of rotating phase which serves to to adjust the frequency positive or negative as needed, as well as gain as needed to provide the cancellation required. The proposed use of the 50 Hz BPF + PLL is that it would be feasible of providing faster corrections to phase changes in the input signal, once the cancellation loop has converged (since the PLL will track the input and can be designed with wider loop bandwidth than the nulling loop)- however if the amplitude change in the interference dominates, then I don't see a strong advantage to having the PLL.

Without the PLL this would reduce to the following block diagram. Here everything is done at complex baseband, so a low pass filter instead of a bandpass filter is used (which I would likely implement as an exponential averaging filter).  k is the complex convergence gain and c is the determined complex correction coefficient that will scale the gain and phase of the generated 50 Hz cancellation signal (which as suggested earlier c would be a changing phase versus time to adjust for the small frequency offsets that would also exist).

The PLL implementation has a clear advantage over a 2nd order IIR notch filter in that it can track and null a tone interference over a relatively much wider bandwidth, within the tracking range of the PLL. Without further study of the implementation below, I am not confident that it would out-perform a 2nd order IIR notch filter, but it is interesting to consider given it is much simpler than the PLL approach.

As a fixed 50 Hz source with a 2KHz sampling rate, the tone generator is trivial to implement, but if you extended it for other applications and used an NCO for the fixed 50 Hz complex tone generator, you end up with a tuneable single tone adaptive canceller that will capture and cancel any single tone within the proximity to any frequency setting of the NCO, within range of the bandwidth of the LPF! The cancelling loop will null the 50 Hz interferor down to the noise floor coming out of the filter- so the tighter the filter, the deeper the null as long as the the signal is within the filter bandwidth.

The amount of cancellation will depend on the time rate of change of the amplitude and phase of the input signal. The approach as shown above has zero delay to the main signal of interest, but is estimating the characteristics of the interference from a lag time T in the past, where T is impacted by the bandwidth (delay) of the low pass filter and the loop BW of the nulling loop (set by k), so will have an upper frequency limit on how much of the noise bandwidth of the interference signal can be rejected (that noise bandwidth is visible in the power spectral density of the interfering signal, basically AM and PM noise side-bands). This nulling loop is a high pass filter to that interference; which is desired as the signal of interest is indistinguishable from noise.

If this approach proves useful, you can improve it further by instead using the method of recursive least squares (RLS adaptive algorithm) to determine the complex weight that is applied to the vector modulator based on the error terms out of the quadrature downconverter. This converges faster at the expense of computation complexity.


What is the difference between undersampling and oversampling in analog to digital conversion ? When converting a signal from analog to digital, we observe undersampling and oversampling.

Does oversampling means that the sampling frequency is greater that the signal's frequency and does undersampling means that the sampling frequency is less than the signal's frequency ?

Good question as you can actually undersample and oversample at the same time! See my "DSP Puzzle" question on that specifically here:
How do you simultaneously undersample and oversample?

To best explain undersampling and oversampling, it is worthwhile understanding the concept of "Nyquist Zones" first.  This was explained in detail recently at this post:

Higher order harmonics during sampling

From which I am copying some graphics for further explanation here.

Sampling in the 1st Nyquist Zone, as in the graphic below, occurs when the real analog signal is located anywhere in the frequency range of DC to half of the sampling rate. Sampling at twice the highest frequency (and actually twice the highest bandwidth, as we will see more clearly when describing "under-sampling") is the minimum necessary condition to capture all information within that bandwidth according to the Nyquist-Shannon theorem. In practice to allow for achievable anti-aliasing filtering in the analog we need to sample a little higher than this, typically 20 to 30%.

The graphic below shows the related processing spectrums when sampling an analog input signal that is in the "first Nyquist Zone".  The colors are related to other discussion that is detailed in the link I gave above but not important here.

Oversampling

The above introduced standard "Sampling".  "Over-Sampling" is when we sample higher than the condition set by Nyquist, and typically much higher for us to be concerned with the possible benefit of doing that.

The main thing that happens when we over-sample is the quantization noise from the conversion gets spread over a much wider frequency range, thus in most cases reducing the noise from quantization in our frequency band of interest and therefore acting as a converter with a higher number of bits!  I have explained this further at this link so will refer to that for further detail on over-sampling:
What are advantages of having higher sampling rate of a signal?

Under-sampling

An analog signal can be higher in frequency than the first Nyquist Zone and still be sampled as long as the analog bandwidths of the system (including the A/D converter itself) are wide enough to pass the signal.  Such is the case in the graphic below copied for the link I first gave above that details this further. Here we are sampling a signal that is in the 3rd Nyquist Zone, and the resulting digital signal is indistinguishable from a similar signal that was in the 1st Nyquist Zone. So assuming a wide enough analog bandwidth, we can select any Nyquist zone to sample through the design of a anti-alias filter prior to the A/D converter. For standard and over-sampling conditions we would use a low-pass filter. And for under-sampling conditions we would use a band-pass filter. With no filter all the Nyquist zones would be sampled, thus we would have significant degradation due to noise and interference folding (aliasing). This is specifically the important reason to always have an anti-alias filter prior to the A/D converter!  "IF-Sampling" and "Direct-Sampling" are implementations of under-sampling. IF-Sampling is when we sample the IF frequency in a super-heterodyne receiver architecture, and Direct-Sampling is when we have an A/D converter with a high enough input analog BW (or a low enough input RF signal) that we can sample the RF signal directly.

Over-sampling implies having many more samples than the highest frequency of interest, and under-sampling implies we are down-converting the bandwidth of interest with a higher harmonic of the sampling clock (effectively). Thus as given in the first link, we can simultaneously under-sample and over-sample such as the following example using GPS:

Consumer GPS approximately occupies 2 MHz of BW at a carrier of 1575.42 GHz.  With an A/D converter that had an analog input bandwidth that exceeded this carrier frequency, we could directly sample the amplified and filtered signal from a GPS antenna with 100 MHz sampling rate (for example), which would result in a (spectrally inverted) copy of the GPS signal at 16* 100 MHz -1575.42 MHz = 24.58 MHz; well positioned for a digital IF in our 100 MHz sampling space. This would be digitally down-converted to baseband and then filtered where we would benefit from the advantages of over-sampling at 100 MHz the GPS signal that occupies significantly less bandwidth.  Thus we under-sampled by down-coverting with the sampling process, and over-sampled by choosing a sampling rate much higher than the minimum necessary.


Best Approach to Rank Complex Magnitude Comparision Problem This is in reference to the responses for an efficient algorithm for the comparison of bounded fixed point complex numbers at this post.

Efficient Magnitude Comparison for Complex Numbers

See the details of that post for the objectives of the problem. I am now determining my approach to ranking the algorithms to determine who best met the objectives I was seeking, and welcome debate on the ranking approach before I dive in.

Key qualifying factors:

I will baseline a standard CORDIC approach (rotate both vectors to real axis and compare absolute magnitude), as well as what can be done with the use of multiplier operations. The algorithm must be more efficient than either of these approaches (using the equivalent score for the multipliers - see below).

Algorithm must be 100% correct for for magnitude differences within $|z_2- z_1| \gt e$ for any e, where Where $z_n = \sqrt{I_n^2 + Q_n^2}$ with I and Q as bounded signed integers and e is any positive real number>0. (it is understood that it will likely take more operations as e decreases; in fact it would be attractive to be more efficient for large e). The CORDIC is a good example of this, you can select any error bound at the expense of the number of iterations needed.

Acceptable answers may return incorrect results for $|z_2- z_1| \le e$, but a bonus score is included for implementations that provide equivalence resolutions given by the following definitions, with a higher score for tight equivalence:

Loose Equivalence:

$|z_1| \gt |z_2| + e$ returns 1

$|z_1| \lt |z_2| -e$ returns -1

$|z_2- z_1| \le e$ returns 0

Tight Binary Equivalence:

$|z_1| > |z_2|$  returns 1

$|z_1| < |z_2|$  returns -1

Tight Ternary Equivalence:

$|z_1| > |z_2|$  returns 1

$|z_1| < |z_2|$  returns -1

$|z_1| = |z_2|$  returns 0

The function prototype is

result = CompareMagntitudes( I1, Q1, I2, Q2 )

With return values of either $-1,0,1$ corresponding to $<, =, > $ of the first compared to the second (and $0, 1$ for $<, \ge$ for binary solutions).

Test cases will be run with bit ranges from $b = 8$ to $b = 32$ bits for I and Q but nothing in the algorithm should prevent implementation with any size b.

Consider the following closely spaced complex points A, B, C, D as depicted in the diagram below ($A=3+j4$, $B=4+j4$, $C=3+j5$, $D4+j5$).

The true radii are A = 5, B =5.66, C = 5.83 and D = 6.403. In this case, the algorithm should provide a solution to resolve all 4 with 100% confidence (setting e to be $e \le 0.17$ corresponding to the minimum distance between B and C), however it would be acceptable and even beneficial if the algorithm allowed for larger e with an associated efficiency gain.

For example if $e=0.5$ then the following results must be returned using the format $f(z_1,z_2)$ with relation to the function prototype given above:

$f(A,B) \rightarrow -1$

$f(C,A) \rightarrow 1$

$f(A,D) \rightarrow -1$

$f(B,D) \rightarrow -1$

Since all those points have a difference in magnitude from the origin > 0.5.

However the following would be acceptable:

$f(B,C) \rightarrow X$

Where X can be 1, 0 or -1 since B and C have a difference in magnitude from teh origin < 0.5.

The algorithm must be implementable with only the equivalent of standard Boolean operations, binary shifts and compares. Look-up tables if used would add to  buffer size considerations in the scoring.

QUESTION: Please suggest / justify alternative metrics (including alternate scoring the starting numbers I list in my answer, or completely different approaches. It is understood that ultimately there is a trade space and can't be a one size fits all simple score, so consideration is weighted toward the objectives in the original question.

Here are some of the latest results:

Algorithm         Correct    Time      Score    Penalties  Eggs
---------------   -------    ------    -------  ---------  ----
  Empty Economy    49.86     2.8104     472849    2378650    0
   Empty Deluxe     0.05     2.8311       1944  474168000  243
Starter Economy    89.75     2.9663     851367     486060    0
 Starter Deluxe    90.68     2.9764    1663118     441920  151

    Dan Beast 4    99.85     3.2622    1750076       7130  151
Cedron Unrolled   100.00     3.2721    1898616          0  243
  Cedron Deluxe   100.00     3.3255    1898616          0  243
 Cedron Revised   100.00     3.2128    1898616          0  243
   Olli Revised    99.50     3.1893    1728065      23880    0
  Olli Original    99.50     3.2464    1728065      23880    0

Cedron Multiply   100.00     3.2042    1898616          0  243
  Matt Multiply   100.00     3.3146    1898616          0  243

The timing for the contenders is too close and too noisy to show a clear favorite.  Benchmarking on the target platform would be much more useful now.

The code has been updated.  It is as it is.

import numpy as np
import timeit

# The passed arguments to a running horse.
#
#   ( I1, Q1 ) First  Complex Value (or Point)
#   ( I2, Q2 ) Second Complex Value (or Point)

# Its return values are
#
#   ( rc ) Comparison Result (Return Code)
#   (  l ) Locus of the Exit

# The return value can be one of these
#
#    -2  The first is for sure less than the second
#    -1  The first is likely less than the second
#     0  The two are equal for sure
#     1  The first is likely greater than the second
#     2  The first is for sure greater than the second
#
# Routines that only return {-1,1} can be called Economy
# Routines that only return {-2,0,2} can be called Deluxe
#

# How Scoring works
#
#   S Score
#   P Penalties
#   E Egg Count
#   W Wrong
#
#                    Correct      Marginal     Wrong
#       {-1,1}        S+=2         S+=1         W+=1,P+=10
#       {-2,0,2}      S+=4(E+=1)   S+=2,P+=10   W+=1,P+=1000
#
#

#====================================================================
#====================================================================
#    W A L K    O N S
#====================================================================
#====================================================================
def WalkOnOne( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
def WalkOnTwo( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
def WalkOnThree( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
#====================================================================
#    S T A R T E R    C O D E
#====================================================================
#====================================================================
def EmptyEconomy( I1, Q1, I2, Q2 ):

        return 1, 0

#====================================================================
def EmptyDeluxe( I1, Q1, I2, Q2 ):

        return 0, 0

#====================================================================
def StarterEconomy( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Return Results

        if x1 < x2:
           return -1, 0

        return 1, 0

#====================================================================
def StarterDeluxe( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              return 2, 0
           thePresumedResult = 1
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              return -2, 0
           thePresumedResult = -1
        else:
           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Return the Presumed Result

        return thePresumedResult, 2

#====================================================================
#====================================================================
#    C E D R O N ' S
#====================================================================
#====================================================================
def CedronRevised( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:

           if x1 + y1 >= x2 + y2:
              return 2, 0

           thePresumedResult = 2
           dx = x1 - x2

        elif x1 < x2:

           if x1 + y1 <= x2 + y2:
              return -2, 0

           thePresumedResult = -2
           dx = x2 - x1

        else:

           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- Bring Factors into 1/2 to 1 Ratio Range

        while dx <  sx:
              dx += dx

              if dy <  sy:
                 dy += dy
              else:
                 sy += sy

        while dy <  sy:
              dy += dy

              if dx <  sx:
                 dx += dx
              else:
                 sx += sx

#---- Use Double Arithmetic Mean as Proxy for Geometric Mean

        cx = sx + dx
        cy = sy + dy

        cx16 = cx << 4
        cy16 = cy << 4

        if cx16 - cx > cy16:
           return thePresumedResult, 2

        if cy16 - cy > cx16:
           return -thePresumedResult, 2

#---- X Multiplication

        px = 0

        while sx > 0:
           if sx & 1:
              px += dx

           dx += dx
           sx >>= 1

#---- Y Multiplication

        py = 0

        while sy > 0:
           if sy & 1:
              py += dy

           dy += dy
           sy >>= 1

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
def CedronUnrolled( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:

           if x1 + y1 >= x2 + y2:
              return 2, 0

           thePresumedResult = 2
           dx = x1 - x2

        elif x1 < x2:

           if x1 + y1 <= x2 + y2:
              return -2, 0

           thePresumedResult = -2
           dx = x2 - x1

        else:

           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Estimate First Multiplied Magnitude

        if y1 < (x1>>1):
           if y1 < (x1>>2):
              m1 = (x1<<8) - (x1<<1) \
                 + (y1<<5) + (y1<<1)
           else:
              m1 = (x1<<8) - (x1<<4) \
                 + (y1<<6) + (y1<<5) - (y1<<2) - (y1<<1)
        else:
           if y1 < (x1>>1) + (x1>>2):
              m1 = (x1<<8) - (x1<<5) - (x1<<2) - (x1<<1) \
                 + (y1<<7) + (y1<<3) - y1
           else:
              m1 = (x1<<7) + (x1<<6) + (x1<<1) \
                 + (y1<<7) + (y1<<5) + (y1<<3)

#---- Estimate Second Multiplied Magnitude

        if y2 < (x2>>1):
           if y2 <   (x2>>2):
              m2 = ( (x2<<7) - x2 \
                   + (y2<<4) + y2 ) << 1
           else:
              m2 = ( (x2<<7) - (x2<<3) \
                   + (y2<<5) + (y2<<4) - (y2<<1) - y2 ) << 1
        else:
           if y2 <   (x2>>1) + (x2>>2):
              m2 = ( (x2<<8) - (x2<<5) - (x2<<2) - (x2<<1) \
                   + (y2<<7) + (y2<<3) - y2 )
           else:
              m2 = ( (x2<<6) + (x2<<5) + x2 \
                   + (y2<<6) + (y2<<4) + (y2<<2) ) << 1

#---- Return Results (1000 is a temp hack value!)

        if m1 > m2 + (m2>>6):
           return 2, 2

        if m2 > m1 + (m1>>6):
           return -2, 2

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- X Multiplication

        px = 0

        while dx > 0:
           if dx & 1:
              px += sx

           sx += sx
           dx >>= 1

#---- Y Multiplication

        py = 0

        while dy > 0:
           if dy & 1:
              py += sy

           sy += sy
           dy >>= 1

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
def CedronDeluxe( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              return 2, 0
           dx = x1 - x2
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              return -2, 0
           dx = x2 - x1
        else:
           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Employ a DanBeast

        L1 = DanBeast_2_8_Level( x1, y1 )
        L2 = DanBeast_2_8_Level( x2, y2 )

#---- Early Out Return

        if L1 > L2 + (L2>>6):
           return 2, 1

        if L2 > L1 + (L1>>6):
           return -2, 1

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- Do the Multiplications

        px = UnsignedBitMultiply( sx, dx )
        py = UnsignedBitMultiply( sy, dy )

#---- Account for Swap

        if x1 > x2:
           thePresumedResult = 2
        else:
           thePresumedResult = -2

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
def DanBeastFour( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with Quick Exit

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              return 2, 0
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              return -2, 0
        else:
           if y1 > y2:
              return 2, 0
           elif y1 < y2:
              return -2, 0
           else:
              return 0, 0

#---- Estimate First Multiplied Magnitude

        if y1 < (x1>>1):
           if y1 < (x1>>2):
              m1 = (x1<<8) - (x1<<1) \
                 + (y1<<5) + (y1<<1)
           else:
              m1 = (x1<<8) - (x1<<4) \
                 + (y1<<6) + (y1<<5) - (y1<<2) - (y1<<1)
        else:
           if y1 < (x1>>1) + (x1>>2):
              m1 = (x1<<8) - (x1<<5) - (x1<<2) - (x1<<1) \
                 + (y1<<7) + (y1<<3) - y1
           else:
              m1 = (x1<<7) + (x1<<6) + (x1<<1) \
                 + (y1<<7) + (y1<<5) + (y1<<3)

#---- Estimate Second Multiplied Magnitude

        if y2 < (x2>>1):
           if y2 < (x2>>2):
              m2 = ( (x2<<7) - x2 \
                   + (y2<<4) + y2 ) << 1
           else:
              m2 = ( (x2<<7) - (x2<<3) \
                   + (y2<<5) + (y2<<4) - (y2<<1) - y2 ) << 1
        else:
           if y2 < (x2>>1) + (x2>>2):
              m2 = ( (x2<<8) - (x2<<5) - (x2<<2) - (x2<<1) \
                   + (y2<<7) + (y2<<3) - y2 )
           else:
              m2 = ( (x2<<6) + (x2<<5) + x2 \
                   + (y2<<6) + (y2<<4) + (y2<<2) ) << 1

#---- Return Results

        if m1 < m2:
           return -1, 2

        return 1, 2

#====================================================================
def CedronMultiply( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        x1 = abs( I1 )
        y1 = abs( Q1 )

        x2 = abs( I2 )
        y2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if y1 > x1:
           x1, y1 = y1, x1

        if y2 > x2:
           x2, y2 = y2, x2

#---- Primary Determination with X Absolute Differences

        if x1 > x2:

           if x1 + y1 >= x2 + y2:
              return 2, 0

           thePresumedResult = 2
           dx = x1 - x2

        elif x1 < x2:

           if x1 + y1 <= x2 + y2:
              return -2, 0

           thePresumedResult = -2
           dx = x2 - x1

        else:

           if y1 > y2:
              return 2, 1
           elif y1 < y2:
              return -2, 1
           else:
              return 0, 1

#---- Sums and Y Absolute Differences

        sx = x1 + x2
        sy = y1 + y2

        dy = abs( y1 - y2 )

#---- X Multiplication

        px = 0

        while dx > 0:
          if dx & 1:
             px += sx

          sx += sx
          dx >>= 1

#---- Y Multiplication

        py = 0

        while dy > 0:
          if dy & 1:
             py += sy

          sy += sy
          dy >>= 1

#---- Return Results

        if px > py:
           return thePresumedResult, 2

        if px < py:
           return -thePresumedResult, 2

        return 0, 2

#====================================================================
#====================================================================
#    O L L I    L I K E
#====================================================================
#====================================================================
def MyVersionOfOllis( I1, Q1, I2, Q2 ):

# Returns  ( c )
#
#    c   Comparison
#
#   -1     | (I1,Q1) |  <  | (I2,Q2) |
#    1     | (I1,Q1) |  >  | (I2,Q2) |
#
#    t   Exit Test
#
#    1     (Partial) Primary Determination
#    2     CORDIC Loop + 1
#    6     Terminating Guess

#---- Set Extent Parameter

        maxIterations = 4

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- (Partial) Primary Determination

        if I1 < I2 and  I1 + Q1 <= I2 + Q2:
           return -2, 1

        if I1 > I2 and  I1 + Q1 >= I2 + Q2:
           return 2, 1

#---- CORDIC Loop

        for n in range ( 1, maxIterations+1 ):
            newI1 = I1 + ( Q1 >> n )
            newQ1 = Q1 - ( I1 >> n )
            newI2 = I2 + ( Q2 >> n )
            newQ2 = Q2 - ( I2 >> n )

            I1 = newI1
            Q1 = abs( newQ1 )
            I2 = newI2
            Q2 = abs( newQ2 )

            s = n + n + 1

            if I1 <= I2 - ( I2 >> s ):
               return -1, 1 + n

            if I2 <= I1 - ( I1 >> s ):
               return 1, 1 + n

#---- Terminating Guess

        if I1 < I2:
           return -1, 7

        return 1, 7

#====================================================================
def MyRevisionOfOllis( I1, Q1, I2, Q2 ):

# Returns  ( rc, l )
#
#    c   Comparison
#
#   -1,-2   | (I1,Q1) |  <  | (I2,Q2) |
#    1, 2   | (I1,Q1) |  >  | (I2,Q2) |
#
#    t   Exit Test
#
#    1     (Partial) Primary Determination
#    2     CORDIC Loop + 1
#    6     Terminating Guess

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- (Partial) Primary Determination

        if I1 < I2 and  I1 + Q1 <= I2 + Q2:
           return -2, 1

        if I1 > I2 and  I1 + Q1 >= I2 + Q2:
           return 2, 1

#---- CORDIC Loop Head

        s = 3

        for n in range ( 1, 5 ):

#---- Apply the Rotation

            newI1 = I1 + ( Q1 >> n )
            newQ1 = Q1 - ( I1 >> n )

            newI2 = I2 + ( Q2 >> n )
            newQ2 = Q2 - ( I2 >> n )

#---- Attempt Comparison

            if newI1 <= newI2 - ( newI2 >> s ):
               return -1, 1 + n

            if newI2 <= newI1 - ( newI1 >> s ):
               return 1, 1 + n

            s += 2

#---- Advance the Values

            I1 = newI1
            I2 = newI2

            Q1 = abs( newQ1 )
            Q2 = abs( newQ2 )

#---- Terminating Guess

        if I1 < I2:
           return -1, 7

        return 1, 7

#====================================================================
#====================================================================
#    M A T T   L    L I K E
#====================================================================
#====================================================================
def MattMultiply( I1, Q1, I2, Q2 ):

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- Ensure the first value is rightmost

        swap = 0;

        if I2 > I1:
           swap = 4
           I1, I2 = I2, I1
           Q1, Q2 = Q2, Q1

#---- Primary determination

        if I1 + Q1 > I2 + Q2:
           return 2 - swap, 2
        else:
           DI = I1 - I2
           if DI < 0:
              tmp1 = -UnsignedBitMultiply( I1 + I2, -DI )
           else:
              tmp1 =  UnsignedBitMultiply( I1 + I2,  DI  )

           DQ = Q2 - Q1
           if DQ < 0:
              tmp2 = -UnsignedBitMultiply( Q1 + Q2, -DQ )
           else:
              tmp2 =  UnsignedBitMultiply( Q1 + Q2,  DQ  )

           if tmp1 == tmp2:
              return  0       , 2
           elif tmp1 > tmp2:
              return  2 - swap, 2
           else:
              return -2 + swap, 2

#====================================================================
#====================================================================
#    U T I L I T I E S
#====================================================================
#====================================================================
def UnsignedBitMultiply( a, b ):  # Smaller value second is faster.

        p = 0

        while b > 0:
           if b & 1:
              p += a

           a  += a
           b >>= 1

        return p

#====================================================================
def DanBeast_2_8_Level( x, y ):

        if         y+y   <  x:               # 2 y < x
           if     (y<<2) <  x:               # 4 y < x
              L = (x<<8)   -x-x \
                + (y<<5)   +y+y                  # y/x = 0.00 to 0.25
           else:
              L = (x<<8) - (x<<4) \
                + (y<<6) + (y<<5) - (y<<2) -y-y  # y/x = 0.25 to 0.50
        else:
            if    (y<<2) <  x+x+x:           # 4 y < 3 x
              L = (x<<8) - (x<<5) - (x<<2) -x-x \
                + (y<<7) + (y<<3) -  y           # y/x = 0.50 to 0.75
            else:
              L = (x<<7) + (x<<6)   +x+x \
                + (y<<7) + (y<<5) + (y<<3)       # y/x = 0.75 to 1.00

        return L

#====================================================================
#====================================================================
#    T E S T I N G    H A R N E S S
#====================================================================
#====================================================================
def Test( ArgLimit, ArgThreshold, ArgLane, ArgTestName ):

#---- Set the Parameters

        t = ArgThreshold

#---- Initialize the Counters

        theCount           = 0
        theWrongCount      = 0

        theEggs            = 0
        theScore           = 0
        thePenalties       = 0

#---- Start Timing

        theStartTime = timeit.default_timer()

#---- Test on a Swept Area

        for i1 in range( -ArgLimit, ArgLimit, 10 ):
          ii1 = i1 * i1
          for q1 in range( -ArgLimit, ArgLimit, 7 ):
            d1 = np.sqrt( ii1 + q1 * q1 )
            for i2 in range( -ArgLimit, ArgLimit, 11 ):
              ii2 = i2 * i2
              for q2 in range( -ArgLimit, ArgLimit, 5 ):
                d2 = np.sqrt( ii2 + q2 * q2 )

                D = d1 - d2    #  = |(I1,Q1)| - |(I2,Q2)|

                theCount += 1

#---- The Fast Side Bench Mark Lanes

                if ArgLane == 0:
                   rc, l = EmptyEconomy( i1, q1, i2, q2 )

                if ArgLane == 1:
                   rc, l = EmptyDeluxe( i1, q1, i2, q2 )

                if ArgLane == 2:
                   rc, l = StarterEconomy( i1, q1, i2, q2 )

                if ArgLane == 3:
                   rc, l = StarterDeluxe( i1, q1, i2, q2 )

#---- The Slower Pace Horses

                if ArgLane == 8:
                   rc, l = TwoMultiply( i1, q1, i2, q2 )

                if ArgLane == 9:
                   rc, l = FourMultiply( i1, q1, i2, q2 )

#---- Walk Ons

                if ArgLane == 11:
                   rc, l = WalkOnOne( i1, q1, i2, q2 )

                if ArgLane == 12:
                   rc, l = WalkOnTwo( i1, q1, i2, q2 )

                if ArgLane == 13:
                   rc, l = WalkOnThree( i1, q1, i2, q2 )

#---- Cedron D.'s Lanes

                if ArgLane == 20:
                   rc, l = CedronRevised( i1, q1, i2, q2 )

                if ArgLane == 21:
                   rc, l = CedronDeluxe( i1, q1, i2, q2 )

                if ArgLane == 22:
                   rc, l = CedronUnrolled( i1, q1, i2, q2 )

                if ArgLane == 23:
                   rc, l = DanBeastFour( i1, q1, i2, q2 )

                if ArgLane == 24:
                   rc, l = CedronMultiply( i1, q1, i2, q2 )

#---- Olli N.'s Lanes

                if ArgLane == 30:
                   rc, l = MyVersionOfOllis( i1, q1, i2, q2 )

                if ArgLane == 31:
                   rc, l = MyRevisionOfOllis( i1, q1, i2, q2 )

#---- Dan B.'s Lanes

#                if ArgLane == 41:
#                   rc, l = Dan1( i1, q1, i2, q2 )

#---- Matt L.'s Lanes

                if ArgLane == 50:
                   rc, l = MattMultiply( i1, q1, i2, q2 )

#---- Assess Scores, Penalties, and Eggs

                if rc == -2:
                   if D < -t:
                      theScore      += 4
                   elif D < 0:
                      theScore      += 2
                      thePenalties  += 10
                   else:
                      theWrongCount += 1
                      thePenalties  += 1000

                elif rc == 2:
                   if D > t:
                      theScore      += 4
                   elif D > 0:
                      theScore      += 2
                      thePenalties  += 10
                   else:
                      theWrongCount += 1
                      thePenalties  += 1000

                elif rc == -1:
                   if D < 0:
                      theScore      += 2
                   elif D <= t:
                      theScore      += 1
                   else:
                      theWrongCount += 1
                      thePenalties  += 10

                elif rc == 1:
                   if D > 0:
                      theScore      += 2
                   elif D >= -t:
                      theScore      += 1
                   else:
                      theWrongCount += 1
                      thePenalties  += 10

                elif rc == 0:
                   if abs( D ) <= t:
                      theScore      += 8
                      if D == 0:
                         theEggs    += 1
                   else:
                      theWrongCount += 1
                      thePenalties  += 1000

                else:
                   print "Disqualification -- Invalid c value:", c, "Lane", ArgLane
                   return

#---- Finish Timing

        theDuration = timeit.default_timer() - theStartTime

#---- Calculate the Results

        theCorrectCount = theCount - theWrongCount

        theCorrectPct   = 100.0 * float( theCorrectCount ) \
                                / float( theCount )

#---- Return Results

        return "%15s  %7.2f %10.4f %10d %10d %4d" % \
               ( ArgTestName, theCorrectPct, theDuration,\
                 theScore,    thePenalties,  theEggs )

#====================================================================
def Main():

#---- Set Run Time Parameters

        L = 101   # The Limit
        T = 0     # Threshold

#---- Print Headers

        print "Algorithm         Correct    Time      Score    Penalties  Eggs"
        print "---------------   -------    ------    -------  ---------  ----"

#---- The Calibrators

        print Test( L, T, 0, "Empty Economy" )
        print Test( L, T, 1, "Empty Deluxe" )
        print Test( L, T, 2, "Starter Economy" )
        print Test( L, T, 3, "Starter Deluxe" )

#---- The Walk Ons

#        print
#        print Test( L, T, 11, "Walk On One" )

#---- The Contenders

        print
        print Test( L, T, 23, "Dan Beast 4" )
        print Test( L, T, 22, "Cedron Unrolled" )
        print Test( L, T, 21, "Cedron Deluxe" )
        print Test( L, T, 20, "Cedron Revised" )
        print Test( L, T, 31, "Olli Revised" )
        print Test( L, T, 30, "Olli Original" )

#---- The Pace Setters

        print
        print Test( L, T, 24, "Cedron Multiply" )
        print Test( L, T, 50, "Matt Multiply" )

#====================================================================
Main()

Earlier, I pledged a 50 point bounty to the best horse (fastest time 99%+ correct) that wasn't one of my entries.  I'm sticking with that, and right now Olli is leading.  (My optimized version is DQ'd)


Omega-K algorithm implementation for synthetic aperture radar I have been trying to develop Omega-k algorithm for SAR image formation. I am using the equations from chapter 8 in digital processing of Synthetic Aperture Radar data by Cumming and Wong. The steps I implemented are:

Range compression
Range IFFT
Two DFFT
Bulk Azimuth compression to a reference range (using eqn. 8.3)
Stolt Mapping (using eq. 8.5)
IFFT

I have confusions about the steps 4 to 6. Is there any norm for choosing the reference range in eqn. 8.3? I substituted eqn. 8.3 with 8.5 for stolt mapping and tried to interpolate the 2dfft spectra. However, the interpolation doesn't work. Do I need to compress the signal again after stolt mapping?

Any relevant answers will be highly appreciated. Thanks.

I am using this link for getting the reference by Richard Balmer :

There are no references 30 and 32


How many types of Bandpass filters are there for image processing? This article shows that:

Ideal Bandpass filter
Butterworth Bandpass filter
Gaussian Bandpass filter

Is that classification correct?

Are there any other types?

That article discusses just 3 types of filters (honestly a 2D Butterworth filter is the first time I hear)  and is clearly lacking any generality. 1D analog and digital filter design is a mature subject. However at some points 2D filter design dramatically differs from a 1D filter design,  as reflected by the active research on the area.

Fundamentally image data is exceedingly sensitive to phase distortion and any operations (including filtering) must therefore preserve the waveform shape as much as possible. (i.e. they have minimum, if not zero, phase distortion.)

Both 1D and 2D FIR filters are prime candidates for this required linear phase property. Their extension to 2D is straightforward.

A Butterworth filter is basicaly an IIR type and hence it has nonlinear phase response especially at the transition band edges. However, butterworth filters do have the smoothest passband characteristics among the basic filter types.

In contrast Chebyshev types I-II and elliptic filters are IIR filter types, with Equiripple characteristics in the passband and/or stopband. They have much sharper transitions than attainable via any FIR or Butterworth filter given the same order N.Elliptic filters are considered as the optimal IIR filter which will have smallest transition width for a given order N, at an accepted maximum passband and stopband deviation. However I have not seen any of these 2D IIR filter types to be preferred over a carefully optimized 2D FIR filter.

For 2D filter design, there are also other methods such as frequency sampling or 1D to 2D circular transforms. A particular class of such mappings can produce practically optimal 2D filters.

Finally a variation of Parks-McClellan minimax filters can also be adopted to 2D , however with more difficulties, since fundemantal theorems of 2D polynomials are more complex and less fruitful than the 1D case.


Shortest geometric distance from surface in 3d dataset? I have a three-dimensional binary image of a collection of discrete, individual voxels ("seeds") contained in a connected 3-dimensional surface ("skin").
(Like a small fruit, with a surface delineated by a one-pixel boundary, that contains seeds.)

The binary matrix was derived from a three-dimensional intensity image of the same size (grayscale).
(It's easy to label the surface voxels "2", the seed voxels "1", and the background "0" using MATLAB's label function).

In MATLAB,
 How do I compute the shortest geometric distance between each individual "seed" and the surface?
In the (crude) 2-dimensional diagram in the link below, each desired distance is represented by the magnitudes of the red lines (in number of voxels).

Do you just need the distance, or do you need the closest point?

For the closest point, the FLANN library can help, and it has Matlab bindings.

If you only need the distance, you can also use a distance transform. Try googling for "distance transform 3d matlab" for implementations.

Which one is faster depends on the number of "seeds" and "skin voxels".


Amplitude estimation of sinusoid in known spiky spectral noise What is the "best" way to estimate amplitude of a known-frequency sinusoid in the presence of known spiky spectral noise (i.e. noise comprising a few spectral peaks at known frequencies)?

By "best", I mean highest accuracy and lowest variance for a given sampling period (assume the sampling rate is greater than the Nyquist frequency).
The phases of the spectral noise peaks are unknown, but the phase of the sinusoid of interest can be provided if it's useful.

The approaches I am aware of are:

Perform a DFT at the frequency of interest, and design the window function such that its nulls/zero-crossings are located at the known noise frequencies.
Use a least squares estimator (which is basically a DFT with a rectangular window?, and which wrongly assumes that the noise samples are uncorrelated), e.g. as described in this survey paper.

I can't think of a better way than method (1), which does not benefit from phase information, but I'm wondering if a better approach could be used which may even take advantage of phase information of the signal of interest.

Just as a note, though I don't quite understand why, I noticed that the nulls of a rectangular window fall at integer multiples of the frequency that completes exactly one cycle in the sampling period.

The best fit time domain solution can be found by constructing two two basis vectors with your known frequency and calculate the coefficients directly.  The magnitude and phase can then be directly determined from these values.

Let C be a vector of cosine values over your frame and S be a vector of sine values.  You then want to find $(a,b)$ so that $aC + bS$ is a close as possible to your function.

$$ F = a C + b S $$

Dot this with the basis vectors.

$$ C \cdot F = a C \cdot C + b C \cdot S $$

$$ S \cdot F = a S \cdot C + b S \cdot S $$

Since the dot products are scalars, this is simply a linear system of two equations with two unknowns.  If you have a whole number of cycles in your frame then $ C \cdot S $ and $ S \cdot C $ become zero and the solution is trivial.

Note, this operation is the equivalent of a single bin in a DFT with the right parameters.  That is, pick a DFT frame with a whole number of cycles of your known frequency, and look at the corresponding bin.

The general case, without cycle alignment, can also be solved in the frequency domain with way fewer calculations.  See my blog article:

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

Do this for each of your peaks.  For better results, subtract out the peaks you have found before processing the others.


MATLAB Hamming window I have been trying to achieve a Hamming window length of 23 using the equation.
I tried several coding, one of the is the one below:
M=23
w = .54 + .46*cos(pi*(-M:M-1)'/M)
plot(w)

If the length of the window should be $23$, $M$ must equal $(23-1)/2=11$:

M = 11;
m = -M:M;
w = .54 + .46*cos(pi*m/M);
plot(m,w)


Acoustic room simulation - Craven's hypothesis I'm reading Udo Zölzer's "Digital Audio Signal Processing" and there's a piece about acoustic room simulation.

This question will be more about maths and acoustics, but I assume anybody who does DSP with sound knows the topic.

I'm trying to understand an excerpt about Craven's hypothesis (conerning human's sound source distance evaluation) which is used at Gerzon's algorithm.

It looks like this:

My understanding: Craven was trying to say that as we perceive signals with our ears, we can easily evaluate two things:

time difference between arrivals of the direct signal and the first reflection
the ratio of their amplitudes

Then, he says, our brain evaluates these two and can determine (more or less) the distance of the sound source.

I don't get the thing with g. It's said, that g stands for relative amplitude of direct signal to first reflection.

Shouldn't it be the other way? I mean, shouldn't it stand for "relative amplitude of first reflection to direct signal"? If we assume same SPL at both sources, the source that is further provides a weaker signal, than the one that is closer to the receiver.

If I would know only the verbal description of g, I would expect it to always be greater than 1, as direct signal is always louder than any reflected.

But how does it happen, that
g = SPL/SPL' = d/d'
?

According to inversed square distance law for sound, if we use SPLs, sound pressure level at distance r:

SPL(r) = (SPL@source)/r

So, if g would stand for the "relative amplitude of direct signal to first reflection" I would calculate:

g = SPL/SPL' = (P/d) / (P/d') = d'/d

SPL is Sound Pressure Level of the direct sound at the receiver,
SPL' is SPL of the first reflection at the receiver, P is sound source's SPL, d and d' are distances from sources.

Question 1.: is there a mistake at my interpretation or in the description in the book?

Question 2.: after all the calculations, when we have the formula for d, that Craven says allows our brain to evaluate the distance of the source, where does brain take the speed of sound from to calculate it? Should it be trained?

Question 1: The equations (6.7), (6.8), and (6.9) are all consistent, so I suggest that the wording for the definition of $g$ is incorrect in the book.

Question 2: The brain just hears the time difference and, by experience, can figure out what that means in terms of distance. So the brain has effectively "calibrated" for whatever the usual speed of sound in air is.  I find that hearing things under water (where the speed of sound is about five times faster) makes things seem much closer.


Explanation of FIR and IIR Filters with Simple Words It seems I am unable to understand how these filters work and what is the point of using them, no matter how many tutorials I watched.

Could someone please explain, in simple words, what FIR and IIR filters do?

Filters modify the frequency content of a signal.

The terms Finite Impulse Response (FIR) filter and Infinite Impulse Response (IIR) filter refer to Linear Time Invariant (LTI) filters.

LTI filters can change the amplitude and phase of specific frequencies in a signal, but they cannot shift any frequency content to other frequencies.

LTI filters are usually used to remove undesired frequencies from an input signal.

That's it.

As an example, your FM radio receives a signal that is the entire FM broadcast band, from 88 MHz to 108 MHz in the electromagnetic frequency spectrum, at its antenna.  In a simplified view of an FM radio receiver, an LTI filter is used to filter out everything but the 200 kHz around 99.5 MHz, so that the radio can demodulate the audio from the station at 99.5 MHz without the frequencies from all the other stations interfering with that process.  (In actuality FM radio receivers are more complicated, using a cascade of LTI filters with some non-linear "mixing" processes, but we're going for the concept here.)

To describe some of the math, at a high level:

When FIR and IIR filters are applied to signal in time domain, their application is equivalent to a mathematical operation called "Convolution".

The mathematical operation of convolution of a signal with an FIR or IIR filter, is equivalent to time-reversing the filter response, shifting it relative to the signal, performing pointwise multiplication, then summing up the values; then shifting the time-reversed filter response over slightly relative to the signal, preforming pointwise multiplication, then summing up the values; then ... keep repeating that for all of time (or until the input signal disappears to 0).


why does the decision boundary change in BPSK? is it true? I am working on modulation and demodulation of BPSK using python. i made a code for modulation and demodulation,while demodulate i tried to compare my actual data bit stream with demodulated bit stream,when no noise is added to the modulated signal i have zero errors after i demodulate.
But when i tried to send the modulated wave(iam working on sound waves) from one PC and i used another PC to receive the signal(using pyaudio) and when i try to apply the same demodulation method on the received sound i am getting 50% error(which i think the total demodulation technique is wrong or decision boundary changed because of noise that get added while recording the sound from transmitting pc).

So i need help to figure out which is causing error in demodulation(error i mean after comparing my sent data and the data that i get after demodulation)

is it possible like there will be shift in decision boundary after adding noise in it? if so how can we find the shifted decision boundary for my BPSK.

thanks in advance, any help is greatly appreciated :).

Zero-mean noise by itself can't modify the decision boundary. However, the number of things that can go wrong in your system is large.

Is your channel flat or frequency-selective?

You can think of the channel as the sound card and its driver, the loudspeaker, the air, the microphone, and the receiver's sound card and its driver. In my experience, there's very large variability in the frequency response of such a system. You can find the channel's response (for example by transmitting noise or an impulse) and use only the frequency band where it's flat. Or, you can design a wide-band system, but then you'll need an equalizer in the receiver.

Is your carrier synchronizer working?

Plot the received signal and see if its envelope is constant. If it's changing over time, you may have a faulty synchronizer. One quick way to get around this problem is to use DSB-LC modulation and demodulate the signal with an envelope detector.

Are you using a matched filter?

You need one to get optimal signal-to-noise ratio. After your filter, plot the signal's eye diagram. If it's open, then you know that your system up to this point is correct.

Is your symbol timing correct?

Since the clocks of the transmitter and receiver are different, you need to do adaptive symbol timing recovery to properly sample the matched filter's output. To test this, plot the sampled symbols over a long time. If they are correct for some time, but then start to get closer and closer, your symbol timing recovery is off.

Is your frame synchronization correct?

Finally, you need to verify your frame synchronization is correct; otherwise it's very hard to align the transmitted and received streams for bit error calculations.

I'd suggest verifying the correct operation of your system block by block, in the order I have listed them.


Matlab FFT not producing symmetric spectrum I am plotting a FFT of a sampled RC pulse but my resulting spectrum isn't symmetric - it's offset.

pulse = rcosdesign(.5, 6, 100)
T = .0001
ft = fft(pulse);
ft = fftshift(ft)
ft_mag = abs(ft);
freq = 1/(T*length(pulse))*(-length(pulse)/2:(length(pulse)/2-1))

plot(freq,ft_mag)

What's going on?

The problem comes from the DFT size. Try the code below that uses power of 2 fft.

T = 1e-4;
fs = 1/T;
pulse = rcosdesign(.5, 6, fs);
L = length(pulse);
Nfft = 2^nextpow2(L);
ft = fft(pulse,Nfft);
ft = fftshift(ft);
ft_mag = abs(ft);
freq=(-Nfft/2:Nfft/2-1)/Nfft*fs;
plot(freq,ft_mag)
axis([-3 3 0 max(ft_mag)])


What are the restrictions on DFT input sequences in order to produce valid spectral results In a recent comment (as I understand it) our colleague OverLordGoldDragon’s stated that the only restrictions on DFT input sequences, to produce valid DFT spectral results, are that (1) the input sample values must be finite in amplitude and (2) finite in length. That sounds sensible to me.
Are there are other restrictions on DFT input sequences necessary to produce valid spectral results?

The only restriction imposed on the input of the discrete finite dimensional transform DFT (Discrete Fourier Transform) is that the input has finite energy ; i.e., $$\sum_{n=0}^{N-1} |x[n]|^2 < \infty$$
And a necessary & sufficient condition for an $N$-point sequence $x[n]$ to have finite energy is that all of its samples are bounded.
There's no more restriction on the discrete-time sequence $x[n]$ to be transformable into $X[k]$ by means of the linear mapping implied by DFT.
However, I guess, from your previous question, that you want to investigate whether periodicity is a restriction on the input or not? The (my) answer is of course not. DFT works on finite length input, and since it's impossible to tell (by anyone) whether a block of data belongs to one period of a periodic waveform, or it's just the whole data itself, hence it's not possible to restrict the input to be periodic (because it cannot be tested)...
But the DFT algorithm has been designed to operate on one period of periodic inputs. If you evaluate the inverse DFT of a given set of DFT coefficents for all $n$, you can see that data periodically extends outside of the base interval $0 \leq n < N$.