Why do Dyadic filterbanks downsample the high pass signal portions?

[Good question, that made me rethink of stuff I believed natural. I shall incorporate them in future lectures]

Downsampling the highpass (and the lowpass) provide you with a critically sampled filterbank, in other words, not redundant. Critical sampling is a price to pay in some applications like compression, and in the context of finding an orthogonal transform. So yes, this is a limitation: one should take care of aliasing, and of the imperfect character of the designed filters. And yes too, filterbank coefficients are subject to folding/mirroring. If one were to analyse these coefficients as if they were a genuine signal, for instance with linear filters or Fourier transforms, that could raise some concerns.

But remember that the DFT can be interpreted as a filter bank as well, downsampled by the length of the signal (plus zero-padding). Each DFT coefficient becomes, in this interpretation (correct me if wrong), a single "sample", rejected to a baseband.
I understand your concern, but it is very rare to apply a linear filter to highpass subbands. On the contrary, one usually employs nonlinear filters and selection tools.
Would you perform linear filtering or Fourier analysis on it without caution? Of course, if you keep this coefficient, set the others to zero, and inverse the DFT, you get a signal back, a sine, and then you could perform traditional processing on it.

The same applies to FIR (even IIR ones) filter banks. Here, the downsampling is much smaller. And the series of coefficients longer. But in many cases, these coefficients are not "signals" per se. Very often, they are not processed "linearly": they are cancelled, shrunk, thresholded. Most of the work focuses on both amplitudes AND  location, due to the relative locally of the filters (contrarily to DFT).

If you want to keep with the linear interpretation, the best is probably to set the other coefficients to $0$ and reconstruct with the inverse FB. Then, you have a real signal back again, and you can safely Fourier, filter, etc.

But not all wavelet representations do that: "I notice in all of the visual representations" and "why does every representation insist on downsampling the high pass subband?". This is not always the case, and you have a wealth of oversampled, shift-invariant, etc. filter banks, not only dyadic, that do a much better job.

For instance, you can avoid the downsampling on the highpass (and also on the lowpass), and you end up with a more redundant representation, that can be useful for some processing cases. Then, you can process the subbands more normally.

For some (personal) background on (oversampled or not) filter banks, which content I can speak of: in 1D Optimization of Synthesis Oversampled Complex Filter Banks, in 2D A Panorama on Multiscale Geometric Representations.


Is discrete cosine transform differentiable?

The Discrete Cosine Transform of a discrete sequence $x$ can be defined as:

$$
{\tt DCT(}x{\tt )} = X_k = \frac{1}{2} (x_0 + (-1)^k x_{N-1})
 + \sum_{n=1}^{N-2} x_n \cos \left[\frac{\pi}{N-1} n k \right] \quad \quad k = 0, \dots, N-1.
$$

If, by differentiable, you mean:
$$
\frac{\partial {\tt DCT(}x{\tt )}}{\partial x}
$$
then isn't it just
$$
\frac{\partial {\tt DCT(}x{\tt )}}{\partial x} = \frac{\partial X_k}{\partial x} =\left [
\begin{array}{ccccc}
\frac{1}{2} & \cdots & \cos \left[\frac{\pi}{N-1} n k \right] & \cdots & (-1)^k\frac{1}{2}
\end{array}
\right ]
$$
where $x = [ x_0, x_1, \ldots, x_{N-2}, x_{N-1} ]^T$ ?

What do you mean by a "mathematical proof"?


I think you can't do Homography on lines.

I don't remember the equations but I don't think lines in one image are kept as lines in the other image.

What you can do is treat a line as set of points and have Homography estimation using many point samples.


Demodulation in Rayleigh Flat Fading Channels

The short answer is: pilot symbols.

The long(er) answer is: The sender periodically transmits known data symbols; this sequence is known to the recipient ahead of time.


Why is the gain of my IIR filter positive?

You have your a and b vectors reversed. For a first-order IIR filter, the a vector should have two elements (and you almost always want the first element to be unity). So, you're plotting the frequency response of the filter with difference equation $0.544062y[n]=x[n]+0.45594x[n−1]$, which is not what you want.

Also, as a note on MATLAB syntax, you'll also need to negate the second element of a (multiply it by -1) in order to get what you want. Read the documentation on functions that take IIR filter coefficients to learn why.


Discrete Time Fourier Transform (DTFT) for an unstable system (Ideal Low Pass Filter)

It is generally not true that the relation $Y(\omega)=X(\omega)H(\omega)$ is independent of the system's stability. For systems with a rational transfer function (i.e., systems that can be described by differential or difference equations), the frequency response, which is the Fourier transform of the impulse response, only exists for stable systems. For marginally stable systems with poles on the imaginary axis (for continuous-time systems), or on the unit circle (for discrete-time systems), the Fourier transform exists if we allow Dirac delta impulses in the expression for the Fourier transform. So for marginally stable systems you can see that the output can become unbounded, even for bounded input signals, due to the Dirac delta impulses in their frequency response. For unstable causal systems with transfer function poles in the right half plane, the frequency response doesn't exist.

For systems for which there exists no transfer function, but which do have a frequency response, such as ideal frequency-selective filters (low pass, high pass, etc.), the relation $Y(\omega)=X(\omega)H(\omega)$ holds, regardless of stability. Ideal frequency-selective filters are unstable (and non-causal), but their frequency response $H(\omega)$ is finite. This means that the frequency response of the bounded input signal which generates an unbounded output (if it exists) must be unbounded for at least one value of $\omega$, and the inverse Fourier transform of $X(\omega)H(\omega)$ must result in an unbounded time-domain signal. Also note that there is no guarantee that the Fourier transform $X(\omega)$ of the bounded input signal that generates an unbounded output signal even exists. And if it exists, it's generally not immediately clear from the expression $X(\omega)H(\omega)$ that the corresponding inverse transform is unbounded.

I think it is a misunderstanding that the relation $Y(\omega)=X(\omega)H(\omega)$ somehow implies that a bounded input must result in a bounded output. It's just not so easy to see whether the output is unbounded because one would have to compute

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)H(\omega)e^{jn\omega}d\omega$$

and check for each value of $n$ if the expression becomes unbounded. Direct time domain analysis is more straightforward in this case.


Shifting the data points (fftshift) and zero-padding the exact center of the FFT aperture has the property that all the even (symmetric) components with reference to the center of the original data set end up in the real part of the complex FFT result, and all the odd components end up in the imaginary part. e.g. the evenness to oddness ratio is preserved, which allows phase (with reference to the center of the window) to be easily interpolated.

Being able to interpolate phase is important in the case of a zero-padded FFT because the zero-padding interpolates the spectral magnitude as well.  Since zero-padding results in an interpolated FFT result, that means that any original non-interpolated FFT result points will need to be interpolated from the zero-padded result.  Without an FFT shift of both the non-padded and zero-padded data, these interpolated results will be different from the non-interpolated FFT results (in phase).

This technique is a simple (homework or quiz level) result of the FT property that a shift in one domain is a complex frequency modulation in the other domain.


"Group delay" isn't the delay between the change on the input and the first effect; it's the delay that a packet of oscillations of different frequencies experience.

In the case of a linear phase filter (and your moving average, its impulse response being symmetric, is linear phase) that is the "center" of the effects of a single impulse. So, it's the center between you seeing the first effect of your input change, and the last effect.


You do get a time-invariant response. Your code produces the same output for all three signals. In particular, it produces the same output for $y(\sigma^T\{x(t)\}))$ as for $\sigma^T\{y(x(t))\}$ (plots 2 and 3 in your code). It is hard to see in your case because you have shifted the signal $3\cdot 2\pi$ in time. Whats a cosine shifted by $6\pi$? The same cosine. Everything will be clearer if you plot the signals on top of each other and shift by a more informative value.

Try this code:

clc
clear all
close all
t=0:.001:10;
x=cos(2*pi*t);
y=x.^2;
plot(t,y)
hold on;

plot(t+0.25,y)

xn=cos(2*pi*(t-0.25));
y2=xn.^2;
plot(t,y2)

legend('y(t)','y(t-3)','y[x(t-0.25)]')

hold off

You can see in the attached image, the orange and red lines, corresponding to shifted output and output on shifted input respectively, align with each other. That is why you can not see the red line up until the orange line ends. The unshifted output of the unshifted signal in blue does not align with the two other.


Help with a product of a sum of infinite series

let me start from the general form
$$ \sum_{ k = M}^{ N } r^{k} = \frac{r^M - r^{N+1}}{ 1 - r}$$

then if you multiply $r^k$ by a unit step $u[n]$ the above series become

$$ \sum_{ k = M}^{ N } r^{k} \cdot u[k] = \sum_{ k = 0}^{ N } r^{k} $$

As you can notice the lower index become zero, because the unit step is equal to one only if $k\geq0$, otherwise is equal to zero.

Therefore in this specific case we have $u[n+4r]$ that is equal to one if the index $n+4r$ is greater or equal to zero, then:
$$ r \geq - \frac{n}{4}$$

At this point I think that you have all the necessary information to resolve the series


First of all, computing the poles of an (ideal) Gaussian filter is an impossible task, because its transfer function is not a rational function, and there are simply no poles. This is in contrast to the classic filter approximations, like Butterworth, Chebyshev, and Cauer filters, all of which have rational transfer functions, and, consequently, all of them have poles.

The frequency response of an ideal Gaussian filter is

$$H(\omega)=e^{-a\omega^2}\tag{1}$$

The corresponding transfer function is a function of $s=j\omega$:

$$G(s)=e^{-a(s/j)^2}=e^{as^2}\tag{2}$$

which is of course not rational, and which doesn't have any finite poles.

What you did is compute a rational approximation by developing $1/G(s)$ into a power series, and using only a finite number of terms. However, since you used $\omega$ as a complex variable instead of $s$, the poles are turned by $90$ degrees. This was already mentioned in the comments. The correct way of doing this is as follows:

$$1/G(s)=e^{-as^2}=\sum_{k=0}^{\infty}\frac{(-a)^k}{k!}s^{2k}\tag{3}$$

Using only a finite number of terms gives the following rational approximation of $G(s)$:

$$G(s)\approx\frac{1}{\sum_{k=0}^{n}\frac{(-a)^k}{k!}s^{2k}}\tag{2}$$

the poles of which are of course the roots of the denominator polynomial.


As you can see in the following image, the image shows the first order 1D derivative.
Now you can write this equation in terms of the previous pixel rather than the following pixel.

For a 1D differentiation, you are only interested in either the x direction (Horizontal) changes of pixel intensity values or the y direction (Vertical).

Therefore the mask will be [0 -1 1] in case of next intensity value.
Or [1 -1 0] for previous intensity value.

You can apply the same idea for vertical differentiation by taking the transpose.


In your question, you postulate the existence of two filters, let's say $p(t)$ and $g(t)$, such that $p(t) \star g(t)=\delta(t)$. This assumption is problematic because it implies that $G(f)=1/P(f)$, which requires $P(f) \neq 0$ for all $f$. Furthermore, for small values of $P(f)$, $G(f)$ will take arbitrarily large values (see Dilip Sarwate's comments above and below this answer).

Assuming (for the sake of discussion) the existance of such a pair of filters, there may not be any advantage to using them. I'll describe a possible application from digital communications. Let's say you want to transmit a number $A$ over an analog channel. You choose an appropriate analog pulse $p(t)$ and transmit $$s(t) = A\delta(t)\star p(t).$$ (The pulse $p(t)$ could be chosen to fit the channel bandwidth, and/or to have a certain energy). The receiver could perform the "inverse filter" operation on the signal $s(t)$: $$s(t) \star g(t)=A\delta(t),$$ where $g(t)$ is your "inverse matched filter". However, in practice it turns out that detecting $A$ in this way is not optimal, because it does not have the best signal-to-noise ratio. In other words, when you add noise to the transmitted signal, the received signal becomes $$r(t)=s(t)+n(t).$$ In this case, filtering $r(t)$ with the "inverse matched filter" $g(t)$ results in a worse signal-to-noise ratio than filtering with the filter matched to $p(t)$.

If you have access to it, I highly recommend "An Introduction to Matched Filters", by G. Turin, IRE Transactions on Information Theory, June 1960. Your question is answered in page 318.


Reproducing paper results about a wavelet transformation using python

With a fake imitation of your data, here is a first try in Matlab, which I may improve later on. It uses the stationary wavelet transform swt.m:

nSample = 256;
timeOnset = 25;

time = linspace(1,nSample,nSample)';
%ramp = max(time,timeOnset)/nSample;
ramp = abs(time-timeOnset)/nSample;

nLevel1 = 4;
nLevel2 = 5;
waveletName = 'sym8';
data = 0.15*((1+rand(nSample,1)/1).^0.5).*sin(2*pi*time/nSample*42)+0.65*(ramp+1);
dataSWT1 = swt(data,nLevel1,waveletName);
dataSWT2 = swt(data,nLevel2,waveletName);
subplot(3,1,1);
plot(time,data);axis tight
subplot(3,1,2);
imagesc(dataSWT1);axis tight
xlabel([num2str(nLevel1),' levels'])
subplot(3,1,3);
imagesc(dataSWT2);axis tight
xlabel([num2str(nLevel2),' levels'])


DFT exercise in the book Understanding digital signal processing 3 Ed

@Henrique Luna. Please forgive me. In Part (b) of the problem the words "sequence time" should "time sequence". Sorry for the confusion!
Years ago when I created that Part (b) question I was thinking about the answer to the question; "What is the time duration of an N-sample time-domain sequence?" Back then I believed the time duration under consideration to be the time difference between the sequence's first and last samples, rather than N times the time interval (1/fs) between the samples. Thinking this way leads to the conclusion that Part (b)'s answer is (N-1)*1/fs = 0.999977 seconds.
But subsequent discussion, years ago, with Dilip Sarwate convinced me that my concluded answer to 3.3 Part (b) may be somewhat valid when thinking about the time-domain but my concluded answer leads to confusion when thinking about signal durations in the context of spectrum analysis.
I now believe the 3.3 Part (b) question leads to more confusion than enlightenment. That Part (b) question was not my most well thought-out homework problem.
By the way Henrique, if you have not done so already, please have a look at the following web page: https://www.dsprelated.com/showarticle/1094.php


What is an octave? It is said to be eight, but I cannot see the eight anywhere

This comes from music terminology. The name "octave" comes from the fact that in the heptatonic musical scales (which are the prevalent scales in western music), the note with a 2:1 frequency ratio is the eighth note in the scale.

For example, in the C major scale (C D E F G A B C) the eighth note is one octave above / has a 2:1 frequency ratio with the tonic.


What does MATLAB display when plotting a sound signal?

What you are observing is the digital representation of the voltage, which in fact represents the acoustic pressure.
Workflow would be something like:

Vibrating larynx is producing Acoustic Pressure [Pa]
Variations of that pressure are converted by the microphone into Voltage [V]
Voltage is being sampled and quantized by the ADC (Analog to Digital Converter) and represented as series of numbers.

Now knowing your microphone's sensitivity [mV/Pa] you can relate the voltage and acoustic pressure. Additionally if ADC's parameters are known (or measured) then you can directly relate sample values to Pascals.


In order to avoid inter symbol interference (ISI) the baseband transmit signal must meet Nyquist's first criterion. In frequency domain it can be formulated as follows.
$$
\omega_\mathrm{N} = \frac{v_\mathrm{s}}{2}
$$
With $v_\mathrm{s}$ the symbol rate and $\omega_\mathrm{N}$ the 3dB frequency of a low pass filter with symmetric flank. The flank must be point symmetric with regard to $\omega_\mathrm{N}$ and zero for $\omega > \omega_\mathrm{N}(1 + \alpha)$. Consequently, the baseband bandwidth $B=\omega_\mathrm{N}(1 + \alpha)$.

After modulation the bandwidth is doubled and the radio frequency bandwidth is $2B = 2\omega_\mathrm{N}(1 + \alpha) = v_\mathrm{s}(1 + \alpha)$. For the ideal case of a rectangular low pass filter $\alpha=0$ and $2B = v_\mathrm{s}$ which corresponds to the statement that the required bandwidth is equal to the symbol rate *). Similarly, if the bandwidth $B$ is given, the symbol rate is limited by $v_\mathrm{s} \leq \frac{2B}{1+\alpha}$

The required bandwidth indirectly depends on the bitrate, because a certain number of bits is transmitted per symbol.

The lowpass mentioned above is also referred to as impulse shaper and $\alpha$ as roll-off factor.

*) As an ideal lowpass isn't feasible, the required RF bandwidth in practice is always larger than the symbol rate.


comparison between frequency offset estimators

Well, since you simulate you could calculate the residual frequency offset which is basically the remaining CFO after correction or the difference between the true CFO and the estimated CFO.

Another thing that is often done is to give the resulting Bit Error Rate with different estimators.

But really, please ask your supervisor what he has in mind. Everything else is guessing and probably wasted time.


How do I find the difference equation from a transfer function?

A transfer function $H(Z)$ can be written as $H(Z)=\frac{Y(Z)}{X(Z)}$. Then, your $H(Z)$ can be written as

$\frac{Y(Z)}{X(Z)}=1-\cos\theta~Z^{-1}+Z^{-2}$ or

$Y(Z)=X(Z)(1-\cos\theta~ Z^{-1}+Z^{-2})$

Now taking the inverse $Z$transform, we get the difference equation as

$y(n)=x(n)-\cos\theta~x(n-1)+x(n-2)$.

I hope this help you.


LMS and delta impulse response- Equalization concepts

You're right. LMS equalizer uses a known input to minimize the error. For communication purposes, this is either provided by a training sequence, or in a decision directed mode, the detector decisions are fed back as known data.

The delta function is also correct. Suppose that the channel impulse response is $h(t)$ and frequency response $H(f)$. Then the equalizer should ideally be the inverse of this function in frequency domain (remember the convolution in time domain and hence multiplication in frequency domain).
$$H(f)\cdot \frac{1}{H(f)} = 1$$
and its inverse Fourier transform is
$$F^{-1}(1) = \delta (t)$$

The fourth tap being 1 implies that there is a delay of 4 symbols with respect to your time base at the Rx. Finally, to further verify your results, plot the frequency responses of the channel and the equalizer response, and you will see an approximately inverse relationship.


That is because envelope centers the signal before computing the standard deviation and afterwards corrects for the mean. Try this instead:
x = randn(50, 1);
xRMS = sqrt(movmean(x.^2, 21));
xRMS_centered = sqrt(movmean((x-mean(x)).^2, 21)) + mean(x);
xRMSref = envelope(x, 21, 'rms');
plot(xRMSref,'x-','DisplayName','xRMSref');hold on;plot(xRMS,'o-','DisplayName','xRMS');plot(xRMS_centered,'s-','DisplayName','xRMS centered');hold off;
legend()

*edit: since you were asking why envelope behaves this way: I think it's because it tries to return an envelope of the signal. It actually comes with two output arguments, a lower and an upper envelope with the idea that the signal should be between the two. Clearly if that's what we want, we need to account for the mean. See example in the following figure:

If you didn't center, you wouldn't capture the local variation around the mean accordingly.


Software Mixing of Multiple Audio Streams From Network

Generally you will need some sort time sync protocol.
Typically you want create locally one circular buffer for each network channel. As soon as a packet comes in, you shove it into circular buffer ort that channel at the write positions.
Playback is clocked by your local D/A clock. As soon as the D/A buffer is running empty, you grab a frame from all the client circular buffers at the current read position, mix them together and put them into the D/A buffer.
There are two things you need to manage.

Drop out: That happens if you need a new playback frame and the circular buffer of a client is empty OR a packet got dropped (they should be numbered to detect this). There isn't a whole lot you can do about. You can play zeros or try to extrapolate from the previous buffer A drop out means that your network can't keep up.
Clock drift. Even if everyone is playing at nominally 44.1kHz, no two clocks are identical and over time they will drift apart. To fix these you need to a) measure the clock drift and b) correct it. a) can be done through either a global time reference (something like NTP), by averaging over packet arrival times, or monitoring buffer levels. b) will require some sort of sample rate conversion: the most crude option would be to occasionally drop a frame insert a zero frame.

This is age old trade off: the lower the latency, the higher your drop out rate will be. So dialing this in "optimally" is the key here. It's easy enough to adjust: simply make your circular buffers smaller/bigger. However finding the optimum is hard since it depends very much on current network conditions so in practice this really should be adaptive.


How to derive the "well-known" solution to Unconstrained Array Gain?

You can solve such a problem using the method of Lagrange multipliers. First note that maximizing the expression in your question is equivalent to minimizing the inverse function:
$$\min_{\mathbf{w}}\frac{\mathbf{w}^H\mathbf{Q}\mathbf{w}}{|\mathbf{w}^H\mathbf{d}|^2}\tag{1}$$
Next note that the solution of $(1)$ is invariant to scaling of $\mathbf{w}$, i.e., replacing $\mathbf{w}$ by $c\cdot\mathbf{w}$ in $(1)$ with an arbitrary scalar constant $c$ will not change the value of the function. So we may as well use a scaling such that $\mathbf{w}^H\mathbf{d}=1$ is satisfied. This scaling corresponds to a unity response for the desired signal. With this constraint, problem $(1)$ can be reformulated as
$$\min_{\mathbf{w}}\mathbf{w}^H\mathbf{Q}\mathbf{w}\qquad\textrm{s.t.}\qquad \mathbf{w}^H\mathbf{d}=1\tag{2}$$
We can solve $(2)$ using the method of Lagrange multipliers by minimizing
$$\mathbf{w}^H\mathbf{Q}\mathbf{w}-\lambda(\mathbf{w}^H\mathbf{d}-1)\tag{3}$$
Formally taking the derivative of $(3)$ with respect to $\mathbf{w}^H$ and setting it to zero gives
$$\mathbf{w}=\lambda\mathbf{Q}^{-1}\mathbf{d}\tag{4}$$
The constraint in $(2)$ is satisfied for
$$\lambda=\frac{1}{\mathbf{d}^H\mathbf{Q}^{-1}\mathbf{d}}\tag{5}$$
From $(4)$ and $(5)$ we finally obtain
$$\mathbf{w}=\frac{\mathbf{Q}^{-1}\mathbf{d}}{\mathbf{d}^H\mathbf{Q}^{-1}\mathbf{d}}\tag{6}$$
Note that the scaling in $(6)$ is optional and the general solution is given by $(4)$.


What should the amplitude be when plotting 1-sided Amplitude Spectrum?

To complement @wrapperapps answer, you can express sinusoids as sum of exponentials: $\cos(2 \pi f_0 t) = (e^{j 2 \pi f_0 t} + e^{-j 2 \pi f_0 t})/2$.

Thus, the fourier transform yields $\delta(f+f_0)/2 + \delta(f-f_0)/2$: each frequency delta has amplitude $1/2$.


There is a mistake in the last step. When you invert both sides of the inequality, the $>$ inverts as well. So, for $K>0$ (as you found that $K>2$, so it must be positive):

$$\frac 1 K > y \iff K<\frac 1 y$$

If you do this when $y\to 0$, then you get the desired result.

Also, note that $\frac 1 K >0$ would hold for any $K>0$, as the only way to get a negative number from $1/K$ is letting $K$ be negative.

Regarding your other question, it would completely depend on the transfer function and its respective Nyquist. But don't worry, the procedure is exactly the same. For example, try doing the Nyquist plot for $-\frac{K}{1-s}$ and see what you get.


Any estimation algorithms that allow a user to indicate that certain Solutions are impossible?

This series of papers tries to enforce such constraints.

The constraints are on the state of the Kalman filter:


The Fourier transform of a triangular window function

Not the perfect fit, but a quick-and-dirty base to play on:

%nBin = 81;
nBin = 2049;

y=(20*log(fftshift(abs(fft(triang(81),nBin)))))';
x = linspace(-40,40,length(y));
basey = min(0,min(y));
ShadingColor = [0 0 1];
h = fill([x x(end) x(1)], [y basey basey], ShadingColor);
set(h, 'EdgeColor','none', 'FaceAlpha', 0.5);


Under what conditions is there a one-to-one mapping between continuous-time and discrete-time signals?

Apart from practical errors related to the truncation of the infinitely-long energy signal when discretized to digital, can we find a UNIQUE mapping between the two signal forms?

Yes.
Say you're starting with $x(t)$, which is perfectly bandlimited to less than $f_B = \frac{1}{2 T_s}$, and you take samples at every $t = T_s k$: $x_k = x(T_s k)$.
All you need to do is to treat the  as a train of Dirac impulses (i.e., let $x_p(t) = \sum x_k \delta(t - T_s k)$), and then filter $x_p(t)$ with a brick-wall low-pass filter.  I.e., let $x_f(t) = h(x_p(t))$, where $h$ is that brick-wall filter.  $h$ will have an infinitely long impulse response, so you cannot realize this system in the real world - but not only will the mapping from $x_k$ to $x_f$ be exact, but if $h$ has zero phase delay then you'll find that $x_f(t) = x(t)$.


Ideal derivative filter

Let $f(x)$ be a signal bandlimited to frequencies $(-\pi,\, \pi)$. Given $f(x)$ as input, the same $f(x)$ is given as output by a system that has as its impulse response the sinc function:

$$\operatorname{sinc}(x) = \left\{\begin{array}{ll}1&\text{if }x = 0,\\
\frac{\sin(\pi x)}{\pi x}&\text{otherwise.}\end{array}\right.\tag{1}$$

Taking the derivative $f'(x)$ of signal $f(x)$ is a linear time-invariant operation. By associativity, $f(x)$ is differentiated by a system that has as its impulse response the derivative of $\operatorname{sinc}(x):$

$$\operatorname{sinc}'(x) = \left\{\begin{array}{ll}0&\text{if }x = 0,\\
\frac{\cos(\pi x)}{x} - \frac{\sin(\pi x)}{\pi x^2}&\text{otherwise.}\end{array}\right.\tag{2}$$

Both of $\operatorname{sinc}(x)$ and $\operatorname{sinc}'(x)$ are bandlimited  to $(-\pi,\, \pi)$ and can thus be sampled at integer $x$ without aliasing. Sampling $\operatorname{sinc}'(x)$ at integer $x$ gives the ideal discrete-time impulse response for you to use in filtering samples of $f(x)$ at integer $x$ to obtain samples of $f'(x)$.

Figure 1. $\operatorname{sinc}(x)$ (red) and its derivative $\operatorname{sinc}'(x)$

Windowed filter

The infinitely long ideal impulse response can be multiplied by a window function to obtain a realizable impulse response. SciPy can calculate several types of window functions and do finite-impulse-response (FIR) filtering with an arbitrary impulse response by scipy.signal.convolve. You need to delay the impulse response to make it that of a causal system.

Least squares derivative filter

A least squares filter impulse response is obtained by sampling $\operatorname{sinc}'(x)$ from Eq. 2 symmetrically at integer $-N\le x\le N$ with $N$ determining the filter order. This includes the samples of the impulse response that have the largest absolute value and consequently have the largest contribution to the mean square error.

Least squares derivative filter for pre-oversampled data

For sampled data that is oversampled by a factor $\beta \ge 1$ compared to the sufficient sampling frequency of $1$, a least squares derivative filter is obtained by minimizing mean square error (MSE) or deviation of the Fourier transform of the impulse response from the Fourier transform $i \omega$ of derivation, over the bandwidth $\omega\in(-\pi,\,\pi)$ of sinc:

$$\begin{align}\mathrm{MSE} &= \frac{1}{2\pi}∫_{-\pi}^{\pi}\left|iω - \sum_{n=1}^N \left(c_n e^{-niω/β} - c_n e^{niω/β}\right)\right|^2\\
&= \frac{1}{\pi}∫_{0}^{\pi}\left(ω + 2 \sum_{n=1}^N c_n \sin\left(\frac{nω}{β}\right)\right)^2,\end{align}\tag{3}$$

where the unorthodox notation $e^{-iω/β}$ represents a one-sample delay at the oversampled sampling frequency or a delay of $1/β$ samples at the critical sampling frequency of the sinc, and $c_n$ is the value of the antisymmetric impulse response at discrete time index $n=1\ldots N$ with a complementary value $-c_{n}$ at time $-n$. The reason for the unorthodox notation is that it gives MS values that are comparable between different $\beta$. If you require coefficients for calculation of the derivative as if the sampling frequency was $1$, then divide each coefficient (solved in the following) by $\beta$. MSE is minimized when partial derivatives of MSE with respect to all $c_n$ are zero. At $\beta = 1$ we get the least squares filter discussed earlier. The least squares solutions and the corresponding impulse responses can be calculated by this Python script:

from sympy import *
for N in [1, 2, 3, 4]:  # <------ number of non-zero coefs / 2, too large is too slow to solve
    omega, beta = symbols('omega beta', real=True)
    c = [Symbol('c_'+str(i + 1), real=True) for i in range(N)]
    MSE = integrate((omega + 2*sum([c[n]*sin((n + 1)*omega/beta) for n in range(N)]))**2, (omega, 0, pi))/pi
    for use_beta in [1, 1.5, 2, 4, 8]:  # <------- Oversampling factor
        use_MSE = MSE.subs(beta, use_beta)
        sol = solve([diff(use_MSE, c[n]) for n in range(N)], [c[n] for n in range(N)])
        for n in range(N):
            print(str(-sol[c[N - n - 1]].evalf()), end=',')
        print('0', end=',')
        for n in range(N-1):
            print(str(sol[c[n]].evalf()), end=',')
        print(str(sol[c[N-1]].evalf()), end='  ')
        print('# N='+str(N)+', beta='+str(use_beta)+', MSE='+str(use_MSE.subs(sol).evalf()))
    print()

SymPy's solver will hang on some inputs, but manages to find the solutions for $N = 1\ldots4$ and $\beta = 1, 1.5, 2, 4, 8$:

-1.00000000000000,0,1.00000000000000  # N=1, beta=1, MSE=1.28986813369645
-1.13548530933771,0,1.13548530933771  # N=1, beta=1.5, MSE=0.178081981619031
-1.27323954473516,0,1.27323954473516  # N=1, beta=2, MSE=0.0475902571416442
-2.12680294939990,0,2.12680294939990  # N=1, beta=4, MSE=0.00251926307592025
-4.06211519814366,0,4.06211519814366  # N=1, beta=8, MSE=0.000151083662191318

0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000  # N=2, beta=1, MSE=0.789868133696453
0.330596796032270,-1.24876549439064,0,1.24876549439064,-0.330596796032270  # N=2, beta=1.5, MSE=0.0130608565500205
0.288948415598137,-1.51850657748724,0,1.51850657748724,-0.288948415598137  # N=2, beta=2, MSE=0.000919718098315458
0.382389721267460,-2.75841263168952,0,2.75841263168952,-0.382389721267460  # N=2, beta=4, MSE=2.58832725686194e-6
0.689919219439685,-5.37907341717520,0,5.37907341717520,-0.689919219439685  # N=2, beta=8, MSE=9.33496011568591e-9

-0.333333333333333,0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000,0.333333333333333  # N=3, beta=1, MSE=0.567645911474231
-0.116742660811357,0.423760333529959,-1.31069003715342,0,1.31069003715342,-0.423760333529959,0.116742660811357  # N=3, beta=1.5, MSE=0.00108046987847814
-0.0790631626589370,0.433017934639585,-1.64079658337694,0,1.64079658337694,-0.433017934639585,0.0790631626589370  # N=3, beta=2, MSE=2.01122717118286e-5
-0.0825922199951234,0.659247420212084,-3.07101401580036,0,3.07101401580036,-0.659247420212084,0.0825922199951234  # N=3, beta=4, MSE=3.01168290348385e-9
-0.140650397261751,1.22875551919474,-6.03556856432330,0,6.03556856432330,-1.22875551919474,0.140650397261751  # N=3, beta=8, MSE=6.53134493440793e-13

0.250000000000000,-0.333333333333333,0.500000000000000,-1.00000000000000,0,1.00000000000000,-0.500000000000000,0.333333333333333,-0.250000000000000  # N=4, beta=1, MSE=0.442645911474231
0.0443010798143126,-0.173517915867411,0.482863594553193,-1.34856865228067,0,1.34856865228067,-0.482863594553193,0.173517915867411,-0.0443010798143126  # N=4, beta=1.5, MSE=9.55482684708157e-5
0.0232132887981686,-0.144357365428612,0.528060960838853,-1.71358998904610,0,1.71358998904610,-0.528060960838853,0.144357365428612,-0.0232132887981686  # N=4, beta=2, MSE=4.70708311492379e-7
0.0191194948759654,-0.179273078273186,0.859556991552367,-3.25775951118075,0,3.25775951118075,-0.859556991552367,0.179273078273186,-0.0191194948759654  # N=4, beta=4, MSE=3.75015864227226e-12
0.0307242274505004,-0.317445566053152,1.62921582434694,-6.42899176483137,0,6.42899176483137,-1.62921582434694,0.317445566053152,-0.0307242274505004  # N=4, beta=8, MSE=4.88983660976190e-17

The filter with N=2, beta=2, MSE=0.000919718098315458 has the same number of taps as Rick Lyons' filter −3/16, 31/32, 0, −31/32, 3/16 that is slightly suboptimal in least squares sense with MSE=0.0009390870.


Measuring how "non-noisy" a sound signal is

If the noise is (close to) stationary, you could use one of the many noise-floor estimation methods to get an estimate of the noise power and, consequently, of the SNR. Search for "noise floor estimation" and "noise floor tracking". You'll find a lot of publications and patents.


Reason that it's clicking is that you are jumping from 1 to 0 in the time that the audio signal reads the last of the 1 control signal and then the 0 control signal.

You have to slow fast signal transitions over multiple samples, for example 1 to 0 should take 5ms or 10ms which is dozens or 100ds of samples.

To do that, you can use a lowpass filter at about 300Hz more or less, which can actually be very cheap for a single monophone LFO its less that 1% cpu useage, or you can use a signal smoothing algorithm, for exampel a linear interpolation function or a SmoothStep which is awesome general knowledge for audio programmers.
http://paulbourke.net/miscellaneous/interpolation/

BTW what is the platform? i can't give you code exampels if i dont know what code you use.

the cool thing about smoothstep is how easy it is, its very good knowldedge in synths:
https://en.wikipedia.org/wiki/Smoothstep

https://documentation.apple.com/en/logicstudio/instruments/index.html#chapter=A%26section=3%26tasks=true

check the filtered sawtooth there.

You can also a smoother that says "if difference between last signal and this signal was larger than 0.02, signal - 0.02... and create a secondary smothed signal... there are fairly simple smoothing algo's based on that it's pretty straightforward.

There is also an entire page on the topic.


Practical cross-spectrum estimation using Blackman-Tukey approach

Let $$\hat{\phi}_{yu}(\omega) = \sum_{k = -M}^{M}  w(k) \ \hat{r}_{yu}(k) e^{-i\omega k}\tag{1}$$
apply a change of variable $$k'=k+M+1\Rightarrow k=k'-M-1$$

and $(1)$ becomes

$$\begin{align}\hat{\phi}_{yu}(\omega) &= \sum_{k' = 1}^{N}  w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega (k'-M-1)}\\[10pt]
&=e^{i\omega (M+1)}\sum_{k' = 1}^{N}  w(k'-M-1) \ \hat{r}_{yu}(k'-M-1) e^{-i\omega k'}\\
&=e^{i\omega (M+1)}\sum_{k' = 1}^{N}  w'(k') \ \hat{r}'_{yu}(k') e^{-i\omega k'}
\end{align}$$

where $N=2M+1$ and $w'$ and $\hat{r}'_{yu}(k')$ are the shifted vectors as you did in your own try.

The above answer is when you want to use linear shift. If you wan to use circular shift, then do it like the following:

 x=w.*rhat;
 phi=fft(circshift(x',M+1)');

where M is even and w and rhat are just the symmetric vectors around M+1 (for instance w=[2 4 6 4 2] is symmetric around M=3).


What technique can I use to remove useless information from an image?

Although you did not spent even a minute in researching the question I will post an answer to it. There are multiple ways; I will try to demonstrate it using wavelets in Mathematica.
So, first of all we need an image.

img = ExampleData[{"TestImage", "Mandrill"}]

Then we apply the DiscreteWaveletTransform using the HaarWavelet

dwd = DiscreteWaveletTransform[img]

We plot the transformation

WaveletImagePlot[dwd]

We threshold it, i.e. removing the noise according to some criteria

wtd = WaveletThreshold[dwd]

Again plot the transformation and the filtering

WaveletImagePlot[wtd]

We now inverse the transformation

iwd = InverseWaveletTransform[wtd]

And visualize the differences between the two images

ImageDifference[iwd, img]


Demodulating wideband FM signal in Octave from recorded baseband IQ samples

This should work if you further filter to 15KHz after your FM discriminator. FM stereo is actually sub-modulated as follows:

L+R: 0-15KHz

Pilot Tone: 19KHz

L-R: 23KHz-53KHz as DSB-SC on a 38 KHz sub-carrier

Where L is the left channel and R is the right channel.

Update: I found a nice graphic showing this from Keysight with further details on FM broadcasting at this link: http://rfmw.em.keysight.com/wireless/helpfiles/n7611b/Content/Main/FM_Broadcasting.htm

If the pilot tone is present, that indicates the station is transmitting stereo but with this approach mono receivers would still work as you would hear the combined left and right channels.

So the composite signal is formed by first AM modulating using Double-Sideband-Supressed-Carrier (DSB-SC) the L-R channels onto a 38 KHz sub-carrier (subtract R from L and then simply multiply this with a 38Khz tone).
The tone is frequency divided by two and this coherent 19Khz pilot tone is added to the DSB-SC L-R waveform. Finally R is added to L and this baseband signal is added.

This is the signal that is then FM modulated, so would be the signal you would expect to get if you are properly FM demodulating.

A stereo demodulator after FM demodulation locks onto the 19KHz pilot to create the AM sub carrier at 38Khz (2x) and with that easily demodulates the AM DSB-SC sub-modulation (by multiplying the two) to achieve L-R at basebnad, and then can extract L and R by summing and subtracting the L+R and L-R results.

If you want to do a simple test of audibly hearing the demodulation of the FM signal, filter out the 19 KHz pilot and higher frequencies after you have FM demodulated this composite signal and you will hear the combined L+R channels.


The total (standard) RMS, continuous time:
$$
x_{RMS}=\sqrt{\frac 1{T_0} \int_0^{T_0} x^2(\tau)d\tau}
$$

And the moving RMS, continuous time:
$$
x_{RMS}(t)=\sqrt{\frac 1{T} \int_{t-T}^t x^2(\tau)d\tau}
$$

Are essentially the same operation, the first one taken over the total number of samples $T_0$ and the second taken over a small number of samples $T$.

Note also, the moving operation can be taken at any reference $t$, hence it behaves like an operator, like a filter.

In discrete form, it is only required to exhange the integrator by a sum, and the time with the indexes.
The total (standard) RMS, discrete time:
$$
x_{RMS}=\sqrt{\frac 1{N_0} \sum_{k=1}^{N_0} x^2[k]}
$$

And the moving RMS, discrete time:
$$
x_{RMS}[i]=\sqrt{\frac 1{N} \sum_{k=i-N+1}^i x^2[k]}
$$

Note that both discrete and continuous tends to be the same quantities, provided that $T=Nd\tau$ with $d\tau$ the sampling time.


DTFT of window function applied to input signal

It seems you are mixing discrete and continuous notations and it's unclear which it is for time and frequency. You use the term DTFT which would indicate both are discrete but you use the frequency index $\omega$ which would be continuous in the frequency domain.

Your time domain equation doesn't make any sense, since your left side depends on $n$ which doesn't show up in the right hand side at all. The way it's written $x[n]$ is simply a constant. I'm not just nit-picky, the exact way you write & define this can make a big difference.

  I was thinking that multiplication in time domain is convolution in
  frequency domain.

Correct

  But X(ω) is 0 except for two frequencies

Maybe. $X(\omega)$ is, but $X[k]$ may or may not be depending on the exact relationships of $N$, $\omega_1$, $\omega_2$ and the sample rate.

  so the convolution integral is 0?

Nope, it's not.

If the frequency domain is continuous  $X(\omega)$ is non-zero in a special way, so that the convolution integral is perfectly well defined and non-zero. If it is discrete, the convolution is a sum, not an integral.


discrete Haar wavelet transform, fast and efficient method?

The Haar wavelet is actually a part of the Daubechies wavelet, for the case D=2. There's some example code on wikipedia that shows the Daubachies transform.

The Haar transform is just a low pass filter combined with a high pass filter, with the coefficients being placed in the first and second halves of the signal. Then, this iteratively (or recursively) keep going.

The low pass filter is just the sum of two components while the high pass filter is their difference.

For the Haar wavelet case, you can check out my vectorized code on Github. Here, idwt2 stands for "Inverse Discrete Wavelet Transform 2D."


Relationship between Carrier Frequency Offset and the Constellation Diagram

I think that you misunderstand the problem. The actual problem is that the frequency offset is unknown to the receiver. I won't give away the answer since this appears to be a homework project, but I'll try to help you by giving a hint. First, imagine what would happen if there was an unknown carrier phase offset at the receiver. How would this affect the constellation? If you know the answer to that question, then think of a carrier frequency offset as a continuing change in phase offset. So how would this affect the constellation?


Window functions and indexes larger than N/2

The window is defined as $-N/2 \le n \le N/2$ i.e. symmetric about $n=0$.

So, for a 1024 window size, $n=-512$ to $n=+511$ should give you the appropriate coefficients.


If you're familiar with Fourier transforms, I think the bridge between the Fourier worlds and the wavelet worlds is the Gabor transform (a Gaussian-windowed STFT) and the complex Morlet wavelet transform.  This is historically how they developed, too.  They are basically the same thing, breaking down a signal into "blips" of complex sinusoids:

But the time-frequency space occupied by the blips are spaced differently:

The wavelet version has more frequency resolution at low frequencies and more time resolution at high frequencies, which is usually a good tradeoff (similar to the way the human ear works).

The Morlet is a continuous wavelet, though, so there is overlap/redundancy in the representation, a discrete version is not a minimal representation of the signal, and does not meet the "admissibility condition", which apparently means it cannot be inverted perfectly back into a signal(?), and Parseval's theorem can't be used on it.  Modifying the wavelet so these things are possible results in other types of wavelets, and you can eventually work back to things like the Haar wavelet (I think).

Also see What's the difference between the Gabor-Morlet wavelet transform and the constant-Q transform?


Can I use SINAD for THD+N measurements? Are they inverse of each other?

So I guess I should have looked around a little better. I've found that SINAD and THD+N are reciprocals.

These references were useful:

https://www.ap.com/technical-library/more-about-thdn-and-thd/

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.476.1903&rep=rep1&type=pdf


Deconvolution of Synthetic 1D Signals - How To?

You cant't recover the original signal through deconvolution.

A Gaussian kernel is in essence a lowpass filter, i.e. it will remove information at higher frequencies from the signal. Once it's gone, it's gone and you can't recover it.

This problem shows up as "divide by zero" or "divide by a very small number", which then amplifies numerically noise of the original convolution to a massive amount.

In order to recover the original signal the system must have "sufficient signal to noise ratio in the bandwidth of interest". What exactly that is, depends heavily on your specific application.

You can de-convolve from a low-shelf or high-shelf filter, but not from a low-pass or anything with a stop band or lots of attenuation at specific frequencies.

You'd also lose some information at the beginning and end of the signal.


Splitting a Complex FIR Filter into Real Filters

I'd say there are two ways to look at this problem.

If you're receiving a conventional $M$-QAM signal with AWGN noise, then the I and Q streams are actually two independent, real, $\sqrt{M}$-PAM signals. You can filter and process them independently. The same is true over the wireless channel, as long as the fading is flat and you do the appropriate detection in the receiver.
If $s(t)=s_I(t)+js_Q(t)$ is a complex signal and $h(t)=h_I(t)+jh_Q(t)$ is a filter's complex impulse response, then

$$s(t)\star h(t)=s_I(t)\star h_I(t)-s_Q(t)\star h_Q(t)+j(s_I(t)\star  h_Q(t)+s_Q(t)\star h_I(t)),$$

where all the filters involved are real. This follows from the convolution's properties.


Finding number of alternations in parks mcclellan

This question is from Oppenheim, Schafer and Buck's textbook Discrete-time Signal Processing. There must be a typo either in the problem or in the solution manual. But the only important thing is to be able to find the number of alternations, and to be able to check whether a given filter can be optimal for a given filter length.
For a type I (even symmetry, odd filter length) linear phase low pass filter, the degree of the corresponding polynomial is $L=(N-1)/2$, where $N$ is the filter length. Note that the filter order $M$ is equal to $N-1$.
According to the alternation theorem, a Chebyshev optimal type I low pass filter must have $L+2$ or $L+3$ alternations. Alternations are points of maximum error, with the only restriction that adjacent alternations must have errors of different signs (that's why we talk about alternations after all). Note that the band edges also count as alternations because at these points the approximation error attains its maximum.
Looking at the figures we see that the first filter has $3$ alternations in the passband (one at $\omega=0$, one at the maximum, and one at the passband edge), and $3$ alternations in the stopband (one at the band edge, one at the minimum, and one at $\omega=\pi$), which makes for a total of $6$ alternations. In the same way you can see that the second filter has $7$ alternations.
If the filter length is $N=11$, we get $L=(N-1)/2=5$, and the consequently, the optimal filter must have $L+2=7$ or $L+3=8$ alternations. In this case, only the second filter satisfies the alternation theorem.
On the other hand, for $N=9$ we have $L=4$, and both filter satisfy the alternation theorem because $L+2=6$ and $L+3=7$.


How to adapt 4th order IIR filter coefficients from matlab to cmsis dsp biquad struct

% design the filter as poles, zeroes and gain
[z,p,k] = butter(2,[445 561]*2/48000);
% turn into biquads
sos = zp2sos(z,p,k);
% eliminate a0
coeffs = sos(:,[1 2 3 5 6]);
% make a linear array
coeffs = [coeffs(1,:) coeffs(2,:)];


MIDI out from freeware sequencer to PD or Max/MSP

Max/MSP has the midiparse object which is incredibly useful for handling midi data. Midi Ox is a program that will help you with routing occasionally. As far as freeware, you get what you pay for. Hydrogen is an open source a drum machine that will send midi out. You'll just to configure Max or Pd to receive the channel you're sending out of your sequencer.


Power spectral density vs. FFT bin magnitude

The fast Fourier transform ($\textrm{FFT}$) algorithms are fast algorithms for computing the discrete Fourier transform ($\textrm{DFT}$). This is achieved by successive decomposition of the $N$-point $\textrm{DFT}$ into smaller-block $\textrm{DFT}$, and taking advantage of periodicity and symmetry.

Now, the $N$-point $\textrm{DFT}$ of a sequence $\{x[0], x[1],\cdots, x[N-1]\}$ is:
\begin{equation}
X(f_k) = \displaystyle \sum_{n = 0}^{N - 1} x[n]\exp\left(-j2\pi f_kn\right) \tag{1}
\end{equation}
For $f_k = k/N$ and $k = 0, 1, \cdots, N - 1$. And the $\textrm{FFT}$ magnitude at bin $k$ is the $\textrm{DFT}$magnitude at bin $k$. For a given $N$ that is:

\begin{equation}
\left|X(f_k)\right| = \left|X\left(\frac{k}{N}\right)\right| = \left|X(k)\right| = \displaystyle \left|\sum_{n = 0}^{N - 1} x[n]\exp\left(-j2\pi nk/N\right)\right| \tag{2}
\end{equation}

You talk about power spectral estimation when the signals being analyzed are characterized as random processes. With random fluctuations in such signals, statistical characteristics and average characteristics are normally adopted. For a wide sense stationary $(\textrm{WSS})$ discrete random process, the PSD is defined as:

\begin{equation}
P(f) = \displaystyle \sum_{m = -\infty}^{\infty} r_{xx}[m]\exp\left(-j2\pi fm\right) \tag{3}
\end{equation}

For reasons you can find in this answer, you see that the squared magnitude of the signal's $\textrm{DFT}$ is taken as the estimate of the PSD in most practical situations. One form, among other variations/methods, is:

\begin{equation}
P(f_k) = \frac{1}{N} \displaystyle \left| \sum_{n = 0}^{N - 1} x[n]\exp\left(-j2\pi f_kn\right) \right|^2 \tag{4}
\end{equation}

  What's the difference between these?

Comparing $(2)$ and $(4)$, you have:

\begin{equation}
P(f_k) = \frac{1}{N} \left| X(f_k)\right|^2
\end{equation}

From the bin number $k$ to frequency in $\textrm{Hz}$, $F = \frac{F_s}{N}k$

For more reading on PSD estimation check this question, that question, and this question.

EDIT:

The power spectral density, $\textrm{PSD}$, describes how the power of your signal is distributed over frequency whilst the $\textrm{DFT}$ shows the spectral content of your signal, the amplitude and phase of harmonics in your signal. You pick one or the other depending on what you want to observe/analyze. And no they're not the same as you can see from the equations above and links given. Their spectra are generally not the same. One is estimated as the squared magnitude of the other.


sampling a continuous time signal using MATLAB?

This information can be found other places as well but I will step through it here using MATLAB.

You say your sampling rate is $100$ Hz and the signal is $10$ Hz. Lets define those along with the sampling period ($1/100$ seconds).

sampleRate = 100;
samplePeriod = 1/sampleRate;
signalFreq = 10;

Now you want $2$ seconds of this, so define your sampled time vector.

nT = 0:samplePeriod:2;

Now we create the signal, and plot the result.

signal = cos(2*pi*signalFreq*nT);

figure
stem(nT, signal)
grid on
xlabel('Time (sec)')
ylabel('Amplitude')

If you only want to plot $0 \leq n \leq 19$, then you may instead set nT = [0:19]*samplePeriod;.


Eliminating inertia effect from the acceleration signal

It sounds like you want to remove the frequency component at about 2.5Hz in the plot, but leave the rest of the plot as untouched as possible.

The usual first approach to that problem is to just use a notch filter:

$$
H(z) = K \displaystyle\frac{1 - 2 \cos \omega_0 z^{-1} + z^{-2}}{1 - 2r \cos\omega_0 z^{-1} + r^2 z^{-2}}
$$

where $\omega_0$ is the (normalized) frequency you want to remove and $r \lt 1$.

There are more complex approaches that use, for example, Kalman filtering... but that might require more assumptions about the dynamics of your system that you are prepared to make given your "black box" comments.


Digital Phase Modulation baseband or passband?

Assuming $g(t)$ is a baseband signal, $s_m(t)$ is passband. You can see that it is by multiplying the two complex exponentials, expanding the result into sine and cosine terms, and then taking the real part.
An interesting thing about PSK signals is that they are an example of quadrature modulation, and in consequence, you can represent them either as real passband signals, or complex baseband signals. If you assume $s_m(t)$ is a PSK signal, and given that it is real, it must also be passband.


For ideal performance, you can DFT your input, pad out the spectrum with zeros, and inverse DFT the result.  This is not a particularly efficient approach, though.

Pretty much the same result is obtained with interpolation by zero-padding and then running that stream through a low pass filter.  An FIR sinc filter is preferred for this for two reasons--(1) it has nearly ideal performance and introduces negligible distortion into the baseband signal, and (2) when zero padding, most of your input samples are zero.  Since 0*x=0, you don't have to do all those multiplies; for your 6:1 example, you can design a 72-point FIR and still only do 12 multiplies per output sample.

Once the filtering is done, you can upconvert (multiply by the exponential).  This should be $e^{j\omega t}$ where $\omega=2\pi f$ and f is your RF carrier frequency.


Removing 'high-frequency' components from an array?

You mention convolution which implies that your arrays represent time domain signals. As this is digital signal processing, the time is discretized and every element of the array corresponds to a sample. Two consecutive samples have a fixed time distance, say $T$.

If you convolve the discrete time signal with the impulse response of a low-pass filter (LPF) the resulting ouput signal will be longer than the input signal, because in general, linear time-invariant (LTI) systems have a memory (length of impulse response). You cannot reduce the array length without distorting the output signal (unless you use block-wise processing with overlap-add/-save, but that's another story). So, no matter what type of filter the impulse response of your LTI system represents, you cannot just remove elements from the output array.

To find out which frequency components are attenuated by a system you should analyze its transfer function (sometimes also refered to as frequency response). The transfer function is obtained by transforming the impulse reponse into the frequency domain. A tool for doing that with DSP is the discrete Fourier transform (DFT) which almost always comes in form of its FFT implementation (Fast Fourier transform). In case that the system in question is a black box and its impulse response is unknown the frequency domain representation of input and output signal can be compared to find out what frequencies are attenuated or amplified.


Are there theories to manage states of signal volatility around a limit

In electronics, a Schmitt trigger converts an analog signal to a digital signal by using two different thresholds (limits as you call them), one for an output state change from 0 to 1 and another for a change from 1 to 0. The first is somewhat higher so that if the input fluctuates between those two thresholds, then there will be no change in the output. The general concept of reluctancy to change state is called hysteresis.

It can also help to lowpass filter the signal first, which is often used in debouncing a mechanical switch, to remove fluctuation very much like in your picture.


Charles E.Radke, "Necessary and sufficient conditions on conditional probabilities to maximize entropy", Information and Control 9, 279-284 (1966), https://doi.org/10.1016/S0019-9958(66)90165-3 gives a proof of Shannon's Theorem 8:

  A necessary and sufficient condition such that $\eta_c$ is maximized and equal to $1$ is that the set of symbol lengths ${l_{ij}}^v$ must satisfy:
  $${p_{ij}}^v = (\beta_j/\beta_i)\omega_r^{-{l_{ij}}^v},\tag{3}$$
  where $\beta_i$ and $\beta_j$ are respectively the $(k, i)$th and $(k, j)$th cofactor of the matrix $A$ for any $k$.

The matrix $A$ is defined earlier as:

  $$A \triangleq \left[\textstyle\sum\omega^{-{l_{ij}}^v} - \delta_{ij}\right];$$

which is something you can recognize in Shannon's Theorem 1.


When is the sum of two (parallel) minimum-phase filters also minimum-phase?

Just looking at it I can see it's related to control theory.  If you make a contrived open-loop gain $$H(z) = \frac{A(z)D(z)}{B(z)C(z)}$$ and then wrap it with unity-gain feedback, you get $$G(z) = \frac{A(z)B(z)C(z)D(z)}{A(z)D(z) + B(z)C(z)}$$.
So any stability test that takes $H(z)$ as an open-loop gain will equally be a test to see if your sum is a minimum-phase system.  If you have some parameter variation then you could use root-locus plots, Bode plot or Nyquist plot analysis, or robust control methods to determine the amount of variation that could be tolerated while keeping the system minimum phase.
As an example, if you let your transfer function be $$k\frac{A(z)}{B(z)} + \frac{C(z)}{D(z)}$$ then you could easily cast this as a root-locus problem, and with a bit of effort you could put it into Bode plot form or feed it into a robust control solution that's designed to work in the $z$ domain.
There's probably a more direct path from here to there than going through control theory, but my hammer says "control theory" on it, and by gosh your problem is looking a lot like a nail to me!


how can we align wto signals with a time lag in between?

I thought you might have to smooth, sort or create an envelope of your signal, which might help. But then i found this. This should help you.

Or if you wish to use MATLAB, you can use this command.


Why integrate over $2\pi$ in inverse DTFT?

There are two reasons:

As you say, the sequence is periodic so all of the information is contained in the $-\pi < w <\pi$ region.  Thus, integrating over that region tells you everything that there is to be learned from the inverse transform.
Integrating from $-\infty < w< \infty$ would only give answers of $0$, $\infty$, or $-\infty$.  The reason for that is that integrating over $-\pi < w <\pi$ will generally give a finite answer which can be something negative (e.g. -4), 0, or something positive (e.g. 2.5).  If you integrate over $-\infty < w< \infty$ each $2\pi$ region will give the exact same answer, so the negative results will accumulate to $-\infty$, the zeros will sum to $0$, and the positive results will sum to $\infty$.  (Note that the results are often complex, but that doesn't change the primary point.  It just increases the dimensionality of the infinities.)


From MFCC to Machine learning. What are the steps?

There are two general things you can do with a machine learning algorithm:
1. Regression-- take input data and return an output numerical value. This is also known as curve-fitting, parameter estimation, etc.
2. Classification-- take input data and return a decision from a discrete set. In traditional statistical terminology, this is similar to hypothesis testing.

(I'm assuming that you are interested in classification.)

If you are doing supervised learning, you will need some training data, i.e. sets of MFCC coefficients that have been labeled with a priori knowledge. E.g. if you are identifying American English accents, you will need some audio data from New Englanders, Southerners, Midwesterners, etc., from which you extract the MFCC coefficients and label them manually. These coeffcients are known as features and the algorithm that distills down the high-dimensional dataset (i.e. the audio clip) down to a few coefficients (i.e. the MFCC coefficients) is known as feature extraction.

Once you have your training data and labels, you can apply a simple algorithm like KNN quite easily using libraries such as the Python sklearn library.


Some more disadvantages:

Channel estimation using pilots is harder, because you do not have orthogonal data on each frequency bin, compared to OFDM. That's why e.g. LTE uplink resorts to a midamble for channel estimation in uplink
Application of non-linear MIMO detection algorithms is harder (i.e. more complex), because you essentially have coupled MxM elements from each transmit antenna, compared to 1x1 elements when using OFDM.
In high SNR regions I think SC-FDMA does not perform as good as OFDM when considering a coded transmission with linear detectors.

Regarding the data loss: Even if $M<N$, there is no real loss, because the used bandwidth of your signal is only M frequency bins (ignoring any guard carriers). So, any other system can use the remaining N-M frequency bins.


How to make color balance of photoshop using opencv

OpenCV don't have any apis to separate data into lowlight,midtones and highlights.But you can categorize pixels in image into shadows,mid-tones and highlights by setting some thresholds.As they depend on the luminance first convert your image into YUV. If its of 8 bit depth then the pixels with luminance less than 70 are lowlight, pixels ranging from 70 -180 are midtones and 180-255 are highlight.

When cyan-red is -20,decreasing the red by -20 is just not enough. Because it will also affect the luminance of the image. So instead of decreasing red by 20 decrease it by 10 and increase green and blue by 5. Decreasing the red is same as increasing blue and green,while it also maintains the luminance. You can do the same for Magenta-Green and Yellow-Blue. Hope this helps you.


inter symbol interference vs digital modulation

ISI is an (typically undersired) effect of the channel – that can be a cable, an optical fiber, free space or the air you use to transport your wave.

The different kinds of constellations are a method of putting data in a waveform.

So, they are connected: The modulations define which symbols your system uses. ISI is what happens to the symbols when you rapidly send symbols one after the other.

Hence, your receiver needs to deal with ISI, and it needs to understand how to get the data back from the symbols it received.


Frequency range for set of all discrete-time sinusoids

assuming $\omega$ is real, if the form of $x[n]$ is how you've defined it above, your thinking is correct.  doesn't matter if it's $+\omega$ or $-\omega$, the sinusoid behaves exactly the same.  if instead, the sinusoid had form

$$ x[n] = e^{j \omega n}  \ ,$$

then your teacher would be correct.


The CWT is not, per se, superior to STFT. Due to its variations in scale, it can be better, for instance, when:

you address natural signals where transients are shorter than more stationary parts,
you do not know the appropriate scale of observation.

which is neither the case from your signal: precise frequency on the same support.

On the one hand, you can work with a continuous wavelet with more oscillations, and better refine the scales and voices.

On the other hand, there are so-called STFT with windows whose width vary with frequency.

Add-ons like hybrid time-scale/time-frequency transformations, reassignement, etc. may be used, but choosing the appropriate tool requires to refine your goal.


Reconstructing/interpolating small regions of a bandlimited signal by taking the fewest possible samples

I have some very short signals in the range of 8 to 16 samples. These represent a bandlimited signal, sampled at or slightly above the Nyquist rate.

Nope. A signal can't be limited in time and in frequency at the same time. If it's very short, than chances are the bandwidth is a lot higher than you think it is and that you've already picked up some significant aliasing in the sampling.

  I've experimented with sinc interpolation with a few different windowing methods, and also tried numerically optimising a kernel to achieve the interpolation, without success.

Can you elaborate what your problem is? If it's sampled correctly, sinc interpolation typically does a good job, although it does destroy any causality and will extent the signal in time. But that goes back to my first point: a band-limited signal can not be limited in time as well, so it can't be causal either.


Writing the optical flow vector as $(u,v)$, then the optical flow equation only states that light is preserved.
Basically, it states that variations in the observed light patterns (measured by the image derivative with respect to time) are explained by the motion of objects (smaller light patterns inside the image).
Object motion is measured through spatial variations of the said light patterns:

$$\frac{\partial I}{\partial t} + u\cdot\frac{\partial I}{\partial x} + v\cdot\frac{\partial I}{\partial y} = 0$$

This problem is ill-posed for any image, be it synthetic or real.
This is easy to see since 1 pixel provides 1 measurement $\frac{\partial I}{\partial t}$ only, while at the same time there are two unknowns $\left(u(x,y), v(x,y)\right)$.

As you said, since this equation involves image spatial gradients $\nabla I = (\frac{\partial I}{\partial x}, \frac{\partial I}{\partial y})$ it can be solved for edges (by tracking them).
Other ways to overcome this issue:

assuming that the optical flow is constant over a square of several pixels and then solving the resulting system. Typically, using 3x3 squares yields 9 equations and a least squares solution to the optical flow can be computed;
more involved approaches rely on regularized problems. Regularized approaches look for a flow that will fit the observed data while having good properties such as smoothness, exhibiting thin discontinuities...


Direction of sound using microphone array

If you connect four discrete USB microphones a challenge will be to synchronize them. Even when you start them simultaneously (from a programming point of view) you cannot be sure that they all start the audio acquisition at exactly the same time so there will be an inherent unknown time delay between the microphones. This can be solved by having a reference signal coming from a known direction. Then, by comparing the expected delays with the actual delays, synchronization can be achieved. With your setup I would expect the synchronization delays to be slightly different every time the system is started so that calibration needs to be done every time.


Confusion related to PN sequence : Terminology and application in statistical signal processing

If the output $x[\cdot]$ of a discrete-time linear system is related to the input via
$$x[n] = \sum_{i=0}^m h[i]u[n-i] \tag{1}$$
where $u[\cdot]$ is a periodic sequence with period $N \gg m$, then $x[\cdot]$ is also a sequence with period $N$ and
the
periodic cross-correlation function $R_{x,u}[\cdot]$ between $x[\cdot]$
and $u[\cdot]$ is the periodic convolution between the sequence $h[\cdot]$
and the periodic autocorrelation function $R_{u,u}[\cdot]$. That is,
$$R_{x,u} = x \star u = (h * u) \star u = h * (u \star u) = h* R_{u,u}
\tag{2}$$
  If
$u[\cdot]$ is a binary pseudonoise (PN) a.k.a. pseudorandom noise (PRN) sequence that takes on values $\pm 1$, then
$$R_{u,u}[n] = \begin{cases}N, & n \equiv 0 \bmod N,\\
-1, &\text{otherwise}\end{cases}\tag{3}$$
which is essentially a periodic digital impulse function,
 and by ignoring that $-1$ as inconsequential compared to $N$, we get from $(2)$ that
$$N h[n] = R_{x,u} \implies h[n] = N^{-1} R_{x,u}[n].\tag{4}$$
If we include additive white Gaussian noise $w[n]$
of variance $\sigma^2$ on the right side of $(1)$,
then $(2)$ becomes
$$R_{x+w,u} = (x + w)\star u = (h * u + w) \star u = h * (u \star u)
+ w \star u= h* R_{u,u} + \hat{w}
\tag{5}$$
where $\hat{w} \sim \mathcal N(0, N\sigma^2)$. Thus, the measured
cross-correlation function is now a Gaussian random variable
with mean $Nh[n]$ and variance $N\sigma^2$.  The estimator for
$h[n]$ is thus just $N^{-1}$ times the measured value of the
cross-correlation function, and its variance is thus
$\displaystyle \frac{\sigma^2}{N}$, This matches up perfectly with
Kay's result
$$\operatorname{var}(\hat{h}[n]) = \frac{\sigma_{w}^2}{Nr_{u,u}[0]}$$
since we have assumed that $u[n] \in \{+1, -1\}$ and hence
$r_{u,u}[0] = 1$: all the $u[n]$, regarded as random instead of
pseudorandom, have variance $1$.


What are $H(i\omega)$, $H(s)$ and $H(z)$ called?

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.


I will try to address @Hilmar answer and answer the general question.

Hilmar Answer

  Deconvolution is in the general case not possible, so it needs to be
  approximated with application specific constraints and requirements.

  Let's look at a simple 1-dimensional example that illustrates the
  problem. Assume you have an impulse response like your kernel, i.e. h
  = [1 1 1]. Then let's look at the output of a signal (which you can think of a line of pixels), x = [2 -1 -1 2 -1 -1 2 -1 -1 2 -1 -1 2 -1
  -1]. The result of the convolution is

    2   1   0   0   0   0   0   0   0   0   0   0   0   0   0  -2  -1

The case above is simply wrong. Specifically if one looks on my comment:

  The result being zero doesn't mean data is lost. It has the same information as any other number.

Let's do the following example:

% Input Signal
vX = [2; -1; -1; 2; -1; -1; 2; -1; -1; 2; -1; -1; 2; -1; -1];
% Convolution Kernel
vH = [1; 1; 1];

% Convolution Result
vC = conv(vX, vH, 'full');

% Convolution Operator - Matrix Form
mH = CreateConvMtx1D(vH, length(vX), 1);

% Deconvolution
vXR = mH \ vC;

Where CreateConvMtx1D() is given in my GitHub Code.

The result is given by:

vXR.'

ans =

    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000    2.0000   -1.0000   -1.0000

Which is a perfect reconstruction of the input signal.

General Case

Indeed the properties of Deconvolution in 1D and 2D are similar.
The ability to reconstruct the signal are basically a function of 2 parameters:

The SNR of the data (Added noise).
The Condition Number of the Convolution Operator of the problem.

Let us assume there is no added noise in the problem.
Then we're only limited by the Condition Number.
In theory, for any bounded Condition Number one would be able to perfectly reconstruct the signal. Yet since we use finite and quantized representation of numbers (Let's say Float64) we are limited with the condition number we're able to tackle.

For more details see - 1D Deconvolution with Gaussian Kernel (MATLAB).


It is called a cardioid pattern. (cardioid because it looks vaguely heart-shaped)

You've shown us a polar frequency magnitude plot with the frequency given in radians. You'll notice by going counter-clockwise around the unit circle that low frequencies are attenuated while the higher frequencies are passed through.


Efficient Magnitude Comparison for Complex Numbers

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.


filtering -> thresholding -> edge detection

apply low pass filter:

averageFilterCoefficient = 0.01
for each sample of x:
  y[i] = y[i-1]*(1-averageFilterCoefficient) + x[i]*averageFilterCoefficient

this simple code should remove all high frequencies and produce more smooth signal. Now apply thresholding:

threshold = 100
for each sample of y:
  if y[i] > threshold then
    z[i] = 1
  else
    z[i] = 0

Now you have nice logical signal, lets find rising edges:

for each sample of z:
  if z[i-1] = 0 and z[i] = 1 then
    print i

You just need to find the right values for exponential average filter
coefficient and thresholding level.