Understanding unexpected i/q calculation on a downsampled signal

Because the decimate function does not try to keep power constant.

Remember that decimation is not inherently an LTI operation; so, probably, the authors of that function simply did not care too much for keeping the amplitude of a passband signal constant.

Note that a decimation step typically includes filtering to $\frac1{\text{ds_factor}}$ of the original bandwidth to avoid aliases. That filter might or might not have a passband gain $\le 1$.

A general remark:

  For each pair in (for i,q in zip(x_real, x_imag)), I calculate h=math.sqrt(i**2 + q**2). I don't expect h to ever be larger than sqrt(2), because the real and imaginary parts are normalized to [-1, 1].

In fact, if coming right out of a receiver with good IQ balance, don't expect $h\ge1$, even!

Remember that I and Q aren't independent; they are the inphase and quadrature component of a signal whose max power, no matter which phase, is limited! Thus, as long as you don't risk clipping in your receiver, $i^2 \le {1-q^2}\iff q^2\le 1-i^2$!

If your I and Q would be independently limited to an amplitude of 1, your $\sqrt 2$ would be right. They are not – imagine, for example, the passband signal that gives you $i\equiv1$: it's the cosine with the carrier frequency $1\cdot\cos(2\pi f_c t + 0)$ (note the $0$ phase). Now, any cosine with larger power would clip, right? Obviously, cosines of the same amplitude would also be at the extreme power limit, no matter whether they have zero phase. The moment the phase is not zero, you'll see some of the cosine's energy appear on Q, and missing on I. However, it's the same power.


How do I model a system to use with a Kalman Filter?

The following article explains the questions you have on modeling the system mathematically, in a very lucid way with examples and equations.

http://www.anuncommonlab.com/articles/how-kalman-filters-work/index.html


One sided frequency spectrum (Matlab vs. Origin)

The indexing is a matter of convention.  "Natural" is zero centered, "Computer implementation" is zero based.

We are having a big discussion right now over several questions on why your $c_{-5}$ should actually be considered as $(c_{-5} + $c_{5})/2$.  You'll find part of the discussion here Convergence of periodic sinc interpolation and the links to the rest are embedded within.

The "fluffy cloud" drawings show what happens when you consider the upper half (from a zero based perspective) as positive frequencies.  How to get Fourier coefficients to draw any shape using DFT?

If I have a circle, and I check your location every minute and it seems you have moved a quarter circle counter-clockwise each time, how do I know you haven't actually moved three quarters of the circle clockwise?

Answer:  You don't.  The convention decides which interpretation to use.


Why does this Fourier series give lower amplitude for max value of signal?

The formula for the Fourier series is

$$x(t)=\sum_{k=-\infty}^{\infty}c_ke^{jk\omega_0 t}\tag{1}$$

You only use the non-negative indices $k$, but you must also use the negative ones. Since your signal $x(t)$ is real-valued, you can in fact compute everything using only the positive indices, but then you must take two times the real part:

$$x(t)=c_0+2\Re\left\{\sum_{k=1}^{\infty}c_ke^{jk\omega_0 t}\right\}\tag{2}$$

Note that you must add the DC term $c_0$ separately.

This is the corresponding Matlab code:

t=-10:0.001:10;
x_t = 0;
T = 10;
w0 = 2*pi/T;
numOfTerms = 100;

for k=1:numOfTerms-1
   a_k = (1/(w0*i*k*T))*(  exp(-1*i*k*w0)-exp(-3*i*k*w0)+ exp(-7*i*k*w0)- exp(-9*i*k*w0));
   x_t = x_t + a_k * exp(i*k*w0*t);
end

x_t = 2*real(x_t) + 2/5;

plot(t,x_t);


Physical interpretation of 4th-order correlations

If you are willing to relax the WSS assumption, and consider cyclostationary signals, you might check out my theory of pure and impure sine-wave components of higher-order probabilistic parameters:

On pure and impure sine-waves, higher-order moments, and cyclic cumulants

and

Oddities concerning pure nth-order sine waves

I've thought about your question a lot, more so in the past than lately, and all I've been able to come up with is that higher-order moments contain characteristic sine-wave components (characteristic of the process, the order, and the number of conjugated factors), and higher-order cumulants characterize what is new in the moment that cannot be accounted for by products of lower-order moments. Unlike the easy physical interpretation of autocorrelation as power (think of the zero-lag autocorrelation), I don't have a simple physical unit to assign here, but maybe my hint will lead you to a more satisfying answer.


Proving the upper bound of cross correlation

The  document cited by the OP actually states that
$$|R_{X,Y}(\tau)| \leq \sqrt{R_X(0)R_Y(0)}\tag{1}$$
which implies what the OP wants to prove, namely that
$R_{X,Y}(\tau) \leq \sqrt{R_X(0)R_Y(0)}$ but also tells us
that $R_{X,Y}(\tau) \geq -\sqrt{R_X(0)R_Y(0)}$.
Recall that for jointly wide-sense stationary processes (which are, of course, individually wide-sense stationary too and thus have constant (possibly nonzero) means $\mu_X$ and $\mu_Y$), the definitions of the functions in $(1)$ are
$$ R_X(\tau)\stackrel{\Delta}{=} E[X(t)X(t+\tau)],\quad R_Y(\tau)\stackrel{\Delta}{=} E[Y(t)Y(t+\tau)],\quad R_{X,Y}(\tau)\stackrel{\Delta}{=} E[X(t)Y(t+\tau)]$$
where, because of the wide-sense stationarity assumption, the functions depend only on $\tau$, and not upon $t$. Note also that the means $\mu_X$ and $\mu_Y$ of the random processes are not involved at all, much to the sorrow of those who prefer to demean their processes before calculating correlations and thus end up with the autocovariance  and cross-covariance functions instead of autocorrelation functions and cross-correlation functions as they are defined and understood in the communications and signal processing literature.
Now consider the quantity $$E\left[(\alpha X(t)\pm\beta Y(t+\tau))^2\right]
= \alpha^2 R_X(0)+\beta^2 R_Y(0) \pm 2\alpha\beta R_{X,Y}(\tau)$$ which, being the expected value of a nonnegative random variable, must be nonnegative for all choices of real numbers $\alpha$ and $\beta$. Thus,
$$\pm 2\alpha\beta R_{X,Y}(\tau) \leq \alpha^2 R_X(0)+\beta^2 R_Y(0)\quad \forall ~\alpha, \beta$$ or equivalently, for every possible choice of positive real numbers $\alpha$ and $\beta$, it must be that
$$2\alpha\beta |R_{X,Y}(\tau)| \leq \alpha^2 R_X(0)+\beta^2 R_Y(0). \tag{2}$$
Choosing $$\alpha^2 = \frac 12 \sqrt{\frac{R_Y(0)}{R_X(0)}},
\quad \beta^2 = \frac 12 \sqrt{\frac{R_X(0)}{R_Y(0)}}$$
so that $\alpha^2\beta^2=\frac 14$, we see that Eq.$(2)$ reduces to
$$|R_{X.Y}(\tau)| \leq \frac 12 \sqrt{R_X(0)R_Y(0)} + \frac 12 \sqrt{R_X(0)R_Y(0)}$$
which is $(1)$.


Calculating spectral coherence coefficient (scc) - implementation help

I am just working through this information myself, but I believe I understand your questions.

  If I put $S_x^0(f-\alpha/2)$ in the denominator then won't the denominator essentially become $S_x^0(f)$, since $\alpha = 0$?

It will become $|S_x^0(f)|$ at $\alpha = 0$, but at other values of $\alpha$, you would need the conjugate multiply at offsets of $f$.

  Will $C_x^\alpha$ be of the same size as $S^\alpha_x$?

Yes, one would think. The coherence SOF function, is just supposed to be some normalized version of the SCF.

I am struggling with the FAM, it appears to output data on a tilted spectral plane, where $\alpha=0$ lies on the matrix diagonal.


$\mathrm{FFT}$ generates non-linear noise

First of all, look at the power of what you call background noise:

It's between -150 dB and -120 dB below your peak's power. Which means that the magnitude of these FFT points is somewhere in $[10^{-7.5};10^{-6}]$ of the main peak magnitude.

That screams "numerical inaccuracy" all over the place:

Assuming your FFT was implemented using 32 bit IEEE754 floating point numbers, you'd have to realize that for each point of the FFT output, there's at least an FFT-length amount of additions taking place with a very finite numerical precision.

Now your FFT has an enormous length; $1323000\approx 1.3\cdot 10^6$ additions bring the most stable algorithm to its knees. For a bit of details on error analysis on commonly used FFT algorithms, refer to [1].

Have a look at the DFT matrix to understand the concept of how these values come to be: For the frequency of your spectral peak, signal-period ($\frac{1}{f_\text{signal}}$) spaced samples are added up; since they always have the same value (definition of "period"), the result is large; for all the others, the should cancel each other out. For example, the bin at $2\cdot f_\text{signal}$ always adds up the same absolute value, but with alternating sign. For other bins, you just take values from the sine and add them up. What happens now is that due to limited floating point math, when you add up a large and a small number, and then subtract the large number again, you won't get the same small number you had – accuracy reduces with magnitude!

Now, the FFT is an algorithm that is pretty resilient with respect to that (because it doesn't try to add up all FFT-length values at once), but there's only so much stability that an algorithm can achieve with a given machine $\epsilon$.

  I would have expected a more even noise, or rather a linear shaped noise.

Why that? If anything, you're observing a process where the imperfections of the FFT are correlated with the input signal. Something like that won't have "white" PSD! Remember, the PSD of a process is the Fourier transform of its autocorrelation; in your case, it seems there's more autocorrelation at higher frequencies.

Do your understanding a favor and add a bit of as-good-as-feasible white, gaussian noise to the input, maybe with a variance of $N=10^{-3}$ of the sine's amplitude, and do your mag squared FFT again. You'll see the noise floor where you expect it, and it will be flat, because the (dominant) source of noise is actually white, and not strongly correlated.

Also, your $1.3\cdot10^6$ FFT really isn't all that helpful, unless you actually need a frequency resolution of $\frac{f_\text{sample}}{N}=\frac{4.41\cdot 10^4}{1.323\cdot10^6}=33\,\text{mHz}$ (which I dare to doubt). Try with shorter pieces of your signal!

[1] Tasche, M., and H. Zeuner, Roundoff error analysis for fast trigonometric transforms, in Handbook of Analytic-Computational Methods
in Applied Mathematics, G. Anastassiou (ed.), Chapman & Hall/CRC, Boca Rota, 2000, 357–406.


What Does 'Zero Order Hold' and 'First Order Hold' Mean?

We need to assume the reader knows some basic stuff to answer that.
Let's give it a try.

Lets understand the sentence - Zero / First Order Hold.
We have the Zero / First Order and the Hold.

Zero / First Order hold means the order of the Taylor Series of the function we use to interpolate. In other words, the degree of the Polynomial we can write the function with.
Zero Order means the function is constant, we interpolate the same value in the missing parts.
First order means we can use linear function to interpolate (Line with a slope).

Hold means we hold the parameters to be the same until the next sample.


This derivation is a tricky one. The approach suggested before has a flaw. Let me demonstrate this first; then I will give the correct solution.
We wish to relate the $\mathcal{Z}$-transform of the downsampled signal, $Y_D(z) = \mathcal{Z}\{x[Mn]\}$, to the $\mathcal{Z}$-transform of the original signal $X(z) = \mathcal{Z}\{x[n]\}$.
The wrong way
One could think of simply plugging the expression for the downsampled signal into the expression of the $\mathcal{Z}$-transform:
$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{x[Mn] z^{-n}}$
A change of variable $n'=Mn$ seems obvious:
$Y_D(z) = \sum\limits_{n' \in M\mathbb{Z}}{x[n'] z^{-n'/M}}$
However, it is important to realize that even though the new summation index $n'$ still runs from $-\infty$ to $\infty$, the sum is now over 1 out of M integer numbers. In other words,
$n' \in M\mathbb{Z}=\{...,-2M,-M,0,M,2M,...\}$,
while the definition of the $\mathcal{Z}$-transform requires
$n \in \{...,-2,-1,0,1,2,...\}$.
Since this is no longer a $\mathcal{Z}$-transform, we cannot write:
$Y_D(z) = X(z^{1/M})$
The right way
Let us first define a 'helper' impulse train signal $t_M[n]$ as:
$\begin{align*}
t_M[n]
& = \sum\limits_{k=-\infty}^{+\infty}{\delta[n-kM]} \\
& =
\left\{
     \begin{array}{lr}
       1 & : n \in M\mathbb{Z}\\
       0 & : n \notin M\mathbb{Z}
     \end{array}
   \right.
\end{align*}
$
This function is $1$ at one out of every $M$ samples, and zero everywhere else.
Equivalently, the pulse train function can be written as:
$ t_M[n] = \frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}} $
Proof: We need to consider separately the cases $n \in M\mathbb{Z}$ and $n \notin M\mathbb{Z}$:
$$\begin{align*}
 t_M[n]
& = \frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}} \\
& =
\left\{
     \begin{array}{lr}
       \frac{1}{M} \sum\limits_{k=0}^{M-1}{1} && : n \in M\mathbb{Z}\\
       \frac{1}{M} \frac{1-e^{j2\pi n}}{1-e^{j2\pi n / M}} && : n \notin M\mathbb{Z}
     \end{array}
   \right. \\
& =
\left\{
     \begin{array}{lr}
       \frac{1}{M} M && : n \in M\mathbb{Z}\\
       \frac{1}{M} \frac{1-1}{1-e^{j2\pi n / M}} && : n \notin M\mathbb{Z}
     \end{array}
   \right. \\
& =
\left\{
     \begin{array}{lr}
       1 && : n \in M\mathbb{Z}\\
       0 && : n \notin M\mathbb{Z}
     \end{array}
   \right.
\end{align*}$$
In the case $n \notin M\mathbb{Z}$, we used the expression for the finite sum of a geometric series.
Now let's go back to our original problem of finding the $\mathcal{Z}$-transform of a downsampler:
$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{x[Mn] z^{-n}}$
We apply the substitution $n'=Mn$, keeping in mind that this makes the summation run only over integer multiples of M:
$Y_D(z) = \sum\limits_{n' \in M\mathbb{Z}}{x[n'] z^{-n'/M}}$
We can now use the above impulse train function to safely rewrite this as a summation over all $n \in \mathbb{Z}$:
$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{t_M[n] x[n] z^{-n/M}}$
Using the above formulation for the impulse train function as a finite sum of exponentials, we get:
$\begin{align*}
Y_D(z)
&= \sum\limits_{n=-\infty}^{+\infty}{(\frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}}) x[n] z^{-n/M}} \\
&= \frac{1}{M} \sum\limits_{k=0}^{M-1}{\sum\limits_{n=-\infty}^{+\infty}{e^{j2\pi k n / M} x[n] z^{-n/M}}} \\
&= \frac{1}{M} \sum\limits_{k=0}^{M-1}{\sum\limits_{n=-\infty}^{+\infty}{x[n] (e^{-j2\pi k / M} z^{1/M})^{-n}}}
\end{align*}$
The summation on the right is a summation over all integers, and is therefore a valid $\mathcal{Z}$-transform in terms of $z'=e^{-j2\pi k / M} z^{1/M}$. Therefore, we can write:
$
Y_D(z) = \frac{1}{M} \sum\limits_{k=0}^{M-1}{X(e^{-j2\pi k / M} z^{1/M})}
$
This is the formula for the $\mathcal{Z}$-transform of a downsampler.


Why/How a Gaussian Channel is both band-limited and Power limited?

What is usually meant by a Gaussian channel is a channel with additive (white) Gaussian noise (AWGN), which simply means that the received signal will equal the transmitted signal plus some unknown noise with a certain power and a Gaussian probability density function. A pure AWGN channel is neither band-limited nor power-limited. However, these restrictions are present in real-world systems, so it is usually assumed that the channel is a band-limited AWGN channel. Furthermore, it is normally assumed that the transmitter can transmit with some given maximum power. This latter fact will result in a bound on the achievable signal-to-noise ratio (SNR).

So, to answer your question, the finite band-width and the power limitation are realistic constraints, but they have nothing to do with the AWGN channel, which only specifies the type of noise (Gaussian) and how it influences the signal (by addition).


Why is random noise assumed to be normally distributed?

Starting at an even more basic level than the other (much smarter) answers, I'd like to pick up on this part of the question:

  This seems contradictory to me as on one side it is random then on the other side their distribution is considered normally distributed.

Perhaps the issue here is what ‘random’ means?

To be clear: ‘random’ and ‘normally-distributed’ do not contradict each other.  ‘Random’ simply means that we can't predict exactly what the next value will be.  But that doesn't mean we can't make probabilistic statements about it.

Consider two experiments:

If you throw a (fair) die, then it could show any number from 1 to 6.  We can't tell which number will come up, but we can say that all numbers are equally likely (i.e. the distribution is uniform).
If you throw two dice and take their sum, that can be any number from 2 to 12.  Again, the sum is still random — we can't predict what it will be — but we can say that those values are not equally likely.  (For example, 7 is six times more likely than 12.)  So in this case it has a non-uniform distribution.  (You can plot all the probabilities; they take on a peaked shape a bit like a normal distribution.)

So there's no contradiction: both cases are random and have a known distribution.

In fact, most things that are random tend to have a non-uniform distribution: electrical noise, weather, the wait for the next bus, voting patterns…  Being able to make general statements about them without being able to predict the exact values is one of the strengths of statistics.

(As for why you often end up with a normal distribution, that's a result of the Central Limit Theorem, which says that when you combine many independent random variables, the result tends towards a Gaussian (normal) distribution.  So you see that crop up a lot.)


No, recovering an "arbitrary" signal from its samples taken below its Nyquist rate (assuming it is bandlimited) is not possible in general.

Compressive sensing is based on the assumption that if you know a-priori that a signal will be sparse in some of its support domains, then those bands will not be carrying information about the signal and can be discarted in principle, therefore, reducing the pre-computed Nyquist sampling rate, which was computed based solely on the bandlimitedness (Shannon-Nyquist sampling theorem) of the signal.

The trick here is on the "a-priori" knowledge which also contains information, and this must be taken into account when computing the total information present in a signal. This a-priori knowledge can be the result of imposed structures on the signal generation mechanism, or well known models of the signal source.

If you can impose no constraints on the signal other than its bandlimit, you have to sample it above its Nyquist rate, so that exact reconstruction (using unimplementable ideal lowpass reconstruction filters) is possible at least in principle.

The topic of sampling is bound by information theory. There is no way to walk around it. Why would you ever sample a DC signal, or a single sine wave, or a periodic signal? It is based what is known and what is not.


What is the z-domain transfer function of variable delay?

A (discrete-time) system with a variable delay must be time-varying, and for this reason it cannot be described by a transfer function of the form

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}\tag{1}$$

where $h[n]$ is the impulse response of the system.
Only linear time-invariant (LTI) systems can be described by a transfer function as given in $(1)$.


Hamming window for filtering unwanted terms

Below is some code (in scilab which is similar to but not the same as matlab) that attempts to do what you're saying.

The output is the following graph which has:

The original signal
The signal multiplied by the window to leave just the left-hand peak.
The signal modulated with a complex exponential to shift the remaining peak to the centre (DC).

// 25709
N = 1024;
dc = 1 + rand(1,N,'norm');
fc = 0.388729874;
a = 1+0.5*rand(1,N,'norm');
t = [0:N-1];
phi = 0.187238947;
x = dc + a.*sin(2*%pi*fc*t + phi);
X= fft(x);

f = fftshift((fftshift([0:N-1]) - N/2)/N);

clf;
subplot(311);
plot(f, fftshift(abs(X)))

subplot(312);
wind = [zeros(1,N-500) ones(1,500) ];
plot(f,fftshift(abs(X.*wind)));

subplot(313);
y = ifft(X.*wind);
yy = y.*exp(%i*2*%pi*fc*t);
plot(f,fftshift(abs(fft(yy))))


In addition to the answer from Matt L., I want to emphasize the fact why the comment mentions the polynomial division. Consider to finite-length sequences $h[n]$ and $x[n]$, with length $N$ (one sequence can be zero-padded to length $N$).

Now, look at the convolution of both:

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

As an example, consider the two sequences $h[n]=[1, 2, 3]$ and $x[n]=[4, 6, 5]$. Their convolution is (in Matlab)

>> conv([1 2 3], [4 6 5])

ans =

     4    14    29    28    15

Now, consider the product of the two polynomials $h(x)=1+2x+3x^2$ and $x(x)=4+6x+5x^2$ (in Mathematica):

Expand[(1 + 2 x + 3 x^2) (4 + 6 x + 5 x^2)]
>>> 4+14 x+29 x^2+28 x^3+15 x^4

As you can see, the coefficients of the polynomial of the product are equal to the values of the convolution of both sequences.

So, now to get $x[n]$ from $y[n]$, you need to perform the polynomial division of $y(x)/h(x)$:

PolynomialQuotient[4 + 14 x + 29 x^2 + 28 x^3 + 15 x^4, (1 + 2*x + 3*x^2), x]
>>> 4+6 x+5 x^2

Hence, the coefficients of the polynomial yield the original sequence.

This technique can be extended to infinite sequences.

To get the connection to the Z-Transform here: The Z-Transform of the finite sequences $x[n]$ and $h[n]$ are given by

$$H(z)=1+2z^{-1}+3z^{-2}\\ X(z)=4+6z^{-1}+5z^{-2}$$

And the Z-Transform of the output $y[n]$ is given by

$$Y(z)=H(z)X(z).$$

Here, you again have the polynomial multiplication. Also, here it becomes even easier to spot, why it also works with infinite sequences: Then both $H(z)$ and $X(z)$ become fractions of polynomials, but the principle of the polynomial multipication remains.


How Do You Save / Convert an Image to Be in the Frequency Domain?

A trick would be to use the TIFF format which supports many channels.
For each channel in your original image create 2 channels, one for the magnitude and the other for the phase.

TIFF supports 32 bit hence data will be well preserved.


Information Theory - units of channel capacity

You are talking of capacity of two different types of channel.

In one case, the channel inputs and outputs are discrete in time.
At the $i$-th time instant, the received signal is $X_i+N_i$ where
$X_i$ is the received symbol of average energy $E$ and $N_i$ is the
noise (typically modeled as a sequence of i.i.d. $\mathcal N(0,\sigma^2)$
random variables).  The channel capacity of this discrete-time
Gaussian channel $~$is
$$C = \frac{1}{2}\log_2\left ( 1 + \frac{E}{\sigma^2}\right) ~\text{bits per channel use}$$
and so
bits per channel use is the natural metric.  If we are told how far
apart in time the discrete time instants are, e.g. one channel use per
microsecond, then a capacity $C$ bits per channel use can be stated
as bits per second, e.g. $C$ Mbps for our one microsecond
example.

In the second case, the inputs and outputs are continuous-time signals which occupy bandwidth and so the natural measure is bits per second per Hertz.  There
are more complications involved in making the transition from the continuous-time
channel to the discrete model, and in connecting the bandwidth $W$, the received
signal $P$ and the noise spectral density $N_0$ to $E$ and $\sigma^2$
(see here
for some details), but when all
this is done we get Shannon's celebrated formula
$$C = W\cdot \log_2\left ( 1 + \frac{P}{N_0W}\right) ~ \text{bits per second}$$
for the capacity of the additive white Gaussian noise (AWGN) channel of
bandwidth $W$.  This capacity can also be expressed as $C/W$ bits per
second per Hertz.


TLDR; The LFO should not be anti-aliased. The final signal, after applying the modulation, should be alias free. In practice, you either ignore this or approximate the latter.

Details:

The correct (although not necessarily always practical) way to think about this is to think of the whole signal processing/synthesis chain as a continuous time process, and then apply an anti-aliasing filter in the end, before discretizing.

So say, for example, that you're modulating the amplitude of an audio oscillator with a square LFO, i.e. $x(t) = (1 + 0.5 v_\mathrm{LFO}(t))v_\mathrm{osc}(t)$, where $v_\mathrm{LFO}$ is the output of the LFO with amplitude $1$ and $v_\mathrm{osc}$ is the output of the oscillator, and assume that the audio oscillator is perfectly anti-aliased. Now, the x(t) will have discontinuous edges at transitions of the LFO, which will generate an infinite number of frequency components, which will alias. However, if you had done the modulation in the continuous time domain and then sampled with appropriate anti-aliasing filters, the transition edges would be smoothed and the signal as a whole would be band-limited.

Another effect is that the spectral components of the LFO itself will be aliased, which may or may not have audible significance. This can be seen for example by the fact that if the period of the LFO is not an integer number of samples, the lengths of the sequences of $1$'s and $-1$'s will vary (by one sample) depending on which side of a sample instant the transition happens to land on on which cycle.

As a further extreme example, say you're modulating the amplitude of a constant DC offset. In this case the output will be that of the LFO itself, and if you bring it to audio frequencies you will certainly hear aliasing if the LFO isn't anti-aliased.

On the other hand, it is not correct to apply an anti-aliased LFO to an anti-aliased oscillator, since modulation is a non-linear effect: ideally you should apply a continuous time LFO to a continuous time oscillator, and anti-alias (i.e. sample with bandlimiting) the result. However, smoothing the LFO transitions a little may restrict the bandwidth of the frequencies produced by the modulation, thus partially remedying the extra aliasing. You should not use a sharp anti-alias filter though, since its step response ripple could cause unwanted modulation, but rather a filter which has a good phase response (simple one-pole or a higher order Bessel).

As far as I know, there is no general method to remove aliasing in a system with modulation, and if there were, it would certainly depend on the specific modulation and what is being modulated (frequency vs. amplitude, oscillator vs. filter, and so on). However, oversampling is a straightforward if computationally costly method which can be used to reduce all aliasing in the audio band, including that caused by modulation.

Finally, as you say, if the LFO is well below audio frequency, it is quite unlikely that any of this will be audible. However, you should be careful in that also modulation depth will have and effect: if you're for example modulating the frequency of an oscillator, even with a low frequency, but such a large modulation depth that the change in the frequency is appreciable over a period of the oscillator, you may get audible effects (at least in principle, I have not made experiments nor do I have a reference to cite).

Addendum:

Why do we care about the continuous time model on a discrete system, when working with audio?

(I'm of course assuming here you mean by plug-in a VST or some other audio processing plugin).

Consider what is your end goal here: to produce a sound to be heard by a person. That is inevitably a continuous time (CT) process. What you have control over is the stream of samples sent to the DAC, which is where the signal is transformed from discrete to continuous time.

In other words, you are actually in the business of producing continuous time signals, and the discrete time (DT) representation is simply a tool to reach that end goal. This is what I mean by talking about discontinuities: they are discontinuities in the imaginary CT signal, samples of which form your DT signal.

The easiest way then to reason about the output CT signal implied by your DT signal is to form the DT signal by properly sampling a mathematically formed CT signal. For example with a sawtooth oscillator, you take an ideal mathematical ramp, you bandlimit it and sample it (by whatever method you like, such as BLIP, BLEP, oversampling etc), and then the sampling theorem will tell you that the output of the DAC will be that same bandlimited sawtooth. If you had instead thought of the sawtooth as a DT process and simply formed a naive discrete ramp, you would have been gravely disappointed, as you would have gotten the nasty sound of an aliased sawtooth, because the DT process you used happens to be equivalent to the naive sampling of a CT non-bandlimited sawtooth. I.e. your ears agree with the conclusions you get from thinking of the DT signal as a sampling of a CT signal!

Note that the above applies just as well to all parts of the synthesis/processing chain, which includes the modulators such as LFO's. The only difference is that the bandlimiting must take place (conceptually, this is often not doable in practice) as the final step in an imaginary processing chain (including all modulation etc etc) which happens in CT, the result of which you then sample to form your DT signal.

To take your S&H example: in a CT world, the S&H clock could well sample its input between the samples, and the resulting discontinuous transition would then happen also between the audio samples. If you model your synthesis chain in CT and then bandlimit the whole system, first of all the CT discontinuous transition in the output signal will be slightly smoothed by the anti-alias filter. Further, the sub-sample position of the transition will be reflected in the exact shape of the anti-aliased edge.

Again, whether all of this makes a difference you can actually hear depends on the circumstances, but the point is that a naively sampled LFO does not give you an ideal modulation, but neither does applying modulation from an anti-aliased LFO.


What's the reasoning behind changing pilot tones in OFDM?

Advantage of pilot pseudo-random sequences in standards (with examples are taken from LTE) :

Interference from other sources with different identities (UE ID, Cell ID, ...) are mitigated if detection by correlation operation is used (cell search for example).
Allow including identity in sequence by changing the value of sequence initialization. For example, pilot sequence in LTE is 31-length Gold code which is reinitialized at the start of each OFDM symbol with a value that depends on cell id.


4th order butterworth filter using frequency response formula

Matlab's fft function returns the two-sided spectrum. To get your graphs to agree, start by only looking at the first half of the spectrum.

% note *2 to double the length
in = ones(1,length(t)*2);
in = ifft(in);
[b,a] = butter(n,lo/(fs/2),'low');
outA = filter(b,a,in);
fftA = fft(outA);

% snip out the first half
fftA = fftA(1:length(fftA)/2);

fourierabsA = abs(fftA);
fourierangA = angle(fftA);


Yes, deconvolution.  This page describes a number of deconvolution methods and methods for estimating the point spread function:

Removing Motion Blur from Astrophotographic Images

They say the deconvolution literature is "extremely extensive".   They choose Lucy-Richardson algorithm for deconvolution and develop their own motion estimation algorithm for determining the point spread function.


Looking for low-noise audio data to benchmark some noise-processing techniques

That is a simple question:

All sound that passes into analogue domain using digital audio converters (DACs) have a SNR inferior to about 113 at best, and more probably 96-100... everything recorded using microphones will have inherant noise, unless the recording has been copy pasted onto a silent track leaving intervals of silence.

So what you need to do to have a very clean signal, is to take a totally digital sound which was generated with i.e. sine waves and has not had reverb and echo effects that raise the noise in the silent time(s) of the recording, and digital synthesizers which return to 0 after the note envelope. You can go though a track that contains elements of silence in Audacity or a more easy to zoom editor, and check that you are having times in the track where the amplitude = 0.0; or program a fast test code to count that the track has amplitudes of 0.0 lasting for more than 500ms or something.

if you take any recording and vary the volume of the track, for example you take an AC-DC rock song, and apply an artificial envelope that fades the track to 0.00 regularly i.e. every couple of seconds, it also makes a test file with totally silent intervals.

you can also write a code that sais : if volume < than 0.05, volume =0; to artifically make any track you want totally silent at times.

find digital recordings on freesound and measure them for 0.00 intervals in between sounds, check digital music artists that boast about 192khz recording equipment for a particular album, get a feel of which artists use state of the art production equipment, check standard commercial tracks that are used to measure monitor speakers, there are many options.


Equalizing frequency response in software

Whilst I am not well versed in Scipy, or how the firwin2 function works that you have quoted, this is a very complicated problem that you propose, and there is a lot more to think about than you may realise.

There is one very important component of how sound is radiated by a loudspeaker that is quite often overlooked (especially in marketing specifications for loudspeakers).  Directivity.  You state that your loudspeaker has a particularly poor response below 1 kHz.  I assume that means a “non-linear” response on-axis.  What about off-axis?  This off-axis energy is what interacts with your room, i.e. reflections, and the frequency response in these off-axis directions most likely have a response different to your perfect on-axis response. If you correct just the on-axis loudspeaker response what has that inverse filtering done to the other positions off-axis?  Will you always listen, with one ear, at the intersection of your loudspeaker on-axis directions?

Some loudspeaker manufacturers do actually use technology that attempts to apply “correction” (Genelec, Bose etc…) to overcome room anomalies and positional errors in loudspeaker placement. Whilst the methods employed are very smart, we have to remember that these are just best approximations.  It is also not uncommon to use some form of linear phase FIR filters to correct for minimum-phase anomalies in loudspeaker response for mid-high frequencies, and alternative approaches can be applied for lower frequencies.  Some of these techniques can be better understood by reading text such as: Immersive Audio Signal Processing

Other things to consider: Just how much time data do you window in creating your corrective filter?  And where should you make your measurement? Will you always listen in that position? What I’m emphasising here is that temporal and spatial variation of the soundfield need to be considered also.

So what you are attempting to do isn’t exactly “doomed to failure”, it’s just very complicated, and at the end of the day, the desired result is going to be that which suits your subjective tastes.  I hope that some of the questions I have raised have been helpful, I have raised them simply to remind you that it’s not as simple as inverting a single omnidirectional measurement of a loudspeaker at some point in a room…  Otherwise, yes, every loudspeaker manufacturer would have already done this. Good luck on your endeavour!

I’ll leave a random selection of some references to technical literature related to research in this topic space that you may find helpful:

Combined Quasi-Anechoic and In-Room Equalization of Loudspeaker Responses

Invertibility of a room impulse response

Equalization of loudspeaker response using balanced model truncation

Equalization in an acoustic reverberant environment: robustness results


Simple problem regarding rectangular pulse shifting

Clearly, in this example you have $h(t)=h(T-t)$. The easiest way to see this is to check for which value of $t$ the argument becomes zero. So the point $t=0$ of the original function is mapped to the point $t=T$, and because you have $-t$ in the argument, the function is inverted on the time axis.

The mistake in your first approach is the way you added the value of $T$. What you should have done is replace $-t$ by $T-t$ (which is what's happening on the left-hand side of the equation), to obtain $h(T-t)=1$ for $0<T-t<T$, which is equivalent to $0<t<T$, corresponding to the correct solution.


Are there preferred patterns for creating texture on objects to aid computer vision?

There are some generally applicable advices, and some application-specific advices.

Shi and Tomasi's paper, Good features to track explains the criteria for choosing patterns: two-dimensional localizability, or "cornerness".

To put it simply, suppose you are trying to find an object at position (x,y), but instead the object appears in the image at (x + dx, y + dy). It is not very useful if our vision system can only tell us that "no, the position is wrong." Instead, we expect the vision system to be able to estimate the amounts dx and dy provided that it is not too far off.

A sharp point (dot) is the most cornerful, but it is also easily buried in random pixel noise. By following through with the mathematics, we learn that there are other patterns that are just as cornerful as a sharp point. (Think about a 1D "edge" being a 1D delta transformed by integration.)

Some applications would call for localizability in fewer, or higher dimensions.

Added 8/25

Two line-like patterns can also be "intersected" to yield a point during calibration, provided that lens distortion is not significant or has been parameterized.

In deblurring applications, a sharp point is often used for recovering the point spread function (psf). However, in theory any arbitrary-shaped objects could be used, provided that the ground truth is available to the calibration software.

In some applications, we would deliberately make the pattern un-sharp. Depth from defocus uses the bluriness to reason about the position of the focal plane relative to the object, which gives an estimate of the object depth.


Are the error probability of orthogonal signal and on/off signal the same?

I found that the on/off binary signal is also $\int s_0  s_1 dt = 0$, so can I say the on/off signal is the orthogonal binary signal too, so their error probability should be the same theoretically?

No! That's a wrong inference: Just because the zero vector (the off signal in your case) is orthogonal to everything doesn't mean everything has the same error probability.

When comparing error probabilities, you first have to agree on a way of making the energy spent on each transmitted bit equal.
Otherwise, you might be comparing an On/off-keying transmitter with 1 kW to a BPSK transmitter with 1 mW output power – you'll notice that this comparison will probably be won by the on/off-keying, no matter whether the other one is better.

So, same average bit energy is what we demand. Since one symbol carries one bit in your example, we can set the average symbol energy to 1 for both cases (without loss of generality, as we can later scale), and compare then.

So, sadly, the book you're learning from doesn't do that energy normalization to begin with.

Let's say for the on-off case, $A_{OOK} = 1$, and $T_b=2$, so that the energy the on-period is $E_{On}=P_{On} T_b = A_{OOK}^2 T_b = 1^2 2 = 2$. Obviously, during the off-period, the power is $P_{Off}=0=E_{Off}$, so that the average energy is $E_{OOK}=1$ (you always assume you send as many 0s as 1s, otherwise your source coding is broken).

For what you call the "orthogonal case", the power is always the same – the sign of the signal doesn't matter to the power, since you square. So, we set $A_{ortho}=\frac1{\sqrt2}$ and get $P_{ortho}=A_{ortho}^2=\frac12$ and thus also $E_{ortho}=1$. So, now we're fair!

That means that the pictures you posted aren't fair: The $A$ from the "orthogonal signalling" picture should be $\sqrt2$ smaller than the one from OOK.

We also see that we need a "stronger" amplifier for the On/off case just to get the same signal power. That alone might be an argument why OOK isn't that good.

Now, the rest of the derivation of what an error probability is depends on whether you assume coherent or incoherent reception – you don't define that, but your textbook will certainly do.

You'll, in any case, notice that the two things have different symbol error rates at the same noise power when you normalize their power.


Link between LTI systems and square-integrable functions eigenfunctions

There are three kinds of Hilbert spaces: finite dimensional ones, separable infinite dimensional ones and inseparable infinite dimensional ones. The last category does not play an important role unless you practice quantum field theory, so you can normally ignore them.

Two finite Hilbert spaces are equivalent if and only if they have the same dimension. Two separable Hilbert spaces are always equivalent. Separable Hilbert spaces have a countable basis.

Your question is about the spectral theorem. The spectral theorem comes in different flavours. Finite dimensional Hilbert spaces allow for a very simple formulation: Every self-adjoint operator $A$ has a real eigen-spectrum $\{ \lambda_k \in \mathbb{R} \vert k=1..K \}$ and decomposes into a finite sum of projectors $P_\lambda$ onto the eigen-subspace $E_\lambda$ associated with the eigenvalue $\lambda$.
$$A=\sum_{k=1}^K \lambda_k P_{\lambda_k}$$
The entire finite dimensional Hilbert decomposes into an orthogonal direct sum of the eigensubspaces.
$$\mathcal{H} = \bigoplus_{k=1}^K E_{\lambda_k} $$

This spectral decomposition only implies a unique orthogonal basis if all eigen-subspaces are one-dimensional in which case the corresponding operator is called non-degenerate. If you have a degenerate operator, you can choose an orthogonal basis for each eigen-subspace and get a orthogonal basis.

The more interesting case of a separable infinite dimensional Hilbert space makes things quite a bit more complicated. First, the concept of the spectrum of an operator needs to be generalised. In the finite dimensional case, the eigen-vector equation $A v = \lambda v$ for $v\neq0$ leads to the equivalent condition $\mathrm{ker}(A-\lambda)\neq\emptyset$. For infinite dimensional Hilbert spaces these two conditions are not equivalent and the requirement of a non-empty kernel is more general. Therefore a complex value $\lambda$ is defined to be part of the (generalised) spectrum of the operator $A$ exactly if $A-\lambda$ is singular. The associated kernel or nullspace is the generalised spectral subspace $E_\lambda:=\mathrm{ker}(A-\lambda)$.

The spectral theorem for separable infinite dimensional Hilbert spaces now states that for a compact self-adjoint operator $A$ we get a countable real spectrum and the spectral subspaces add up to the entire space. This is practically the situation we had before, only with infinite sums:
$$A=\sum_{k=1}^\infty \lambda_k P_{\lambda_k}$$
$$\mathcal{H}=\bigoplus_{k=1}^\infty E_{\lambda_k} $$
And again, an orthogonal basis of the Hilbert space is only implied if the operator is non-degenerate.

Compared to the finite dimensional case the big restriction is that the operator has to be compact. Unfortunately, many interesting operators are not compact, because they are not bounded. The spectral theory for non-compact operators is non-trivial and best avoided if possible. This is where Sturm-Liouville (SL) comes in. SL theory is a more specific version of the spectral theorem that works for a certain set of unbounded operators and also provides a method for calculating the spectral vectors. So if you can map an operator to a SL problem, you can make the same statements about the spectrum as for compact operators, even though the operator is unbounded.

Now, how does all this related to linear time-invariant (LTI) systems? Linear operators describe linear systems. Time-invariance can be added as an additional constraint. If you have a linear system operator $A$ and a unitary time-shift operator $S(\tau)$, then $A$ is time-invariant if and only if $AS(\tau)=S(\tau)A$, that is $A$ commutes with $S(\tau)$, for all real $\tau$. This is often written as $0=[A,S(\tau)]:=AS(\tau)-S(\tau)A$ with the commutator $[\cdot,\cdot]$. There always exists a common spectral basis for two commuting operators, that means both operators can be made diagonal in that basis. We will also need that an operator commutes with any function of that operator $[A,f(A)]=0$.

The time-shift operator can be expressed using a Taylor series generated by the exponential function $S(\tau)=\exp(\tau \frac{d}{dt})$. With what has been said above, it follows that LTI operators commute with $\frac{d}{dt}$ and therefore also with $\frac{d^2}{dt^2}$ and they share a spectral basis.

This should answer your first question. I will add more later to answer the two remaining ones.


Calculating poles with complex numbers and quadratic equation

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

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

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

agreeing with Wolfram Alpha.


Compute Hann window without cos function

i think this generates a sinusoid:
$$ y[n] = 2 \cos(\omega_0) y[n-1] - y[n-2] $$
if you initialize
$$ y[-2] = \cos(2 \omega_0) $$
$$ y[-1] = \cos(\omega_0) $$
then the result is
$$ y[n] = \cos(\omega_0 n) \qquad n \ge 0$$
the period of this cosine function is $N=\frac{2\pi}{\omega_0}$ and you want $N-1$ to be the non-zero width of the window.
finally you construct your window with
$$ w[n] = \tfrac{1}{2} (1 - y[n]) \qquad  0 \le n \le N$$


get poles and zeros of frequency response

This is a very complicated problems and I don't think there exists a one-size-fits-all solution. You can try Matlab's $invfreqz()$ and see if it works for your purposes https://www.mathworks.com/help/signal/ref/invfreqz.html

In general this is a error minimization problem but the actual data and the way you set the your error function and the search strategy makes a huge difference in the outcome. Things to watch out for

How good is your data? What's the frequency coverage range, density and shape of frequency grid, signal to noise ratio, etc.
What do you care about in the result. Dips vs peaks ? Phase vs amplitude, ?
Stability and causality constraints. How do you deal with non-causalities in the data (real or noise induced)
Search space representation, search strategy & algorithms.
Error surface shaping and optimization.
Convergence tracking, stop criteria, metrics for success (or lack thereof)

I developed one of these that is optimized around the specific constraints and requirements of my application and found invfreqz() to be totally useless. Your mileage may vary.


Sox, mixing mono inputs into mono output (efficiently)

As far as I can understand, you want to add two signals, together?
I guess that you are missing -m switch.

Generate two test files (sweep and sine):

sox -Dr 16000 -c 1 -b 16 -n sine.wav synth 5 sin 1000
sox -Dr 16000 -c 1 -b 16 -n sweep.wav synth 10 sine 100:8000

Merge them together:

sox -m sweep.wav sine.wav -Dr 16000 out.wav

Using soxi out.wav you can see that produced file is mono:

Input File     : 'out.wav'
Channels       : 1
Sample Rate    : 16000
Precision      : 16-bit
Duration       : 00:00:10.00 = 160000 samples ~ 750 CDDA sectors
File Size      : 320k
Bit Rate       : 256k
Sample Encoding: 16-bit Signed Integer PCM

Now you can observe the spectrograms:

sox sine.wav -n spectrogram -o sine.png
sox sweep.wav -n spectrogram -o sweep.png
sox out.wav -n spectrogram -o out.png

Sine:

Sweep:

Merged together:


Swift vDSP rFFT not same as Python np.fft.rfft()

Sanity check your data layouts.
It is very common for FFT libraries to pack the (real-valued) Nyquist-bin into the unused imaginary component of the DC-bin and it happens almost naturally when you implement a real-transform as a half-size complex-transform.
As far as I can tell, vDSP wants such "packed" layout, where as Python apparently stores it as a separate "unpacked" real-component in a larger array.


Is the typical implementation of low pass filter in C code actually not a typical low pass filter?

The "simplistic" implementation that you referenced above is often called a leaky integrator, a special case of a first-order IIR lowpass filter. They are discussed in more detail in several previous questions:

Is a leaky integrator the same thing as a low pass filter?
Is there a technical term for this simple method of smoothing out a signal?
Designing a practical integrator
"Low-pass filter" in non-EE, software API contexts

To answer your question specifically, yes, it is a lowpass filter. You can certainly construct more complicated filters with a better lowpass frequency response. The main advantages of the leaky integrator construct are its very simple implementation, and the fact that you can adjust its cutoff frequency by simply changing the value of $a$ (as $a \to 0$, the cutoff of the filter decreases).


How to set the "duration" parameter for a flanger with LFO oscillator

The LFO rate is the speed of the delay modulation. One full cycle is going from "down" to up" back to "down" again. It's in Hz, i.e. "per second". Let's says the music is 120 bpm and your LFO rate is 0.5Hz. 120 pbm is 2 beats per second or 4 beats per flanger cycle.
Beats Per Cycle = Music_speed/60/flanger_rate


Here's a short snippet of code that returns the domain of the FFT depending on if you have the onesided or twosided FFT. This is from some code I posted on a somewhat related question here.

% Nfft - FFT size
% Fs - Sampling frequency in Hz
% oneside - flag indicating FFT results will be onesided or two sided

% Compute the frequency vector
if oneside == true
  f = fs * (0:Nfft/2).'./Nfft;
else
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end


Deblurring algorithm to precede thresholding - speed over accuracy

You only need to find and report positions of nine specific symbols.  You can do this by computing correlations (done quickly by multiplying in Fourier space) with samples of the nine digits snipped from a typical puzzle image, blurred as it is.  All blurred '9' will look the same, and likewise all other digits.  I'm assuming puzzle images are taken with the same camera the same way.

It might work better working on high-pass filtered versions of the input puzzle image and image snips, by subtracting from each a more-blurred copy.


Functionally PCA, FFT could be interchanged without changing anything.

e.g.

[U,~,~] = svd(A);
[FU,~,~] = svd(fft(A));

fft(U) and FU are the same (up to an arbitrary phase rotation)

However, it may help speed to do BPF first, then FFT, and discard the stopband to reduce dimensionality before doing PCA.


Is there a better Fourier Transform-based algorithm to use in Image recognition

The Fourier Transform is good for picking up some similarities between images:

Shifts in $x$ and $y$ are just phase changes in the Fourier domain.
Rotations in the $x-y$ plane cause rotations in the Fourier domain.

So, if this is what you are trying to achieve, then it makes sense.

See this paper for a performance analysis of doing this using the following algorithm


Generating Music Tones, How to handle when tones overlap? (16bit samples)

If you can't tolerate any distortion, then divide by 2 the amplitude of the signals you want to mix.

Please check this question


Can disparity range from negative to positive values for parallel cameras?

It was the way I was setting principal points during rectification, to center the images back for the warping. If I have same principal points in both cameras I only get positive disparities, as expected. If cameras are parallel, it is not possible to have both negative and positive disparities. That only happens with converging stereo cameras (negative disparities for objects in front of point of convergence and positive for objects in the back).


How to filter objects on basis of orientation of regionprops

I don't have MATLAB, but here are a few ideas:

Since you use regionprops, you could use the Orientation parameter - which is pretty much exactly what you want to have.

Another idea

Create an matrix of the same size as your image. The matrix has entries that give a linear variation in the perpendicular orientation that you desire (it basically is a gradient field in the perpendicular orientation, e.g. with values ranging from 0 to 1).
Use bwlabel to label all your candidate lines
Loop over all labels and create a binary mask for each label. Multiply the mask with the gradient matrix.
Calculate the variance inside the mask. Since the gradient is perpendicular to your direction, the variance should be 0 for a perfect line. For a not-so-perfect-line, there is a certain variation. You could use

an empirical threshold
check the math for an more elaborate threshold, because one can certainly find a function $\sigma^2(\alpha)$, with $\alpha$ being the angle of the region of interest in the mask.


A zero coefficient at DC simply means the mean of the waveform is 0; bin 0 of the DFT is identical to calculating the mean value scaled by $N$.
Consider the general expression for the DFT:
$$X(k) = \sum_{n=0}^{N-1}x(n)e^{-j k \omega_n n}$$
with $\omega_n = 2\pi/N$ and $k = 0,1 \ldots N-1$
(This is a correlation of the waveform x(n) with every integer multiple of the tone with frequency $\omega_n$)
When $k=0$ this reduces to:
$$X(0) = \sum_{n=0}^{N-1}x(n)$$
While the mean of $x(n)$  is
$$\bar x = \frac{1}{N}\sum_{n=0}^{N-1}x(n) = \frac{1}{N}X(0)$$
This is how the DFT works fundamentally; for each index k in the DFT, the signal $x(n)$ is frequency translated by $-k \omega_n$, and then we take the average of that result (without dividing by $N$, so a scaled average). If any frequency component was at $k \omega_n$ prior to the frequency translation, the resulting average would be large, indicating the presence of that frequency component.
Note for just one bin of the DFT how this is the basic structure of a Digital Downconverter as shown in the block diagram below for $x(n)$ real, where the low pass filters (LPF) are basic moving average filters (such as a CIC filter). The DFT computes one time sample, and replicates this for every single bin, and supports real or complex input signals where with a complex input there would be 2 additional multipliers and 2 adders for a full complex multiplication. So to put that in perspective, a 1023 point DFT is 1023 digital down-converters all running in parallel, which makes the FFT algorithm very impressive! Every computation of the DFT is one time sample of the result, so to process as a streaming receiver the input would be shifted one sample and then a new output computed (overlap-add and overlap-save methods improve the efficiency further for such streaming applications but shifting by one sample and computing each output makes one bin of the DFT and the digital-downconverter receiver shown below equivalent).
One Bin of the DFT as a DDC (for case of real input)


Efficient FFT computation of a zero-padded vector

You can decompose the FFT that you describe into some smaller transforms. Take a look at what you're calculating when you zero-pad a vector $x[n]$ of length $N \over 2$ to length $N$ and calculate a DFT:

$$
X[k] = \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi nk \over N}
$$

To split this up into two transforms, first look at just the even-indexed terms in the $N$-point sequence $X[k]$:

$$
\begin{align}
X[2k] &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n2k \over N} \\
&= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi nk \over \frac{N}{2}}
\end{align}
$$

Via inspection of the above, we can conclude that the even-indexed terms of the zero-padded DFT are just equal to the $\frac{N}{2}$-point DFT of the original input sequence (with no zero-padding). What about the odd-indexed terms?

$$
\begin{align}
X[2k+1] &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n(2k+1) \over N} \\
&= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n \over N} e^{-j2\pi nk \over \frac{N}{2}}
\end{align}
$$

This is a similar expression to the even-indexed terms $X[2k]$. We take the original $N\over 2$-length sequence, multiply it by a complex exponential sequence $e^{-j2\pi n \over N}$, then calculate the $N \over 2$-length DFT of the result. This makes intuitive sense if we write the exponential as $e^{-j2\pi 0.5n \over \frac{N}{2}}$: multiplying by such a complex exponential yields a shift in the frequency domain of half of a $N \over 2$-length DFT bin width. This gives the interpolation-by-2 effect that the zero-padding is known to yield.

In summary: To decompose an $N$-point DFT that consists of an $N \over 2$-point signal $x[n]$ followed by $N \over 2$ zeros, do the following:

Calculate an $N \over 2$-point DFT of the original signal $x[n]$. This make up the even-indexed values in the zero-padded DFT result.
Multiply the original signal by the complex exponential function $e^{-j2\pi n \over N}$. Calculate an $N \over 2$-point DFT of the product. This makes up the odd-indexed values in the zero-padded DFT result.

The caveat: While it depends on your DFT size $N$, this may not effect any tangible reduction in complexity beyond just calculating the zero-padded DFT directly. I have an application where I've tried to make a similar optimization myself in the past, although it is on a very different platform from what you described (mine runs on high-performance server machines using highly-optimized FFT libraries like Intel Math Kernel Library). While you may be able to achieve a theoretical reduction in the total number of arithmetic operations that you need, due to other effects like memory accesses, caching, and so on, it's common that trying to do 2 $N \over 2$-point FFTs is slower than just doing one $N$-point FFT directly. In my case, I've never been able to beat the existing library's $N$-point transform.

To estimate the complexity reduction by doing this, we can estimate the number of arithmetic operations that each method needs. Using the rule of thumb that an $N$-point FFT requires approximately $5 N \log N$ arithmetic operations, we can estimate the savings in operations:

$$
\text{savings} = 5 N \log N - \left(2\left(5 \frac{N}{2} \log \frac{N}{2} \right)+5 \frac{N}{2}\right)
$$

$$
\text{savings} = 5 N \log N - \left(5 N \log \frac{N}{2} +5 \frac{N}{2}\right)
$$

The first term in the parentheses counts the approximate number of operations required by the 2 $\frac{N}{2}$-point DFTs, while the second term counts the number of operations needed for the $\frac{N}{2}$-point complex vector multiplication (there are 5 scalar operations per complex multiply).

The problem is, as $N$ increases, the amount of savings diminishes because of the logarithmic scaling of the DFT. For typical DFT sizes, you will likely find that this approach yields little complexity reduction.


Make your block diagram
At each point where significant quantization can happen, add noise
Analyze your system's behavior with that added noise

If you know that the quantization effects will be essentially random, and if the following stages tend to low-pass or band-pass filter, then model the quantization noise as Gaussian with $x_n \sim N(0, q/12)$, where $q$ is one LSB.

The justification for this is that the quantization noise will be uniformly distributed (hence $\sigma = q/12$), and the following filter's output will contain a weighted sum of a bunch of samples of the quantization noise, which will tend to Gaussian by the central limit theorem.

If you want to be Maximally Paranoid, then figure out what the worst-case behavior of the quantization noise would be for your application (stuck high, stuck low, some specific signal, etc.), and model the quantization noise as having that shape, with a magnitude of $q$.  In the case of a linear system that has any frequency selectivity, this works out to a square wave at a frequency equal to the systems highest sensitivity to noise at the injection point (or a sine wave with an amplitude of $1.09 q/2$.

The justification for this is really just paranoia -- but if you have an absolute upper bound on the effects of quantization, and the system still works well enough, then you just know you're done.

To design data path widths, you can turn this around, and determine the acceptable level of quantization at each step, and make sure that your data paths are wide enough (and scaled correctly) so that quantization is smaller than your acceptable level.

It's dangerous to make general statements about nonlinear systems, but I would hazard a guess that if you can analyze the algorithm on paper at all, there's a good chance that you can analyze the algorithm plus quantization on paper.


I can't figure out why spectral features are "displaced" on a spectrogram

After subtracting Fs/2 from F to get a frequency axis from -Fs/2 -> Fs/2 rather than from 0->Fs

You can't center the frequency by changing the labels on the x-axis. Use 'fftshift' instead.


Adaptive filter to scale and phase shift two sensors output

It seems that you want to equalize one time series to the other or vice versa or some combination. The general problem is then (in shorthand)
$$
h_1{\Large *}y_1 \approx h_2{\Large *}y_2
$$
where $\Large *$ denotes convolution, $h_i$ are filters, and $y_i$ are your time series. If you simply equalize one time series to the other, this is
$$
y_1 \approx h_2{\Large *} y_2
$$

In matrix form, this becomes
$$
\left[\begin{array}{c} y_1[p]\\y_1[p+1]\\\vdots\\y_1[p+N-1]\end{array}\right] \approx
\left[\begin{array}{cccc}y_2[2p] & y_2[2p-1] & \cdots & y_2[0]\\
y_2[2p+1] & y_2[2p] & \cdots & y_2[1]\\
\vdots & \ddots & \ddots & \vdots \\
y_2[2p+N-1] & y_2[2p+N-2] & \cdots & y_2[N-1]
\end{array}\right]
\left[\begin{array}{c}h_2[0]\\h_2[1]\\\vdots\\h_2[2p] \end{array}\right]
$$
or more concisely
$$
\underline{y}_1 \approx {\bf Y_2}\underline{h}_2
$$
The least squares solution for the equalizing filter coefficient vector is
$$
\underline{h}_{2,LS} = \left({\bf Y_2^HY_2}\right)^{-1}{\bf Y_2^H}\underline{y}_1
$$

You can alternatively use LMS or RLS to do recursive updates to find/track $\underline{h}_2$.

One additional comment on equalizing one time series to the other: by choosing to equalize $\underline{y}_2$ to $\underline{y}_1$, we choose one spectral weighting - implicitly. If we instead equalized $\underline{y}_1$ to $\underline{y}_2$, we would have a slightly different spectral weight for the overall least squares fit. If the filters are relatively minor in magnitude variation over the signal band, the spectral weight difference is similarly minor, but if the equalizer filter response is quite different for the two cases (can just as well solve and analyze both), you will need to assess which is preferred/better and why.


Automatic detection of noise-only segments in audio

The way I read it, the algorithm works by figuring out where the signal part lives in the frequency domain by using the statistics of the noise part. This means you don't need a noise-only signal for each signal you want to process as long as the statistics of the noise are the same for each signal realization.
For example, if you had 1000 signal realizations that you want to perform this noise reduction on, you would use a noise-only signal to compute the statistics of the noise. Then using those statistics, the algorithm will figure out signal versus part by comparing to a threshold (average noise power or something I would guess).


Converting mel spectrogram to spectrogram

Both taking a magnitude spectrogram and a Mel filter bank are lossy processes.  Important information needed to reconstruct the original will have been lost.  Thus you need to go back and use the original audio samples to do the reconstruction by determining a time or frequency domain filter equivalent to your dimensionality reduction.

You can make assumptions about the lost information, but those assumptions themselves usually sound inaccurate, artificial and/or robotic.  Or you can use only specially synthesized input, where the assumptions will be correct by design of that input.


How wide is the word width of the processor receiving these 10-bit values from the A/D converter?  Wider than 8 bits?

What is there to gain by adding another layer of quantization error?

What you have here is a highly dithered signal that is dithered naturally.  Sometimes in audio we actually go out of our way to dither a signal so that further processing and quantization is better.

If your processor is 8 bits wide and you absolutely do not want to do multi-byte arithmetic (with an Add-with-carry instruction) or do not have the computational bandwidth to do higher precision arithmetic, then perhaps it makes sense to quantize your noisy 10-bit value to an 8-bit value.  Otherwise it does not.


Cascade Transfer Function; Odd or Even Transfer Function

If the order $N$ of the system is even, you get $N/2$ second-order sections, where each second-order section combines either a pair of complex conjugate poles and zeros, or two real-valued poles and zeros. If $N$ is odd, you get $(N-1)/2$ second-order sections and one additional first-order section for the remaining real-valued pole and zero.


Discrete Fourier transform of a finite length signal which saturates at non-zero value

You can & should window your signal. Choice of window depends on the specific requirements of your application. Wikipedia has a really good overview on properties and trade-offs: https://en.wikipedia.org/wiki/Window_function
Keep in mind that you already ARE windowing: your finite length window is simply an infinite length signal multiplied with a rectangular window.


You are right. DFT and FFT are the same things; they produce the same outputs and they operate on the same samples. FFT is just a fast implementation of DFT with much less artihmetic operations. They only have roundoff error difference which is due to the different number of arithmetic operations involved. Therefore you either misunderstood the article or what they say is wrong.


Difference between compressive sensing and DCT-based compression?

Comparing compressive sensing and traditional compression is a complicated task. Even  comparing two compression algorithms.

One the one hand you have compression ratios, making sure the bits sent can be decompressed without other hidden bits. On the other hand, error measurements, rarely perfectly matched to the actual meaningful distorsions. [EDIT] Standard SSIM or NMSE are not well adapted to sparse signal.

Then, knowing the range over which one algorithm is better than the other is important, because one can be better at low-rates, the other at high-rates. There is rarely a free lunch in compression.

From a very mundane point of view, if you are in a compression range where the signal is not "sparse" enough (with respect to some basis) and the CS quantization is not really taken into account, yet sufficiently harmonic, it is well understandable that a DCT-like engineered algorithm can do better than a CS method if the CS theoretical conditions are not fulfilled.


Your question (as expanded in the comments) is asking if we start with
$$
x(t)  = x_{\rm true}(t) + n_1(t)
$$
and filtering it using a filter $H(\omega)$ to get
$$
x_{\rm hp}(t) = h(t) * x(t) = x_{\rm true}(t) + n_2(t)
$$
where the variance of $n_1$ is $v_1$, then what is $v_2$, the variance of $n_2$?

That seems ill-posed because the filter $H$ will change $x_{\rm true}$ as well as the noise.

This page looks at the problem when $x_{\rm true}(t) = 0,\ \forall t$.

Here it is assumed that $n_1$ is white noise.

Where they get it is that
$$
r_{x_{\rm hp}}(l) = v_1 \cdot (h * h)(l)
$$
and $v_2$ is just this evaluated at $l=0$ i.e.
$$
v_2 = r_{x_{\rm hp}}(0)
$$

So you need to find the impulse response of your filter, convolve it with itself and find the zero-lag value of this.

As I said above, this doesn't quite answer your question but it might be a start.


If you filter complex data $x[n]$ with a filter with a real-valued impulse response $h[n]$, your complex-valued output is given by the convolution sum

$$y[n]=\sum_kx[k]h[n-k]\tag{1}$$

The real and imaginary parts are

$$y_R[n]=\sum_kx_R[k]h[n-k]\tag{2}$$

and

$$y_I[n]=\sum_kx_I[k]h[n-k]\tag{3}$$

So whatever you do, in order to compute the complex-valued output $y[n]$ you need to perform two real-valued convolutions / filtering operations.


It's been some time since I asked this question and after some work done on this subject, I think it's time to revisit it. It should be noted that
while the spectral centroid, pretty much like any other spectral feature, is calculated using the FFT magnitudes, and not raw FFT data, we don't
have to manipulate the magnitudes directly. Instead, we modify FFT data, which in turn, affects the magnitudes. This makes the process
a lot easier than having to manipulate the magnitudes.

The spectral centroid is a single value (per FFT frame) so in the case of STFT we have as many centroid as we have STFT frames. The process of manipulating
the spectral centroid is pretty straightforward. If, for instance, the spectral centroid is 600Hz and we wish to increase its value by 10% so that the new value is
660Hz, we would simply rearrange our frequency bins by a certain amount (a bin step). The actual calculation is qute simple

const binStep = frameLenHalf - frameLenHalf/(1+(float)freqDelta/100);
   //frequency delta % (e.g. 10%, 20%, -10%, -20% etc)

    //for increasing frequency
    float *pCfftLReal = new float[frameLen];
    float *pCfftLImg =  new float[frameLen];
    //for increasing frequency
    for(j=frameLenHalf-1;j>binStep;j--){
            pCfftLReal[j] = pFFTReal0[j-binStep];
            pCfftLImg[j] = pFFTImg0[j-binStep];
                }
        }
    //for decreasing frequency
    for(int j=0; j<frameLenHalf-1-binStep; j++){
            pCfftLReal[j] = pFFTReal0[j+binStep];
            pCfftLImg[j] = pFFTImg0[j+binStep];
                }

Needless to say, this process suffers from rounding errors (primarily calculating the binStep), but in general, the larger the bin resolution, the larger the errors.
This can be alleviated by having a greater frame length (but for STFT this is usually undesirable). Other more complex approaches are also possible but I haven't investigated any of these.

Also, sometimes the results can be skewed if the amount of frequency shift exceeds the Nyqiust frequecy and we keep the original sample rate (upsampling would help here).

Before we can actually listen to our manipulated sound, an inverse FFT shoud be carried out on our modified data.

What I've described here is just a rough, naive approach to spectral centroid manipulation. It should also be said that this kind of manipulation
doesn't retain the original harmonics ratio, given that the frequency delta shift is applied equally on every bin. Keeping the harmonics ratio intact would require
calculating a new bin step for every frequency, which is only slightly more complex that the algo given above.

I've tried both approaches (with a constant bin step and a changing bin step) and perceptually, the first doesn't really change the original sound ( I guess
this has to do with higher frequencies affecting timbre more than lower ones and the first approach changes the higher frequencies less than it does the lower ones).

I'm not done with this yet, so I will update this post as soon as I have something interesting to share.


I have an acoustic signal of a test and I want to remove the "thin"
  peaks, because they don't really affect the structure that I'm
  studying. For that I've thought of a Low pass filter in order to
  retain only frequencies below 300 Hz.

If this is accurate, that is, you have an acoustic signal that exists from baseband to 300 Hz, then yes, what you want to do is a digital low-pass filtering operation to retain all energy below 300 Hz, and discard energies above that.

  My aim is to count the peaks of the filtered signal above 3 sigma
  (Standard Deviation). I'm programming in C# and I don't know if there
  is a direct (and simpler way) of doing this but what I'm doing now is

Your intuition is correct - you want to remove frequencies below 300 Hz, so naturally you think to transform your signal into the Fourier domain, clobber unwanted frequencies, and then transform back. There are some issues with this method:

Firstly, this method is taboo in the DSP community, (and for decent reasons), although there are cases where you can get away with it, depending on your application. Some background: Nulling out unwanted co-efficients in a transform domain and inverse transforming back into the original domain is used in denoising especially in wavelets, and it also makes intuitive sense, (and correctly so). The problem with doing this in the Fourier domain is that the basis function of the DFT are global, and not local like they are in wavelets. In other words, you can clobber our wavelet co-efficients and affect only local areas in your signal (ie, small parts of your signal) - but doing so in the Fourier domain affects all your signal.
Expanding on the above background, when you "boxcar" your frequency domain, (ie, mask it with 1s and 0s that either allow or disallow certain frequencies), this is the same as convolving your time-domain signal with a sinc function, that carries on forever. This is why simply culling out unwanted frequencies - while acceptable in some applications - will introduce 'ringing/smearing' in the time domain.
Instead what is usually done, is that people create filters (masks in the frequency domain that are not sharp 1s and 0s, but more graceful), such that their equivalent transform back into the time-domain does not have the effective infinite extent of sinc function, (that is, the decay of the sidelobes is much faster than that of a sinc). This then provides a nice trade-off between removing unwanted frequencies, and not having to pay the price of extreme ringing/smearing in the time domain.

Armed with this task at hand, what you need to do is the following:

Create a digital low pass filter in some other software such as MATLAB or Python. You can create a simple FIR filter for your purposes. (Do not use an IIR because you are working with an audio application and default IIR filters will warp your phase, but thats another story). When you make your filter, you will simply have a vector of numbers, of a length you chose. The larger the length of the filter, the sharper the transitions in the frequency domain, but also the larger the delay, so this is a trade-off you can experiment with. At the end of the say though, you will $1$ by $N$ vector of numbers. This is your filter. Let us call it $h[n]$
You have your audio signal, let us call that $s[n]$. The planets are not aligned, so most likely your audio signal is of a different length than your filter, and let us say it is of length $M$. So your signal $s[n]$ is simply a $1$ by $M$ vector of numbers.
Your objective is to filter $s[n]$ by $h[n]$. There are a number of ways to do this. Since you already seem to have access to the FFT in C#, we will use the FFT method. In the FFT method, you want to multiply the Fourier transform of your filter, (call it $H(f)$), with the Fourier Transform of your signal, (call it $S(f)$).
We doing a linear convolution, (as opposed to a circular one), but in the frequency domain, so we have to take extra care when it comes to FFT sizes. If you did the linear convolution in the time domain, the resulting filtered signal would be of length $N + M -1$. This is then the FFT size you need to use.
You are almost done. All you need to do now is take the $N+M-1$ length FFT of your filter, $h[n]$. Take the $N+M-1$ length FFT of your signal, $s[n]$. Multiply them together. Then take the $N+M-1$ length IFFT of this element-by-element product, and viola, you have filtered your signal.

Hope this helped.


Fourier transform of dirac comb with function: The scaling factor

Your first formula is wrong. If you use angular frequency ($\omega$), then multiplication in the time domain corresponds to $1/(2\pi)$ times convolution in the frequency domain, just as in the last formula in your question:

$$\mathcal{F}\{x(t)\cdot y(t)\}=\frac{1}{2\pi}X(j\omega)\star Y(j\omega)\tag{1}$$

If you use frequency $f=\omega/(2\pi)$, then you get

$$\mathcal{F}\{u(t)\cdot v(t)\}=U(f)\star V(f)\tag{2}$$

where I've used the definitions

$$X(j\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{3}$$

$$U(f)=\int_{-\infty}^{\infty}u(t)e^{-j2\pi f t}dt\tag{4}$$


Calculating an in-loop signal as part of a hierarchical control loop

UPDATE: My original answer incorrectly assumed the outer loop would not impact the result, but as pointed by Ben, this transfer function is dependent on the fast feedback.

The resulting equations are found by solving for $A/X$ from the block diagram below, with X representing the disturbance input:

Solving this:

$$Y = P(SY+FY+X)$$
$$A/S = P(A + FA/S + X)$$
$$A = PAS + PFA + PXS$$
$$A(1- PS - PF) = PXS$$

$$\frac{A}{X} = \frac{PS}{1-PS-PF}$$

Note that this follows the general equation for the closed loop gain using:

$$G_{CL} = \frac{G_F}{1+ G_{OL}}$$

Where:

$G_{CL}$: Closed Loop Gain

$G_{F}$: Forward Gain (In this case $G_{F}  = PS$)

$G_{OL}$: Open Loop Gain (In this case $G_{OL} = -P(S+F)$

Note how negative feedback is assumed in the generic expression, but since the block diagram does not show negation at the summations, the open loop gain is then negative as given above. Using the general equation and recognizing that the parallel result of S and F as a gain given as S+F allows us to solve this from inspection.


I've heard rumors of covert acoustic room impulse response measurement being done during concerts. You would need the people in the room if you want to use the impulse response for artificial reverberation and you don't want the room to sound overly echoey like there's no-one in it. As received at the microphone, the signal would be a rather quiet noise-like test signal, and the music, the sounds from the people, and equipment noise would be the noise. Because the impulse response is in practice only a few seconds long, a repetitive test signal can be used. The recording would be accumulated into a circular buffer. Because noise has on average zero correlation with what has already been recorded in the buffer, its root mean square amplitude will grow on average as $\sqrt{N}$ for $N$ recording cycles, while the signal's root mean square amplitude will grow as $N$. At the end of the session, deconvolution is used to get from the test signal response to the impulse response. For a spectrally flat test signal such as a maximum length sequence (MLS), it suffices to convolve with the reverse of the test signal. For a buffer of $M$ sample points the deconvolution gain is $M$ for the signal and $\sqrt{M}$ for the noise. Combining the gains from repetition and deconvolution and assuming a sampling frequency of 44.1 kHz, we get a recording time $t$ dependent dB-scale gain difference of $20\log_{10}\left(\frac{44.1\text{ kHz }\times\ t}{\sqrt{44.1\text{ kHz }\times\ t}}\right)$ dB between the two:
Figure 1. Signal-to-noise ratio (SNR) gain for room impulse response measurement as a function of recording time.
Guessing some numbers for an "easy" scenario, if we start with the signal at a sound pressure level (SPL) of 30 dB and the noise at 60 dB SPL which is a starting SNR of -30 dB, in 1 hour of recording we get a SNR of 52 dB, growing very slowly after that, by 10 dB for every 10-times increase of the recording time. 52 dB is already useful for artificial reverberation.


why wavelet coefficients are zeros after decomposition

If you recall, db1 is the Haar Wavelet.  The Haar Wavelet takes either an average of pixels within local windows for the approximation coefficients (or the LL band), or a difference of pixels within local windows for the detail coefficients (or the LH, HL and HH bands).

Be advised that the input image that you specified only consists of two intensities: 0 and 255.  Also, you set a square grid within this image to be 255 and it is uniformly shaped.

For self-containment, this is what your test image looks like:

This uniformly shaped object within a square image is important as part of the reasoning why you are not getting any output for the detail images (HL, LH, and HH).

The best way to describe why you're seeing an output for db2 and not db1 can be shown visually.

This slide is from the University of Toronto's CS 320: Introduction to Visual Computing course, specifically, the Discrete Wavelet Transform lecture:

You are well aware that when you take the 2D DWT, you produce 4 sub-images that are half the resolution of the original image.  The first output of dwt2 are the approximation coefficients where each output pixel is an average of a 2 x 2 window.  The other outputs (second, third and fourth) are detail windows that take two pixels within the window and subtracts these with two other pixels in the window.

As such, the reason why you are not getting an output with db1 is because all of your calculations for the detail images will cancel out.  Specifically, you will get 2 x 2 windows of either completely 0, or completely 255, and when you calculate the detail images, you will get 0s for the output regardless.  You would take add two 0 values, or two 255 values and you would subtract these with two other 0 values, or two other 255 values respectively, thus causing the output to be 0 regardless.

The db2 wavelet is a more complicated transform which is a weighted sum of non-uniform coefficients and so you will certainly get outputs for the detail images rather than just a simple differencing of the 2 x 2 windows.

I would like to stress that if you have a more complicated shape that is non-uniform, db1 will certainly not give you a zero output.  Try it on any test image that comes with MATLAB, like with cameraman.tif.

Hope this helps!


I prefer to solve such problems using the $\mathcal{Z}$-transform. First, write down the time domain equations:

$$\begin{align}\tag{1}
y[n]&=r[n]+x[n]\\
r[n]&=r[n-1]+t[n]+t[n-1]\\
t[n]&=t[n-1]+y[n-1]+y[n-2]
\end{align}$$

Taking the $\mathcal{Z}$-transform of these equations gives

$$\begin{align}\tag{2}
Y(z)&=R(z)+X(z)\\
R(z)&=z^{-1}R(z)+T(z)+z^{-1}T(z)\\
T(z)&=z^{-1}T(z)+z^{-1}Y(z)+z^{-2}Y(z)
\end{align}$$

After simplification you get

$$\begin{align}
R(z)&=T(z)\frac{1+z^{-1}}{1-z^{-1}}\\
T(z)&=Y(z)z^{-1}\frac{1+z^{-1}}{1-z^{-1}}
\end{align}$$

from which

$$R(z)=Y(z)z^{-1}\frac{(1+z^{-1})^2}{(1-z^{-1})^2}\tag{3}$$

follows.
Plugging (3) into the first equation of (2) gives

$$Y(z)=Y(z)z^{-1}\frac{(1+z^{-1})^2}{(1-z^{-1})^2}+X(z)\tag{4}$$

from which you finally obtain afters some algebra

$$H(z)=\frac{Y(z)}{X(z)}=\frac{(1-z^{-1})^2}{1-3z^{-1}-z^{-2}-z^{-3}}\tag{5}$$

In the time domain, Eq. (5) is equivalent to

$$y[n]-3y[n-1]-y[n-2]-y[n-3]=x[n]-2x[n-1]+x[n-2]\tag{6}$$


Determine minimum distance between two symbols in 16-QAM using the average symbol energy

A symbol with coordinates $(x,y)$ has energy $x^2+y^2$. In 16-QAM, minimum distance $2A$ implies that the values of $x$ and $y$ are restricted to the set $\lbrace\pm A,\pm3A\rbrace$, and in consequence the possible symbol energies are $2A^2$, $10A^2$ and $18A^2$.

Furthermore, there are 4 symbols with energy $2A^2$, 8 symbols with energy $10A^2$, and 4 symbols with energy $18A^2$. Calculating the averagey symbol energy, we conclude that $E_s=160A^2/16=10A^2$.

Note that a similar procedure can be used to calculate the average symbol energy of any quadrature modulation, even if they're not square (or even rectangular).


What are the 'phases' in this 2D scaler filter

My best guess is that they limit the available interpolation points to 16 discreet points in x and the same for y, and so any fractional scaling where a sample point falls between 2 of these 1/16th sampling points, will just use the closest point to it?

Yes, that's it.  Probably on the assumption that within 1/16 of a pixel shift, and with the further interpolation on the output, that filter will be good enough.


Yes the signal is perfectly reconstructed. Consider the process at each stage as I show using the block diagram below:

Consider each sample of the signal at each node in the diagram (each sample is shown using the sample index at the node for each row):

(Note: You see the same form of reconstruction in the FFT algorithm).

I will attempt to illustrate how the aliasing is cancelled (addressing MBaz's comment):

First consider the frequency response of $z^{-1}$ (and $z^{+1}$).  $z^{-1}$ is the z-transform of a DELAY of 1 sample in the time domain (see How/why are the $\mathcal Z$-transform and unit delays related?). The frequency response of a delay is constant magnitude, and linear phase (such as we would see with a cable of fixed length; very low frequencies would only be delayed a fraction of a cycle, while higher frequencies will be delayed several cycles:

Specifically for $z^{-1}$, since the delay is 1 sample per cycle, the amplitude will be constant (1) and the phase will be linear from 0 to $-2\pi$ for the frequency axis going from 0 to our sampling rate. (And $z^{+1}$ which is an advance instead of delay will be similar with positive phase versus frequency:

Phase only shown below (the magnitude is 1 for all frequencies)

Consider the following digital spectrum to see what happens as we pass through the system presented by the OP. (Note, when dealing with multi-rate systems, as well as mixed signal analog-digital systems, it has helped me to view the frequency axis from $-\infty$ to $+\infty$, which I depict below).

Now observe the beautiful thing that happens when we compare our original spectrum to that after a unit delay. I am emphasizing the that phase at $F_s/2$ has rotated 180°, and then fully 360° which is back to 0° at $F_s$.  This rotate is continuous from the linear phase shift from 0 to the sampling rate (which I have trouble showing in this picture-- I would be eternally grateful if someone could make me a 3d plot showing the actual spiral effect with the red and blue spectral artifacts since I use this plot in my class and as is takes too long to explain)

Now review what happens to the spectrum when we decimate by 2 (the aliasing that MBaz was astutely concerned about). For me, the analogy of the aliasing that occurs when you do analog to digital conversion has helped me intuitively see this process (since decimation is a "digital to digital conversion". That comment was just for those already familiar with aliasing in the A/D process.):

Finally, review what happens when we do a "Zero-Insert" interpolation. Zero-Insert interpolation perfectly maintains the spectrum with no distortion other than the image which now becomes part of our primary digital spectrum. This is because a zero-insert is done by convolving our waveform with the unit-sample response (impulse response).

With that important background, let's follow the system given by the OP. First we decimate our input signal into even and odd samples. Note the rotation caused by the delay and then the aliasing caused by the decimation (as described above, to really visualize this properly the spectrum after the delay would be a spiral where the blue is exactly 180° out of phase at the $f_s/2$ point only):

The even and odd samples are then interpolated by two by doing a zero-insert, and the bottom path is advanced one sample (so note that doing that is not a causal process, we would end up with a "parasitic delay" of one sample to be able to actually implement this, following what was done in the chart). Note that the $z^{+1}$ creates a spiral in the opposite direction. The combined effect of [$z^{-1}$, decimate by 2, interpolate by 2, $z^{+1}$] will shift the phase 180° for every frequency component above $fs/2$ in our original spectrum, while passing every frequency component below $fs/2$ with 0° phase shift (as well as create the images in the center of the new spectrum).  Thus by combining (adding) the two paths, the original spectrum is recovered!

Note: You can also subtract the two spectrums to get a reversed spectrum! Notice that the adding and subtracting is exactly what you do with a 2pt DFT: $$\begin{bmatrix}1 & 1\\1 & -1\end{bmatrix}\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}y_1\\y_2\end{bmatrix}$$

We learn that aliasing can catastrophically ruin our spectrum and think that once the spectrums are contaminated we cannot separate the interference. That would certainly be the case here if we only had one of either path, but what is occurring is that our second path is keeping track of our aliasing in such a way that as long as we have both paths, we can constructively cancel out the interference.  This also gives insight into Multi-User Detection schemes.


Relationship between z-transform and DFT

For a) you're correct. For b), $x_1$ is a length $2N$ signal, and its DFT is given by

$$X_1[k]=\sum_{n=0}^{2N-1}x_1[n]e^{-j2\pi kn/2N}=\sum_{n=0}^{2N-1}x_1[n]e^{-j\pi kn/N}\tag{1}$$

With $x_1=x[n]+x[n-N]$ you get

$$\begin{align}X_1[k]&=\sum_{n=0}^{N-1}x[n]e^{-j\pi kn/N}+\sum_{n=N}^{2N-1}x[n-N]e^{-j\pi kn/N}\\&=X\left(e^{j\pi k/N}\right)+\sum_{n=0}^{N-1}x[n]e^{-jk\pi(n+N)/N}\\&=X\left(e^{j\pi k/N}\right)+e^{-jk\pi}X\left(e^{j\pi k/N}\right)\\&=X\left(e^{j\pi k/N}\right)(1+(-1)^k)=\begin{cases}2X\left(e^{j\pi k/N}\right),&k \text{ even}\\
0&k\text{ odd}\end{cases}\tag{2}\end{align}$$

Eq. $(2)$ shows that decimation in the frequency domain corresponds to aliasing in the time domain.


Vector Quantization in Hidden Markov Models

do multiple HMM classifiers (one for each class) share common Codebook
  (i.e common set of symbols) in Vector Quantization?

Yes, for example in large vocabulary speech recognition there is a big codebook of about 4000-10000 gaussians which are shared across 20-30k HMMs

  To which Codebook the test input vector belongs to (we need this to
  determine the symbol, in turn which is fed to the HMM for
  classification)?

In modern ASR systems codebook defines a set of GMM probability distributions, not the discrete symbols. So "belongs" doesn't really have sense, it's more a probability value. But during the classification a test vector is usually evaluated for the subset of all codebooks which is active for this stage.

Active codebooks are tracked across test vectors based on active HMMs and preceeding history. Whole HMM set is pruned in order to get active ones. Usually there are about 10k active HMMs on each step and about 1k codebooks.


This is just a 1D example that shows the idea (I can't seem to get scilab to cooperate with real images; at least not yet).

The top plot in blue shows the "image".  The "patch" is indicated by the red dots.  I've selected the location of the "patch" randomly throughout the "image".

The bottom plot shows the cross-correlation between the patch and the image (in green), and the red dot indicates the actual location (start position) of the patch.

The code below does no prewhitening because the "image" is generated using white noise.  Apart from the mean-correction, that is.

// 25803
clf;
NIM = 100;
IM = rand(1,NIM,"uniform");
IM = (IM - min(IM))/(max(IM) - min(IM))*255;
subplot(211)
plot(IM)

NPT = 24;
idx = ceil(rand(1,1,"uniform")*75);
PT = IM(idx + [0:NPT-1]);
plot(idx+[0:NPT-1],PT,'r.')
title('Original and patch')
NFFT = NIM + NPT - 1;

IMFFT = fft([IM-mean(IM),zeros(1,NPT-1)]);
PTFFT = fft([PT-mean(PT),zeros(1,NIM-1)]);

subplot(212)
CORR = ifft(IMFFT.*conj(PTFFT));
plot(CORR,'g')
plot(idx,CORR(idx),'r.')
title('Cross correlation and actual patch position')


Input and Output SNR Calculation in AM receiver input and output

yes, you are right. But there is misconception that figure of merit is the ratio of output to input SNR. It is actually the ratio of output SNR of a receiver to the output SNR of a baseband system (without modulation). Here, the output SNR of a baseband system is taken as a benchmark for judging noise performance of receiver. Hence,
              FoM = (output SNR)/(Si/NoW)


Books that explain DSP well to those not directly in engineering?

I am surprised no one has mentioned Richard Lyon's book - by far one of the BEST books out there on understanding digital signal processing in a very clear, concise, and methodological way.

Its excellence comes from that fact that he explains concepts to you in a very easy way to grasp, without loss of rigor or detail needed to get to the heart of various DSP concepts. He does not shy away from math, (which is good), but you are guaranteed to understand all of it.


Eye pattern construction and interpretation

It is the objective of the receiver to make the best estimate for each symbol as to what was transmitted. This is often done by ultimately determining a decision time in each sample (through timing recovery) on the waveform after it has been processed by the receiver (equalization and matched filtering) in which to sample the waveform and make a decision as to what symbol was transmitted.

The eye diagram superimposes the waveform for each symbol at the point in the receiver just prior to decision. The eye diagram readily shows the ability to make a decision between signal levels; so the ideal eye would show that there is a lot of margin both in the vertical and horizontal axis to allow for a minimum error rate (what we would call a wide eye opening). As the opening narrows, it would show that there is an increased sensitivity to noise since any added noise would cause a value to pass a decision threshold resulting in an error (this could be phase noise or jitter in the vertical direction and amplitude noise in the horizontal direction, or intersymbol interference which effects both directions).

I have this graphic pasted below that I created for one of my courses that may be useful in answering your question. This particular eye diagram is over two symbol durations. The waveform directly below the eye-diagram plot is replicated in the eye diagram as the darker blue lines, amidst the rest of the trajectories of the waveform from symbol to symbol over a much longer duration as given in the lighter blue lines:

Other examples of eye diagrams are shown below. The first is for QPSK, where we see the constellation on the left with the ideal sample locations in red, and the remaining trajectories from symbol to symbol everywhere else. On the right is the eye diagram formed from the real and imaginary components of the waveform.

And below is an eye diagram for a 16 QAM signal before and after the root-raised-cosine filter in the receiver. Here we clearly see the utility of the eye-diagram in that we see due to the proper filtering we have significantly increased our margin against noise. This would be similar to using the eye diagram to assess the effectiveness of equalization or to determine if certain filtering in our receiver is degrading our ability to obtain a low error rate due to the eye closing.


I believe the best strategy is to filter prior to calculating the magnitude. To see this easily, consider the low pass filtering as an averaging process and consider the noise as a zero mean Gaussian white noise. The distribution of the noise prior to the magnitude computation as zero-mean Gaussian will average to zero if the noise is white. The distribution of the magnitude of a Gaussian process is Rayleigh which has a maximum at the variance $\sigma$ of the noise distribution with an average of $\sigma \sqrt{\pi/2}$.  With a signal present the distribution is Ricean with similar considerations.

Filtering after will help in removing the variability of your result to the extent the system is stationary but will not reduce the error due to the noise; simple case in point is as above with noise only (and stationary), the result converges to zero; while after the magnitude the result will converge to a constant proportional to the noise level.

I suspect for both filters you will be limited by the dynamics involved in how much you can filter, which is ultimately the trade you will make with desired SNR of your result.

For more details on this specific to filtering before or after computing the magnitude (or equivalently before or after decision) see this great app note by Agilent:

http://my.ece.ucsb.edu/York/Yorklab/Useful%20Stuff/Tutorials/SpecAnalyzerMeas-Noise.pdf

There are applications where it is theoretically better, such as when the phase of the information is not preserved.

Here is a decent reference on both cases:

http://www.acfr.usyd.edu.au/pdfs/training/sensorSystems/10%20Detection%20of%20Signals%20in%20Noise.pdf

If you don't want to read or understand the reference materials in detail, or have little time- I suspect filtering prior to taking the magnitude will be better but you can easily test this by establishing SNR metrics on your result; just pay attention to rate of change as well with regard to how fast your system will be changing and how fast you want to be able to track that versus averaging out those changes.


Yes. Your code does compute the same thing as the formula from the paper.


Why are all real frequencies located on the Imaginary axis of the s-plane. I'd like to understand intuitively why this is true?

Take the complex (Laplace transform) variable $s$, and look at the basic exponential signal

$$x(t)=e^{st}\tag{1}$$

With $s=\alpha+j\omega$ we have

$$x(t)=e^{(\alpha+j\omega)t}=e^{\alpha t}e^{j\omega t}=e^{\alpha t}(\cos(\omega t)+j\sin(\omega t))\tag{2}$$

From $(2)$ you see that the real part $\alpha$ takes care of the exponential amplitude, whereas the imaginary part $\omega$ determines the frequency of the signal. If $\alpha=0$ (i.e., $s=j\omega$ is on the imaginary axis) you have an undamped sinusoidal signal with frequency $\omega$.