Correlation between slightly misaligned images I am matching camera images of projection screens with rendered pages of PDF documents using Pearson's and Spearman's correlation coefficient. When the images (3, 4, 5) are aligned, this works fine in spite of the clear shortcomings (optical aberration in the camera images, the random sample is not independent, uneven illumination, etc.).

In the figure below, Pearson's r between the left and middle images is 0.8395 and Spearman's r is 0.6499. Pearson's r between the left and right images is 0.7139 and Spearman's r is 0.5066.

However, when the content of the rendered PDF page differs slightly from the document in the camera image or when the images (8, 9, 10) are more than slightly misaligned, the correlation coefficient ceases to be useful.

In the figure below, Pearson's r between the left and middle images is 0.771 and Spearman's r is 0.4719. Pearson's r between the left and right images is 0.779 and Spearman's r is 0.4963.

Reducing the images to density estimates (i.e. discarding the pixel positions and comparing histograms) would likely remove valuable information. Non-feature-based image registration seems like a more promising direction. What technique would you recommend to attack this problem?

As suggested by Dan Boschen, I maximized the Spearman's correlation coefficient over all possible “small” translations – where small is less than 15 pixels in L1-distance – of the rendered PDF page images. The results were underwhelming, but I suppose this method would work much better if the camera images weren't optically distorted.

In the figure below, Pearson's r between the left and middle images is 0.789 and Spearman's r is 0.479. Pearson's r between the left and right images is 0.8008 and Spearman's r is 0.5094.

As suggested by Marcus Müller, I used a feature detector (ORB, not SURF) to find keypoints in the images. Rather than directly use the feature descriptors for matching, I used them to find a homography that registers the two images.

In the figure below, Pearson's r between the left and middle images is 0.6244 and Spearman's r is 0.4921. Pearson's r between the left and right images is 0.0213 and Spearman's r is 0.036. For my test case, this works remarkably well.


Effect of windowing on FFT magnitude I have 2 sin waves:

Fs = 1500;
N = 100;
t = (0:N-1)'/Fs;

sinwave = sin(2*pi*15*t);
sinwaveshifted = sin(2*pi*15*t + 1.3);

If I get the magnitude at 15Hz using an fft (need to use an FFT that is a power of 2 for reasons beyond my control):

fftsin = fft(sinwave,128);
fftsinshift = fft(sinwaveshifted,128);

% 2 corresponds to 15Hz

abs(fftsin(2)) = 51.819;
abs(fftsinshift(2)) = 41.6151;

Why the difference?  I'm assuming because my bins do not line up perfectly there is some leakage occurring. Also, I haven't done any windowing.

If I do window my signals (I used a flat top because I read that was best in terms of magnitude accuracy), these are my results:

w = window(@flattopwin, 100);
sinwin = sinwave .* w;
sinshiftwin = sinwaveshift .* w;

fftsin = fft(sinwin ,128);
fftsinshift = fft(sinshiftwin ,128);

% 2 corresponds to 15Hz

abs(fftsin(2)) = 2.7283;
abs(fftsinshift(2)) = 17.822;

Why the huge difference?

I presume you are correct in that the small difference in the original signal fft's is due to leakage. However, the huge difference in the windowed results arises because you simply do not have enough data. Your original signals only have data for one period of the sinusoid. When you window this, you are clamping the edges of your signal to zero, effectively leaving you with only the middle portion of the signal. So you are left with signals that have only a fraction of a single period worth of information, definitely not enough to have any accuracy in the extracted frequency analysis. To illustrate this, here are matlab images of your original functions before and after windowing:

and you can easily see the loss of much of the meaningful frequency information from the signals. A better candidate for windowing would be something like below, using 10 periods of your sinusoid signals:

which demonstrates the balance between retaining frequency information and clamping the edges to reduce leakage. This paper provides a more in depth look at different windowing functions and their performance on different classes of signals if you are interested.


Plotting the magnitude response of a filter? I have designed a second order butterworth filter.

Sampling frequency: $4\textrm{ kHz}$,   Cut-off frequency: $500\textrm{ Hz}$

$$e = \tan\left(\frac{\pi \times f_c}{f_s}\right) = \tan\left(\frac{500\pi}{4000}\right) = 0.41421$$

Difference equation:

$$y[n] = 0.0976\cdot x[n] + 0.1952\cdot x[n-1] + 0.0976\cdot x[n-2] + 0.9429\cdot y[n-1] - 0.3334\cdot y[n-2]$$

How can I plot the magnitude response of this filter to prove that cut-off occurs at $500\textrm{ Hz}$?

If you want to code your own application that plot the magnitude response, you first need to extract the poles and zeros from your transfer function in the $Z$ domain. The process that follows can either be analytic or graphical. I will try to cover both, starting with the analytic approach, then graphical.

Extracting the poles and zeros

Taking your time domain equation. :
$$y[n] = 0.0976⋅x[n]+0.1952⋅x[n−1]+0.0976⋅x[n−2]+0.9429⋅y[n−1]−0.3334⋅y[n−2]$$
Transfer function in the Z domain is
$$H(z)=\frac{Y}{X}=\frac{0.0976+0.1952z^{-1}+0.0976z^{-2}}{1-0.9429z^{-1}+0.3334z^{-2}}=\frac{0.0976(1+2z^{-1}+1z^{-2})}{1-0.9429z^{-1}+0.3334z^{-2}}$$

Solving the numerator and denominator for $z=0$ will yield some poles,zeros and gains. I will not detail this step as this subject is widely covered. In your case, you will find:
$$zeros = \{-1, -1\}$$
$$poles = \{(0.4746 + 0.3289j), (0.4746 - 0.3289j)\}$$
$$K_{num}=0.0976$$
$$K_{denom}=1$$
The analytic approach

Once you have the poles and zeros, you can rewrite your transfer function in this different form :

$$H(z)=\frac{K_{num}}{K_{denom}}\cdot\frac{(z-zero_{1})(z-zero_{2})...(z-zero_{n})}{(z-pole_{1})(z-pole_{2})...(z-pole_{n})}$$

Or in a condensed form :
$$H(z) = \frac{K_{num}}{K_{denom}}\cdot\frac{\prod_{n=0}^{N_{zeros}}{(z-q_{n})} }{\prod_{n=0}^{N_{poles}}{(z-p_{n})}}$$

The magnitude response of your filter is basically the magnitude of your transfer function when $z=e^{j\omega}$. We can define  $|H(z)|\biggr\rvert_{z=e^{j\omega}}$

You then get :
$$|H(z)|\biggr\rvert_{z=e^{j\omega}} = \biggr\rvert \frac{K_{num}}{K_{denom}}\cdot \frac{\prod_{n=0}^{N_{zeros}}{(e^{j\omega}-q_{n})} }{\prod_{n=0}^{N_{poles}}{(e^{j\omega}-p_{n})}}\biggr\rvert$$

Translate that into a code, and you get something like : (matlab example)

h = ones(1,n);
knum=0.0976
kdenom=1
%Zeros
for k=1:length(q)
    v = exp(1i*w)-q(k);
    h = h .* v;
end
%Poles
for k=1:length(p)
    v = exp(1i*w)-p(k);
    h = h ./ v;
end
h = h * knum/kdeom
plot(abs(h))

The graphical approach

What we will see here, is exactly what we just saw in the analytic approach, but we will try to visualize it a little bit. Let's plot you poles and zeros into the $Z$ plane:

The unit circle, or $z=e^{j\omega}$, contains all the frequencies from $\omega=0$ to $\omega=Nyquist=\frac{2\pi f_{s}}{2}$

In order to know the frequency response of your filter at a specific value of $\omega$. Draw a line from each poles/zeros to the corresponding point on the unit circle.

Take the products of the line length originating from a zero and divide by the product of the line length originating from a poles. You'll get the magnitude response of your filter.

In other word:

$$ |H(z)|\biggr\rvert_{z=e^{j\omega}} = \frac{d_{3}d_{4}}{d_{1}d_{2}}\cdot\frac{k_{num}}{k_{denom}}$$

Take few seconds to understand what we just did there, and you will see that it is exactly what the analytical approach suggests.

Finally

I wrote some matlab code to plot the frequency response of a filter not so long ago and I posted this question to get help getting this done. It might also help you.

Good luck


Why are LTI systems defined by convolution, why not in any other way? I am curious about why are LTI systems defined by convolution, why not in any other way?
This Might be silly question but I want to know that.

Also why polynomial multiplication is convvolution?

it takes more work to do this for continuous-time systems, so i'll just do it for discrete-time and we can all wave our arms to extend to continuous-time.

a general discrete-time input signal can be expressed as

$$ x[n] = \sum\limits_{i=-\infty}^{\infty} x[i] \delta[n-i] $$

where $\delta[n]$ is the Kronecker delta or "unit impulse function"

$$ \delta[n] \triangleq
\begin{cases}
1, & \text{if }n = 0 \\
0, & \text{if }n \ne 0
\end{cases} $$

so think of your input $x[n]$ as the sum of a bunch of individual impulse functions, each delayed to time $i$ and weighted by a constant which is equal to $x[i]$.

so, let's say we have a Linear and Time-Invariant discrete-time system, represented with operator, $\mathbb{T}$, and output:

$$ y[n] = \mathbb{T}\left\{ x[n] \right\} $$

and we input into that system a single unit impulse, $\delta[n]$.  by definition, the "impulse response" of the LTI system is the output coming out when a unit impulse, $\delta[n]$, is input.

$$ h[n] \triangleq \mathbb{T}\left\{ \delta[n] \right\} $$

now, imagine a general input, $x[n]$ is input to $\mathbb{T}\{ \}$.  then, because it's Linear:

$$ \begin{align}
y[n] & = \mathbb{T}\left\{ x[n] \right\} \\
& = \mathbb{T}\left\{ \sum\limits_{i=-\infty}^{\infty} x[i] \delta[n-i] \right\} \\
& = \sum\limits_{i=-\infty}^{\infty} \mathbb{T}\left\{ x[i] \delta[n-i] \right\} \\
& = \sum\limits_{i=-\infty}^{\infty} x[i]  \mathbb{T}\left\{ \delta[n-i] \right\} \\
\end{align} $$

and because it's Time Invariant,

$$ h[n-i] = \mathbb{T}\left\{ \delta[n-i] \right\} $$

then, putting the two together, you get

$$ \begin{align}
y[n] & = \sum\limits_{i=-\infty}^{\infty} x[i]  \mathbb{T}\left\{ \delta[n-i] \right\} \\
& = \sum\limits_{i=-\infty}^{\infty} x[i]  h[n-i] \\
\end{align} $$


Coefficients of Farrow structure? What is the method for generating coefficients for the Farrow structure in Matlab? I am designing a low pass filter using firpm, but the output of the Farrow structure is not delayed.

The farrow structure effectively interpolates between different sets of filter coefficients, so if you need an adjustable delay ($\delta$) of, say, 0 to 1 samples, you need to create a number (K) of FIR filters covering that delay interval. The matlab command then creates a $ K-1^{th} $ order polynomial in $ \delta $ to represent each filter coefficient.

As Matt L said, if you don't need an adjustable delay, or to interpolate between filter coefficients based on some other parameter, there's no need for the farrow structure. If, however, you're just interested in how to get from FIR coefficients to farrow structure coefficients in general, I can post some rough code as I was messing about with this the other day. (hint: it's just simple least squares fitting).


Generating Kawase Blur Kernels to Approximate a Gaussian Blur on an Image Is there an algorithm to generate the Kawase Blur Kernels to approximate the actual Gaussian Blur (of a specific kernel size, and if possible, a sigma)?

Looking at Intel - An Investigation of Fast Real Time GPU Based Image Blur Algorithms By Filip Strugar it looks like the Kawase kernel is just a way of implementing a linear kernel quickly, but in a way that constrains the kernel somewhat.

This means that you could make such an algorithm.  Either choose a set of spreads and adjust their weights (if that is possible) for a best fit, or search a number of possible spreads (with the best weights) and find the ones that fit best.

Given the ad-hoc nature of the kernel I suspect you'll end up with a table of best fits for each set of spreads, and then you (or anyone using the table) would want to make a judgement call about effect vs. clock ticks.

The kernel is designed with some specific hardware in mind, so you'll need to check to see if you can adjust the weights at all, or if you can only adjust the spacing and the repetition at each spacing.  At any rate, because the spacing is discrete, its a fairly highly nonlinear optimization problem.  It may be best to just do a brute-force search at each size of Gaussian you want to emulate.  Just try out a bunch of different spacing combinations, and see which ones give you the closest match to a Gaussian of a certain size for the number of iterations.


What are Autoregressive Coefficients? Can anyone explain what are Autoregressive Coefficients? What is their meaning that is.

Consider a method:

public double[] calculateARCoefficients(double[] inputseries, int order)

When this method is called on 256 values, lets say these values represent some kind of signal, I will get an array of n numbers, where n=order. What do these numbers represent?

I mean for example variance tells me how spread out the values of the signal are. The higher variance value is, the more spread out values are. The smaller variance is, the more unified signal is.

Are you able to say, in a similar fashion to above paragraph about variance, what value of AR coefficients says about the signal?

The coefficients $a_i$ of an autoregressive (AR) model of a signal $x(n)$ enter the model definition in the following way:

$$x(n)=\sum_{i=1}^Na_ix(n-i)+e(n)$$

where $e(n)$ is zero-mean white noise and $N$ is the model order. I.e. the signal is modeled as the output of a linear time-invariant filter with transfer function

$$H(z)=\frac{1}{1-\sum_{i=1}^Na_iz^{-1}}$$

and with input signal $e(n)$. There is no simple interpretation of specific values of the coefficients $a_i$. You have to look at the whole polynomial $1-\sum_{i=1}^Na_iz^{-1}$, e.g. by computing the power spectrum of the AR model of $x(n)$:

$$S_x(z)=\sigma_e^2\left|\frac{1}{1-\sum_{i=1}^Na_iz^{-1}}\right|^2,\quad z=e^{j\theta}$$

where $\sigma_e^2$ is the variance of $e(n)$. So the AR coefficients are the parameters of a specific signal model, and they do not have a simple interpretation comparable to mean or variance of a random process. However, in general they represent a much more detailed description of the signal than mean and variance alone.


Relation between Power Spectral Density and RMS of the signal I have a pressure signal at two locations (r=1 and r=1.3). I found the Power spectral density (PSD) and the Root Mean Square (RMS) of the siganl. My question is why at r=1 the PSD is higher than the PSD at r=1.3 while the RMS at r=1 is less than RMS at r=1.3.

I mean is there any relation between PSD and RMS of the signal?

instantaneous power is:

$$ p_x(t) = |x(t)|^2 $$

mean power is:

$$\begin{align}
 P_x &= \lim_{T \to \infty} \quad \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} p_x(t) \, dt \\
\\
&= \lim_{T \to \infty} \quad \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} |x(t)|^2 \, dt \\
\end{align}$$

the relation between Power Spectral Density to mean power is:

$$ P_x = \int\limits_{-\infty}^{\infty} S_x(f) \, df  $$

the relation between Root Mean Square to mean power is:

$$ P_x = \lVert x(t) \rVert^2 $$


Finding the maximum of a real valued trigonometric polynomial I have been posting all over SE looking for a good method of finding the maximum of a trigonometric polynomial. A trigonometric polynomial is the sum of a set of sinusoids,

$$a_0 + \sum_{k = 1}^{N} a_k \cos(k \theta + \phi_k)$$

and can also be written in the form

$$ \sum_{k=-N}^{N} c_k \, e^{i k \theta} $$
where $c_k$ = $\bar{c}_{-k}$

Substituting $z = e^{i \theta}$ turns this into the typical form of a complex valued polynomial.

The method I am currently using is to take the find the roots of the derivative of the polynomial, then evaluate the original polynomial at these locations to find the maximum. While quite fast with MATLABs roots command, it is not quite fast enough.

I need a method where I can trade accuracy for speed. Does anyone know of any other fast methods? This is a similar problem to estimating the crest factor of a signal. In my problem $N$ ranges from 3 to 15.

This function is periodic with 2*pi and it it's highest frequency component is N times the fundamental frequency. As such it can be sampled at a rate of 1/2*N without losing information. However, there is no guarantee that sampling will hit the maximum or that maximum sampled value is close to the analog maximum. The likelihood of this increases drastically by oversampling, e.g sample 4*N or 8*N values.

So one thing to try would be to basically calculate the function values on a equidistant grid of, say 8*N, values, locate the maximum and then do a binary search around the maximum until the desired accuracy is reached.

Again there is no guarantee that this will ALWAYS work, however you can try it with a few actual data sets and see how much over sampling is required and whether this is actually faster than the root search.


Efficient Magnitude Comparison for Complex Numbers Is there a more efficient algorithm (or what is the most efficient known algorithm) to select the larger of two complex numbers given as $I + jQ$ without having to compute the squared magnitude as

$$I^2+Q^2$$

In particular binary arithmetic solutions that do not require multipliers?  This would be for a binary arithmetic solution using AND, NAND, OR, NOR, XOR, XNOR, INV, and shifts and adds without simply replacing required multiplication steps with their shift and add equivalents (or what would be equivalent in processing steps). Also assume the number is represented with rectangular fixed point (bounded integers) coordinates (I, Q).

I am aware of magnitude estimators for complex numbers such as $max(I,Q) + min(I,Q)/2$, and more accurate variants at the expense of multiplying coefficients but they all have a finite error.

I have considered using the CORDIC rotator to rotate each to the real axis and then being able to compare real numbers. This solution also has finite error but I can choose the number of iterations in the CORDIC such that the error is less than $e$ for any $e$ that I choose within my available numeric precision. For that reason this solution would be acceptable.

Are there other solutions that would be more efficient than the CORDIC (which requires time via iterations to achieve accuracy)?

Determining Best Answer

There was incredible work done by the participants (and we still welcome competing options if anyone has other other ideas). I posted my proposed approach to ranking the algorithms and welcome debate on the ranking approach before I dive in:

Best Approach to Rank Complex Magnitude Comparision Problem

PROLOGUE

My answer to this question is in two parts since it is so long and there is a natural cleavage.  This answer can be seen as the main body and the other answer as appendices.  Consider it a rough draft for a future blog article.

Answer 1
  * Prologue (you are here)
  * Latest Results
  * Source code listing
  * Mathematical justification for preliminary checks

Answer 2
  * Primary determination probability analysis
  * Explanation of the lossless adaptive CORDIC algorithm
  * Small angle solution

This turned out to be a way more in depth and interesting problem than it first appeared.  The answer given here is original and novel.  I, too, am very interested if it, or parts of it, exist in any canon.

The process can be summarized like this:

An initial primary determination is made.  This catches about 80% of case very efficiently.
The points are moved to difference/sum space and a one pass series of conditions tested.  This catches all but about 1% of cases.  Still quite efficient.
The resultant difference/sum pair are moved back to IQ space, and a custom CORDIC approach is attempted

So far, the results are 100% accurate.  There is one possible condition which may be a weakness in which I am now testing candidates to turn into a strength.  When it is ready, it will be incorporated in the code in this answer, and an explanation added to the other answer.

The code has been updated.  It now reports exit location counts.  The location points are commented in the function definition.  The latest results:

   Count: 1048576

    Sure: 100.0
 Correct: 100.0

Presumed: 0.0
Actually: -1

  Faulty: 0.0

    High: 1.0
     Low: 1.0

1   904736   86.28
2     1192   86.40
3     7236   87.09
4    13032   88.33
5   108024   98.63
6     8460   99.44

Here are the results if my algorithm is not allowed to go into the CORDIC routine, but assumes the answer is zero (8.4% correct assumption).  The overall correct rate is 99.49%  (100 - 0.51 faulty).

   Count: 1048576

    Sure: 99.437713623
 Correct: 100.0

Presumed: 0.562286376953
Actually: 8.41248303935

  Faulty: 0.514984130859

    High: 1.05125
     Low: 0.951248513674

1   904736   86.28
2     1192   86.40
3     7236   87.09
4    13032   88.33
5   108024   98.63
6     8460   99.44

Okay, I've added an integer interpretation of Olli's algorithm.  I would really appreciate somebody double checking my translation into Python.  It is located at the end of the source code.

Here are the results:

   Count: 1048576

 Correct: 94.8579788208

1   841216   80.22  (Partial) Primary Determination
2    78184   87.68  1st CORDIC
3   105432   97.74  2nd
4    10812   98.77  3rd
5        0   98.77  4th
6    12932  100.00  Terminating Guess

Next, I added the "=" to the primary slope line comparisons.  This is the version I left in my code.

The results improved.  To test it yourself, simply change the function being called in the main() routine.

 Correct: 95.8566665649

1   861056   82.12
2   103920   92.03
3    83600  100.00

Here is a Python listing for what I have so far.  You can play around with it to your heart's content.  If anybody notices any bugs, please let me know.

import array as arr

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

#---- Initialize the Counters

        theCount      = 0
        theWrongCount = 0

        thePresumedZeroCount    = 0
        theCloseButNotZeroCount = 0

        theTestExits = arr.array( "i", [ 0, 0, 0, 0, 0, 0, 0 ] )

#---- Test on a Swept Area

        theLimit = 32
        theLimitSquared = theLimit * theLimit

        theWorstHigh = 1.0
        theWorstLow  = 1.0

        for i1 in range( theLimit ):
          ii1 = i1 * i1
          for q1 in range( theLimit ):
            m1 = ii1 + q1 * q1
            for i2 in range( theLimit ):
              ii2 = i2 * i2
              for q2 in range( theLimit ):
                m2 = ii2 + q2 * q2
                D  = m1 - m2

                theCount += 1

                c, t = CompareMags( i1, q1, i2, q2 )

                if t <= 6:
                   theTestExits[t] += 1

                if c == 2:

                   thePresumedZeroCount += 1
                   if D != 0:
                      theCloseButNotZeroCount += 1

                      Q = float( m1 ) / float( m2 )
                      if Q > 1.0:
                         if theWorstHigh < Q:
                            theWorstHigh = Q
                      else:
                         if theWorstLow > Q:
                            theWorstLow = Q

                      print "%2d  %2d   %2d  %2d   %10.6f" % ( i1, q1, i2, q2, Q )

                elif c == 1:
                   if D <= 0:
                      theWrongCount += 1
                      print "Wrong Less ", i1, q1, i2, q2, D, c
                elif c == 0:
                   if D != 0:
                      theWrongCount += 1
                      print "Wrong Equal", i1, q1, i2, q2, D, c
                elif c == -1:
                   if D >= 0:
                      theWrongCount += 1
                      print "Wrong Great", i1, q1, i2, q2, D, c
                else:
                   theWrongCount += 1
                   print "Invalid c value:", i1, q1, i2, q2, D, c

#---- Calculate the Results

        theSureCount   = ( theCount - thePresumedZeroCount )
        theSurePercent = 100.0 * theSureCount / theCount

        theCorrectPercent = 100.0 * ( theSureCount - theWrongCount ) \
                          / theSureCount

        if thePresumedZeroCount > 0:
           theCorrectPresumptionPercent = 100.0 * ( thePresumedZeroCount - theCloseButNotZeroCount ) \
                                        / thePresumedZeroCount
        else:
           theCorrectPresumptionPercent = -1

        theFaultyPercent = 100.0 * theCloseButNotZeroCount / theCount

#---- Report the Results

        print
        print "   Count:", theCount
        print
        print "    Sure:", theSurePercent
        print " Correct:", theCorrectPercent
        print
        print "Presumed:", 100 - theSurePercent
        print "Actually:", theCorrectPresumptionPercent
        print
        print "  Faulty:", theFaultyPercent
        print
        print "    High:", theWorstHigh
        print "     Low:", theWorstLow
        print

#---- Report The Cutoff Values

        pct = 0.0
        f   = 100.0 / theCount

        for t in range( 1, 7 ):
          pct += f * theTestExits[t]
          print "%d %8d  %6.2f" % ( t, theTestExits[t], pct )

        print

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

# This function compares the magnitudes of two
# integer points and returns a comparison result value
#
# Returns  ( c, t )
#
#    c   Comparison
#
#   -1     | (I1,Q1) |  <  | (I2,Q2) |
#    0     | (I1,Q1) |  =  | (I2,Q2) |
#    1     | (I1,Q1) |  >  | (I2,Q2) |
#    2     | (I1,Q1) | ~=~ | (I2,Q2) |
#
#    t   Exit Test
#
#    1     Primary Determination
#    2     D/S Centers are aligned
#    3     Obvious Answers
#    4     Trivial Matching Gaps
#    5     Opposite Gap Sign Cases
#    6     Same Gap Sign Cases
#   10     Small Angle + Count
#   20     CORDIC + Count
#
# It does not matter if the arguments represent vectors
# or complex numbers.  Nor does it matter if the calling
# routine considers the integers as fixed point values.

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

        a1 = abs( I1 )
        b1 = abs( Q1 )

        a2 = abs( I2 )
        b2 = abs( Q2 )

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

        if b1 > a1:
           a1, b1 = b1, a1

        if b2 > a2:
           a2, b2 = b2, a2

#---- Primary Determination

        if a1 > a2:
           if a1 + b1 >= a2 + b2:
              return 1, 1
           else:
              thePresumedResult = 1
              da = a1 - a2
              sa = a1 + a2
              db = b2 - b1
              sb = b2 + b1
        elif a1 < a2:
           if a1 + b1 <= a2 + b2:
              return -1, 1
           else:
              thePresumedResult = -1
              da = a2 - a1
              sa = a2 + a1
              db = b1 - b2
              sb = b1 + b2
        else:
           if b1 > b2:
              return 1, 1
           elif b1 < b2:
              return -1, 1
           else:
              return 0, 1

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

        db, sb = sb, db

        while da < sa:
            da += da
            sb += sb
            if sb > db:
               db, sb = sb, db

#---- Ensure the [b] Factors are Both Even or Odd

        if ( ( sb + db ) & 1 ) > 0:
           da += da
           sa += sa
           db += db
           sb += sb

#---- Calculate Arithmetic Mean and Radius of [b] Factors

        p = ( db + sb ) >> 1
        r = ( db - sb ) >> 1

#---- Calculate the Gaps from the [b] mean and [a] values

        g = da - p
        h = p - sa

#---- If the mean of [b] is centered in (the mean of) [a]

        if g == h:
           if g == r:
              return 0, 2;
           elif g > r:
              return -thePresumedResult, 2
           else:
              return thePresumedResult, 2

#---- Weed Out the Obvious Answers

        if g > h:
           if r > g and r > h:
              return thePresumedResult, 3
        else:
           if r < g and r < h:
              return -thePresumedResult, 3

#---- Calculate Relative Gaps

        vg = g - r
        vh = h - r

#---- Handle the Trivial Matching Gaps

        if vg == 0:
           if vh > 0:
              return -thePresumedResult, 4
           else:
              return thePresumedResult, 4

        if vh == 0:
           if vg > 0:
              return thePresumedResult, 4
           else:
              return -thePresumedResult, 4

#---- Handle the Gaps with Opposite Sign Cases

        if vg < 0:
           if vh > 0:
              return -thePresumedResult, 5
        else:
           if vh < 0:
              return thePresumedResult, 5

#---- Handle the Gaps with the Same Sign (using numerators)

        theSum = da + sa

        if g < h:
           theBound = ( p << 4 ) - p
           theMid   = theSum << 3

           if theBound > theMid:
              return -thePresumedResult, 6
        else:
           theBound = ( theSum << 4 ) - theSum
           theMid   = p << 5

           if theBound > theMid:
              return thePresumedResult, 6

#---- Return to IQ Space under XY Names

        x1 = theSum
        x2 = da - sa

        y2 = db + sb
        y1 = db - sb

#---- Ensure Points are in Lower First Quadrant (First Octant)

        if x1 < y1:
           x1, y1 = y1, x1

        if x2 < y2:
           x2, y2 = y2, x2

#---- Variation of Olli's CORDIC to Finish

        for theTryLimit in range( 10 ):
            c, x1, y1, x2, y2 = Iteration( x1, y1, x2, y2, thePresumedResult )
            if c != 2:
               break

            if theTryLimit > 3:
               print "Many tries needed!", theTryLimit, x1, y1, x2, y2

        return c, 20

#================================================
def Iteration( x1, y1, x2, y2, argPresumedResult ):

#---- Try to reduce the Magnitudes

        while ( x1 & 1 ) == 0 and \
              ( y1 & 1 ) == 0 and \
              ( x2 & 1 ) == 0 and \
              ( y2 & 1 ) == 0:
            x1 >>= 1
            y1 >>= 1
            x2 >>= 1
            y2 >>= 1

#---- Set the Perpendicular Values (clockwise to downward)

        dx1 =  y1
        dy1 = -x1

        dx2 =  y2
        dy2 = -x2

        sdy = dy1 + dy2

#---- Allocate the Arrays for Length Storage

        wx1 = arr.array( "i" )
        wy1 = arr.array( "i" )
        wx2 = arr.array( "i" )
        wy2 = arr.array( "i" )

#---- Locate the Search Range

        thePreviousValue = x1 + x2  # Guaranteed Big Enough

        for theTries in range( 10 ):
            wx1.append( x1 )
            wy1.append( y1 )
            wx2.append( x2 )
            wy2.append( y2 )

            if x1 > 0x10000000 or x2 > 0x10000000:
               print "Danger, Will Robinson!"
               break

            theValue = abs( y1 + y2 + sdy )

            if theValue > thePreviousValue:
               break

            thePreviousValue = theValue

            x1 += x1
            y1 += y1
            x2 += x2
            y2 += y2

#---- Prepare for the Search

        theTop = len( wx1 ) - 1

        thePivot = theTop - 1

        x1 = wx1[thePivot]
        y1 = wy1[thePivot]
        x2 = wx2[thePivot]
        y2 = wy2[thePivot]

        theValue = abs( y1 + y2 + sdy )

#---- Binary Search

        while thePivot > 0:
            thePivot -= 1

            uy1 = y1 + wy1[thePivot]
            uy2 = y2 + wy2[thePivot]

            theUpperValue = abs( uy1 + uy2 + sdy )

            ly1 = y1 - wy1[thePivot]
            ly2 = y2 - wy2[thePivot]

            theLowerValue = abs( ly1 + ly2 + sdy )

            if theUpperValue < theLowerValue:
               if theUpperValue < theValue:
                  x1 += wx1[thePivot]
                  x2 += wx2[thePivot]
                  y1  = uy1
                  y2  = uy2

                  theValue = theUpperValue

            else:
               if theLowerValue < theValue:
                  x1 -= wx1[thePivot]
                  x2 -= wx2[thePivot]
                  y1  = ly1
                  y2  = ly2

                  theValue = theLowerValue

#---- Apply the Rotation

        x1 += dx1
        y1 += dy1

        x2 += dx2
        y2 += dy2

#---- Bounce Points Below the Axis to Above

        if y1 < 0:
           y1 = -y1

        if y2 < 0:
           y2 = -y2

#---- Comparison Determination

        c = 2

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              c = argPresumedResult
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              c = -argPresumedResult
        else:
           if y1 > y2:
              c = argPresumedResult
           elif y1 < y2:
              c = -argPresumedResult
           else:
              c =  0

#---- Exit

        return c, x1, y1, x2, y2

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

# Returns  ( c, t )
#
#    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 -1, 1

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

#---- CORDIC Loop

        Q1pow1 = Q1 >> 1
        I1pow1 = I1 >> 1
        Q2pow1 = Q2 >> 1
        I2pow1 = I2 >> 1

        Q1pow2 = Q1 >> 3
        I1pow2 = I1 >> 3
        Q2pow2 = Q2 >> 3
        I2pow2 = I2 >> 3

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

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

            if I1 <= I2 - I2pow2:
               return -1, 1 + n

            if I2 <= I1 - I1pow2:
               return 1, 1 + n

            Q1pow1 >>= 1
            I1pow1 >>= 1
            Q2pow1 >>= 1
            I2pow1 >>= 1

            Q1pow2 >>= 2
            I1pow2 >>= 2
            Q2pow2 >>= 2
            I2pow2 >>= 2

#---- Terminating Guess

        Q1pow1 <<= 1
        Q2pow1 <<= 1

        if I1 + Q1pow1 < I2 + Q2pow1:
           return -1, 6
        else:
           return 1, 6

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

You want to avoid multiplications.

For comparison purposes, not only do you not have to take the square roots, but you can also work with the absolute values.

Let

$$
\begin{aligned}
a_1 &= | I_1 | \\
b_1 &= | Q_1 | \\
a_2 &= | I_2 | \\
b_2 &= | Q_2 | \\
\end{aligned}
$$

Note that for $a,b \ge 0$:

$$ (a+b)^2 \ge a^2 + b^2  $$

Therefore
$$ a_1 > a_2 + b_2 $$
means that

$$ a_1^2 + b_1^2 \ge  a_1^2 > ( a_2 + b_2 )^2 \ge a_2^2 + b_2^2  $$

$$ a_1^2 + b_1^2 > a_2^2 + b_2^2  $$

This is true for $b_1$ as well.  Also in the other direction, which leads to this logic:

(The previous pseudo-code has been functionally replaced by the Python listing below.)

Depending on your distribution of values, this may save a lot.  However, if all the values are expected to be close, you are better off buckling down and evaluating the Else clause from the get go.  You can optimize slightly by not calculating s1 unless it is needed.

This is off the top of my head so I can't tell you if it is the best.

Depending on the range of values, a lookup table might also work, but the memory fetches might be more expensive than the calculations.

This should run more efficiently:

(The previous pseudo-code has been functionally replaced by the Python listing below.)

A little more logic:

$$
\begin{aligned}
( a_1^2 + b_1^2 ) - ( a_2^2 + b_2^2 ) &= ( a_1^2 - a_2^2 ) + ( b_1^2 - b_2^2 ) \\
&= (a_1-a_2)(a_1+a_2) + (b_1-b_2)(b_1+b_2) \\
\end{aligned}
$$

When $a_1 > a_2 $ ( and $a_1 \ge b_1 $ and $a_2 \ge b_2 $ as in the code):

$$ (a_1-a_2)(a_1+a_2) + (b_1-b_2)(b_1+b_2) >= (a1-a2)(b1+b2) + (b1-b2)(b1+b2) = [(a1+b1)-(a2+b2)](b1+b2) $$

So if $a_1+b_1 > a_2+b_2$ then

$$ ( a_1^2 + b_1^2 ) - ( a_2^2 + b_2^2 ) > 0 $$

Meaning 1 is bigger.

The reverse is true for $a_1 < a_2 $

The code has been modified.  This leaves the Needs Determining cases really small.  Still tinkering....

Giving up for tonight.  Notice that the comparison of $b$ values after the comparison of $a$ values are actually incorporated in the sum comparisons that follow.  I left them in the code as they save two sums.  So, you are gambling an if to maybe save an if and two sums.  Assembly language thinking.

I'm not seeing how to do it without a "multiply".  I put that in quotes because I am now trying to come up with some sort of partial multiplication scheme that only has to go far enough to make a determination.  It will be iterative for sure.  Perhaps CORDIC equivalent.

Okay, I think I got it mostly.

I'm going to show the $ a_1 > a_2 $ case.  The less than case works the same, only your conclusion is opposite.

Let

$$
\begin{aligned}
d_a &= a_1 - a_2 \\
s_a &= a_1 + a_2 \\
d_b &= b_2 - b_1 \\
s_b &= b_2 + b_1 \\
\end{aligned}
$$

All these values will be greater than zero in the "Needs Determining" case.

Observe:

$$
\begin{aligned}
D &= (a_1^2 + b_1^2) - (a_2^2 + b_2^2) \\
&= (a_1^2 - a_2^2) + ( b_1^2 - b_2^2) \\
&= (a_1 - a_2)(a_1 + a_2) + (b_1 - b_2)(b_1 + b_2) \\
&= (a_1 - a_2)(a_1 + a_2) - (b_2 - b_1)(b_1 + b_2) \\
&= d_a s_a - d_b s_b
\end{aligned}
$$

Now, if $D=0$ then 1 and 2 are equal.  If $D>0$ then 1 is bigger.  Otherwise, 2 is bigger.

Here is the "CORDIC" portion of the trick:

Swap db, sb   # d is now the larger quantity

While da < sa
  da =<< 1
  sb =<< 1
  If sb > db Then Swap db, sb # s is the smaller quantity
EndWhile

When this loop is complete, the following has is true:

$D$ has been multiplied by some power of 2, leaving the sign indication preserved.

$$ 2 s_a > d_a \ge s_a > d_a / 2 $$

$$ 2 s_b > d_b \ge s_b > d_b / 2 $$

In words, the $d$ will be larger than the $s$, and they will be within a factor of two of each other.

Since we are working with integers, the next step requires that both $d_b$ and $s_b$ be even or odd.

If ( (db+sb) & 1 ) > 0 Then
   da =<< 1
   sa =<< 1
   db =<< 1
   sb =<< 1
EndIf

This will multiply the $D$ value by 4, but again, the sign indication is preserved.

Let
$$
\begin{aligned}
p &= (d_b + s_b) >> 1 \\
r &= (d_b - s_b) >> 1 \\
\end{aligned}
$$

A little thinking shows:

$$ 0 \le r < p/3 $$

The $p/3$ would be if $ d_b = 2 s_b $.

Let
$$
\begin{aligned}
g &= d_a - p \\
h &= p - s_a \\
\end{aligned}
$$

Plug these in to the $D$ equation that may have been doubled a few times.

$$
\begin{aligned}
D 2^k &= (p+g)(p-h) - (p+r)(p-r) \\
&= [p^2 + (g-h)p - gh] - [p^2-r^2] \\
&= (g-h)p + [r^2- gh] \\
\end{aligned}
$$

If $g=h$ then it is a simple determination:  If $r=g$ they are equal.  If $r>g$ then 1 is bigger, otherwise 2 is bigger.

Let
$$
\begin{aligned}
v_g &= g - r \\
v_h &= h - r \\
\end{aligned}
$$

Evaluate the two terms on the RHS of the $D2^k$ equation.

$$
\begin{aligned}
r^2 - gh &= r^2 - (r+v_g)(r+v_h) \\
&= -v_g v_h - ( v_g + v_h ) r \\
\end{aligned}
$$

and

$$ g - h = v_g - v_h $$

Plug them in.

$$
\begin{aligned}
D 2^k &= (g-h)p + [r^2- gh] \\
&= (v_g - v_h)p - v_g v_h - ( v_g + v_h ) r \\
&= v_g(p-r) - v_h(p+r) - v_g v_h \\
&= v_g s_b - v_h d_b - \left( \frac{v_h v_g}{2} + \frac{v_h v_g}{2} \right) \\
&= v_g(s_b-\frac{v_h}{2}) - v_h(d_b+\frac{v_g}{2}) \\
\end{aligned}
$$

Multiply both sides by 2 to get rid of the fraction.

$$
\begin{aligned}
D 2^{k+1} &= v_g(2s_b-v_h) - v_h(2d_b+v_g) \\
\end{aligned}
$$

If either $v_g$ or $v_h$ is zero, the sign determination of D becomes trivial.

Likewise, if $v_g$ and $v_h$ have opposite signs the sign determination of D is also trivial.

Still working on the last sliver......

So very close.

theHandledPercent 98.6582716049

theCorrectPercent 100.0

Will continue later.......Anybody is welcome to find the correct handling logic for the same sign case.

Another day, another big step.

The original sign determining equation can be factored like this:

$$
\begin{aligned}
D &= d_a s_a - d_b s_b \\
&= \left( \sqrt{d_a s_a} - \sqrt{d_b s_b} \right)\left( \sqrt{d_a s_a} + \sqrt{d_b s_b} \right) \\
\end{aligned}
$$

The sum factor will always be positive, so it doesn't influence the sign of D.  The difference factor is the difference of the two geometric means.

A geometric mean can be approximated by the arithmetic mean.  This is the working principle behind the "alpha max plus beta min algorithm".  The arithmetic mean is also the upper bound of the geometric mean.

Because the $s$ values are bounded below by $d/2$, a rough lower bound can be established for the geometric mean without much calculation.

$$
\begin{aligned}
s &\ge \frac{d}{2} \\
ds &\ge \frac{d^2}{2} \\
\sqrt{ds} &\ge \frac{d}{\sqrt{2}} \\
&= \frac{\frac{d}{\sqrt{2}}}{(d+s)/2} \cdot \frac{d+s}{2} \\
&= \frac{\sqrt{2}}{1+s/d} \cdot \frac{d+s}{2} \\
&\ge \frac{\sqrt{2}}{1+1/2} \cdot \frac{d+s}{2} \\
&= \frac{2}{3} \sqrt{2} \cdot \frac{d+s}{2} \\
&\approx  0.9428 \cdot \frac{d+s}{2} \\
&>  \frac{15}{16} \cdot \frac{d+s}{2} \\
\end{aligned}
$$
If the arithmetic mean of a is greater than b's, then if the upper bound of the geometric mean of b is less than the lower bound of the geometric mean of a it means b must be smaller than a.  And vice versa for a.

This takes care of a lot of the previously unhandled cases.  The results are now:

theHandledPercent 99.52

theCorrectPercent 100.0

The source code above has been updated.

Those that remain unhandled are "too close to call".  They will likely require a lookup table to resolve.  Stay tuned.....

Hey Dan,

Well, I would shorten it, but none of it is superfluous.  Except maybe the first part, but that is what got the ball rolling.  So, a top posted summary would be nearly as long.  I do intend to write a blog article instead.  This has been a fascinating exercise and much deeper than I initially thought.

I did trim my side note about Olli's slope bound.

You should really be studying the code to understand how few operations actually have to be done.  The math in the narrative is simply justification for the operations.

The true "winner" should be the algorithm that is most efficient.  A true test would be both approaches programmed on the same platform and tested there.  As it is right now, I can tell you that mine (coded in C) will leave his in the dust simply due to I am prototyping with integers and he is using floats with a lot of expensive operations.

My thoughts at this point are that the remaining 0.5% cases I'm not handling are best approached with a CORDIC iterative approach.  I am going to try to implement a variation of Olli's approach, including his ingenius varying slope, in integers.  The "too close to call" category should be very close indeed.

I agree with you that Olli does excellent work.  I've learned a lot from him.

Finally, at the end, aren't we all winners?

Dan,

Your faith in the CORDIC is inspiring.  I have implemented a lossless CORDIC different than Olli's, yet might be equivalent.  So far, I have not found your assertion that it is the ultimate solution true.  I am not going to post the code yet because there ought to be one more test that cinches it.

I've changed the testing a little bit to be more comparable to Olli.  I am limiting the region to a quarter circle (equivalent to a full circle) and scoring things differently.

Return       Meaning
 Code
  -1     First Smaller For Sure
   0     Equal For Sure
   1     First Larger For Sure
   2     Presumed Zero

The last category could also be called "CORDIC Inconclusive".  I recommend for Olli to count that the same.

Here are my current results:

   Count: 538756

    Sure: 99.9161030225
 Correct: 100.0

Presumed: 0.0838969774815
    Zero: 87.610619469

  Faulty: 0.0103943157942

    High: 1.00950118765
     Low: 0.990588235294

Out of all the cases 99.92% were determined for sure and all the determinations were correct.

Out of the 0.08% cases that where presumed zero, 87.6% actually were.  This means that only 0.01% of the answers were faulty, that is presumed zero erroneously.  For those that were the quotient (I1^2+Q1^2)/(I2^2+Q2^2) was calculated.  The high and low values are shown.  Take the square root to get the actual ratio of magnitudes.

Roughly 83% of cases are caught by the primary determination and don't need any further processing.  That saves a lot of work.  The CORDIC is needed in about 0.7% of the cases.  (Was 0.5% in the previous testing.)

***********************************************************
*                                                         *
*   C O M P L E T E   A N D   U T T E R   S U C C E S S   *
*                                                         *
*   H A S   B E E N   A C H I E V E D  !!!!!!!!!!!        *
*                                                         *
***********************************************************

   Count: 8300161

    Sure: 100.0
 Correct: 100.0

Presumed: 0.0
    Zero: -1

  Faulty: 0.0

    High: 1.0
     Low: 1.0

You can't do better than that and I am pretty sure you can't do it any faster.  Or not by much anyway.  I have changed all the "X <<= 1" statements to "X += X" because this is way faster on an 8088.  (Sly grin).

The code will stay in this answer and has been updated.

Further explanations are forthcoming in my other answer.  This one is long enough as it is and I can't end it on a nicer note than this.


What is 'antenna-chain permutation' and its advantages? I am working with Wi-Fi channel state information (CSI, a sort of channel matrix with OFDM subcarriers). I am using Intel 5300 Wi-Fi card and this tool provides a way to play around with extracted CSI information via MATLAB.

Interesting thing is, there is perm, a permutation vector in a figure below:

It says a permutation vector is [2, 1, 3]. In FAQ page of the tool, it says:

  perm tells us how the NIC permuted the signals from the 3 receive antennas into the 3 RF chains that process the measurements. The sample value of [3 2 1] implies that Antenna C was sent to RF Chain A, Antenna B to Chain B, and Antenna A to Chain C. This operation is performed by an antenna selection module in the NIC and generally corresponds to ordering the antennas in decreasing order of RSSI.

What I am curious about is, why does this Intel 5300 Wi-Fi card do antenna-chain permutation and what are advantages of performing it?

I'll highlight terms that are worth googling in cursive below:

What you're looking at is a diversity receiver.

In this case, the idea is that by combining the signals of multiple antennas, the SNR can be maximized, minimizing the BER.

There's multiple methods of doing so – simply selecting the antenna with the highest RSSI (being probably the one with the best SNR) is called selection combining. When thinking about outage probabilities, this leads to a Bernoulli distribution, ie. the more antennas you have, the less like it is that your communication fails totally (but the gain of having more antennas quickly decreases after the second).

Another method is just adding up the signals – you can do that, and you'd get what's called equal gain combining. You'd get some sort of "oversampling" gain, supressing noise that is independent in each receiver branch, but a strong interferer might distort the result.

Then, you can add the receive signals up, but weigh them by some function that is proportional to their "quality". That way, you get maximum ratio combining, at least if the weights are basically proportional to SNR.

In a real-world receiver with a limited amount of ressources dedicated to assessing the signals, a simple ranking in three tiers of signal quality does sound like a feasible approach. For example, spending much power on synchronizing the strongest signal (i.e. estimating the phase of the channel impulse response) makes sense, since an error here has the gravest effect, whereas simpler approaches for the signals that are used to "push down" uncorrelated noise power (maybe aided by the strongest' chains info) might still make sense. Notice that modern WiFi receivers are pretty mighty beasts – not what you'd need to feel bad about if you don't understand them at once.

Another explanation might be much more in the analog than in the digital domain: For the strongest signal, you might want to use a different receive amplifier than for the weakest one – amplifiers are either a) low-noise, or b) high dynamic range or c) power-efficient (not to mention d) through z) cheap).


Extract jump location in data I would like to extract the location of data jump (see figure below), the difficulty is that the data are noisy as you can see. I try with a median filter or a other regularization and it could move the location of my jump.

Perhaps i could compute differences between consecutive data and then extract the most important one. Nonetheless my data could contain two jumps... And for information I m working with Python

Use a low-pass filter with linear phase. Since group delay is the derivative of phase over frequency, linear phase filters have constant group delay.

That way, yes, you shift that jump, but simply by subtracting the group delay, you get the original position.

Group delay for a discrete causal filter is half its length.


Is it possible to use k-means algorithm with just one vector Suppose I have a vector $X = (x_1, x_2 , . . . ,x_n)$, $x_i$ is the maximum of $X$ and $x_k$ is the minimum. Is it possible to use k-means algorithm to cluster the values in vector $X$ into two cluster. First, the values which are near to maximum value and the other cluster includes the values which are near to minimum value?
If not, Is there another algorithm can do that task?

Thank  you

Certainly possible. You can think of it in the way you described but I would say it just becomes an adaptive thresholding assuming a bimodal distribution: We admit that there are two main modes in the data and we try to identify those by grouping the values around the means of the clusters. These groups form clusters. The means are positioned such that each element contributes to the average of the cluster it belongs and thus creating a mode nearby.
It is almost sure that the minimum and maximum will go to separate clusters.


Fourier Transform of a PSD and response of a PSD input Does a Fourier Transform of a white noise exist?
If so, what is its general form?

It is possible to compute the fourier transform of a Power Spectral Density?

My problem in detail. I have to compute the response of a dynamic system given the PSD of an input ($ \Phi_{vv} $).

I thought to solve the problem in three step:

Build a shape filter (system which, given a white noise input ($ n $), gives a colored noise as output)
Express both the shape filter and the system dynamic with a fourier transform
If i can manage to express also the input (the white noise) in the fourier domain I should have solved the problem

The Fourier Transform of stationary white noise exists: in fact "white" here refers to the characteristics of the Fourier Transform, in that the power spectral density will be constant over the entire spectrum, and also implies that each sample in the time domain is independent of all other samples.
The distribution of the magnitude of the samples in time is also usually defined, such as being "Gaussian" or "Uniform" depending on what distribution is actually used. For white noise samples the phase is also distributed uniformly  $\in [0, 2\pi)$. Due to the central limit theorem, the Fourier Transform of any distribution in time that is white will tend toward Gaussian in frequency. The resulting Fourier Transform will be Gaussian distributed in magnitude, uniformly distributed in phase, and each sample in the frequency domain is independent of all other samples.
The Power Spectral Density is computed from the Fourier Transform as the complex conjugate product, and not the other way around. The reason is once computed all phase information is lost so there is a many to one relationship in the computation.


Autocorrelation function between received and training signal used for channel estimation I have read a paper that performs Channel Estimation of Wireless Channel as follows.

A training sequence of length $N$, lets call it $a_N(i)$ for $i\in[0:N-1]$  is repeated twice then sent out over the channel.

Assume the transmit sequence formed of this two sequeses is denoted by $s_{CE}(t)$.

The received signal is given by (assuming channel of length $T_{CH}$) is
$$r_{CE}(t) = \sum_{t'=0}^{T_{CH}-1}h(t')s_{CE}(t-t')+n(t)$$

The authors claim that if one takes the autocorrelation between a sequence $a_N(i)$ and the received signal then we can mathematically write the autorcorrelation as

$$R(t)= \frac{1}{N}\sum_{d=0}^{N-1}r_{CE}(t+d\underbrace{-N+1}_{????}) a_{N}^*(d) $$

My question is why is there a need for the term I have underbrace $-N+1$. I thought that the autocorrelation is in general expressed as

$$R(t)= \frac{1}{N}\sum_{d=0}^{N-1}r_{CE}(t+d) a_{N}^*(d) $$

Thanks looking forward for your view!

Update: The following is the reference

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=4699819&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D4699819

I particular (25) and (26) are my main concerns.

I think it's just a trick which makes $R(t)=0, \mathrm{when}\quad t<0$.

First, $-N+1$ is just a time offset and not affect the essence of function.

$r_{CE}(t)$ is bounded to $[0, 2N+T_{CH}-2]$, and $a_N(t)$ to $[0, N-1]$. In the form you gave, $R(t)$ is bounded to $[-N+1, 2N+T_{CH}-2]$. May be the author want the non-zero value of $R(t)$ also starts at 0, so he just move it right with $N-1$ length, which causes the $-N+1$ term in the equation.

--------- Why $R(t)$ is bounded as that -------------

This is just a explanation to give some intuition, not a formal proof. Consider two function $f(t)$ and $g(t)$, which $f(t)=0, t \notin [a, b]$ and $g(t)=0, t\notin [c, d]$. Now we want to calculate the sum of products, i.e., $S=\sum_t f(t)g(t)$, $[a,b]$ and $[c,d]$ must be overlapped, right? Thus we get $a \le d $ and $c\le b$.

Get back to the $R(t)=\sum _d r_{CE}(t+d)a^*(d)$(note here the argument is $d$), $a^*(d)$ is bounded to $[0, N-1]$ and $r_{CE}(t+d)$ to $[-t, L-1-t]$, in which $L=2N+T_{CH}-1$ is the length of $r_{CE}(d)$. Here we get $a=0, b=N-1, c=-t, d=L-1-t$. Using $a \le d $ and $c \le b$, we get $-(N-1) \le t \le L-1$.

As to why $r_{CE}$ is bounded as that, you can make a same process to understand it.


What is the Roll-Off factor of a Root Raised Cosine Filter? When designing a root raised cosine filter on Matlab, there is an option to set its "roll-off factor". What does this mean? I have been looking around for a while but I can't seem to find a clear explanation. Any help will be appreciated.

Think of the root-raised cosine as a family of functions. The roll-off factor is a parameter that chooses one of the possible functions.

All root-raised cosine pulses look similar to the sinc pulse. The interesting thing about the roll-off factor is that it controls two features of the function:

The rate at which the function's lobes (or ripples) decrease. At one extreme, with roll-off set to 0, they decrease slowly (as $1/t$, with $t$ equal to time). At the other end, when the roll-off factor is 1, they decrease as $1/t^3$.
As the roll-off factor increases, so does the function's bandwith. Basically, if the function's bandwidth is $B$ when the roll-off is 0, then it becomes $2B$ when the roll-off is 1.

You can choose any value for the roll-off factor between 0 and 1, so you can select the function with the ripple and bandwidth that are best suited to your needs.


Are location in BRIEF feature descriptor reused? In the BRIEF feature descriptor ("BRIEF: Binary Robust IndependentElementary Features", Calonder et al.,  2011) they investigate five different methods of choosing point pairs to compute the binary values of the descriptor vectors:

Among these methods there are the methods named GI up to GIV that sample the point pairs randomly, according to various distributions.
Are these point pair sampled for every keypoint that we want to compute the descriptor for separately, or are these pairs sampled once and then remain the same for every keypoint?

When using a randomized pattern in BRIEF, this means that you computed random positions inside the patch once in an offline procedure, then used these random locations every time you computed the descriptors. This makes sense, as it means that when comparing descriptors you will actually compare the same locations, it's simply that the sampling pattern was not obtained by some regular or geometric scheme (which is probably what makes it a bit more generic than hand-crafted patterns).
For example, here is the sampling pattern used by OpenCV.


In a JPEG file, how much is noise and how much is the image? I have an IP webcam that generates very low light images of a plain background.  The camera spits out JPEG only files.  The following is one such frame...

and this is with graphic equalisation applied...

You can see that there is a lot of noise, but also a discernable image of the stepped lens mount.  These JPEG files are approximately 20KB in size.

I want to determine what proportion of the 20KB is due to the fixed image (of the lens mount), and what proportion is entirely due to noise.  It's important to reiterate that I'm only concerned with the proportions of the JPEG file, and not the expanded raster images.

So noise.jpg(part) + mount.jpg(part) = 20KB.

My instinct suggested repetitively sampling the static image and frame stacking to isolate the lens mount.  That fundamentally doesn't work.

How can I determine these proportions?

Your problem here is twofold:

Image Noise / low  SNR due to darkness
problems to average things out due to JPEG compression

The first point is something we cannot really change per se; your averaging approach surely is a good one.

The second problem now is due to how the JPEG compression works in your camera:

Gamma correction is applied to the RGB image first. That is a nonlinear operation, "pushing together" high values and stretching the lower values. That way, you don't get "pitch black" everywhere for dark images. Sadly, this also makes our quantization noise problem harder.
The resulting R'G'B' pixel data, is converted to the Y'CbCr colorspace, meaning you convert the threedimensional red,green,blue signal to the threedimensional brightness,Blue-difference,Red-difference signal.
The Cb/r channels are low-pass filtered and downsampled
The original-resolution Y', and the reduced-resolution Cb and Cr images are split into 8x8 pixel blocks separately
These blocks are separately DCT'ed
A quantization matrix is chosen according to pick a quality/compression trade-off. This results in different DCT coefficients being omitted, saved with fewer or more bits.

As shown in the other answer, averaging on the reconstructed RGB image is a bad idea, since the random noise must lead to you basically seeing an effect of the IDCT of the quantization matrix in each block – meh.

So, first of all, notice that it's pretty certain that the for the DCT coefficient that defines the constant intensity offset ("DC component", an EE would say), the highest quantization level was chosen.

Conversely, at least for most  higher-frequency components in the 8x8 blocks, quantization error energy is higher.

Thus, inherently, the DCT's constant bin is the most reliable – also, it happens to be the average of the 64 pixels in the block, and thus reduce noise already (if noise's frequency distribution was white, that would be the case for all components, but neither is it reliably white in cameras nor is the gamma correction linear...).

So, what you probably want to do as a first step is

work on the  Y'CbCr, not the RGB reconstruction
only look at the Y component first
ignore everything that is not the DC component in each 8x8 block

You'd get an image of $\frac18$ the width and height of your original photo. It's probably still going to be noisy – with the shot noise that is probably relevant in a low-light setting contributing to the average in each frame. However, averaging out that over time will reduce that effect.

You can probably get an noise reduction, too, by summing up the non-zero-frequency coefficient, weigh them with a factor < 1, and subtracting them from the DC component. The idea behind that is that when you set a single pixel in a 8x8 frame, you  get energy in all coefficients (the cosine transform of a dirac has energy at all frequencies, if you think about it). Since you do not get all coefficients, most of them being thrown away or quantized heavily, you'll never be able to fully remove the effect of shot noise from the average coefficient – but you can, at least on average, reduce the effect. Choice of that factor would depend on camera sensor, gamma factor, quantization matrix, the way the DCT is implemented (truncating fixed point?) and thus, too many unknowns, but you could probably just empirically try a couple of values and minimize the variance of the result over multiple frames.


What is an Adaptive Mean Filter? $$\hat f(x, y) = g(x, y)-\frac{\sigma_n^2}{\sigma_L^2}\left[g(x, y)-m_L\right]$$

What are the meanings of the following terms:

$m_L$
${\sigma_\eta}^2$
${\sigma_L}^2$

Here we see that $m_L$ is subtracted from the Image and then the whole term is multiplied by $\frac{{\sigma_\eta}^2}{{\sigma_L}^2} $. Then, the whole term is again subtracted from $g(x,y)$.

What does that actually mean?

A wild guess. $L$ denotes a region $R$ around pixel $g(x,y)$, potentially a symmetric one of size $(2L+1)\times (2L+1)$. $m_L$ and  $\sigma_L^2$ are the average and variance in $R$. And $\sigma_n^2$? Probably the estimated noise variance over the whole image, supposed stationary.

In a flat region with noise $\sigma_L^2 \sim \sigma_n^2$, the result could be close to the average. Around edges, one expect an higher response, so sharper. Potential problems arise when  $\sigma_L \ll \sigma_n $. Another interpretation:

$$\hat{f}(x,y)=(1-\frac{\sigma_n^2}{\sigma_L^2})g(x,y)+ \frac{\sigma_n^2}{\sigma_L^2}m_L(x,y)$$
is a sort of weighted mean (since weights sum to one) of the image and its smoothed version, whose dependency on $(x,y)$ is made more explicit.


What are $H(i\omega)$, $H(s)$ and $H(z)$ called? What are the Fourier-, Laplace- and Z-transforms of the impulse response of a filter called? I've seen $H(s)$ referred to as the system function or transfer function, what about $H(i\omega)$ and the discretized $H(z)$?

I would consider them to all be 'transfer' or 'system' functions.
Those two names are interchangable, but in control theory books I see then listed more as transfer functions rather than system functions.

The domain you are in, whether it be fourier, laplace, Z or time, isn't entirely relevant to the name in that respect.


How to test the FIR filter in C++? I am almost done with design and implementation of my first FIR filter using C++. Its a Sinc-Windowed FIR filter designed for a cut off frequency of 100 Hz with sampling rate of 500, with a filter order M = 100. All the operations have been carried out in discrete time domain. So basically my question is regarding testing to see if my filter is generating the desired response .

I presume one way to test the output would be to apply DFT to it and check the response spectrum of the filter in the frequency domain, but for this I have to develop a function that performs the DFT on the output .

Question :To avoid this, I am looking,

---->Are there any methods that make it possible to test the filter response in the (without the need for applying DFT) time domain itself?.

----> If the above is not possible, are there are any alternative methods like some C++ DSP library with FT function's you can suggest?

Or Any suggestions and advice regarding this with your experience would be appreciated. Thanks in advance!

The first basic test could be to use a unit impulse as an input signal and see if the output signal equals the impulse response (i.e. the filter coefficients). Another simple test signal is a unit step. The corresponding output should be the cumulative sum of the filter's impulse response, i.e. for $x[n]=u[n]$, the output must be

$$y[n]=\sum_{k=0}^nh[k],\quad n\ge 0$$

If you obtain the desired response for these two simple test signals, that's a good start, but I wouldn't stop there.
In general it is very helpful to use a tool like Matlab or Octave (free), or anything similar, which has built-in filtering routines against which you can check your routine for arbitrary input signals.


Accuracy of dense optical flow Currently I am learning dense optical flow by myself. To understand it, I conduct one experiment. I produce one image using Matlab. One box with a given grays value is placed under one uniform background and the box is translated two pixels in x and y directions in another image. The two images are input into the implementation of the algorithm called TV-L1. The generated motion vector outer of the box is not zero.

Is the reason that the gradient outer of the box is zero?
Is the values filled in from the values with large gradient value?

In Horn and Schunck's paper, it reads

  In parts of the image where the brightness gradient is zero, the velocity estimates will simply be averages of the neighboring velocity estimates. There is no local information to constrain the apparent velocity of motion of the brightness pattern in these areas.

and

  The progress of this filling-in phenomena is similar to the propagation effects in the solution of the heat equation for a uniform flat plate, where the time rate of change of temperature is proportional to the Laplacian.

Is it not possible to obtain correct motion vectors for pixels with small gradients? Or the experiment is not practical. In practical applications, this doesn't happen.

It's normal and due to your test setup.

Let's try to understand a bit. The algorithm is trying to minimize an objective function defined as:

$$ F(u, I_0, I) = \int \| I - I_0(u(x)) \|_1  + \lambda \| \nabla u \|_ dx, $$

where $I_0$ and $I$ are the images, $x$ is an image point (ie a pixel) and $I_0(u)$ a shortcut for "$I_0$ as remapped using the optical flow $u$". The first term is the data term (because it binds the result $u$ to the observed images) and the second is the regularization constraint on the shape of $u$ (the total variation TV in this case).

Let's ignore the regularization term for now. The registration in the data term can be obviously minimized by an optical flow field $u$ which is null everywhere, except for the pixels that belong to the edges of your box. For these pixels, the best solution is a translation that moves them to the nearest box edge in the new image.

So far, so good? Well, not anymore if we consider now the regularization term. The TV constraint will penalize solutions where neighbouring pixels have different optical flow values. This means that our current solution is going to have a high score (bad!) because the optical flow will vary around the box edge pixels. A much better solution with respect to TV is to get a constannt flow field ($\nabla u$ becomes null, and  the regularization term evaluates to 0).

Joining both constraints (registration and regularization), we see that we need a solution flow field that:

brings the box edge pixels to the new box position
is as-constant-as-possible everywhere.

A solution fulfilling both constraints is a constant flow corresponding to the box translation: inside and outside of the box the values are constant so any flow vector works, and for box edges this is the best optical flow vector. In that case, the objective function does even evaluate to 0 (no registration error, TV is 0 too).

How to avoid this behaviour?

You can avoid this behaviour by using some random noise texture to fill your background and box. In that case, the registration constraint will force the flow to be near 0 on the background (the actual solution will depend on the $\lambda$ parameter). And obviously, this won't happen in "real life" applications because of sensor noise, small scene movements, textures...

In practice, you should generate 2 images filled with random noise, for example Gaussian noise, and preferably using different statistics (mean and variance):

the first image will have the size of your test image. Let's call it the background. It is the texture of the immovable part of your test image;
the second image should have the size of your box. Let's call it the object. It is the texture of your moving target.

Note that both images shouldn't be modified after creation: the textures are random, but they are not dynamic. If they are dynamic then there would be additional reasons (other than object motion) for their appearance to change.

You can create your test images by pasting the object layer on top of the background at the desired position. if you test only with integer displacements, you don't even have to take care of interpolation of any kind.


Output of a stable LTI system Let $\mathcal{L}$ be a stable LTI system. Is it true that if input has finite energy then output also has finite energy? I'm not sure about that. We know that $$\int_{-\infty}^{+\infty}|h(t)|dt\lt\infty \tag{1}$$Where $h(t)$ is the impulse response. Also we have $$\int_{-\infty}^{+\infty}|y(t)|dt = \int_{-\infty}^{+\infty}|Y(s)|ds = \int_{-\infty}^{+\infty}|X(s)||H(s)|ds \tag{2}$$Since $y(t) = \mathcal{L}(x(t)) = x(t)\star h(t)$ which implies $Y(s) = X(s)H(s)$. Applying Cauchy–Schwarz inequality to $(2)$, $$\left(\int_{-\infty}^{+\infty}|X(s)||H(s)|ds\right)^2 \le  \left(\int_{-\infty}^{+\infty}|X(s)|^2ds\right)\left(\int_{-\infty}^{+\infty}|H(s)|^2ds\right) \tag{3}$$We know that $$\int_{-\infty}^{+\infty}|X(s)|^2ds = \int_{-\infty}^{+\infty}|x(t)|^2dt<\infty$$Since input is an energy signal but $$\int_{-\infty}^{+\infty}|H(s)|^2ds$$doesn't necessarily exists. So is this indicates that we can find a counterexample to the statement or we can prove that by other methods?
Edit: Here is a counterexample which shows $$\int_{-\infty}^{\infty}|h(t)|dt\lt \infty \nRightarrow \int_{-\infty}^{\infty}|h(t)|^2dt\lt \infty$$

I think I've found the answer. Please correct me if I'm wrong. First of all, I've made a silly mistake  $$\int_{-\infty}^{+\infty}|y(t)|dt = \int_{-\infty}^{+\infty}|Y(s)|ds$$which is clearly false. Let $y(t) = x(t)\star h(t)$. We have $$E_y = \int_{-\infty}^{+\infty}|y(t)|^2dt = \int_{-\infty}^{+\infty}|Y(s)|^2ds = \int_{-\infty}^{+\infty}|H(s)X(s)|^2ds = \int_{-\infty}^{+\infty}|H(s)|^2|X(s)|^2ds$$Also we have $$|H(s)| = \left|\int_{-\infty}^{+\infty}e^{-2\pi ist}h(t)dt \right | \le \int_{-\infty}^{+\infty}|e^{-2\pi ist}h(t)|dt =  \int_{-\infty}^{+\infty}|h(t)|dt \lt\infty$$So  $\exists M\in\mathbb{R}:\ \ |H(s)|\le M$ for all $s$. This means that $|H(s)|^2\le M^2$ and then $$\int_{-\infty}^{+\infty}|H(s)|^2|X(s)|^2ds\le M^2\int_{-\infty}^{+\infty}|X(s)|^2ds$$By assumption $$E_x = \int_{-\infty}^{+\infty}|x(t)|^2dt=\int_{-\infty}^{+\infty}|X(s)|^2ds$$The result is $$E_y \le M^2 E_g$$


Half Window Functions and Formulas I am looking at implementing different window functions in an application I am making.  I have found the formulas for all major symmetrical windows, however I would like to implement half windows.

This may be a stupid question, but does the window literally jump straight from 0 to 1? or is there a small curvature to the start. ie, do I use the symmetrical formula and just use the second half? Or are there specific formulas for calculating half windows?

A use of such window functions would be to window data that "goes straight to business" starting at the first sample and then gradually dies off, like an impulse response of some system. See What is the meaning of half window functions? A smooth rise in the window function would be detrimental as it would attenuate the important peak, unless the rise in the window is quick enough and the peak delayed enough. However, I think it would be a stretch to call something like that a half-window, if it is not a half of some window. So start directly at 1 at the peak of a symmetrical window I think.

Then again, beware that a sine window a.k.a. cosine window may also be called half-sine window or half-cosine window, as it represents half a cycle of a sinusoid (sine or cosine).


ROC of Z Transform of $x(n) = 2(3)^nu(-n)$ Using definition, I got its Z transform as $X(z) = \dfrac{2}{1-\dfrac{z}{3}}$ and the summation converges only when $|z|<\frac{1}{3}$. So its ROC is $|z|<\frac{1}{3}$.

But my question is: for such a left sided signal $x(n)$, its ROC should be inner to the innermost pole. But the pole of $\dfrac{2}{1-\dfrac{z}{3}}$ is $z=3$ and not $z= \frac{1}{3}$. So the ROC I found is incorrect, right?

Using the definition, and letting $k=-n$:

$X(z)=2 \sum_{n=-\infty}^03^nz^{-n}=2\sum_{k=0}^{\infty}\big(\frac{1}{3}z\big)^k$

Now we can use the geometric series (https://en.wikipedia.org/wiki/Geometric_series) formula only if we require $|\frac{1}{3}z|<1$, or equivalently $|z|<3$


Calculating poles with complex numbers and quadratic equation given: $z^2 + 0.8 \sqrt2 z + 0.64 = 0 $

Then, I am using the quadratic equation:

$ z_{1,2} = \frac{-0.8 \sqrt2 \pm \sqrt{(0.8\sqrt2)^2-4 \cdot 0.64}}{2} $

Wolfram Alpha says it the end there should be: $-0.4 \sqrt2 \pm j 0.4 \sqrt2$, but I get something different:

$$ = \frac{-0.8\sqrt2\pm0.8\sqrt2-2\cdot0.8j}{2} $$
$$ = -0.4\sqrt2\pm0.4\sqrt2-1.6j$$

You seem to be computing, on the square root side, something like, $a$ and $b$ being positive, $\sqrt{a-b}=\sqrt{a}-j\sqrt{b}$, which is wrong. The square root is not linear.

Here, the discriminant $\Delta = (0.8\sqrt{2})^2-4\times 0.8^2 = -2\times 0.8^2$ is negative. Hence, the roots are:

$$ z_r= \frac{-0.8\sqrt{2}\pm j\sqrt{-\Delta}}{2} = 0.4\sqrt{2}(-1 \pm j)$$

agreeing with Wolfram Alpha.


Reconstructing a signal from FFT by adding individual signal components I'm attempting to reconstruct a signal from the DFT of the signal. I tried to do it by extracting the individual sinusoids and adding them up, but the answer I get is incorrect.

s = [3 1 4 1 5 9 2 6]#Some random signal(first few digits of pi in this case)
t = 0:(1/7):1;#Sampling at each of these equally spaced time intervals
f = fft(s);
a = abs(f);#Magnitude of each DFT component
an = arg(f);#Phase of each DFT component

recs = (a(1)/8)*cos(0*pi*t + arg(1)) + (a(2)/8)*cos(2*pi*t + arg(2)) + (a(3)/8)*cos(4*pi*t + arg(3)) + (a(4)/8)*cos(6*pi*t + arg(4)) + (a(1)/8)*cos(8*pi*t + arg(5));
#The above equation is derived(hopefully correctly) from the theory given in the question description
recs
s

This is a program I wrote in GNU Octave(should work on Matlab too). And the basis for using that particular expression for recs is based on some info I got from http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf.

I don't understand if my method for getting the individual components was wrong or is it some basic step I'm missing here.

There are a couple of mistakes in your calculations:

The t array is incorrectly specified. It should be 0:1/8:7/8 as the time positions at which the signal is assumed to be sampled for the default FFT calculation performed by fft() are 0, 1, ..., (N-1) where N = 8 in your case, the signal length. Since in your formula you are using t as the ratio n/N (where n is the array of sample times) you get the aforementioned expression to use for t.
You wrote a(1) instead of a(5) at the last term of your recs expression.
You are missing 2* multiplying each term except the first (k=0) and last (k=N/2) ones, where k is the index of the frequency axis of the Fourier transform F (recall that this 2* comes from leveraging the symmetry of the Fourier transform associated to a real signal).

Therefore, your reconstruction expression should read:

recs =  (1*a(1)/8)*cos(0*pi*t + an(1)) + ...
        (2*a(2)/8)*cos(2*pi*t + an(2)) + ...
        (2*a(3)/8)*cos(4*pi*t + an(3)) + ...
        (2*a(4)/8)*cos(6*pi*t + an(4)) + ...
        (1*a(5)/8)*cos(8*pi*t + an(5));

which gives the original signal s as verified by computing:

recs - s

ans =

   1.0e-14 *

         0   -0.1332         0   -0.0444   -0.2665    0.3553    0.1332         0


Sum of Diracs in Frequency Domain I am attempting to compute the sum of dirac functions in the frequency domain via the following formulas:

$$ ReLu(f(x)) = f(x) * \sum_{i} \delta(x-x_i) $$
$$ x_i : f(x_i) > 0 $$
$$ F(ReLu(f(x)))(k) = (\sum_{i} e^{2\pi jk x_i})F(f(x)) $$

From what I understand I need to find the values where x is larger than zero, then sum $e^{2\pi jk x_i}$ over these values. I made myself a trivial example to determine if I am doing this correctly and there is no difference to the image after the convolution by the sum of the exponentials. What am I interpreting wrong here?

Note this idea is from the paper: http://cs231n.stanford.edu/reports/2015/pdfs/tema8_final.pdf

import numpy as np
import matplotlib.pyplot as plt

# Make image array with a single negative value
img = np.array([[-1.0,2.3],[5,7.8]])

# Dirac values
N = 2
x = np.arange(0, N, 1)
y = np.arange(0, N, 1)
xm, ym = np.meshgrid(x, y)
x0 = xm
y0 = ym
shiftMat = np.exp(1j * ((2.0 * np.pi)) * (xm*x0 + ym*y0))

# Set equal to shift mat
# In this trival example I know that [0,0] is only negative position
# So set to 0 and compute sum of all positions in which f(x) > 0 as detailed in paper
freqRelu = shiftMat
freqRelu[0,0] = 0
freqRelu = np.sum(freqRelu)

# Fourier Convolution and IFFT
imgFFT = np.fft.fft2(img)
freqR = np.multiply(imgFFT,freqRelu)
# freqR = np.outer(freqRelu,imgFFT)
reluedFreq = np.real(np.fft.ifft2(freqR))

# Spatial Relu For Comparision
reluedImg = img
reluedImg[0,0] = 0

plt.subplot(121)
plt.imshow(reluedImg)
plt.subplot(122)
plt.imshow(reluedFreq)
plt.show()

# print(np.allclose(reluedFreq,reluedImg))
print(reluedFreq)
print(reluedImg)

ReLu in the frequency domain becomes convolution, NOT multiplication.
See chapter 4.5 of the paper: "...in frequency domain, ReLu() acts as a convolution with the function of a known form. "
But you don't need to calculate it at every point because of the pooling.


is any filter whose impulse response has a discontinuity in the first derivative fundamentally high-pass? a triangular (bartlett) window has a kink (discontinuity in first derivative).  thinking in the continuous (not discrete) case, my intuition is that the fourier transform of any window with a kink must asymptote to a non-zero value (rather than exponentially decay), since it should take infinitely high frequencies to be able to sum to a precisely sharp "corner."  treated as the kernel (impulse response) of a filter, this would correspond to passing arbitrarily high frequencies.

is this right?  in general, how do i calculate the asymptotic value, and why isn't this a standard quoted feature of a window?

This actually isn't the case. Consider the continuous-time triangle function:

$$
\text{tri}(t) = \begin{cases}
1 - |t|, && |t| < 1 \\
0, &&\text{otherwise}
\end{cases}
$$

This function has the Fourier transform:

$$
\text{Tri}(f) = \text{sinc}^2(f)
$$

This follows readily from the fact that the triangle function can be constructed by convolving a rectangle function, which has a Fourier transform of $\text{sinc}(f)$, with itself. The $\text{sinc}$ function decays to zero as its argument approaches infinity (and its square decays even faster yet), so the triangle window does not have a highpass characteristic.


How to interpret these different Fourier analysis of this audio signal? This is my first dive in DSP. I would like to familiarize myself with frequency analysis. I have two audio tracks which should be digitized at 16bit-44.1kHz and 24bit-192kHz (music, presented as a 24bit-192kHz sample) respectively.

I wanted to identify the effect of the low-pass filter around the Nyquist frequency (22.05kHz and 96kHz respectively).

Edit: I completely reworked the question.

Software used:

I basically estimated the power spectral density using Welch's method
as implemented by scipy.signal.welch in the Scipy library of the Python programming language.

Basically, I used a script equivalent to:

import numpy as np
from matplotlib import pyplot as plt
from scipy import signal
from waveio import readwav

# Load data from one channel (#0) for each sample file
wav192 = readwav("24b-192khz.wav")[:,0]
wav44 = readwav("16b-44khz.wav")[:,0]

# DoE: 2 sample size and two windows types
chunks = [256, 4096]
windows = ["hanning", "boxcar"] # boxcar is rectangular

# Prepare a figure
plt.figure()
# Calculate density spectra and plot
for N in chunks:
    for w in windows:
        f, Pxx44 = signal.welch(wav44, fs=44100, window=w, nperseg=N, nfft=2*N, scaling="density")
        plt.semilogy(f, Pxx44)
plt.legend(["chunk=%d; window=%s"%(c, w) for c in chunks for w in windows])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Density (I$^2$/Hz)")

The power spectral density of the 44.1kHz audio sample:

Which is basically just as expected:

The chunk size, i.e. the number of samples per fft-transform segment in the real domain, does not change the density a lot if a Hann window is used.
The chunk size effect is clearly visible with the boxcar (rectangular) window. From what I understand, this is because of spectral leakage which diminishes as the chunk size increases. Is that correct?
The low-pass filter effect at the Nyquist frequency (22.05kHz)

So far, so good.

The power spectral density of the 192kHz audio sample:

Good point:

Same behaviour in regard to the chunk size and the window. Is spectral leakage really that strong? That's pretty impressive.

Oddities:

What the heck is happening?
Where is the low-pass filter near the Nyquist frequency?
Why are very-high frequencies even increasing? Could that be related to the choice of the windowing function?

From my interpretation, there is no low-pass filter visible because basically no audio system would go above 192kHz and generally, the software/hardware creator are smart enough to apply a low-pass filter designed with regard to the actual output bandwidth of the audio system.

As for the increasing audio signal above 57kHz, I really can't explain it: the original audio sample is some classical music. I wouldn't expect any instrument to generate louder sounds in that range or frequencies. Any idea? Could this be an example of upsampling?

If you look at the Rectangular window the best its rejection gets is about 40 dB. So that behavior, especially obvious in the bottom plot, for the rectangular window is to be expected.

I don't know for sure if this explains everything but look at the level between your peak signal and the high-frequency components. There is almost 60 dB of rejection there. I've always heard that a good rule of thumb is to get 60 dB of rejection from your filters. I know that's not the full 96 dB offered by 16-bit music, but I bet you'd be hard pressed to actually hear that. Of course, there's the silliness of having music at that sample rate. Humans just can't hear anything above around 20 kHz, give or take. This article gives a good summary of the issues. It also brings up a good point, that that sample rate has the potential to pick up harmonics and other high-frequency distortion caused by equipment and electronics, despite the fact that we can't hear it. Perhaps something like that is going on?


I can't make sense of the amplitude spectrum of a signal. Help please I have a sinusoidal signal namely:
signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);

so the segment of the code becomes:
freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
% Let the max value of t be 4 times the largest period
t_max = 4 * periods(1);
% Let us have 500 samples over the interval [0, t_max].
t = linspace(0, t_max, 500);
% The synthesized signal
signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);

Now I wanna view the amplitude spectrum for the first 50 samples.
So I did this:
signal_first_50_samples= signal (1:50);
% Display its magnitude spectrum.
N = 64;
signal_spect = abs (fft(signal_first_50_samples,N));
signal_spect = fftshift(signal_spect);
F = [-N/2:N/2-1]/N;
plot (F, signal_spect)
set(gca, 'XTick',[-0.5:0.1:1])
grid on;

I was expecting spikes at each of the frequencies, however, I got this:

I re-did the sampling part as follows, to include more samples:
%% Synthesizing a sampled signal that consists of 4 sine waves
% We are given 4 different frequencies
freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
% Let the max value of t be 4 times the largest period
t_max = 4 * periods(1);
% Let us have 500 samples over the interval [0, t_max].
t = linspace(0, t_max, 500);
% The synthesized signal
signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);
% Visualising the synthesized signal
plot(t, signal, 'b-', 'LineWidth', 1.5);
grid on;
title('The Synthesized Signal', 'FontSize', fontSize);
xlabel('t', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
% Display its magnitude spectrum.
% i was expecting spikes at each freq
N = 64;
signal_spect = abs (fft(signal,N));
signal_spect = fftshift(signal_spect);
F = [-N/2:N/2-1]/N;
plot (F, signal_spect)
%set(gca, 'XTick',[-0.5:0.05:1])
grid on;

but i still got approx the same result, not what i expected

another edit:
i did something else i had only 50 samples for the singal and took all of them
%% 50 samples taking 50
clc;
clear;
% We are given 4 different frequencies
freqs = [0.08, 0.2, 0.32, 0.4];
periods = 1./ freqs;
% Let the max value of t be 4 times the largest period
t_max = 4 * periods(1);
% Let us have 50 samples over the interval [0, t_max].
t = linspace(0, t_max, 50);
% The synthesized signal
signal = sin(2*pi*0.08*t) + sin(2*pi*0.2*t) + sin(2*pi*0.32*t) + sin(2*pi*0.4*t);
% Determine the first 50 samples of this signal
% signal_first_50_samples= signal (1:50);
% Display its magnitude spectrum.
% ??????????????????????????????????????????????????????????????????????
% I was expecting spikes at each freq
N = 64;
signal_spect = abs (fft(signal,N));
signal_spect = fftshift(signal_spect);
F = [-N/2:N/2-1]/N;
figure
plot (F, signal_spect)
%set(gca, 'XTick',[-0.5:0.05:1])
grid on;

I got this

which looks about right
but I wonder why it didn't work when i had 500 samples and i took all 500 of them
any ideas?

A few suggestions and comments:

You calculate the frequency vector with F = [-N/2:N/2-1]/N;. This is wrong. You need to take the sampling frequency into account.
The way you organized your code, you don't have direct control over the sampling frequency. I recommend specifying both the time length and the sampling frequency explicitly in your code.
You need to set up your analysis so that you get high resolution in the frequency domain (so that you can resolve peaks that are very near to each other). The only way to achieve this is to sample for a long time. The more you oversample the signal (that is, the higher the Nyquist frequency is), the longer time you need to observe the signal. This is the reason you want to have explicit control over sampling rate and signal duration.
Do not limit the FFT size by choosing an explicit N, since that may lower your frequency resolution. If you want to set the FFT to a power of 2, then use nextpow2 or similar.


Using Matched Filter To Predict the Existence of An Event I have a stream of signal, and I'm looking for the presence of a certain event (the signal), also I have a "model"(which is a an approximation of the original signal) for that signal, and I will adopt it as my matched filter. So I tried to choose the matched filter, without any mathematical calculation. When I was looking for the signal. I saw a pattern that appears repeatedly. And I expect that pattern is the signal.

Could I cross-correlate the matched filter every unit of time (for example, every 2 seconds) with the signal and see if the value of cross correlation approaches to 1, and based on that decide if that event is existing or not?
To make things more clearly, the following is an algorithm of my approach

while (Signal s is Streaming)
{
  S = epoch_the_signal(s)
  x = cross_correlation(S,M)
  if(x > Threshold)
  {
    print "The signal appears"}
  }
}

What if I choose the matched filter as what is given in the literature (taking the matched filter as to increase the SNR. So I take what I expect to be the signal, and I used the formula, which is shown below to build the matched filter). does that will increase the noise that my model is including it in addition to the signal?

I think the hardest work is how to get the signal, not how to find the matched filter. Especially if it's biomedical signal.

Yes. However cross correlating is matching so you got a mix-up in terminology. You can detect the presence of the signal this way.
Matched filter will always maximize the output SNR so it will not increase your noise. Given Gaussian noise it is also a Most Powerful test (I.e optimal).


Digital filters and the Z-transform Why is it that we Z-transform a difference equation to get a the transfer function of an digital filter?

How come a digital filter is given in the Z-domain, and what is the Z-domain?

And for that sake, why do analog filters operate in the S-domain, and what is the S-domain?

The S-transform allows you to deal with differential equations in an algebraic manner - so they become easier to solve. Since continuous/analog filters consist of integrators and differentiators the S-transform is therefore a natural way to deal with these systems.

The z-transform provides an algebraic way of dealing with finite difference systems and therefore it is a natural way to deal to discrete-time systems i.e. digital.

In the S-transform, setting $s=j\omega$ results in the Fourier transform. For continuous time systems we are interested in whether poles are on the right or left side of the $j\omega$ axis, because that determines the stability.

For the Z-transform, setting $z=e^{j\omega}$ results in the Discrete-Time Fourier transform. To determine stability, we are interested in whether the poles are inside or outside the unit circle.

So the S and Z domains are similar - they allows you solve continuous and discrete time systems, respectively,  using algebra.


Is there an intensity estimation algorithm working solely in time domain Given a second of sampled sound pressure signal, speed of sound and medium density, can the average sound intensity be estimated without performing DFT, that is in linear time by simply traversing all samples?

Could you provide\link a proof that it is correctly evaluating\approximating this definition from wikipedia:

A lot of examples deal with pure tonal sine waves, where particle velocity member can be transformed to pressure, and amplitude is clearly defined, but I'm struggling to find an example for random signal.

Particle velocity and sound pressure are related by the field impedance. The field impedance is a function of the media and the boundary conditions of the field.

For a reasonably shaped wave in a reasonably well behaved medium, that impedance is simply $\rho \cdot c$, i.e. density of the medium times speed of sound. For air that is about 400 Pa/(m/s). There is no phase shift and not a lot of frequency dependency.

In this case the intensity signal can be calculated as

$$I(t) = \frac{p^2(t)}{\rho \cdot c}$$

This will only get you the magnitude. Keep in mind that intensity is a vector with the same direction as the particle velocity. If you want to integrate Intensity to get actual power you need to take into account the angle between the Intensity vector and the normal vector of your integration surface.

If you only have the pressure signal, you have to make some assumptions around the direction of the vector. Pressure is a scalar and there is no spatial information contained.


Ambient noise levels in desert and recording studio Are there mathematical models which suggest which of the two - desert or recording studio - is likely to be quieter? I would think a recording studio would be quieter than a desert. Is there a mathematical model which suggests this?

Is there a mathematical model which suggests this?

No. Background noise levels in recording studios and outside vary a lot from situation to situation. Outside is typically quite loud but if there is no wind, no air movement, no animal noises, and no traffic/planes/machinery/power lines levels can drop quite low.

Studios are equally variable. Some are designed to be very quiet inside, others are not.


Demodulation of a BFSK frequency hopping signal using FFT I am not very familiar with DSP techniques, and I am currently working on my research project to develop the firmware of a DSP processor, using C/Linux. It involves the modulation and demodulation of a dataset using Binary FSK and frequency hopping Tx and Rx algorithms. I am able to modulate the signals and have verified the output using MATLAB. However I an a little stuck at demodulation.

The technique that I am using for demodulation is to compare the amplitude of the FFT of the signal received at every hopping sequence and based on that, determine whether the bit received is 1 or 0. I have a predefined function which is used calculate the FFT of the signals. The function is as follows:

void fft(complex_float* in, complex_float* out, int N, int inv)
{
int i, j;
float a, f;
complex_float s, w;
f = inv ? 1.0/N : 1.0;
for (i = 0; i < N; i++) {
s.re = 0;
s.im = 0;
for (j = 0; j < N; j++) {
  a = -2*PI*i*j/N;
  if (inv) a = -a;
  w.re = cos(a);
  w.im = sin(a);
  s.re += in[j].re * w.re - in[j].im * w.im;
  s.im += in[j].im * w.re + in[j].re * w.im;
}
out[i].re = s.re*f;
out[i].im = s.im*f;
}
}

complex_float is a predefined struct :

typedef struct {
float re;
float im;
}complex_float;

In the function fft() N defines the number of FFT points.

The problem that I am facing is how to determine the frequency points of the individual elements in the array returned by fft() function. In MATLAB, this was very easy to do the FFT function in MATLAB took the sampling frequency in account and plotted the power spectrum accordingly. However, In the fft() function defined above there is no frequency involved and hence I am confused as to how to find the power of the signal at any particular frequency, or how to determine the frequency point of any element of the retuned array.

At some point, before the fft() function is called, the sampling frequency must be known (from ADC? or other resampling functions?)... use that to find the frequencies of the bins as you would in Matlab.


Root raised cosine pulse shaping peak value I am using a FIR filter to generate Root Raised Cosine pulse for an IQ modulator.

How can I predict the peak value at the output, depending on the roll-off factor of the pulse ?

Knowing that :

FIR taps are root raised cosine pulse, with a peak value of 1.0
Input of the filter is I or Q channel, valued +1 or -1, upsampled by adding 0s between each symbol

By experimenting, I can see that depending on the input sequence and the roll-off factor, the peak output varies. But I can't predict it theoretically.

This plot illustrates my question (while this particular sequence may not generates the maximum peak that I want to predict) :

You can derive an upper bound for the peak value of the output signal depending on the filter's impulse response. The output is given by

$$y[n]=\sum_{k=0}^{N-1}h[k]x[n-k]\tag{1}$$

where $h[k]$ is the length $N$ impulse response, and $x[k]$ is the input signal. From $(1)$ we have

$$|y[n]|=\left|\sum_{k=0}^{N-1}h[k]x[n-k]\right|\le\sum_{k=0}^{N-1}|h[k]||x[n-k]|\le |x[n]|_{\text{max}}\sum_{k=0}^{N-1}|h[k]|\tag{2}$$

Since in your case $|x[n]|_{\text{max}}=1$, the magnitude of the peak value of the output signal is guaranteed to be smaller than the sum of the absolute value of all filter coefficients.

However, since your input signal is upsampled by a factor $L$ (i.e., it has $L-1$ zeros between two signal values), we can come up with an improved (i.e., tighter) bound:

$$|y[n]|\le\max_{l=0,1,\ldots L-1}\sum_{k=0}^{\lfloor \frac{N-1}{L}\rfloor}|h[l+kL]|$$

because in the sums $(1)$ and $(2)$ only every $L^{th}$ filter coefficient is multiplied by a non-zero input value.

The following Matlab/Octave script shows the computation of the bound:

h = randn(50,1);    % some impulse response
x = randi([-1,1],1,20); % input signal
L = 10;         % interpolation factor
x = [x;zeros(L-1,20)]; x=x(:);  % interpolate x
y = filter(h,1,x);  % compute output signal
p = max(abs(y));    % peak value

% compute bound
b = zeros(L,1);
for i = 1:L, b(i) = sum(abs(h(i:L:end))); end
b = max(b);

% show peak value and bound
[p,b]

A few runs of the script show that the bound is pretty tight.


Analog butterworth to digital - bilinear transform In my previous question I've designed analog butterworth filter (poles own calculated). But now I would like to transform it to digital domain. I'm using bilinear transformation but not all is clear for me. Please look at my script below:

clc, clear all, close all;

%ANALOG BUTTERWORTH FILTER DESING
%https://engineering.purdue.edu/~ee538/AnalogFilterDesigns.pdf

%Butterworth LP with -3dB at fc and N order and G gain

wc = 1000;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain

%analog filter design with prewarped frequency
%poles and zeros

%frequency normalizing and prewarping
wc = 2*tan(wc/2);

%poles searching
for k = 0:N-1
   sk(k+1) = wc*exp(j*(pi/2))*exp(j*(2*k+1)*(pi/(2*N)));
end

[NUM,DEN] = zp2tf([],sk,G);
b0 = DEN(end)*G      %gain normalization

Hs = tf(b0,DEN)

figure();
freqs(b0,DEN)
title('lowpass filter with pre-warped frequency visualization')
figure();
plot(sk./10^2, 'x')
xlim([-10,10]);
title('Poles placement in s-plane')

%continuous to digital conversion
[NUMD,DEND] = bilinear(b0,DEN,1);

fvtool(NUMD,DEND)

Results looks like below:

I'm sure that is not what I want, because I have nowhere applied frequency sampling that is interesting me. I'm not also sure if I've done prewarping well.

EDIT: frequency response for ws = 2000 in bilinear transformation:

look at MATLAB bilinear().  you are specifying fs=1 in your call to it.

if you're gonna do any digital filtering, sometime, somewhere you need to commit to a sampling frequency and you haven't yet.