Reconstruction of BandPass Signal from Complex Envelope

Let’s get one thing straight. You cannot recover the original real and imaginary parts of a complex signal based only on the magnitude values of the complex signal. That’s equivalent to saying, “I have the magnitude of a complex number. How do I compute that complex number’s real and imaginary parts?” There are an infinite number of possible correct answers to that question!

However, as MBaz rightly said, you don’t have any “envelope” signal in your block diagrams. And by the way, there are no such things as “complex envelope” signals. All “envelope” signals are real-only. Having said all of that, you can learn everything about your two block diagrams by writing down the algebraic equations for all their various signals. (Your textbook should have done that for you.) Then pay attention to all the various trigonometric identities in your math reference book. It’s a mild pain to do all of that, but it’s terrifically educational.

And another thing, your two block diagrams are not the normal “complex modulator” and “complex demodulator” used in modern digital communications systems. In a normal digital comms system the real-valued gI(t)/2 and gQ(t)/2 would be added (or subtracted, I forget) to produce a new real-valued signal that’s transmitted via an antenna. Have a look at Figure 1 of the following PDF file: http://www.testequity.com/documents/pdf/keysight/complex-modulation-generation-wp.pdf


MATLAB: Applying Inverse Fast Fourier Transform to an array of real numbers

Since $\vert X_1 \vert$ is a real even function, its FFT (IFFT) is a real function. This is a basic property of Fourier transforms.

The power spectrum can be real but not even, in which case the ifft will give complex coefficients. In the example you linked, they start with a non-even real signal and then take the FFT. This will give complex values coefficients for the FFT. The answer then says that

  The problem is that we have only our log-energies, no phase
  information, so after applying the ifft we get complex-valued
  coefficients

So, after taking the FFT, the phase is discarded, and this loss of information makes it impossible to get back to the original real signal. (The Fourier transform is bijective). Once the phase is discarded you have a real but non-even function, the IFFT of which will yield complex coefficients.


What's the fastest way to get the phase angle and amplitude of a COSINUS 50hz signal?

Depending on how the signal will look like you can get away with only two or three samples:

$$
y(t) = A + B \cos(\omega t + \theta)
$$

If you know what the main value, $A_0$, would be, then only two samples would be required. Namely the cosine, shifted by a phase angle $\theta$, can be rewritten as a linear combination of a unshifted sine and cosine:

$$
y(t) = A_0 + A_1 \cos(\omega t) + B_1 \cos(\omega t) = A_0 + \sqrt{A_1^2 + B_1^2} \cos\left(\omega t - \tan^{-1}{\frac{A_1}{B_1}}\right)
$$

The gains can be solved for by solving the following linear system of equations:

$$
\begin{bmatrix}
1 & \sin \theta_0 & \cos \theta_0 \\
1 & \sin \theta_1 & \cos \theta_1 \\
1 & \sin \theta_2 & \cos \theta_2
\end{bmatrix}
\begin{bmatrix}
A_0 \\
A_1 \\
B_1
\end{bmatrix}
=
\begin{bmatrix}
y_0 \\
y_1 \\
y_2
\end{bmatrix}
$$

where $\theta_i=2\pi f(t_0+i\Delta t)$, with $f=50\ Hz$ and the $t_0+i\Delta t$ are the consecutive moments in time used for for the three measurements and $y_i$ is the measured signal at those moments in time. The values found for $A_0$, $A_1$ and $B_1$ can then be used to find the constants of the initial equation, namely:

$$
A = A_0
$$

$$
B = \sqrt{A_1^2 + B_1^2}
$$

$$
\theta = -\tan^{-1}{\frac{A_1}{B_1}}
$$

This should work for any different points in time, which are not a multiple of the period of the signal spaced apart, however this does assume that there is no noise present in the signal. If there is noise then you might want to consider using more points and performing a least squares method. Also a shorter sample time/higher sample rate would make these calculations more sensitive to the noise, but you will be able to take more samples in the same amount of time.


Determining the average probability of error

Note that
$$
\mathbb{P}(n_m+i_m<-1\cap i_m = \frac{-1}{2})
$$
is equal to
$$
\mathbb{P}(n_m < -0.5)\mathbb{P}(i_m = -0.5) = 0.25\mathbb{P}(n_m < -0.5)
$$
since $n_m$ and $i_m$ are independent.
Another way to obtain the solution is to add the the probability of these three disjoint events:

$n_m < -1.5$,
$i_m = 0$ and $-1.5 < n_m < -1$,
$i_m = -0.5$ and $-1 < n_m < -0.5$

In this case, we can add the three probabilities because the events are disjoint. To keep them disjoint, it is important to keep the noise ranges disjoint too; otherwise, you would be adding them multiple times.


How to prove that Dilution of a time domain signal results in duplication in Frequncy domain?

Unfortunately I can't come up with an intuitive explanation for this, but lemme try to prove this. During reading please be very careful with indexes as they can be confusing! I use zero-based indexing. And please correct me if I'm wrong.

Let  be our undiluted time domain of N elements.

So, diluted vector looks like this:

Frequency spectrum of x is given by .

Frequency spectrum of y, teh diluted x is . Note the length of Y equals 2N.

Now I'd like to prove that the right half of the diluted frequency domain is "copied" from the left half, that is . I guess it's sufficent to prove only the last equation, but the first two are good for grasping the idea.

By expanding a sum for Y[0] and Y[N] we get:

The fact that Y[0] and Y[N] are equal can easily be seen by noticing that  because of it's periodic nature- see Euler's relation: . Later we will exploit this property of exp(j*phi) to "move" exponent back or forward by 2nPi.

The next (unnecessary, but good to examine) step is to prove that Y[1] = Y[N+1]. Let's expand:

Because Y[1] = Y[N+1] must hold regardless of x values, the only way that it can do so is when exponents next to proper xn's are equal. Let's prove that . If we can do this, then when a=0, we prove that exponents next to x0 in Y[1] and Y[N+1] are equal. When a=1, exponents near x1 are proven to be equal and so on, until a=N-1.

As you can see, it was not hard. It's more like showing than proving. If you are wondering why did 2pi*a dissapear, look again at the Euler's relation and maybe rewrite the middle step of the proof in cos - isin form. It should be clearer now that you can add or subtract 2pi from argument of the sine or cosine function without affecting it's value.

Okay, let's move on to the next, (the only necessary) part. Y[k]=Y[k+N].

Like the last time, let's show equality of exponentials near xn's:

The idea is the same, only the last point has the least numbers and the most variables, making it the most general case. It's general enough that I think it's the end of the proof.

Hope it helps :)


The answers to this are now on my site, which I'm not going to keep updating into this answer when I modify it


SineWave generation and BandPass filter : Audio

As someone suggested in the comments, filtering a sine wave is just a change in amplitude and possibly phase, although you will have some edge effects at the beginning. I assume this is not simply an excersize in changing amplitudes and that you really need a real filter.

I don't know of an open source library in Java, although this book contains some simple source for a few filters (not sure it has a bandpass, but lowpass + highpass = bandpass).

I don't recommend the suggestion in the comments of using an FFT, although if you can perform a transform on your entire input that is an option, and one that many people find very intuitive. Keep in mind it's an O(n log n) operation, so it may not be feasible for large or even moderately-sized audio files. I explain here more on why EQ is usually not done in the frequency domain. Some people design their filters in the frequency domain and use the IFFT to create a filter, but this results in a design that is inferior in many ways. There are some rules of thumb that make this better, but there are better techniques, so this is not a well developed art.

You haven't specified your filter in any meaningful way, but assuming a basic biquad will suffice, checkout RBJ's audio EQ cookbook. If you need help implementing those filters, I've written a tutorial.

Finally, someone else on StackOverflow just asked about a bandpass filter in Java, and offered some not yet working code. Their design is based on a blackman windowed sinc filter. This is a high-quality FIR filter.


A short answer to your question is "NO". To my best knowledge, there is no such a smart end-to-end OCR system that can handle various OCR tasks. However, there are OCR products could reasonably handle some tasks under controlled environments.

You seem to be confused of different OCR related tasks and you erroneously believe these tasks are quite similar ( and thus something works well on document could also work well on reading signs ).

Here lists several levels of OCR tasks:

level 0: single character recognition
level 1: word recognition level
level 2: sentence recognition
level 3: paragraph recognition
level 4: document recognition
level 5: OCR in the wild

Depending on the source of characters, you may classify texts into 1) machine printed texts ; 2) handwritten texts; and 3) scene text. All these are OCR related tasks, and thus can all be called "OCR".

These tasks are clearly of different difficulties. Basically, the higher the level the more difficult it is. More importantly, all these tasks are NOT just about recognition. Sometimes pre- and post- processing are both important. For example, it is well known that the line segmentation problem in handwritten document recognition is very challenging when text lines are cursive, or touched, or mixed of texts from different directions, etc. If you cannot first segment text line correctly, how to recognize them? For another example, it is also well known that deciding the reading order of multi-column document mixed with pictures is challenging. Even if you recognize every sentence correctly, you may still make mistakes in organizing different paragraphs. In your interested scene text recognition tasks, text localization itself is challenging in the sense that texts can be any font in anywhere with any direction and perspective change.

For these reasons, a good end-to-end OCR engine is at least a mix of image processing, pattern recognition, and natural language processing, and thus this is a quite big system requiring knowledge from all these parts.

I put my comments of 'prior knowledge' below.
To discuss 'prior knowledge', it is better to clarify its definition.
- For recognition, it often means there is prior distribution to tell you how likely a character is expected to see. If you said you have no prior knowledge, then you believe it is equally likely to see all possible characters. Think about how many valid characters in all languages. This number is big, and the number of valid words in all languages is even bigger.
- For image processing, it often means you have (explicit or implicit) assumptions about a given image. Such assumptions could be 'one-column document', 'car plate', 'taken by iphone', etc. If you say I have no prior knowledge about where a car plate is, then this simply means you treat the possibility of seeing a car plate is equally likely for any image patch size and location. Think about how many possible image patches you will have.

It is clear that doing OCR without prior knowledge is the worst thing we can try. In other words, we could do better when we have more prior knowledge for our task. In practice, people are using various prior knowledge, but often in an implicit way, and thus they even donot notice it.

One typical example is training data. Because your training data is not the entire universe, any rules, distributions, characteristics etc you learnt from this dataset may not necessarily apply to the rest of the world. However, you are still use your trained models to do various tasks. Your prediction is based on the prior knowledge that you learnt from data.

Back to your 'car plate' OCR problem. We know a car plate is of a rectangular shape under some perspective transformation. We know it is very likely to locate at the lower part of a car. We know it cannot be to too thin nor too small. Such prior knowledge helps quickly locate a candidate plate while avoiding a lot false alarms. When you recognize a plate, there is more noticeable prior knowledge. If you know only digits are allowed in your plate, then you cut your solution space to the digit space. Even if you see something more similar to letter 'g' than digit '9', you can safely say "this must be '9'" because of your prior knowledge. For the same reason, if you know plate are valid, you can safely thrown away any candidate solution that is invalid. Even 'not that useful prior knowledge' helps. For example, if you know no plate in our world has more than N characters, you know a plate candidate you are recognizing is NOT a valid plate because its length is above N.

If you donot have any prior knowledge, what happens? For me, it is similar to ask me find some object named '%#$$' on an image, which I completely have no idea what it is. And this is an impossible task. However, if I have the knowledge that '%#$#' is the ciphertext of 'ball', then I can easily do this task.

I hope that I convince you that prior knowledge is useful and closely associated with each processing step, and we cannot do OCR without prior knowledge. Indeed, this is the exact reason why we need to combine models from different aspects for OCR.

Donot be fooled by the name of OCR "optical character recognition". It is not just character recognition: we are actually reading words, sentences, paragraphs from an image that is taken by someone using some particular device... and thus we need prior knowledge at many different levels.


Considering the FFT of Real & Complex Signals

There's nothing wrong here -  complex sinusoids like your signal really have only one peak in frequency domain! This is the fundamental idea of why we use the Fourier transform for periodic (even complex) signals.

You can think of it this way:

the cosine has two peaks, one at +f, the other at -f. That's because Euler's formula actually says $\cos x = \frac12\left(e^{ix}+e^{-ix}\right)$. You can see the two complex sinusoids that lead to your two peaks.

Now, doing the same for the sine, we see that $i \cdot\sin x =i\frac1{2i}\left(e^{ix}-e^{-ix}\right)$. That means, the sine kills one of the two peaks of the cosine, if you will so.

Everything's alright :)


how a signal for which calculate RMS can be filtered and have a fast settling time?

I believe you are running into stability issues in an output power control loop design.
See below a diagram of similar power control loops that I have implemented, where for stability reasons any filtering in the loop is minimized and only done with the loop filter itself which is designed for stability. The noise you are trying to filter gets filtered by the loop bandwidth itself.
The diagram below is very simple to implement (first order type 1 loop), where the accumulation of the derived error sets the loop bandwidth by adjusting K.  Make K very small for a long averaging time and better filtering but slower response.

Ultimately if the trade of filtering versus response time with a simple implementation like shown above is insufficient, you can change the loop filter to be a PID loop. This would more complicated and require more understanding of control loop design and a more detailed loop model but give you better performance in noise filtering versus response time. Essentially the accumulator (I) shown above would also be summed with a proportional path (P) and a derivative path (D) all with different gain weights and followed by a second accumulator for a second order type 2 PID loop.
Once a loop model is created, the modeling of the response time and effective filtering is very straightforward with tools such as Matlab, Octave or Python. As an example of the design process, below shows the loop model for the Power Control Loop given by the block diagram above along with the modeled response times versus K:

The above case results in an open loop gain given by:
$$G_{OL}(z) = \frac{0.83155K}{z(z-1)}$$
Below demonstrates given the Matlab and Octave code snippets how simple it is once the model is established to create the step response for the loop. Note the overshoot indicative of a higher order loop- this is due to the parasitic $z^{-1}$ in the digital loop implementation, and with a high enough K the loop will be unstable (unlike a true first order type 1 loop). Similarly the noise bandwidth, etc is easily established and optimized.


Image processing: using PCA for 3d object recognition

Principle Component Analysis by itself is not going to provide rotational or scale invariance. You still have to provide templates for rotational and scaled models for image/template matching. What PCA does do is reduce the dimensionality of the search space  by converting a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components.

Now since number of principle components <= number of components, you save on the search space.

What the lecture tries to tell you is with the eigenvectors, you get an efficient representation of the scales of the template which are easier to manipulate(a single scale factor) then manipulating the entire image everytime.


How would I spectrally rotate/invert an audio signal (Matlab or Python)?

$$ y[n] = (-1)^n \cdot x[n] $$
This will flip the spectrum in linear frequency. Easy to do, easy to recover.


Signal integrity is a full area of electrical engineering and signal engineering.
A few examples :

To protect a propagating signal in free-space you can emit it louder (=more transmitting power, hence increasing signal to noise ratio (SNR)).
Numerical transmissions often include error correcting codes.
To protect propagating signals in cables, most of them are electromagnetically shielded.

I may detail my answer according to what you want to know exactly.


Hermitian FFT vs discrete Fourier transform. Python's $\tt hfft$ vs $\tt fft$

You have to understand that fft is the general function, which always works. hfft, fftr and their derivatives are optimized for special signal constellations and should either be faster or more accurate (if this is at all necessary with modern floating point accuracy).

The hfft function is used, when you expect the spectrum to be real-valued. This happens, if your input signal has even symmetry. I doubt that your audio signal does have this property. Have you been looking for fftr (i.e. fft for a real-valued signal)?

In case your input signal has the even symmetry property, the output of fft and hfft should be exactly the same (up to numerical residual errors). Also, hfft should be a bit faster in this case (otherwise, it would make no sense to use it).


Significance of outer product in structure tensor / second moment matrix

I have been looking again into this question. The tensor product doesn't really have such a special meaning here.
According to Nonlinear Structure Tensors:
"Although this tensor product contains no more information than the gradient itself, it has the advantage that it can be smoothed without cancellation effects in areas where gradients have opposite signs."
I think it is then just similar to a covariance matrix i.e. $M_{11}, M_{22}$ is the variance of the $x$ and $y$ direction. $M_{12}$ and $M_{21}$ is the covariance of $x,y$ and $y,x$. It is just important to see that we sum each component of the matrix (see Wikipedia). The outer product is then $E[X]E[X]^T$. So we consider an image $I : \Omega \to \mathbb{R}^n$ as a vector field and look for local variations $||dI||^2 = dI^TdI$. The outer product appears because all $4$ directions are compared.
Another explanation comes from diffusion equations (i.e. movement from a region of higher concentration to a region of lower concentration). Without the outer product, we would only look at the gradient at a point. The diffusion is the same in every direction (isotropic). This is called (nonlinear) isotropic diffusion:
$$\frac{\partial u}{\partial t} = \text{div}\left(g\left(|\nabla u|^2\right)\nabla u\right)$$
When $g(\cdot) = 1$ ("linear"), this is the regular heat equation because  $\text{div}\left(\nabla u\right) = \nabla^2 u$. The solution would be a Gauss filter, which blurs everything. Wikipedia calls nonlinear isotropic diffusion anisotropic, even though the scalar product $|\nabla u|^2$ is applied (instead of tensor product).
The second possibility is anisotropic diffusion. The PDE is
$$\frac{\partial u}{\partial t} = \text{div}\left(g\left(\nabla u\nabla u^T\right)\nabla u\right)$$
$\nabla u\nabla u^T$ is the structure/diffusion tensor. The diffusion is different in every direction. Diffusion across edges is prevented. When we only look at the diagonal of the matrix, the diffusion is isotropic again. Some more information can be found here and here.


why 2 vectors are not the same in matlab?

Neither 0.11 nor 1E-4 can be exactly represented in floating point arithmetic. And the product of 2 numbers will introduce an extra error.

Testing for equality in floating point arithmetic is therefore more complex than it seems.


How to verify and plot expression for sampled signal is correct?

The best way is to let the program to sample the function for you.
You can define a function. And you don't even need a separate program, you can use an anonymous function, @(t) expr for MATLAB/Octave (That seems to be the program you are using).
  x_a = @(t) A1*cos(2*pi*f1.*t + th1) + A2*cos(2*pi*f2.*t + th2);

Then you could plot with
  plot(t, x_a(t), '-', n/fs, x_a(n/fs), 'o');

since you are ploting x, x_a(x) there is nothing that could go wrong in the sampling, the time values you pass will be replaced automatically in your expression.
If you want the expression you can use symbolic computation, however when I need to do so I use wxmaxima or sympy


Disjoint ROCs of input and system function

If the ROCs don't overlap, the corresponding convolution sum doesn't converge, hence the output is undefined. Note that in your question it is assumed that the $\mathcal{Z}$-transforms of the input as well as of the system's impulse response both exist.
You can construct examples for which neither the $\mathcal{Z}$-transform of the input nor of the output exist, yet the input and output signals are perfectly defined. E.g., the input sequence $x[n]=A\cos(\omega_0n)$ has no $\mathcal{Z}$-transform, neither has the corresponding output sequence $y[n]=B\cos(\omega_0n+\phi)$, where $B$ and $\phi$ depend on the system's frequency response. [In this example it is assumed that the ROC of the system's transfer function $H(z)$ includes the unit circle, i.e. the system is stable.]


Finding Correlation response in fourier domain

From my understanding, the more spread out, or the less dense the power spectrum is over a set number of frequency bins, the less correlated the two signals are, and vice versa.

I'm afraid you have this backward. Narrow peaks in the original (spatial) domain correspond to strong cross-correlations. Narrow peaks in the spatial domain will be relatively spread out in the Fourier domain. The results of the cross-correlation are usually interpreted in the spatial domain (rather than the Fourier domain).

Transform the image and the (zero-padded) template, multiply in the Fourier domain, and then transform back to the spatial domain. Peaks in the spatial domain correspond to (the possibility of) matches between the template and the image.


Estimating velocity from known position and acceleration

You would not make a constant-acceleration assumption in this case. You would typically do that if you had no means of measuring the system's acceleration, but you say that it is observable for your case. The most obvious way to model this system would be using the state vector

$$
\mathbf{x_k} = \left [
\begin{array}{c}x_k \\
\dot{x}_k \\
\ddot{x}_k
\end{array} \right ] = \left [
\begin{array}{c}x_k \\
v_k \\
a_k
\end{array} \right ]
$$

where $x_k$ is the position, $v_k$ is the velocity, and $a_k$ is the acceleration of the system at time instant $k$. Since you say that you can measure the position and acceleration of the system, the measurement vector $\mathbf{z_k}$ would be:

$$
\mathbf{z_k} = \mathbf{Hx_k} + \mathbf{v_k}
$$

with

$$
\mathbf{H} = \left[
\begin{array}{cc}
1 && 0 && 0 \\
0 && 0 && 1
\end{array}
\right]
$$

resulting in the measurement model:

$$
\mathbf{z_k} = \left[
\begin{array}{c}
x_k \\
a_k
\end{array}
\right] + \mathbf{v_k}
$$

where $\mathbf{v_k}$ is the (Gaussian) noise inherent in the measurement. Now, unless you are adding some known force to the system that would affect its dynamics, then the input vector $\mathbf{u_k} = \mathbf{0}$. See this recent question for more detailed discussion of the various terms in the Kalman model, which you need to know, and which you don't.

Edit: Regarding the below comment on how to handle the state transition matrix for the acceleration term: there are a couple different ways to handle that. You could assume a zero-jerk (jerk is the time derivative of acceleration) model; this is equivalent to assuming that $a_{k+1}=a_k$. Or, you could add a fourth element to your state vector:

$$
\mathbf{x_k} = \left [
\begin{array}{c}x_k \\
\dot{x}_k \\
\ddot{x}_k \\
\dddot{x_k}
\end{array} \right ] = \left [
\begin{array}{c}x_k \\
v_k \\
a_k \\
j_k
\end{array} \right ]
$$

where $j_k$ is the system's jerk, or time rate of change of acceleration, at time instant $k$. Since you aren't measuring jerk directly, you would have a measurement matrix of:

$$
\mathbf{H} = \left[
\begin{array}{cc}
1 && 0 && 0 && 0\\
0 && 0 && 1 && 0
\end{array}
\right]
$$

And the state transition matrix would become:

$$
\mathbf{F} = \left[
\begin{array}{cccc}
1 && dt && \frac{1}{2}dt^2 && \frac{1}{6} dt^3 \\
0 && 1 && dt && \frac{1}{2}dt^2 \\
0 && 0 && 1 && dt \\
0 && 0 && 0 && 1
\end{array}
\right]
$$

This is a "constant-jerk" model (similar to constant-velocity or constant-acceleration models you might see in simple examples); your state transition matrix makes an implicit assumption that the jerk is constant for all values of k. This is not likely to be true. To handle that aspect of the problem, you would introduce a process noise term in the jerk component of the state transition equation.

Qualitatively, this allows you to express that you aren't sure how the jerk term is going to change from time step to time step, but that you expect the changes to be random with some Gaussian distribution. This is a tool often used to "tune" adaptive filter models like this one until you find a set of parameters that works well for your application.


Given measurements
$$\begin{align}
Z_1 &= x + N_1\\
Z_2 &= x + N_2
\end{align}$$
where $N_1$ and $N_2$ are independent zero-mean Gaussian random variables
with variances $\sigma_1^2$ and $\sigma_2^2$ respectively, it can be shown
that the minimum-mean-square-error (MMSE) estimate of $x$ in terms of $Z_1$ and
$Z_2$ is a linear function of $Z_1$ and $Z_2$.  Note that this linearity
is a consequence of the Gaussian assumption; for other kinds of distributions
that might be assumed for the noise variables $N_1$ and $N_2$, the MMSE
estimate is usually not a linear function.

So, let us consider the estimate $\hat{X} = k_0 + k_1Z_1 + k_2Z_2$ for the
unknown parameter $x$.  The mean-square error is thus
$$\begin{align}
E[(\hat{X}-x)^2] &=  E[(k_0 + k_1Z_1 + k_2Z_2 - x)^2]\\
&= E[(k_0 + (k_1+k_2-1)x + k_1N_1+k_2N_2)^2]\\
&= \operatorname{var}(k_1N_1+k_2N_2) + (k_0 + (k_1+k_2-1)x)^2
&{\scriptstyle{\text{because}~E[(Y+a)^2=\operatorname{var}(Y)+(\mu_Y+a)^2}}\\
&\geq \operatorname{var}(k_1N_1+k_2N_2)
\end{align}$$
where equality holds in that last inequality if we choose $k_0=0$ and $k_1+k_2=1$. With a slight change of notation to incorporate these
constraints into our calculations,
let us write
$$\begin{align}
E[(\hat{X}-x)^2] &= \operatorname{var}(\alpha N_1+(1-\alpha)N_2)\\
&= \alpha^2\sigma_1^2 + (1-\alpha)^2 \sigma_2^2
&{\scriptstyle{\text{because}~N_1, N_2 ~\text{are independent random variables}}}\\
&= \alpha^2(\sigma_1^2 + \sigma_2^2) -2\alpha\sigma_2^2 +  \sigma_2^2
\end{align}$$
where it is easily shown that the minimum value of that quadratic function
of $\alpha$ occurs when we choose
$$\alpha = \frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2} = k_1, \qquad
1-\alpha = \frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2} = k_2,$$
that is, when the two measurements are afflicted with Gaussian
noise of different variances (say with $\sigma_2^2 > \sigma_1^2$),
then the measurement with the smaller variance ($N_1$ in this
example) is given greater weight
$k_1=\frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2}$
than the measurement with the larger variance ($N_2$ in this case)
which has smaller weight $k_2 = \frac{\sigma_1^2}{\sigma_1^2 + \sigma_2^2}$.

When $\sigma_2^2 = \sigma_1^2$, the measurements are given equal weight
$\frac 12$ which is, according to the OP's insistence, the weight that
must be
assigned in all circumstances so as to minimize the mean-square error
because that is the only weighting that makes sense to the OP.


Instantaneous phase calculation problems

One way to overcome this problem is phase unwrapping, which is also implemented by the Matlab function unwrap. But I believe that the best option for obtaining the derivative of the phase is to avoid a direct computation, but instead start out from the complex function itself. To simplify notation, let's use $U$ (without arguments):

$$U=re^{j\phi}=a+ib$$
where $r$ is the magnitude and $\phi$ is the phase. The derivative of $U$ is given by

$$U'=r'e^{j\phi}+rj\phi'e^{j\phi}\\
\frac{U'}{U}=\frac{r'}{r}+j\phi'$$

Since $r'/r$ is real-valued, we get

$$\phi'=\text{Im}\left\{\frac{U'}{U}\right\}=\frac{ab'-a'b}{a^2+b^2}\tag{1}$$

where $a$ and $b$ are the real and imaginary parts of $U$, respectively. Using (1) instead of a direct derivative of the phase should solve the problem.


Basically, multipath propagation causes time dispersion (hence causing a delay spread ($\sigma_{\tau}$)) and that dispersion influences the fading characteristics of the wireless channel.So delay spread is related to the fading characteristics of a channel and not really to the path loss. Particularly, if $\sigma_{\tau}$ is greater than the symbol period ($T_s$) then the wireless channel is said to be frequency selective, meaning that a signal will interfere with itself causing additional fading.

You should also note that the rms delay spread is not the maximum delay that your signal experiences as I think you may asssume, $\sigma_{\tau}$ can be calculated as follows:

$\sigma_{\tau} = \sqrt{\overline{\tau ^2} - \overline{\tau} ^2}$

where

$$\overline{\tau} =  \frac{\Sigma P(\tau _i) \tau _i}{\Sigma P(\tau _i)},$$
$$\overline{\tau ^2} =  \frac{\Sigma P(\tau _i) \tau _i ^2}{\Sigma P(\tau _i)}$$
and $P(\tau _i)$ is the power of the component with delay $\tau _i$


Are there frequency detection methods that are subject to less latency than the Goertzel algorithm?

I might be mentioning this too often, but if you need to detect frequencies finer than $f_\text{sample}/N$, then superresolution frequency estimators might be the way to go.

In your case, an algorithm called MUSIC would sound promising. $N$ would in that case be the amount of samples you use to estimate the autocorrelation matrix that algorithm works on, and can be as small as $n_\text{tones}+1$, but realistically, don't go too low, and average the estimate a bit.


How to convert my transfer function to the frequency domain

Since this is not an LTI system, we cannot have $H(\omega) = K\Theta_1(\omega)$, where $K$ is a scalar complex number. But as OP mentioned in the comment he has knowledge of $\Theta_1(\omega)$, we can compute $\Theta_2(\omega)$ and $\Theta_3(\omega)$. As $x^2(t)=x(t)\times x(t)$,
$$
\Theta_2(\omega)= \Theta_1(\omega)*\Theta_1(\omega) =\int_{-\infty}^{+\infty} \Theta_1(\alpha)\Theta_1(\omega-\alpha)d\alpha
$$
where $*$ is the convolution operation. Similarly, $x^3(t)=x(t)\times x(t)\times x(t)$
$$
\Theta_3(\omega)= \Theta_1(\omega)*\Theta_1(\omega)*\Theta_1(\omega)=\Theta_1(\omega)*\Theta_2(\omega)
$$
So $H(\omega)=B\Theta_1(\omega)+a\Theta_2(\omega)+b\Theta_3(\omega)$, where $\Theta_2$ and $\Theta_3$ are computed from $\Theta_1$ above.


Is the discrete Gaussian kernel an eigenfunction of the DFT?

Since the DFT is representable by multiplication with the Fourier matrix, your question is equivalent to asking what the eigenvectors of the Fourier matrix are.

Actually, Wikipedia provides the answer (http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Eigenvalues_and_eigenvectors).

However, since the eigenvalues ($1, -1, i, -i$) are not simple, the eigenvectors are not unique (i.e. linear combinations are also eigenvectors). Also no simple closed formula exists.

One formula for an eigenvector close to what you ask is provided by Wikipedia

$$F_m = \sum_{k = -\infty}^\infty \exp\left(-\frac{\pi\cdot(m+N\cdot k)^2}{N}\right) \quad\quad\quad m = 0,\ldots,N-1$$

So concluding, the Gaussian function itself is not an eigenvector, but an infinite sum of Gaussians.
The infinite sum can probably be interpreted as equivalent to the discretization of the frequency and time domain when we going from the FT to the DFT. So it is not as simple, as just truncating the discrete Gaussian.


In that range it is guaranteed to converge. It doesn't mean it will necesseraly won't converge for higher values. If you want deeper understanding you can read about the step size in Convex Optimization context where there the step size related to the Lipschitz Constant of the function (Which matches the eigen value for Quadratic functions).

If you share the problem itself from the book, we'll be able to solve the problem.


This reference could be useful for blind equalization of CiC Signals


Determining The Sparsity or Density of a Chroma Vector

I am considering a similar problem as well and I found Gini index as one of the solution: Gini Index as Sparsity Measure for Signal Reconstruction from Compressive Samples. Please have a look at the definition in page 3. It is normalized in the range from 0 to 1 as the least sparse or the most sparse the vector is, respectively.

The code below is an implementation in Python, though I am not sure about returning 0 if the input is a zero vector:

def gini(v):
    v = np.array(v)
    sort_v = sort(abs(v))
    norm_v = sum(sort_v)
    if norm_v==0:
        return(0)
    else:
        sig = 0
        for k in arange(len(sort_v))+1:
            sig+=(sort_v[k-1]/(norm_v))*((len(sort_v)-k+0.5)/len(sort_v))
        return(1-2*sig)


Since you already know your background image beforehand, it should be simple. I done many background subtraction before. Here is an example of removal of background by comparison.

    if (newBitmap.Width == backgroundBitmap.Width && newBitmap.Height == backgroundBitmap.Height)
    {
        for(int x = 0; x < newBitmap.Width; x++)
        {
            for (int y = 0; y < newBitmap.Height; y++)
            {
                Color firstPixel = newBitmap.GetPixel(x, y);
                Color SecondPixel = backgroundBitmap.GetPixel(x, y);

                int Totalfirst = (int)((firstPixel.R * .3) + (firstPixel.G * .59) + (firstPixel.B * .11));
                int Totalsecond = (int)((SecondPixel.R * .3) + (SecondPixel.G * .59) + (SecondPixel.B * .11));
                if (Totalfirst < Totalsecond )
                {

                    Color WhiteColor = Color.FromArgb(0, 0, 0);
                    newBitmap.SetPixel(x, y, WhiteColor);
                }
                else
                {
                    Color SameColor = Color.FromArgb(firstPixel.R, firstPixel.G, firstPixel.B);
                    newBitmap.SetPixel(x, y, SameColor);
                }

            }
        }
    }
    FinalImage = newBitmap; // after substaction of back ground.
}

This code is just a copy and paste from one of the program I did back then. Note that newBitmap is the image of the object and the background, your second picture. While the backgroundBitmap is well, the background.

Use this as your reference on how your code should look like for this segment. Should be good enough.

Also, I am new to this website(my first post here), but I think for image processing questions, you can post it to https://stackoverflow.com/questions next time, there will be more people to help you.

If this method doesn't work, let me know, cause I still got some other methods. But from you question, this should be sufficient.


What is the benefit of FFT when we're already sampling at an acceptable rate?

Not all periodic signals are sinusoids.  Items 1 and 3 should not assume such.  Not all sampling is done just for reconstruction.  Look up other needs or uses for sampled signals.
A linear time-invariant system can behave differently to different frequency sinusoidal inputs (as provided by common solutions to certain low order differential equations), so it may be easier to figure how that system responds to other types of far more complicated inputs by breaking them down into the sum of simpler sinusoidal ones.
A DFT or FFT won't break down a pure (unmodulated) sinusoid that is integer periodic in aperture into other sinusoids (other than that a pure sinusoid of an arbitrary phase can be broken down into component cosine and sine sinusoids of that same frequency).  But an FFT can "break down" other types of signals.


Combining Bandpass filters on 1D signal to obtain better SNR?

If you know that the interference is only in a given frequency range, you might use a notch filter to filter out this interference. If it is in several frequency ranges, you can either filter with different filters in parallel and use the output that works best for you, or filter the signal serially through all filters (but this might distort the signal shape a lot, depending on the Q-factors and number of filters). Here I show an example with a single filter and a known interferer frequency.

clear all;

Fs = 1000; % 100Hz sampling frequency
t = -5:1/Fs:5;

f_sig = 0.5;

f_intf = 15;
A_intf = 0.05;

% Generate two exemplary signals and pollute them with interference
signal1 = double(sin(2*pi*f_sig*t) > 0);
signal2 = cumsum(signal1 - 0.5) / (Fs/f_sig);
noise1 = A_intf * sin(2*pi*(f_intf)*t);
noise2 = A_intf * sin(2*pi*(f_intf)*t);

received1 = signal1 + noise1;
received2 = signal2 + noise2;

% create a notch filter to filter out this intreference.
w0 = f_intf/(Fs/2);
bw = w0/1; % 1 is the Q-factor here, you can play with it to see the effect
[num, den] = iirnotch(w0, bw);

filtered1 = filter(num, den, received1);
filtered2 = filter(num, den, received2);

subplot(2,2,1);
plot(t, received1)
title('received 1');
grid;

subplot(2,2,2);
plot(t, received2);
title('received2');

subplot(2,2,3);
plot(filtered1);
title('filtered 1');

subplot(2,2,4);
plot(filtered2);
title('filtered2');

Edit: I have ported your code to python, using the correct filter. Here's the filter design function that was taken from scipy v0.19 (not yet officially released), you can find the source here: Scipy Source

def _design_notch_peak_filter(w0, Q, ftype):
    """
    Design notch or peak digital filter.
    Parameters
    ----------
    w0 : float
        Normalized frequency to remove from a signal. It is a
        scalar that must satisfy  ``0 < w0 < 1``, with ``w0 = 1``
        corresponding to half of the sampling frequency.
    Q : float
        Quality factor. Dimensionless parameter that characterizes
        notch filter -3 dB bandwidth ``bw`` relative to its center
        frequency, ``Q = w0/bw``.
    ftype : str
        The type of IIR filter to design:
            - notch filter : ``notch``
            - peak filter  : ``peak``
    Returns
    -------
    b, a : ndarray, ndarray
        Numerator (``b``) and denominator (``a``) polynomials
        of the IIR filter.
    """

    # Guarantee that the inputs are floats
    w0 = float(w0)
    Q = float(Q)

    # Checks if w0 is within the range
    if w0 > 1.0 or w0 < 0.0:
        raise ValueError("w0 should be such that 0 < w0 < 1")

    # Get bandwidth
    bw = w0/Q

    # Normalize inputs
    bw = bw*np.pi
    w0 = w0*np.pi

    # Compute -3dB atenuation
    gb = 1/np.sqrt(2)

    if ftype == "notch":
        # Compute beta: formula 11.3.4 (p.575) from reference [1]
        beta = (np.sqrt(1.0-gb**2.0)/gb)*np.tan(bw/2.0)
    elif ftype == "peak":
        # Compute beta: formula 11.3.19 (p.579) from reference [1]
        beta = (gb/np.sqrt(1.0-gb**2.0))*np.tan(bw/2.0)
    else:
        raise ValueError("Unknown ftype.")

    # Compute gain: formula 11.3.6 (p.575) from reference [1]
    gain = 1.0/(1.0+beta)

    # Compute numerator b and denominator a
    # formulas 11.3.7 (p.575) and 11.3.21 (p.579)
    # from reference [1]
    if ftype == "notch":
        b = gain*np.array([1.0, -2.0*np.cos(w0), 1.0])
    else:
        b = (1.0-gain)*np.array([1.0, 0.0, -1.0])
    a = np.array([1.0, -2.0*gain*np.cos(w0), (2.0*gain-1.0)])

    return b, a

Here's your corrected code using this function:

length = 1200
frequency = .00259
interference = 0.0537
interference_amplitude = 0.2

def iirfilter_filter(data):
        nyq = 1/2.
        bottom = interference/nyq
        top = interference*10/nyq
        #b, a = iirfilter(10, [bottom,top], rs=60, btype="bandstop" , analog=False , ftype="cheby2" )
        b, a = _design_notch_peak_filter(interference*2, 1, 'notch')
        return lfilter(b,a, data)

time_series = np.linspace(1, length, length )

sine_wave = np.sin(2 * np.pi * frequency * time_series )
square_wave = (sine_wave / abs(sine_wave))

sine_wave = sine_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) - 1
square_wave = square_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) + 1

noisy_sine_wave = sine_wave + ( numpy.random.randn(length) * interference_amplitude)
noisy_square_wave = square_wave + ( numpy.random.randn(length) * interference_amplitude)

a = plt.subplot(2,2,1)
b = plt.subplot(2,2,2)
c = plt.subplot(2,2,3)
d = plt.subplot(2,2,4)

a.plot( sine_wave )
a.plot( square_wave )
b.plot( iirfilter_filter( sine_wave ) )
b.plot( iirfilter_filter( square_wave ) )
c.plot( noisy_sine_wave )
c.plot( noisy_square_wave )
d.plot( iirfilter_filter( noisy_sine_wave ) )
d.plot( iirfilter_filter( noisy_square_wave ) )

plt.show()

And the output:

Note that you cannot filter the white noise from randn (corrected from rand), as it contains all frequencies and the notch filter only removes a single frequency.


Your FIR state space representation seems to be doing too much.

The way I would write it is to have the process noise is $\epsilon(t)$ as your input, and your state as two time-delayed copies of it:

$$x(t+1) = \left[ \begin{array}{c}
\epsilon(t+1)\\
\epsilon(t) \\
\epsilon(t-1)
\end{array} \right] = \left[ \begin{array}{ccc}
0 & 0 & 0\\
1 & 0 & 0\\
0 & 1 & 0
\end{array} \right] x(t) + \left[\begin{array}{c}
1\\
0\\
0
\end{array} \right] \epsilon(t + 1)
$$

then your output equation is just:
$$
y(t) = \left[ \begin{array}{ccc}
1 & h_1 & h_2 \end{array} \right] x(t) + v(t)
$$

Your IIR representation is too confused for me to make head or tail out of. Can you clarify that you have the right information there?

Other questions:

What is $w(t)$? It does not appear to take part in any of your equations?
What are the parameters for $u(t)$?


Here I expected $y(n)$ is to be computed by convolving $x(n)$ with $h(n)$, but in the equation given by Wikipedia it is shown as a matrix multiplication
  $y(n) = h^H(n).x(n)$. Are these two operations(convolution and matrix multiplication) same here?.

The system is an FIR system, so the vector multiplication here is equivalent to convolution --- for one time instant output. Note that $x(n)$ is a vector o the current and last $N$ samples of $x(n)$.

  Why do we need a Hermitian transpose rather than normal transpose in the above equation?

Because the signal is assumed to be complex.

  Cost function is defined as $E\{|e(n)|^2\}$. Why do we need expectation operator here? Why not normal Mean Square Error(MSE)?

Because the error is assumed to be stochastic. We don't want to minimize for one realization, we want to minimize for all realizations.

  Does LMS filter converge to all type of input signal?

You don't believe the proof? :-)

Provided the assumptions made by the proof are met, then the LMS algorithm converges for all signals.

The convergence is only convergence in the mean: what can happen is the filter gets close to the optimal value, but bounces around a bit. It's not guaranteed to stay precisely on the optimal value.

The main constraint is on the value of the step size, $\mu$. It needs to be big enough to ensure convergence and small enough to not diverge.

The picture below shows a couple of examples. One where the convergence is slow but steady (low, but high enough $\mu$) and the other where $\mu$ is a little higher where the convergence bounces around a little.

  I understand LMS make use of Stochastic Gradient Descent and the update equation for the same is $\hat{h}(n+1)=\hat{h}(n)+\mu x(n)e^*(n)$, what would be the update equation if we would have used a normal Gradient Descent?

The main difference is that stochastic gradient descent uses sample-by-sample updates, rather than a larger data set. See this discussion on Quora for more details.


First you have to clearly define what you mean by "noise". Is simply the electronic noise in your measurement system or is it also brain activity that's not related to your stimulus.

One easy way to determine SNR is to compare the signal level during a stimulus with the signal level without stimulus. If noise and signal are uncorrelated you can determine the signal level through energy subtraction and thus determine the SNR.

Depending on your noise definition, this may not be practical cause you can't easily shut off the brain of the subject :-).


Why is a first/second derivative useful in spectroscopy?

Derivative spectroscopy was at its peak in the 70-80s. It is still very useful. The only problem with derivatives is that they enhance noise like crazy i.e., they suppress low frequency signals and enhance high frequency signals.

Before you do anything with the derivatives, you have to apply a smoothing process on them, specially on the second derivative. This "problem" of noise with derivatives was very well known then and smoothing was necessary before anyone used them.

For smoothing, you have apply a filter in such a way that it does not shift the peak positions, and the simplest one is the centered moving average. Choose an odd number of data points, say 5 (or higher), so the data at point x_3 will have two data points above and two data points below. Take the average of all 5, at position number x_3. Continue doing this for the rest of the points. You will start to see features in the second derivative.

Basically, derivatives are used to reveal the fine structure of a band. Shoulders become clear. The second derivative can "act" like a high resolution signature of a given molecule.

A very nice paper on higher order derivatives is
High‐Resolution, Higher‐Order UV/VIS Derivative Spectrophotometry

And there are some other tricks with second derivatives to resolve overlapping signals by my co-author.

https://terpconnect.umd.edu/~toh/spectrum/Differentiation.html


Gain error when performing STFT spectogram of increasing frequency signal

A sinusoidal signal is

$$y(t) = sin(\phi(t))$$

The phase of a linear chirp, is

$$\phi(t) = \phi_0 + 2\pi\left(f_0t+\frac{k}{2}t^2\right)$$
where
$$ k = \frac{f_1-f_0}{t_1-t_0}$$

So your chirp signal is missing a 0.5 correction for the frequency variation.

Use the scipy.signal.chirp:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp

def stft(x, L, overlap, window):
    opoints = int(L*overlap)
    X = np.array([np.abs(np.fft.fft(window*x[i:i+L])[:L//2])
                 for i in range(0,len(x)-L,opoints)])
    return X

tmax = 5
fs = 11025
N = fs*tmax

t = np.linspace(0, tmax, N)
x = chirp(t,0,tmax,2000)
windowsize = 1024

f = np.linspace(0, fs, windowsize)
window = np.hanning(windowsize)
stft = stft(x,windowsize,0.5,window)

plt.imshow(stft.T, cmap = plt.get_cmap('gist_yarg'), origin='lower', aspect='auto',
           interpolation='nearest', extent=(0,t[-1], 0,fs/2))
plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]")
plt.show()

Sources:

1: https://en.wikipedia.org/wiki/Chirp

2: Basic question about chirp signal


Low order sinc interpolation vs. polynomial interpolation for variable fractional delay

If you don't need matching points to come thru the resample interpolation unchanged, then you can use a different low-pass filter kernel.  A minimum phase equivalent to a linear phase windowed Sinc will have an identical frequency response, but less delay, especially when using a longer window.


Your instincts are right. You won't increase SNR by increasing sampling rate because you also have more noise samples. The overall signal energy and noise energy does not depend on sampling rate.

But, I believe you are mixing up two problems. In your case the frequency is not constant, therefore Fourier based analysis yield less than optimal results. Your problem is in the signal itself, not the signal to noise ratio. You can have a look at the high-order ambiguity function (try this or that) or the Wigner-Ville transform. Those are more suitable for linearly changing frequency.


MFCC's are related strictly to spectrum shape and what's more - these are relaying on the mel scale. Additionally each coefficient is a DCT value from fitted cosinusoids to log-energies. Thus interpretation of pitch would be troublesome for you.

Chord recognition is usually based on chromagram calculation or so called chroma features. Check out work and publications of guys from QMUL, i.e: Real-time Chord Recognition for live performance - in here you will find another references.

You might also find their software useful: Sonic Visualizer together with lot's of so called Vamp plugins. You can check out the code and even use command line tool to extract features you want.


doing sine regression to recover the paramters from a given signal

If your samples are equidistant you can use a model function like
$$\hat{x}(n)=a\sin(bn+c)+d$$ where $n$ is your sample index. You can compute the parameters $a,b,c,d$ by setting up an overdetermined nonlinear system of equations.

Let $x(n)$ be your signal. Then you'll have

$$a\sin(b\cdot 0 + c) + d = x(0)\\
a\sin(b\cdot 1 + c) + d = x(1)\\
a\sin(b\cdot 2 + c) + d = x(2)\\\vdots
$$

(Of course, the equalities are only approximate equalities.) You can use a nonlinear method such as the Newton method to solve this nonlinear overdetermined system. The only problem is that you need a sufficiently good initial solution. Try to find a simple initial guess by inspection (if your data looks at all like a sinusoid ...). For $a$ and $d$ you could simply use

$$a_0 = \max(|x(n)|),\quad d_0 = \frac{1}{N}\sum_{n=0}^{N-1}x(n)$$

as initial guess ($N$ is the number of data points).


How to do convolution in Fourier Series?

Following your notation, if we define the coefficients of the Fourier series of $z(t)=x(t)y(t)$ as $c_k$:

$$c_k=\sum_{l=-\infty}^{\infty}a_lb_{k-l}$$

Notice that the $k$ in which you evaluate the left side of the equality has to be the same as the $k$ in the summation. Therefore, for example, if you want to find $c_0$, you should calculate:

$$c_0=\sum_{l=-\infty}^{\infty}a_lb_{-l}$$

Then you have to calculate that summation. Given that $a_l$ and $b_{-l}$ are non-zero only if $l$ equals $1$ or $-1$:

$$\begin{align}
c_0
&=\sum_{l=-\infty}^{\infty}a_lb_{-l} \\
&=a_{-1}b_1 + a_1b_{-1} \\
&=\frac12 \frac{1}{2j} - \frac12 \frac{1}{2j}=0
\end{align}$$

I think this might have clarified the notation for you. You can go on and calculate the rest of the coefficients $c_k$ for all values of $k$.


How to remove 50 Hz frequency noise from a PCM file?

You should just use a notch filter.

Below shows two sinusoids added together (top) and the result of filtering the lower frequency one (bottom).

The filter is just:
$$
H(z) = \frac{1 - 2 \cos(2\pi \omega_{b}) z^{-1} + z^{-2}}
{1 - 2 \alpha \cos(2\pi \omega_{b}) z^{-1} + \alpha^2z^{-2}}
$$
where $\omega_b$ is the bad frequency and $\alpha<1.0$ is related to the notch width / depth.

R Code Below

#29936

fs <- 8000
fbad <- 50
phibad <- 2*pi*runif(1,0,1)
fgood <- 101.98340234
phigood <- 2*pi*runif(1,0,1)

T <- 1000
t <- seq(1,T)

x_noisy <-  sin(2*pi*fgood/fs*t + phigood) + sin(2*pi*fbad/fs*t + phibad)

num <- c(1, -2*cos(2*pi*fbad/fs), 1)
alpha <- 0.99
den <- c(1, -2*alpha*cos(2*pi*fbad/fs), alpha*alpha)

x_filtered <- filter(num, den, x_noisy)

par(mfrow = c(2,1))
plot(t,x_noisy, type="l")

plot(t,x_filtered, type="l")


What is the phase of the Fourier Transform of $h(n)={\frac13,\frac13,\frac13}$ (zero position is at 2nd element)

The result on the phase comes from writing $H(\omega)$ as a modulus/phase product: $H(\omega) = |H(\omega)|e^{-\imath\phi(\omega)}$. The value of $\phi(\omega)$ is consistent with the sign of $(1+2\cos \omega)$.

There is a little inconsistency is your notation (maybe from the book!): the phase seems to takes both values $0$ and $\pi$ at $w=\frac{2\pi}{3}$. This is not vital as the modulus vanishes at this point.

What you might have forgotten in your reasoning is the fact that $t\to \tan t$ is $\pi$-periodic, and that $\tan \pi  =0$ too.


LDPC codes are a subset of linear channel codes. In terms of functionality, they are similar to other linear channel codes. However, they are constructed using a sparse bipartite graph.

Fountain codes are rateless codes. These type of codes are more useful in broadcast/multicast settings. If we use fixed-rate codes in such settings, we need to make sure that each and every packet is correctly received at all receivers. This can become inefficient since each broadcast can potentially be re-iterated  (depending on the re-transmission protocol) until all receivers receive all packets. Consider for example a case with 10 receivers, all have correctly received the packet #2 except one. This specific packet will be broadcast to all until the last receiver gets the packet. However, with a rateless code, the packet index is not a matter anymore (i.e. whether it is packet #2 or not). This is because a group of packets are selected at first and new encoded packets are generated randomly (hence, all randomly encoded packets look just like a new piece of information about the group and the index does not matter anymore). Any linearly independent subset of valid packets is enough to decode (acquire) the original group of packets.

Also I can see you have difficulty understanding acknowledge policy. We can have many different types of such policy. So based on the Ack protocol, the transmitter stops and listens for the ack from the receivers to find out their status in terms of the number or the index of missing/collected packets.


Symbol synchronization with M-QAM. M < 16

Are there any superior ways of doing symbol sync?

Yes! As you've noticed, your initial estimation of timing isn't "valid" forever.
This might be caused by clock jitter, or simply remaining sampling rate offset.

It's a very common problem in communication receivers.

What you can of course do is that you simply do the same estimate you did at the beginning periodically (with the already recovered data instead of a preamble, even!), or even simply for every symbol.

A popular way of doing so is basically using the derivative of the pulse shape to track the "flattest" point in the pulse-shaped receive signal even before applying the (typically matched) filter in the receive chain. You might want to watch this video (slides) of GNU Radio's synchronization expert Andy Walls – I like that talk especially since not only is it clear enough to be used as instructional material, it actually also addresses "how it is actually done".


Does Reed-Muller codes exist when $r = 0$ and $m = 0$?

The codewords in a Reed-Muller code of degree $d$ and length $2^n$ are binary vectors of the form
$$\big(f(0,0,\ldots, 0, 0), f(0,0, \ldots, 0, 1),
f(0,0, \ldots, 1, 0), f(0,0, \ldots, 1, 1), \cdots, f(1,1,\ldots, 1,1)
\big)$$
where $f$ is the corresponding binary polynomial of degree at most
$d$ in $n$ variables.

Suppose first that $n > 0$. The dimension of the code is
$$k = \sum_{i=0}^d \binom{n}{i}\tag{1}$$
If $d=0$, then $k = \binom{n}{0}=1.$ There are only two "polynomials"
of degree $0$ in $n$ variables, namely the constants $0$ and $1$, and consequently
the codewords of length $2^n$ are $000\cdots 0$ and $111\cdots 1$.

If $n=0$ also, the codewords (if any) are of length $2^0 = 1$.
Now, you might want to claim from $(1)$ that $k = \binom{0}{0} = \frac{0!}{0!(0-0)!} = 1$ and so there are $2^k = 2$ codewords
(which, of course, are $0$ and $1$).
Alternatively, the $k$ in Equation $(1)$ is the total number of
subsets of cardinality $d$ or less of a set of $n$ elements.
There is exactly one subset (of cardinality $0$)
of the set of $0$ elements (namely,the empty set $\emptyset$).
and that subset is $\emptyset$ itself. Thus, the Reed-Muller code
of order $0$ and length $2^0$ is the $(1,1)$ code consisting
of two codewords $0$ and $1$. As a linear code, this is the
identity map.


Conversion between DTFT in radians/sample to DTFT in cycles/sample

You can use the following :

Consider the discrete time sequence $x[n]$ and the DTFT $X(\omega)$ related by:
$$ X(\omega) = \sum_k x[n] e^{-j \omega n} \longleftrightarrow x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(\omega) e^{j\omega n} d\omega$$

Then define $X'(\mu)$ as the normalized DTFT, which is related to $X(\omega)$ as
$$X'(\mu) = X(2\pi \mu) = \sum_k x[n] e^{-j 2\pi \mu n} $$ and then
$$x[n] = \int_{-0.5}^{0.5} X'(\mu) e^{j 2\pi \mu n} d\mu $$

Applying this to your example case:
$$x[n] = 1 \longleftrightarrow X(\omega) = 2\pi \sum_r \delta(\omega - 2\pi r)$$ then since $X'(\mu) = X(2\pi \mu)$ we have

$$X'(\mu) = 2\pi \sum_r \delta(2\pi \mu - 2\pi r) = 2\pi \sum_r \delta(2\pi [\mu - r]) = 2\pi \sum_r \frac{1}{2\pi} \delta(\mu - r) $$
$$X'(\mu) = \sum_r \delta(\mu - r) $$

as expected, where we've used the property that $\delta(ax) = \frac{1}{|a|} \delta(x)$.


One straightforward approach is to do blur detection (see here and here) and prune the frames that are detected.


Calculating cross-correlation using Walsh-Hadamard transform

The link provided as a reference is broken. But, regardless of what the OP's source says, the Walsh-Hadamard transform method (multiply the data vector with a $2^n\times 2^n$ Hadamard matrix $H_n$ in Sylvester form $$ H_n = \left[\begin{matrix} H_{n-1} & H_{n-1}\\ H_{n-1} & -H_{n-1}\end{matrix}\right]$$ does not give the cross-correlation of the data vector with a binary maximum-length linear feedback shift register sequence (presumably the MLS referred to by the OP) of period $2^n-1$. It is possible to permute the data vector  entries to get a close approximation of the result, indeed even the exact result, of the desired cross-correlation, but just applying the Walsh-Hadamard transform to the unpermuted data vector results in garbage, as the OP has discovered.
A description of how the MLS method works can be found in this answer of mine. The necessary permutations to make the Walsh-Hadamard transform method work are probably there somewhere on this forum.


Your notation is a bit odd, as you are mixing the time-domain with the $z$-domain in the equation you wrote. If you refer to the system

$$y[n]=x[n-1] \quad \mathrm{or} \quad Y(z)=z^{-1}X(z)$$

then you are correct. To put it mathematically, the transfer function of the filter would be:

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1}{z}$$

The frequency response is then

$$H(z) \Bigg|_{z=e^{j\omega}} =\frac{1}{e^{j\omega}}$$

You can see that the magnitude of $H(e^{j\omega})$ is $1$ for all frequencies. Therefore, it is an all-pass system.


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

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

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


Minimum length of windows and separable signals

I don't know if you can have a formula in terms of window length, number of DFT points and difference between two angular frequencies you are trying to resolve. But you can do some experiments with different types of the windows and determine this in MATLAB. Note you need to consider the following.

For a DFT of size N you get two frequency bins separated by $\frac{2*\pi}{N}$. So you should have N large enough resolve $\omega_2 - \omega_1$
Note that before you take DFT, when you window it the frequency response will be convolution of two frequency response - frequency response of the window and the frequency response of the signal which in your case pretty much a pair of delta function.

Now Hanning window will have a wide main lobe and suppressed side lobes. So if the second signal's frequency falls within the main lobe then you won't be able to detect it unless amplitude of the second signal is large enough compared to the smeared component of the first signal at the location of the second signal. Or if the second signal falls in some of the side lobe location and second signal's amplitude is large enough compared to the first signal's component in side lobe region you should be able to detect it.

If you chose a rectangular window the main lobe will be substantially shorter compared to rectangular window. But the side lobes will be substantially larger compared to the Hanning window case.

There exist the trade-off between using different window sizes.


Can IR emitter signal be distorted by curved glass housing around receiver?

Short answer: No.

Long answer: you can of course create shadowing that way, and that would disturb operation. And of course, a glass wall will refract infrared just as it refracts any other light. But glass is usually IR-transparent, just as it is transparent visible light, and an IR remote isn't exactly a source of a focused beam, so chances are that, no, there's no difference in operation.

What will happen is that the IR will be a tiny bit weaker, due to absorption in the glass. But that won't be really much, I guess.

  Could the jam jar's reflection and refraction be impacting the length of the pulses?

How so? It's the same signal entering and exiting the glass. Admittedly, at a different position, maybe a bit attenuated, but still, if you wanted to change the duration of the pulses, then you'd need a frequency converter (which is a nonlinear device) or a time machine. You have neither.

What could be happening is that due to reflections in the glass, at the receiver, you'd get one received signal that took "the shortest" path, overlaying with a "delayed" path. For that to have any negative effect, the path length difference would need to be in the order of the duration of your pulses.

Infrared remotes work with longer pulses modulated to something like 38 kHz, if I remember correctly. Ignoring that wavelengths are different in glass (factor of maybe 1.5, and we're just here to get a feeling for the orders of magnitudes involved): a signal travelling a speed of light travels, in one period of a 38 kHz signal,
$$\frac{3\cdot 10^8\frac{\text m}{\text s}}{3.8\cdot 10^4\,\text{Hz}} \approx 0.7 \cdot 10^4\,\text m = 7\,\text{km}$$

Since any distance your "detouring" beam might take will be way, way smaller than that, we can rule that out.

I think for the typical roughness of a marmelade jar glass, we can also rule out thin-film effects on the surface.


Phase plot of $H(e^{j\omega})=(1-re^{-j\omega})\left(1-\dfrac{1}{r}e^{-j\omega}\right)$, $0<r<1$

In fact you've already proved that the two expressions must be identical. But if you want to show it the hard way, you can indeed use the formula you suggested. The argument of the combined arctangent function becomes

$$\begin{align}\frac{\frac{r\sin(\omega)}{1-r\cos(\omega)}+\frac{\sin(\omega)}{r-\cos(\omega)}}{1-\frac{rsin^2(\omega)}{(1-r\cos(\omega))(r-\cos(\omega))}}&=\frac{r\sin(\omega)(r-\cos(\omega))+\sin(\omega)(1-r\cos(\omega))}{(1-r\cos(\omega))(r-\cos(\omega))-r\sin^2(\omega)}\\&=\frac{\sin(\omega)(1-2r\cos(\omega)+r^2)}{-\cos(\omega)(1-r\cos(\omega)+r^2)+r(\underbrace{1-\sin^2(\omega)}_{\cos^2(\omega)})}\\&=\frac{\sin(\omega)(1-2r\cos(\omega)+r^2)}{-\cos(\omega)(1-2r\cos(\omega)+r^2)}=-\frac{\sin(\omega)}{\cos(\omega)}=-\tan(\omega)\end{align}$$

Now you just have to be careful with the principal value and you'll obtain the same result as with your first method (which is of course much smarter).


Adaptive equalization vs inverse of transfer function

Inverting a channel can only be done when the channel is a minimum phase system (trailing echos only). A minimum phase system is characterized as having all zeros in the left half plane (for the s plane, or equivalently in a sampled system and the z plane all zeros inside the unit circle). Inverting such a channel results in poles where every zero exists, and a causal system that has any poles in the right half plane (outside the unit circle) is not stable.  So a minimum phase system has a stable causal inverse, while a mixed phase or maximum phase system does not.


Condition for modelling noise as white noise

I will try to answer your three questions as a whole, with the following insights:

"Whiteness" means: flat power spectral density, thus zero autocorrelation for lag different than zero (uncorrelated samples), thus wide sense stationary (time-independent autocorrelation and mean).

White Gaussian noise passing through a linear filter always produces Gaussian noise, only changing its parameters (mean and variance), since any linear combination of uncorrelated Gaussian random variables produces a Gaussian random variable. However, the whiteness is lost in general, and thus samples become correlated after passing a linear filter.

If your system's transfer function has a bandwidth in $\omega$ significantly lower than the noise's bandwidth, you may approximate the input noise as white, since what the system "sees" is a flat spectrum.

Anyway, you must take into account that white noise is always an approximation for real-life systems. There is no such thing as a noise with infinite bandwidth in practice. But you may be in presence of noise with very large bandwidth (compared to your system), and thus successive samples of such a noise are (nearly) uncorrelated.


Methods for calculating signal quality metric

From your question it sounds like you are looking for a SNDR metric.  In that case you are measuring a distortion level by comparing the distorted signal to a reference.
For a discrete-time signal $y$, comparable to reference $r$, the SNDR would be:
$$ SNDR = \frac{E[r^2]}{E[(y-r)^2]} $$


We do not define it in that manner. Note the difference between the time and Laplace domains.

Generally, we characterize a system by the response produced when the system is excited by an impulse response $\delta(t)$. Given the output $y(t)$ for such signal we ask which function $h(t)$ did we convolve with the $\delta(t)$ to produce it. Formulation:
$$\delta(t)*h(t)=y(t)$$

Here $*$ represent convolution and not a product. When we find $h(t)$ we can produce the output for any other input $x(t)$ by means of convolution:
$$x(t)*h(t)=y(t)$$

If we transform into the Laplace domain, a convolution is translated into multiplication:

$$\mathcal{F}\{x(t)*h(t)\}=\mathcal{F}\{y(t)\}$$
$$\mathcal{F}\{x(t)\}\cdot \mathcal{F}\{h(t)\}=\mathcal{F}\{y(t)\}$$
$$X(s)\cdot H(s)=Y(s)$$
$$H(s)=\frac{Y(s)}{X(s)}$$

It is a convention that the capital letters represent the transformed signals.  Also, $t$ has changed into $s$ because we are now in the Laplace domain. Since the convolution transformed into a product we may now divide both sides by $X(s)$. We could not do this on the time domain because there we had a convolution operator and not a product.


The effects of frequency resolution in FFT on the frequency spectrum

What may explain why your spectrum is noisy is that you are computing it using a single burst of data. You will have to smooth it by averaging successive spectrums applied on your measurements. The number of samples that shall be used for the averaging have to be tuned depending on the speed of variation of your phenomena. After having done this, the other questions can be addressed more easily.
Best regards,
Mourad


The problem is that your input vector is not symmetric when considering the definition of the DFT:

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

Note that from (1) the time indices start at $n=0$. If you periodically continue the signal $x[n]$ in your example you get

n:    ... -5 -4 -3 -2 -1 | 0 |  1  2  3  4  5  6  7  8  9  10

x[n]: ...  1  0  0  0  0 | 0 |  0  0  0  1  1  1  0  0  0  0

If $x[n]$ were symmetric, the values for negative $n$ would need to be the same as the values for positive $n$. For $n=\pm 4$ this is not the case. If you defined your signal instead as

x[n] = 0 0 0 0 1 1 1 0 0 0

i.e. a signal of length $N=10$, then it would be symmetric and, consequently, its DFT (FFT) would be real-valued, apart from numerical errors.


How does a Huffman encoded message gets unified before channel encoding?

They don't get unified.

Think of the transmitter pipeline (data source, source encoder, channel coder, modulator, etc) as a sequence of independent blocks. Blocks don't assign any particular meaning or order to their input: they regard the input as just a stream of bits.

So, the output of the Huffman encoder can be regarded as a stream of 0s and 1s. The channel coder will subdivide the stream in groups of $k$ bits and convert them to groups of $n$ bits. Again, the channel coder output is a stream of bits that is processed by the modulator.

This approach has several benefits. One is that each block can be studied independently of the others. Also, you can connect, disconnect and replace blocks, and the system will continue to work, because there is no dependency between the blocks.

The main disadvantage is that there is a little bit of extra performance that can be obtained by making the blocks dependant of each other -- see for instance "joint source-channel encoding".

The receiver is more complicated because a block may need to know the boundaries of the code blocks. This is solved by introducing a frame structure to the data.


This is because intersymbol interference is caused by dispersion in the channel (assuming properly designed pulse shaping filters such that there is no "native" ISI in the waveform itself). This dispersion can be caused by analog filtering or over the air multi-path as illustrated in my graphic below:

For the case of multi-path, the multiple transmission paths result in a "delay spread".  According the the reference in the next graphic, for the macrocellular systems the typical delay spreads are:

Open areas: 0.2 uS
Suburban: 0.5 uS
Urban: 3 uS

When the symbol duration is significantly longer than the delay spread, there is minimum interference from one symbol to the next, but if the symbol duration approaches the delay spread, the interference will be significant! For a multipath case resulting from multiple replicas of the same waveform arriving at different times, this can be corrected with equalization. For the case of OFDM, we are actually converting the waveform into independent channels of significantly longer duration, thus solving this issue. The cyclic prefix allows a guard interval to eliminate the intersymbol interference from a previous OFDM symbol. To achieve this, the cyclic prefix must be at least as long as the delay spread of the channel.


Nyquist noiseless channel capacity; how can bit-rate be two times the bandwidth?

I think you're confusing two different (but related) terms.

Nyquist says that in a channel of bandwidth $B$ you can transmit up to $2B$ orthogonal pulses per second. So, $R_p \leq 2B$, where $R_p$ is the pulse rate.

To achieve $R_p = 2B$, the pulses need to be sinc-shaped. Other, more practical pulses achieve slightly less than that. For example, raised cosine pulses with 50% excess bandwidth achieve $R_p = 1.5B$.

If you send pulses with $L=2^k$ amplitude levels, then the bit rate is $R_b = kR_p$.

So far for Nyquist, who was only interested in maintaining orthogonality. Shannon wanted to find a limit for the bit rate of reliable communication, with probability of error $P_b \rightarrow 0$. This limit is the capacity $C = B \log_2(1+\text{SNR})$, assuming an AWGN channel and orthogonal pulses. Note that $L$ is not mentioned explicitly.

  what is an example of a signal with capacity twice the bandwidth?

The capacity is a property of the channel, not the signal. A signal has a bandwidth, and implicitly a pulse rate and a bit rate. Regardless of that, the channel imposes a limit on the bit rate (its capacity), if you want reliable communication.


How to change the start point of resampling process?

since @Community resurrected this question, i'll pick up a little where hotpaw left off.  i'm not gonna run this through MATLAB or Octave (i can never remember which are columns and which are rows).

let $r$ be the resampling ratio.  $r>1$ means that your sample rate has increased and you're making the sequence longer because there are more samples.  so if you pad a delay of $\tau$ at the beginning of the input $x[n]$, that should result in a delay of $r \cdot \tau$ at the output $y[n]$ (whether $\tau$ is an integer number of samples or not).

now suppose you pad (using concatenation) the input with $P$ samples of zero before $x[0]$.  $P$ is a positive integer. that will result in a pad in the resampled output $y[n]$ of $rP$ samples where $rP$ is not necessarily an integer.  but

$$ \lfloor rP \rfloor \triangleq \operatorname{floor}(rP) $$

is an integer, and you can strip those samples (which will be mostly zero) off of the beginning of $y[n]$.  $\lfloor u \rfloor$ is the most positive integer that is no greater than $u$.  this is always the case:

$$ \lfloor u \rfloor \le u < \lfloor u \rfloor + 1 $$

at least for $u \ge 0$, we might say that $\lfloor u \rfloor$ is the "integer part" of $u$.  the "fractional part" of $u$ is whatever is left:

$$ 0 \le u - \lfloor u \rfloor < 1 $$

so, if what we want left is a fractional delay of $\alpha$ samples where $0 \le \alpha < 1$, then we want to choose an integer $P$ to pad the input such that $\alpha$ is the fractional part remaining in the output delay after stripping away the integer part, $\lfloor rP \rfloor$.  that is, given the sample rate ratio $r$, and fractional delay $\alpha$, we want to choose $P$ such that:

$$ rP = \lfloor rP \rfloor + \alpha $$

or

$$ \alpha = rP - \lfloor rP \rfloor \ . $$

after thinking about this a while, the only algorithm i can suggest is to start with $P=1$ and iteratively increment $P$ so that the above equation is close enough to be accurate.  that is, given an error constraint on $\alpha$, (call it $0 < \epsilon \ll 1$), increment $P$ until

$$ \left|\alpha - \left( rP - \lfloor rP \rfloor \right) \right| < \epsilon $$

pad the input with $P$ samples of zero, and strip the output of $\lfloor rP \rfloor$ samples, and you will be left with $\left( rP - \lfloor rP \rfloor \right)$ samples of delay.


The various weighting systems are a domain-speific attempt to adapt the raw edit distance metric into something more "perceptual". The needs & error tolerances for one domain, say large vocab open speech, may be different from the command-and-control, "assistant" systems.

Humans differentially weigh different types of ASR errors and are less forgiving with errors (like say insertions) which may distort the meaning of the transcription.

Have a look at this paper: Predicting Human Perceived Accuracy of ASR Systems


How can I measure the harmonicity of an audio spectrum?

Please take a look at this previous answer.  Using the semantics and mathematical definitions from that answer...

Let $x[n]$ be the input signal from which you want to measure harmonicity or periodicity.  Assume $x[n]$ has no DC component (like it is the output of a DC blocking filter).  And $n_0$ is the sample index of the neighborhood of signal around where you want the harmonicity to be measured.  $N$ is a large integer defining how wide your rectangular window is.

Average Squared Difference Function, ASDF:

$$ Q_x[k, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 $$

$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $.

From that define an alternative form of Autocorrelation:

$$ R_x[k,n_0] = R_x[0,n_0] - \frac12 Q_x[k, n_0]  $$

where

$$ R_x[0, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} (x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor])^2 $$

Since $Q_x[0, n_0] = 0$ and $Q_x[k, n_0] \ge 0$ for all lags $k$, that means that $ R_x[k, n_0] \le R_x[0, n_0] $ for all lags $k$.

You will find a period $P$ such that $P > 0$ and $R_x[P, n_0] > R_x[k, n_0]$ for all lags $k>k_0$, other than the highest peak at $k=0$.

$k_0$ is the smallest value of $k>0$ such that $ R_x[k, n_0]<0 $.

Once that peak location $P$ is found, the harmonicity is

$$ r_x(n_0) \triangleq \frac{R_x[P, n_0]}{R_x[0, n_0]} $$

If $r_x(n_0) = 1$, then $x[n]$ is perfectly periodic (or "harmonic") in the vicinity of sample index $n_0$. If $r_x(n_0) = 0$, it's noise, no periodicity to be found.


How should calibration points be distributed?

For that I have some guidelines, which I'd like to share. Even though some of these might not be valid for certain calibration software, for typical cases, they hold:

1) The entire pattern should be visible in all sequences.
2) At least 1/4th of the image should be occupied by the pattern.
3) Preferably, the lighting should be uniform everywhere on the pattern.
(As little shadows as possible)
4) The light should reflect from the pattern to the camera (as much as possible).
Thus, objects made up of metal are not suggested.
5) There should be minimal blur. If focus/depth-of-field is an issue,
use a smaller plate or less tilts.
6) If you can't fix the camera, it shouldn't move at all.
For extrinsic calibration, moving the camera after calibrating will ruin the entire
calibration. Be careful about this as even small shakes or movements distort the accuracy.
7) Don't pre-process the calibration images.
Even if you do very smart things, you might perturb the accuracy.
8) Ideally, have 15-50 images fitting the above criterion.
This should give you 450-1500 points to optimize in total.
9) If you are printing the pattern on a piece of paper, then make sure that
contrast between black-white squares is significant and no strange shape is appears on the paper.
10) Don't bend the calibration pattern. Make sure that it is as rigid (and planar) as possible
11) If you are printing on an A4 paper, make sure you have some white space around the pattern,
remember you will want to hold it somehow.
12) If you have distortion, use small circles, they are less affected from distortion.
If you have low resolution but less distortion use larger circles.
13) Start with low tilts and fronto-parallel images.
Move to the tilts and rotations in all axes, gradually, covering the entire visible space.
If you have bias in the space coverage or distribution of points,
then it means that some regions of the space will be better calibrated.
Keep this as equal as possible. Some variation can be tolerated.

I hope these help.


The FFT calculates the spectrum at N points in the range from $[0, 2\pi)$ where due to the cyclic property the relative frequency 0 must be equal to the frequency $2\pi$. This is best shown via an example. Note that $2\pi$ corresponds to your sampling frequency.

Example:
For a length 3 impulse response, say $h=[0.5, 1, 0.5]$, the FFT calculates the spectrum at the relative frequencies $[0, \frac{1}{3}2\pi, \frac{2}{3}2\pi]$, instead of the frequencies $[0, \frac{1}{2}2\pi, 2\pi]$. One way to calculate these frequencies in Matlab for an $N$ points FFT is as follows:

N = 3;
h = [0.5 1 0.5];
fq = linspace(0, 2*pi, N+1);
fq = fq(1:end-1);
H = fft(h);
figure;
plot(fq, abs(H));

The function fftshift takes the right part of the spectrum and places it on the left such that you obtain the frequencies $[\frac{2}{3}2\pi, 0, \frac{1}{3}2\pi]$.

When you are not using the fftshift function and just plot the range $[0, 2\pi]$, then you can use positive time indices for t.


Fourier transform of moving image sensor

I feel that the question needs a little bit more clarifying but there probably is enough material here to provide some sort of response and adjust it later if needed.

To cut a long story short: Rolling Shutter distortion. This concept is key here (irrespectively of whether or not you have a rolling or global shutter camera that is unsynchronised to the interferometer).

Now, I will treat $v, R$ as independent from time here for a minute. Also, $F_{FR}$ will stand for frame rate.

Your "problem" seems to be that due to $v$, each pixel row might actually be sampling spectral lines from two radically different wavelengths. If this is true then each one of your frames will acquire a slant whose angle will be proportional to $\frac{v}{F_{FR}}$.

In the time dependent version, the frames are capturing a "curve" whose curvature is proportional to $\frac{v(t)}{F_{FR}}$ and...good luck with that.

Since you are dealing with interferometry, you are only really imaging ONE single focal plane, therefore, you could apply a simple (anti)-skew transform to counteract the skew caused by the rolling shutter. This, of course, is within limits. In any case, the rolling shutter effect of the camera inserts a phase error in the image which, to an extent, can be counteracted but in the end it all comes down to controlling $v$.

Assuming that your camera can maintain a constant $F_{FR}$ for the duration of the scan, then essentially, you don't have a 2D camera.

I would like us to imagine it as one long linear $M \times N$ camera (because of $v$ and interferometry) only it is wrapped around with a "wrap" factor that depends on $v$.

So, theoretically, you could find a $v$ which scans a spectrum range in one frame. Maybe $v$ is such that the spectral range is contained in two frames or three frames or four...Obviously, the more, the better. Therefore, each pixel row now has a $\lambda_{start}$ and a $\lambda_{end}$. So, I want you to imagine a chirp that is wrapped across the camera sensor or in fact across many different frames.

Therefore, you could convert your single $M \times N$ image, into one long vector and then do a spectrogram on this vector. This vector is essentially encoding wavelength to space (across the camera sensor).

Does this have any distortions? Of course it has, but it would be like trying to obtain the frequency response of a system by feeding it a chirp. Within limits it is perfectly possible. If you scan it too fast, you might skip a spectral line

Perhaps it is easier to see why $v(t)$ would insert a phase error to the "chirp" when considering the camera as one long image sensor (?).

Hope this helps.


How does MATLAB recover picture from magnitude spectrum alone?

Your code uses the phase for the reconstruction. Have a look at the output of fft2(x); they are complex numbers, i.e. the contain phase and magnitude. Have a look at this code:

%%
[x,map] = imread('http://www.cs.cmu.edu/~chuck/lennapg/lena_std.tif'); %import image
%%
i = rgb2gray(x);   % to greyscale

I = fft2(i); % 2D FFT

I = fftshift(I); % centre
mag = abs(I);
phase = angle(I);
I_magonly = abs(I);
I_phaseonly = I ./ abs(I);

figure(1);
subplot(321); imshow(abs(i),[]); title('Original');
subplot(322); imshow(log(double(abs(I))),[]); title('spectrum magnitude');
subplot(323); imshow(phase,[]); title('spectrum phase');
subplot(324); imshow(abs(ifft2(I)),[]); title('reconstructed');
subplot(325); imshow(abs(ifft2(I_magonly)),[]); title('reconstructed magnitude only');
subplot(326); imshow(abs(ifft2(I_phaseonly)),[]); title('reconstructed phase only');

It loads the Lena image, creates the FFT, and displays phase and magnitude information. Then it tries to reconstruct only based on phase or magnitude. Apparently, phase contains more information (for our eyes) than the magnitude. But you need to both exactly reproduce the image from its Fourier transform.


Power-Spectral Density is the distribution of power along the frequency axis.
It is generally used for non-finite energy signals (mostly not limited in time signals), who aren't square-summable.
The signal's PSD is the autocorrelation of the signal's Fourier Transform, as stated by the Wiener–Khinchin theorem.
In Matlab:

N = length(S);
F = fft(S);
F = F(1:N/2+1);
PSD = (1/(2*pi*N)) * abs(F).^2;
PSD(2:end-1) = 2*PSD(2:end-1);
freq = 0:(2*pi)/N:pi;

see: https://de.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

Cross-Spectral Density is the same, but using cross-correlation, so you can find the power shared by a given frequency for the two signals using its squared module, and the phase shift between the two signals at that frequency using its argument.

Cross-Spectral Density can be used to identify the frequency response of a noisy LTI system : if the noise is not correlated to the input or output of the system, its frequency response can be found from the CSD of the input and output.


How would I calculate the state-variables of digital biquad filter (direct form II transposed) at steady-state?

Consider a discrete-time LTI system (that's your second order stage) described by a LCCDE in the form of a recursion equation:

$$y[n] = a y[n-1] + b y[n-2] + c x[n] + d x[n-1] $$

where $a,b,c,d$ are real or complex constants (coefficients of your IIR filter) and the values $y[n-1],y[n-2]$ are previous values of the output, where as $y[n]$ is the current value of it. And $x[n],x[n-1]$ are the current and previous values of the input respectively.

In your case, you assume a constant input $x[n] = C$ for all $n \ge 0$, hence you have $x[n] = C$ and $x[n-1] = C$ as $n \to \infty$. Now assuming that the filter has reached its steady-state output implies that $$y[n] = y[n-1] = y[n-2] = K$$ as $n \to \infty$

Hence putting those values into the recursion equation yeilds a relation between output $K$ and input $C$ as:

$$K = a*K + b*K + c*C + d*C$$ which can be arranged to yield:

$$K(1-a-b) = C(c+d)$$

At this point you tell us that output is $K=20$ mm, then you can find the corresponding input as $x = C = 20*(1-a-b)/(c+d)$, and your states will be:

$$ y[n-1] = y[n-2]= K , x[n] = x[n-1] =x[n-2] = C$$


2D Haar Transform Computational Complexity

This is just to address the lack of multiplications:  why would there be?

The coefficients of the Haar wavelet (in 1D or 2D) are always +1 or -1.  There is no point in doing a multiplication to find the value or its negative.

Check out these notes, particularly this page:


How to tell if a high-res flac file has been upsampled from a CD-quality file?

Let's get some basics out of the way first.
An audio format is basically defined by three different things: sample rate, bit depth and compression type. The compression type can be either lossless (WAV, PCM, FLAC, ALAC) or "lossy" (AAC, MP3, OGG, MQA, etc). Lossy compression works by shaping the quantization noise so it's masked by the remaining content and inaudible. Lossy formats are characterized by their bit-rate: the higher the bit-rate, the less quantization is done.
The typical life cycle is the following: recording, mixing and happens at 96kHz floating point. When that's done, there are typically two masters made: a studio master at 96kHz/24-bit and a CD master at 44.1kHz/16-bit.
The label will send one of these two to the distribution channel, e.g. Spotify. Which one depends on the contract. Spotify will then create the formats that they need for their streaming platform, which almost always includes lossy compression. For example Spotify uses OGG at 320 kb/s. Apple Music uses AAC at 256 kb/s. Which format gets to the customer depends on the contract between the customer and the distribution.

Is there anyway I can tell apart a "fake" high-res file (converted from CD-quality or mp3/aac) apart from a "real" high-res file

Here are a few things to look for in a spectral analysis

Each format has it's own rolloff: CD drops like a rock at 20 kHz and MP3 drops steeply at 16 kHz. If you have a 96 kHz file with these sharp drop offs, it's likely been up-sampled.
Inspect the content above 20 kHz. If there is random "noise like" features in there, it's probably genuine. If it has very little content and/or the content looks like a low-pass filtered mirror image of the content below 20 kHz, it's been up-sampled.
You can look at correlation at high frequencies. For a genuine recording this will mostly be uncorrelated. If there is significant correlation, it's a potential sign of "joint stereo coding" which could hint at lossy compression.
Look at the recording date: If it's been made before 1990 it's almost guaranteed to be up-sampled. There never was a digital studio master and the best they can do is to sample a tape master.

Some formats like MQA, include authentication, i.e. they put a "seal" on the content and can detect whether the "seal" has been tampered with. It's also possible to put a watermark on the content.
Update
You may be able to detect bit padding as well: Any "decent" CD quality file have been properly noise shaped and dithered. It's possible that signs of this will survive the up sampling process can be detected. I've neve tried that, so I don't know how feasible it is, but it sure would be interesting to try.


It works ... bitwise. Assume the value of a pixel is $200_{10} = 11001000_2$; then the bitwise not of that simply is $00110111_b=55$. Do that for all pixels.

It's not a useful transform in the context of intensities, if you ask me.

If the mask was actuall monochromatic (i.e. exactly black or exactly white in each pixel), it would make more sense.

So, don't "come across code" like that :)


Matching guitar signal with tablature, near real-time algorithms?

Assuming that the tablature is well defined; let's say it's a MIDI file complete with pitch bend commands.  you can make a synthesizer where each note is actually a tracking comb filter tuned to the midi note (and offset by the pitch bend command) with notches at each integer harmonic (including the fundamental).  if those tracking comb filters processes a sound file playing that MIDI file and the two are properly aligned in time, the tracking comb filters should null out most of the input audio, even polyphonic.

Then magnitude-square-LPF the result and a minimum value means a pretty good match.

Just an idea.  No equations, but I have posted about pitch detection before. And I thought I did about tuned comb filters, too (you need a precision delay line to tune the comb filter to whatever fine-tuned pitch that is specified).


How to implement an RF Mixer for chirp signal in python?

Wouldn't mixing them give you a sum and difference frequency and you have to lowpass filter to get the difference frequency only?  Also your chirps are the same frequency, just out of phase?  So the difference would be zero and the sum would be double.
So:

Actually create a chirp that is a different frequency from the original, not just out of phase.
After mixing, lowpass filter to keep only the difference frequency


I think this is simply a property of your signal and you indeed have two significant frequency components. I don't really know the nature of your signal, but for me, it makes perfect sense, you can distinguish two periods in your signal, as shown on the plot below:

As you can see there is period $\color{green}{T_1}$ and $\color{brown}{T_2}$ between consecutive "peaks". So by doing very rough approximations from your plot:

Corresponding frequency for $\color{green}{T_1}$ can be estimated as: $\dfrac{1}{7.1-6.3}=\dfrac{1}{0.8}=1.25 \mathtt{Hz} $
In case of $\color{brown}{T_2}$ you get: $\dfrac{1}{6.7-6.3}=\dfrac{1}{0.4}=2.5 \mathtt{Hz} $

Some quick and dirty plot of such case with two sinusoids of a given frequencies and amplitudes:

So indeed you have two frequencies present in your signal. What's more, component around $2.5 \mathtt{Hz} $ can be stronger. Please consider that your signal also has very prominent DC component (also visible at your spectrum).


Upsampling by 100 can be factored into upsampling by 2X, 2X, 5X and 5X for each output block (output blocks sized to be data cache resident), thus requiring a smaller phase table for each FIR interpolation kernel.

Some output points may either need to be overlapped across blocks or re-calculated, depending on computation speed vs cache load latency.


Understanding audio downsampling code (JavaScript) and adapting to a different ratio

The vector of coefficients that you gave implement a finite impulse response (FIR) filter with a cutoff of 8 kHz (for an input rate of 48 kHz). As you surmised, these coefficients specify weights that are applied to a sliding window of input samples to generate each output sample; this is the basic operating principle of an FIR filter.

An illustrative way to look at the filter's behavior is via its frequency response. This indicates what frequencies are passed through the filter unchanged and which are attenuated. Generally, if you want to downsample a signal, you first need to bandlimit it so that any frequency content above the post-downsampling Nyquist rate are sufficiently attenuated, to avoid aliasing as much as needed for a particular application.

jojek provided the frequency response in a comment above, which I will shamelessly reproduce here:

Some observations:

This shows a lowpass characteristic, with a band of low frequencies passed through the filter with little to no attenuation, followed by a transition region where the frequency response falls off. It then enters a stopband where all frequencies are significantly more attenuated than in the passband. There are no hard-and-fast objective terms that can be used in all cases to identify exactly where these regions start and end, but you can get an idea by looking at the plot.
This is not a particularly good filter in terms of stopband attenuation (it only attenuates frequencies above the desired band by ~20 dB at minimum), but better specifications would require more coefficients (and therefore more computations per output sample).

Since the filter only attenuates the stopband content by ~20 dB, this data will alias into the signal that you are left with after downsampling. Any content that was in the stopband region in the original signal will be aliased down to some other frequency in the downsampled signal, with its power level decreased by ~20 dB (or more, depending on exactly where in the stopband it is; notice that there are some narrow regions, often called nulls or zeros of the filter, where the filter generates very high attenuation).

This is an equiripple filter, in that the ripples in the stopband have an equal peak height (at just under -20 dB). Based on some quick playing around with Matlab's Filter Design and Analysis tool, it looks like this was designed for a passband from 0-7 kHz and a stopband starting at 9 kHz, with equal weights to the passband and stopband approximation error (these are parameters for the Remez exchange algorithm used to design equiripple filters).

If you wanted to enhance or duplicate this design, you would need to use a digital filter design procedure. As far as implementation goes, these range from pretty simple (e.g. the window method) to quite complex (Remez exchange as I described before).

However, there are a lot of tools already out there to do this for you, both in paid tools like MATLAB and in free open-source tools like SciPy. Here's an example of how you could use SciPy to design a filter for the 4 kHz output rate that you mentioned (you'll also need matplotlib installed to plot the result):

from scipy import signal
import numpy as np

# Sample rate of the input signal (Hz).
fs = 48e3
# Desired number of coefficients in the lowpass filter.
N = 21
# Passband edge (Hz).
fpass = 3e3
# Stopband edge (Hz).
fstop = 5e3
# Generate the desired filter coefficients. The arguments can be interpreted as follows:
#   - the desired length of the filter
#   - list of band edges; this describes a passband from 0 Hz to `fpass` Hz
#     and a stopband from `fstop` Hz to `fs/2` Hz
#   - a response of 1 in the passband is desired (0 dB attenuation), and
#     a response of 0 in the stopband is desired
#   - since we passed frequencies in Hz to the filter design routine, it needs
#     to know the value of the sample rate
h = signal.remez(N, [0, fpass, fstop, fs/2], [1, 0], Hz = fs)

# Normalize the filter to 0 dB gain at zero frequency (which is just equal to
# the sum of the filter coefficients).
h = h / np.sum(h)

# Calculate the frequency response of the filter, then calculate its magnitude
# response, measured in dB.
freq, response = signal.freqz(h)
ampl = 20*np.log10(np.abs(response))

import matplotlib.pyplot as plt
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.plot(freq/(2*np.pi)*fs, ampl, 'b-')
plt.title('Magnitude response of designed filter')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.grid()
plt.show()

Note that digital filters don't have any inherent knowledge of their sample rate. Instead, their behavior is normalized to whatever sample rate you push through the filter. For example, if you took your original filter and passed data into it at 24 kHz instead of 48 kHz, it would have a cutoff commensurate with a 4 kHz output rate instead of the 8 kHz of the original.

With that said, I'm not saying that this is a great filter design. There is noticeable ripple in the passband and the stopband attenuation is relatively poor. There are a lot of parameters that you can adjust with this method, and there are lots of other filter design techniques as well; this should get you a start. As far as a text goes, I would recommend Understanding DSP by Lyons; it's a good introductory book.