Why does circular symmetry of complex envelopes mean that $\mathbb{E} \Big[\mathbf{x}[i] \mathbf{x}^T[i] \Big] = \mathbf{0}$?

A complex random variable with dimension being 2 $Z=X+jY$ can be defined as a Random vector with components $X$ and $Y$ which take values along 2 different dimensions. The Probability Density Function of $Z$ is defined using the PDFs of $X$ and $Y$. The Expectation is defined as $\mathrm{E}(Z)=\mathrm{E}(X)+j\mathrm{E}(Y)$ and the variance is defined the usual way as $Z$ being a complex function of $X$ and $Y$ and the Variance is $\mathrm{VAR}(Z)=\mathrm{E}(|Z^2|)-|\mathrm{E}(Z)|^2$. Intuitively speaking, $Z$ is a complex function of two random variables which just assumes value $X+jY, \forall x,\forall y$. The PDF is written as follows: Write the 2-dimensional Random variable as $$\tilde{Z}=\left[\begin{array}{ccc}X \\Y\end{array}\right]$$where $\tilde{Z}$ is a Real random vector. The PDF of $Z$ is the PDF of $\tilde{Z}$.
Now with Complex Random vectors, $\mathbf{Z}=\mathbf{X}+j\mathbf{Y}$ is a vector of Complex Random variables which has a lot of interesting properties in a sense that it has more to offer than a Vector of Real Random Variables. Particularly of importance are the Covariance, Pseudocovariance(this is something new) and the Expectation. In our situation, we have a special Complex Random vector called the Symmetric  Random Vector $\mathbf{Z}$.
Define the $2n$-dimensional Complex Random Vector as $$\tilde{\mathbf{Z}}=\left[\begin{array}{ccc}\mathbf{X} \\\mathbf{Y}\end{array}\right]$$ where $\mathbf{X}$ and $\mathbf{Y}$ are $n$-dimensional real random vectors. The PDF definition follows the same from Complex Random Variable. You can look up what the Pseudocovariance and Variance mean for a Complex Random Vector. I'm referring to Digital Communications by Proakis here in my definitions and notations.

Now as in the comments, if $\mathbf{Z}$ is a Symmetric i.e the PDF doesn't change on rotation, we can say that $\mathrm{E}(\mathbf{ZZ^\mathrm{T}})=0$. This is a consequence of the following. Since $\mathbf{Z}$ and $\mathbf{Z}e^{j\theta}$ have the same PDF, $$\mathrm{E}(\mathbf{Z})=\mathrm{E}(e^{j\theta}\mathbf{Z})=e^{j\theta}\mathrm{E}(\mathbf{Z}), \forall \theta$$ If we had $\theta=\pi$ (the condition has to be satisfied for all $\pi$), we can see that $$\mathrm{E}(\mathbf{Z})=\mathrm{E}(-\mathbf{Z})=-\mathrm{E}(\mathbf{Z})$$ which clearly says that $\mathrm{E}(\mathbf{Z})=0$. Going further, we can prove that $\mathrm{E}(\mathbf{Z}\mathbf{Z}^\mathrm{T})=0$ the same way assuming a special case of $\theta=\frac{\pi}{2}$. And hence Circularly Symmetric Complex Random Vectors are Zero-mean and Proper. The last equation left unproven here is the definition of being proper.


Filter response to sinusoidal input signal

The fact that a sinusoidal signal's amplitude is multiplied by the magnitude of the system's frequency response (evaluated at the input frequency) is a general property of LTI systems, regardless of whether they are continuous-time or discrete-time systems.

This is a consequence of the complex exponential being an eigenfunction of LTI systems, and the value of the frequency response evaluated at the input frequency is the corresponding eigenvalue.

For more details on eigenfunctions of LTI system check out this question and its answers. And I'm sure that there's even more information here on this site.


how we can prove that the Variance of (white Gussian noise) is average power of the noise?

Since this is clearly a homework question:
Simply by applying the definition of variance. Remember that white noise has a mean of 0.


Derivative of signal with missing samples

See Smooth noise-robust differentiators. While Pavel doesn't address your specific problem there, he and I corresponded about a similar "missing samples" problem earlier this year, and he has some good ideas that you can try. In particular, he recommends polynomial approximation in the region of the missing samples. In addition I'll mention that if you use Neville's Algorithm you can guarantee that the polynomial passes through your "good" sample values. You can find Pavel's contact info in various places on his Website.

Greg Berchin


Phase offset and timing offset are two different things with two independent tracking loops to resolve. A Phase Lock loop eliminates carrier phase offset (and a frequency offset is a carrier offset changing with time). We can have phase offset with no timing offset, and we can have timing offset with no phase offset.
Once carrier offset is removed, I recommend creating an upsampled eye diagram to evaluate the signal quality and see if further timing recovery is all that is needed. The diagrams below specific to 16 QAM illustrate the difference between carrier offset and timing offset:
Here is the eye diagram of the real portion of the 16 QAM signal, where we see the four levels and we see the signal is sampled at 4 samples per symbol. The signal quality is excellent as given by the wide eye openings and small variation at the ideal sampling locations. However given the present sampling locations there is a clear timing offset where further timing recovery is needed. There is no phase offset in this example.

For comparison, after timing recovery the eye diagram looks like the following:

If we plot only the samples at the correct sampling locations (one sample per symbol) on an an IQ diagram to get the constellation, the above case including both real and imaginary values would be at the orange dots in the following graphic:

If there was a phase error in contrast to timing offset, the entire constellation would be rotated by the phase error, demonstrating phase error with no time error:

Note that if there is phase error even with no timing error the eye diagram will appear very corrupted still (often there is a phase rotation due to frequency offset which makes the eye unrecognizable, these must be corrected first before evaluating signal fidelity with an  eye diagram). Below shows the same eye as was done above for 16 QAM with no timing error, but with the phase offset as shown in the previous constellation:

A Gardner Loop is my recommended go-to for timing offset correction in that it can converge and remove timing errors even in the presence of very large carrier offsets. Then the carrier offset can be determined by simply measuring phase rotation from symbol to symbol using the correct sampling points as determined by the Gardner timing recovery loop.
More details on the Gardner Loop are provided here:
Gardner Timing Recovery for Repeated Symbols
And more details on Carrier Recovery are provided here:
High modulation index PSK - carrier recovery


Crash course on turning engineering formulae into computer code?

I don't usually recommend this series, but Linear Algebra For Dummies might be right for you. Chapter 19, for example, covers all the Greek symbols.

You should also complement this with an introductory book on DSP, like Applied digital signal processing : theory and practice and DSP for audio An Introduction to Audio Content Analysis: Applications in Signal Processing and Music Informatics


Why is it possible to frequency sweep with just one sinusoid?

Thanks to geometrikal and Jason R, the sum of exponentials is equal to zero when $\omega_0 = \frac{2\pi}{N} \hat{k}$ and $0 \leq \hat{k} \leq N-1$, which is obvious when you realize that it's just the sum over $2 \hat{k}$ periods of $e^{-j2\omega_0n}$, where $T = \frac{2\pi}{2\omega_0} = \frac{N \pi}{2\pi \hat{k}} = \frac{N}{2 \hat{k}}$.

For $\hat{k} > N-1$, I suppose you would have to zero pad, which is where the digital sinc Jason R. mentions comes in.


Detecting tennis courts in aerial images?

Since you are trying to discover structure in the images it's better to work in grayscale.

This is a very nice case where the court appears to be a nice rectangle, but in general, courts might come in different sizes and orientations. Also, white lines on green background is not a general rule, consider for example Roland Garos.

Having said this, you can try thresholding your image (you might have to use adaptive thresholding for this) to get something like this:

You can then go through the image creating chains of pixels that are 4 or 8 connected

Since you are looking for lines, you can get each of the chains and fit a line to it.

At this point, the "interesting pixels" in your image have been transformed to a set of line models and some of them will be intersecting. You can now apply simple rules to discover which sets of lines are intersecting (here is a good example ).

To discover the set of intersecting lines that make up a tennis court you would search for (this is an example) "A set of 9 lines where 4 are intersecting at right angles, two more parallel to the long lines intersecting them at 1/8 of the width of the field, two more parallel to the short sides and lying about 1/5 from the beginning of the court and a third running in the middle of the court and intersecting the two shorter ones"...In other words, you will have to express the tennis court as a set of intersecting lines with relative positioning to each other (so that your expression is free of the absolute dimensions).

For a very good application of the above technique you can have a look at this paper that outlines the recognition of shotcodes in an image.

For the theoretical background and more information about these techniques you can refer to the excellent "Digital Image Processing With Matlab" and specifically, chapters 10-13 which are the most relevant to your problem

I hope this helps.


Fourier Transform of image convoluting with kernel

The function
$$
\frac{\sin(x)}{x}
$$
has a well-known Fourier transform:
$$
\int_{-\infty}^{\infty}  \: \frac{\sin{x}}{x} e^{-j 2 \pi f x} dx = \begin{cases} \pi & |f| < \frac{1}{2 \pi} \\ 0 & |f| > \frac{1}{2 \pi} \\ \end{cases}
$$
which is derived here.

Your function:

$$
\frac{\sin(x)\sin(y)}{xy}
$$

is a simple separable two dimensional version of it. The 2D Fourier transform will just be:

$$
\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}  \: \frac{\sin{x}}{x} \frac{\sin{y}}{y} e^{-j 2 \pi f_1 x} e^{-j 2 \pi f_2 y} dx dy =
\int_{-\infty}^{\infty} \: \frac{\sin{x}}{x} e^{-j 2 \pi f_1 x}  dx \int_{-\infty}^{\infty}  \:\frac{\sin{y}}{y} e^{-j 2 \pi f_2 y}  dy =\begin{cases} \pi & |f_1| \mbox{ and } |f_2|  < \frac{1}{2 \pi} \\ 0 & |f_1| \mbox{ or }  |f_2| > \frac{1}{2 \pi} \\ \end{cases}
$$

Below is a plot that takes a random image with a white square in the middle (top left), FFTs it, multiplies the FFT by a brick filter (top right), inverse FFTs it and plots the result (bottom right).

R Code Below

# 26383

img <- runif(128*128, 0, 1)
dim(img) <- c(128, 128)

img[32:96,32:96] <- 1

N <- 50

imgFFT <- fft(img)
imgFFT[N:(128-N+2),] <- 0.0000000000000
imgFFT[,N:(128-N+2)] <- 0.0000000000000
img2 <- fft(imgFFT, inverse = TRUE) / 128 / 128

img2[Mod(img2) > 1] <- 1

mx <- max(log(Mod(imgFFT + 0.00001)))
mn <- min(log(Mod(imgFFT + 0.00001)))

scaled_log_imgFFT <- (log(Mod(imgFFT + 0.00001)) - mn)/(mx - mn)

par(pty="s") # make the plit square
plot(0:256, type='n')
rasterImage(as.raster(img),0,128,128,256)
rasterImage(as.raster(scaled_log_imgFFT),128,128,256,256)
rasterImage(as.raster(Mod(img2)),128,0,256,128)


modulate/demodulate audio file as FSK signal mixed with music signal

So, your unmixed secret signal has energy content in a frequency range of ca -6 kHz to +6 kHz.

As far as I understand, you mix that with a real cosine of ca 19 kHz, so that it now has signal content in -25 kHz to -13 kHz, and +13 kHz to +25 kHz. All things are still nice and symmetric – so if you mix it back down to 0 Hz, you'd still have the original spectrum.

Only that your 48 kHz sampling rate doesn't even allow your to unambigously represent a signal above 24 kHz – it instantly aliases into the lower end, and vice versa, so you have distortion in your signal of interest (still shouldn't be hearable).

But then you high-pass filter at 500 kHz above the frequency you shifted to – but that cuts off 5.5 kHz of your signal! So, afterwards, your signal is no longer symmetrical in spectrum, and mixing it back down will produce non-constructive aliases.

So, to understand what's happening here, plot your secret signal's spectrum after being shifted up. Don't use the logarithmic frequency axis, but simply linear frequency.


Computing N-point inverse DFT with N/2-point forward DFT

Just follow your recipe as written

Conjugate your frequency domain vector
Apply DFT on even samples
Apply DFT on odd samples
Recombine the two $N/2$ DFTs to create the $N$ point DFT result
Conjugate the DFT result and divide by N.

You can swap the last two steps, if you want.

Keep in mind that there is rarely a good reason to split an $N$ point DFT into two $N/2$ point DFTs. Because of the recombination it's typically more as expensive a single full size DFT. In essence you just take one stage out of the FFT algorithm and compute it manually. The only reason to go there if you have severe memory constraints or cache management issues.

That's different from using a $N/2$ point complex FFTs to compute an $N$ point FFT on real value data. That's generally a considerable savings.


Don't discount the contribution of physics to the blur in your pictures: Especially in low light, you need to have long exposure times, to get a picture that isn't completely dominated by grain. Conversely, things become less sharp if they're moving.

In fact, security cameras tend to have low frame rates (that might actually be for historical reasons – because nobody actually cared about the lawn growing in 25 fps, or about getting 100 grainy shots of someone smashing a window in 4 seconds, instead of 10 shots that have better exposure). So, that might actually be the limiting factor here.

  HVEC was developed for entertainment applications mostly, wasn’t it?

Certainly, the highest money flow would be from the movies industry, but they have the same problem, namely that action scenes are harder to compress than still live shots of landscape. So, HEVC tries very hard to be better at accomodating different types of video content, so, I'd assume it's the best standard solution.

All in all, HEVC does a lot of things to make action scenes work best – motion prediction on a sub-pixel level, high dynamic range, you name it.

Especially, it compresses still-standing parts very efficiently, leaving more headroom on storage space when there's movement. Of course, the encoder has to be configured to actually spend the bitrate exactly when necessary. Maybe your vendors misconfigure it. Impossible to tell. Maybe look at the instantaneous bitrate over time; if it doesn't change much when there's something appearing in the picture, the encoder was configured to aim for a constant bitrate rather than a constant quality.

Of course, there are numerous approaches to making surveillance video capture more efficient. (I wouldn't be surprised if a lot of that wasn't available for civilian purposes, though. Automated video surveillance is very much something that states do.) I know that several "smart" home cameras draw rectangles around areas they consider noteworthy based on some more or less artificially intelligent algorithm.

Source: Pyimagesearch/Adrian Rosebrock

Of course, if you attach enough neural network processing power to your camera, and train accordingly, you can make sure your system reliably doesn't even store background, but recognizes e.g. cats and people in the picture.

Source: Micah Scott (Scanlime): Tucoflyer Smart DNN-based cat tracking indoor camera gimbal system

If you know what you don't care about and that you care about the rest, and only record the latter, you can encode more efficiently by not even trying to record your grass growing. But: this takes a solid amount of processing power and hence electricity.

Taking that thought further:

Video compression is literally nothing but taking a high-dimensional input and mapping it to a lower-dimensional description and quantizing that cleverly.

That is also a pretty good description of what a video object classifier neural network does. Thus, there's plenty of approaches to encode video through neural networks.

You might then store video regions of interest not even as "collection of picture information", but simply as coefficients necessary to recreate e.g. the same face from a neural network optimized at representing faces by feature vectors with little error. This works pretty well. In fact, so well, that by synthesizing such feature vectors, we can create convincingly real-looking faces. For example This Person Does Not Exist:

It seems at first that the quality as evidence of a face picture recreated from a camera neural network, which you cannot prove to be right, is very limited. Or at least, should be, in a court of law.

  My algorithm says this is the visual representation of the lowest-error feature vector of a face that would match what the camera saw at 3 a.m. when someone broke into my car. I cannot formally prove an upper bound for the error, but it looks just like the defendant. So, it's the word of my trained-on-thousands-of-pictures neural network against the word of the defendant! If that doesn't prove him guilty, what will?

shouldn't be an all too successful argument when in dubito pro res is a principle in criminal justice.

Then again, all video compression claims the same; got a JPEG picture of the thief breaking into your store done with your iphone? These JPEG coefficients just are the lucky coincidence of image content, shot noise and the JPEG encoder.


Fourier transform of discrete time unit step function

The DFT of a unit step response is
$$U(\omega) = \frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega)$$
Applying the shift property as you did will give:
$$\mathcal{F}(u[n] - u[n-1]) = U(\omega) - U(\omega)e^{-j \omega} = \frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega) - [\frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega)]e^{-j \omega}$$
that is
$$\mathcal{F}(u[n] - u[n-1]) =  \frac{1- e^{-j \omega}}{1 - e^{-j \omega}} +\pi\delta(\omega)( 1 - e^{-j\omega}) = 1+\pi\delta(\omega)( 1 - e^{-j\omega})$$
The second term is always zero because for $\omega = 0$, $1 - e^{-j\omega} = 0 $ and it's zero on any other point, so you get,
$$\mathcal{F}(u[n] - u[n-1]) = 1 = \mathcal{F}(\delta[n])$$


Template matching in the link below will work as a solution.

https://www.docs.opencv.org/2.4/doc/tutorials/imgproc/histograms/template_matching/template_matching.html

The downside is you will need templates of all your valves handles(coloured or not), also the templates and images will have to be the same resolution.
Based on whether the template or the template rotated by 90 degrees matches you can come to a decision on the valves state.


Why is the Z-transform so important in digital filters analysis and design?

The Z Transform is to discrete-time (digital) signals precisely the same role that the Laplace Transform is to continuous-time (analog) signals.

Linear Time-Invariant (LTI) Systems (a.k.a. "filters"), are made up of signal-processing elements that fall into 3 fundamental classes:

adders (devices that add two signals).
scalers (devices that scale a signal by a constant).
"reactive" elements (devices that are able to discriminate w.r.t. frequency).

Element classes 1. and 2. are essentially the same for analog or digital filters.  The are sometimes called "memoryless" devices or elements.

For an analog filter (or "analog LTI system"), those reactive elements would be capacitors or inductors.  They integrate ($s^{-1}$) one signal to become another.  That turns a sine signal into a cosine signal or shifts the phase by $\pm$ 90°.

$$\cos(\Omega t) = \sin(\Omega t + \tfrac{\pi}{2})$$

For digital filters, the reactive elements are delay elements.  A unit delay (a delay of exactly one sample period $T$) will delay any signal, including a sinusoid, by 1 sample or $T$ units of time ($z^{-1}$).  That shifts the phase by an amount that is dependent on frequency

$$\sin(\Omega (t-T) ) = \sin(\Omega t - \Omega T)$$

or

$$\sin(\omega (n-1) ) = \sin(\omega n - \omega )$$

Any LTI system that acts as a "filter", a device to filter out some frequency components and leave others, must have reactive elements (or "non-memoryless" elements or components having memory) in order to discriminate one frequency from another.  And such a filter will shift phase which will normally be different for different frequencies.  But a memoryless LTI system (which is just a scaler) will not discriminate between frequencies nor will shift phase, except for possibly by 180°, which is just a polarity reversal or scaling by a negative constant.


The need for all these coefficients may not be obvious if you are looking at only real number input to the DFT.

But a complex vector input of length N will have twice the degrees of freedom compared to a strictly real input of that length.  The DFT, as an invertible linear transform therefore must therefore retain that same number of degrees of freedom, which for complex input ends up in the relationship between the positive and "negative" frequencies, or upper half of the transform, to produce complex exponential frequency results, instead of just strictly real cosine and sine functions.

For strictly real valued input, you don't really need to "care" about the other half, as it is duplicated information in the form of a complex conjugate (where the imaginary components of positive and "negative" frequency bins all cancel out to zero).  Thus the "mirroring" reduces the number of degrees of freedom to match that of the strictly real-valued input.  But one might still use the other half for computational purposes, as it makes the matrix math more regular and symmetric, and thus optimizable in the form of some FFT algorithms.

Thus the negative frequency coefficients don't need to be computed because they are redundant (unless the input is complex and not strictly real valued). In that case, the frequency coefficients above N/2 don't need to be computed because they are identical to the negative coefficients that don't need to be computed


Is it possible to hack camera calibration without having access to the camera?

If I understand correctly, you don't need the intrinsics or extrinsics to achieve that, if a top-down view is all you want. You could basically define 4 points on your parallel lines and then warp the entire image into a canonical view (say $\{\{0,0\}, \{480,960\}\}$).

To do that in OpenCV, all you need to do is compute the homography using findHomography as described here. Then, use warpPerspective (look here) to warp your images.

The only thing you need to care about is the selection of the dimensions of the fronto-parallel view that you are after. It should more or less match the aspect ratio of the original image region.

Regarding the second part of your question: Yes. First of all, let me remind you that many works such as this choose to ignore the distortion estimation. But, I can think of two ways to accomplish that:

Use the road lanes itself to correct the distortion: In the fronto-parallel view, create another perspective warp so that the lanes themselves would be vertical. In this setting, think about the lane lines as the vanishing lines in the image (as in a perspective view, they always do). From that, you can even calibrate your camera. Look here, here and here for example.
The cars should be following more or less a straight/smooth trajectory. If you compute the motion of all the cars, they should give you similar motions, parallel to the lanes. That information is again useful in constraining and generating a transformation, which would warp your image to a reduced distortion one.


WOLA vs Simple Hamming Window on signal before FFT

I'll try to briefly answer your questions.

No, WOLA (weighted overlap and add) is applied directly on the result(s) of an inverse FFT operation. The IFFT results ( time segments ) are first multiplied by a particular window function (e.g Hamming ) and then overlap-added to produce the final time domain signal ( waveform ).
As explained in #1, WOLA is used to generate a time domain signal from multiple, overlapping ifft segments. If you have only a single fft segment (from your FFT analysis), then there's no need for (W)OLA at all.
Yes


My impulse response does not tend to zero

Note that the DFT transforms a periodic sequence into a periodic sequence. Hence, the result of your IFFT is intrinsically periodic and you should view it as such. The parts you are seeing at the end of it are thus just as well behind the start as they are before it.

Hence, this part may as well be some form of pre-ringing, stemming from a nonzero group delay of your filter.

Note also that the FFT implicitly periodifies your spectrum. It looks a little cut, so you may be missing some essential parts.


use the quadratic formula.  you know:

$$ a \, x^2 \ + \ b \, x \ + \ c \ = \ 0 $$

$\implies$

$$ x = \begin{cases}
 \frac{-b \, \pm \sqrt{b^2 - 4 a c}}{2\,a}  \quad \text{if } b^2 \ge 4 a c \\
  \\
 \frac{-b \, \pm j \sqrt{4 a c -b^2}}{2\,a}  \quad \text{if } b^2 < 4 a c\\
\end{cases}$$


To make it clear, you padded your image in what domain ?

Time domain >> No problem, IFFT and then truncate will bring you the original image. It's still intact.
Freq domain >> Padding in the Freq domain means image resizing. You cant get your original image with 100% anymore. If you really want the original size one, the only thing you can do is to resize it back with some resampling methods.


Your understanding is not correct. The baseband processing maps the message to a complex-valued signal

$$x(t)=R(t)e^{j\phi(t)}\tag{1}$$

or, equivalently, to an in-phase component $I(t)=\text{Re}\{x(t)\}$, and a quadrature component $Q(t)=\text{Im}\{x(t)\}$. It is not clear from the figure how $R(t)$ and $\phi(t)$ are related to the message signal $m(t)$. The modulated signal is constructed by multiplying $(1)$ by a (complex) carrier, and taking the real part:

$$s(t)=\text{Re}\{x(t)e^{j\omega_ct}\}=R(t)\cos(\omega_ct+\phi(t))\tag{2}$$


Extracting motion vectors from H.264 bitstream

I think the following code can help you:
https://github.com/vadimkantorov/mpegflow
The main idea is that you need to decode frames first and query motion frames by calling av_frame_get_side_data function.


"Ensemble averaging ... cannot track dynamic changes"?

You're misinterpreting ensemble averages.... this is really in "batch" like after recording a few runs. You average functions on a time interval, e.g. that's where those doctor's charts for kids' height/weight mean & spread over time arise. Time-averaging is totally different...


How can I convert a step response that is described by the sum of two exponents to an IIR filter analytically?

You can use the step-invariant transformation, which makes sure that the step response of the discrete-time system matches the step response of the continuous-time system at the sampling instants. I suppose this is what you mean by "exact solution".

Sampling the given continuous-time step response at $t=nT$ gives the desired step response of the discrete-time system:

$$g[n]=y(nT)=\left[1-Ae^{-nT/T_1}-(1-A)e^{-nT/T_2}\right]u[n]\tag{1}$$

where $u[n]$ is the unit step.
The $\mathcal{Z}$-transform of $(1)$ is

$$G(z)=\frac{1}{1-z^{-1}}-\frac{A}{1-e^{-T/T_1}z^{-1}}-\frac{1-A}{1-e^{-T/T_2}z^{-1}}\tag{2}$$

The frequency response of the desired discrete-time system is then given by

$$H(z)=(1-z^{-1})G(z)\tag{3}$$

After performing some basic algebra and using $a=e^{-T/T_1}$ and $b=e^{-T/T_2}$ we obtain

$$H(z)=\frac{[1-b+A(b-a)]z^{-1}+[a(b-1)+A(a-b)]z^{-2}}{1-(a+b)z^{-1}+abz^{-2}}\tag{4}$$

As a result, the IIR filter with transfer function $(4)$ has a step response that equals the given continuous-time step response $y(t)$ at the sample instants $t=nT$:


So the problem is that when simulating FMCW signals, there is some ambiguity in the convention on the mixing order. MATLAB's documentation makes reference to it here:

"The mix operation reverses the Doppler shift embedded in x, because
of the mixing order of xref and x. The mixing order affects the sign
of the imaginary part of the output argument, y. There is no
consistent convention in the literature about the mixing order. This
function and the beat2range function use the same convention. If your
program processes the output of dechirp in other ways, take the mixing
order into account."

So you're amost there. MATLAB's functions account for this so you need to account for the change yourself when doing the manual range-Doppler map. You can do this by applying a complex exponential to each pulse that shifts the phase in the other direction so that it will produce the right Doppler shift.
For example, in the simulation loop right before you assign the dechirped signal to xr do something like the code below. I currently don't have the Phased Array Toolbox installed so I did some of the calculations manually:
for m = 1:Nsweep

    % After all the target sim stuff do your own Doppler calculations
    rangeToTarget = norm(tgt_pos - radar_pos);
    radialVel = dot(tgt_vel - radar_vel, tgt_pos - radar_pos)/rangeToTarget;

    % Negate Doppler so that closing targets produce positive frequencies
    fd = -2*radialVel/lambda;

    % De-chirp and apply phase shift for this pulse to reverse Doppler
    % to the correct value
    xr(:,m) = dechirpsig*exp(1i*2*pi*(2*fd)*(m - 1)*tm);

end

This uses the simple "stop and hop" model where the target moves in discrete time steps which are equal to the sweep time tm.
Addressing Range Issue
I totally to forgot to address the range part of your question. The reason this is happening is because when you use MATLAB's functions to generate the range-Doppler map, you are giving it a sampling rate $f_s = 4 \space GHz$, which in turn yields a range-bin (or range-cell) size $\delta R $ of
$$\delta R = \frac{c}{2B} = \frac{c}{2(4e9)} = 0.0375 \space m$$
If you had set RangeFFTLength = numel(waveform) when doing your own 2D FFT, you would have a HUGE number of samples in the range dimension, but it would be correct given the achieved range-bin resolution $\delta R$ and will give you the correct target ranges.
The problem now comes that when you use less range samples than what you started with by setting RangeFFTLength = 2048. This decimates the effective sampling rate, so your new range-bin size is now going to be artificially larger. If you take the ratio of the number of samples produced by using $f_s = 4 \space GHz$ and the size of the range FFT of 2048, your effective sampling rate is now 98 times worse! So your new range bin size becomes
$$\delta R_{new} = 98\frac{c}{2B} = 98\frac{c}{2(4e9)} = 3.675 \space m$$
Applying this range-bin size to the peak range index should yield the correct range in your range-Doppler map. However, note that this is still wrong. I say that because the whole point of putting that much bandwidth into a chirp is so that each range-bin in the range-Doppler map is equal to to 0.0375 m and NOT 3.675 m.
If you did it correctly, you would be able to sample at a much lower rate since your beat frequencies after mixing have lower bandwidth requirements, giving you a smaller and reasonable number of samples. Each range-bin size is still 0.0375 m, preserving the range resolution desired.


Determining Power of White Noise Component of an Incomplete Time Signal

I assume you have information about the signal itself (Not the noise).

What you can do is use LS to estimate the signal (Easily compensate for the missing samples in the Model Matrix).
Once you have it estimated you can do 2 things (Which are the same):

Recreate the samples using the model. Subtract the measurements. You'd be left with IID Samples of the noise -> Use the empirical STD / Variance estimator.
If it is Linear Least Squares use its internal estimation of the noise Variance.

Good Luck!


Interpolation formula for two dimensional signal reconstruction in the frequency domain from polar samples

I haven't worked through the algebra, but this page just before that equation:

suggests that
$$
\operatorname{sinc}\left[2A(\rho-\frac{n}{2A})\right]
$$
might be the correct form of the $\operatorname{sinc}$ term.


Is it necessary to scale a signal for quantization noise analysis?

No, you don't need to scale the stimulus up to the full bit depth. That would actually end up underestimating the effect that roundoff noise would have on your actual system, as you would be artificially increasing the signal power (thus lowering the predicted signal to quantization noise ratio).

It makes sense to simulate the filter with a range of inputs that covers the space that you're designing the system to operate with. If the inputs will always be in the range $[-1, 1]$, then simulate using that range. On the other hand, if it's very likely that the data will be in that range, but possible that it might be larger, then it would be appropriate to test at higher levels. This requirement would likely be driven by higher-level considerations in your design.


Is there any reason that cross correlation would perform well or poorly on sparse binary arrays?

If you upsample your array to a resolution more than the time uncertainty, or measurement time quantization plus jitter noise, by just zero-filling between events in your event count data, you are adding high-frequency jitter noise that can interfere with correlation results.  You could try upsampling your arrays by using a Gaussian distribution of the appropriate width for each event, instead of just zero-padding between events.  Or, if you have a better statistical distribution for the measurement time offsets/errors, use that.


Wavelet transform and FFT using to extract feature power bands with EEG signals

Alpha rhythm
Generally, alpha-band oscillatory activations (8-10Hz) relate to relaxation and in principle accompanied with closure of eyes. This is the prime marker that is used to detect drowsiness, but surely not the only one (see alpha dropout, NREM1, eye-rolling EEG artefacts).
Event detection
To detect alpha, electrodes from occipital or adjacent regions are used. Due to low SNR in EEG signals and its variance across EEG systems, you may want to try different detection methods to decide for your classifier.


Questions about autocorrelation, correlation, and randomness

1) g1 assumes that there is a non-zero mean and that always needs to be subtracted off for the correlation calc to make sense.  g2 assumes that the mean is zero. g2 is incorrect if the mean is not zero.

2) the mean is taken over the whole sample because the autocorrelation calc has the underlying assumption that the process is stationary so the mean doesn't change during the process.

3) As the lag gets larger, then, if the process is stationary, then the autocorrelations should doe down approach zero. Note that there are "long memory" processes that don't have this characteristic but I don't think that they are
considered stationary processes. Someone feel free to correct me if I'm wrong about long memory processes not being considered stationary.


Convolution of delta function in discrete time with parameter

You can use the following argumentation to find the result. The discrete time unit-sample function $\delta[n]$ has the following property for integer $M$:
$$ \delta[Mn] = \delta[n] $$ and more generally you can conlcude that for integer $M$ and $d$ we have
$$ \delta[M(n-d)] = \delta[n-d] $$

Therefore you can replace $\delta[5n-20] = \delta[5(n-4)]$ with $\delta[n-4]$ and proceed as usual to find the result of the convolution as:

$$y[n] = n\delta(5n-20) \star 3^{5n}u[n] = n\delta(n-4) \star 3^{5n}u[n] = 4 \cdot 3^{5(n-4)}u[n-4] $$

where in this last line we have used the sifting and shifting properties of the impulse function.


Introduction to IP : Laplacian of Gaussian?

You are right about the $x$ and $y$ values. The center of the matrix is $(0,0)$ and the corner points are $(\pm 4,\pm 4)$. But they obviously wanted integer values in the matrix, so they simply scaled the LoG function by a factor of $482.75$ (just to get a decent range). Evaluating the function with this scale factor gives you (for the lower right quarter of the total matrix):

LoG =

  -40.000000  -23.086993    0.294243    5.218347    2.080704
  -23.086993  -11.762411    3.077873    4.839310    1.745668
    0.294243    3.077873    5.409023    3.361998    0.998325
    5.218347    4.839310    3.361998    1.456020    0.365518
    2.080704    1.745668    0.998325    0.365518    0.081641

Rounding should give you the final result. However, if you look closely, you'll see that it doesn't (at least not for all $(x,y)$ values). If you check the matrix on this page, you'll see that for the same $\sigma$ it is also different. So people make mistakes and they just copy from each other.


The usual way to do this is to interpolate using surrounding points.  There is a simple 6-point mechanism that uses:

This question on the Mathematics SE site has some more details. Note the answer's correction to the equation in the question.


Filtering 50Hz using a notch filter in MATLAB

I'm not sure what iirnotch does, but this is how to design the notch filter by hand.

fs = 20000;             % sampling rate
f0 = 50;                % notch frequency
fn = fs/2;              % Nyquist frequency
freqRatio = f0/fn;      % ratio of notch freq. to Nyquist freq.

notchWidth = 0.1;       % width of the notch

% Compute zeros
notchZeros = [exp( sqrt(-1)*pi*freqRatio ), exp( -sqrt(-1)*pi*freqRatio )];

% Compute poles
notchPoles = (1-notchWidth) * notchZeros;

figure;
zplane(notchZeros.', notchPoles.');

b = poly( notchZeros ); %  Get moving average filter coefficients
a = poly( notchPoles ); %  Get autoregressive filter coefficients

figure;
freqz(b,a,32000,fs)

% filter signal x
y = filter(b,a,x);


If you are trying to remove the gaps during hough transform, you can set minLineLength(Minimum length of line. Line segments shorter than this are rejected) and maxLineGap (Maximum allowed gap between line segments to treat them as single line) in cv2.HoughLinesP()

EDIT

Some results from the morphological operations:


Effect of Adding the cyclic prefix on the toeplitz matrix in OFDM

What you mean might be circulant matrix instead of toeplitz matrix. See section 3.4.4 in https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter3.pdf about how the circular convolution in OFDM is represented by matrix operations (eq 3.130 onwards).

First, in almost all standard OFDM systems, you can assume $D \le L$. The cyclic prefix will be less than or equal to maximum multi-path delay, so as long as $D \le L$, the linear convolution $x * h$ gets converted circular convolution of $x$ and $h$.

When $y_c = H x_{cp}(k)$, you only need to take $N$ original elements for $x$.  $H$ is size $N \times N$. This is because you have already re-written the circular convolution in matrix form. Each row of $H$ will do dot-product with $x$ to generate $y_c[n]$. Like that there will be $N$ values of $y_c$ corresponding to each row of $H$.

$y$ is of length $N+L +D-1$. As you correctly mentioned $y_c = y(D+1:D+N)$.

Equivalent $H$ (size $N \times N$):

h(0)  0    0 ... h(L-1)  h(L-2) .. h(1)

h(1)  h(0) 0     0  ... h(L-1) h(L-2).. h(2)

h(2)  h(1)  h(0)  0 ..  0 0 h(L-1) h(L-2) .. h(3)

.

.

.

0 0 0 ...    h(L-1) h(L-2) .. h(2) h(1) h(0)

MATLAB command to generate above $H$:

Assuming $h$ is the vector representing channel having $L$ non-zero values.

H = toeplitz([h zeros(1,N-L)][h(1) zeros(1,N-L) flip(h(2:L))])

For a small example where $N=4, L=2, D=1$

$x = [x0, \,x1,\, x2,\, x3]$

$h = [h0\,,h1]$

$y = [x3h0\;, x3h1+x0h0\;, x0h1+x1h0\;, x1h1+x2h0\;, x2h1+x3h0\;, x3h1]$

$H$

h0 0 0 h1

h1 h0 0 0

0 h1 h0 0

0 0 h1 h0


assuming you have a relative high frame rate video (the object is moving relatively slow between two consecutive images) some standard tracking algorithms can be applied.

Mean shift [Wikepedia - Mean shift]
Use some template of the object you are tracking to create a likelihood function for the object location and track the object.

if you need more advanced method that can be more robust, e.g. object occlusions you can use more advanced methods:

kalman filter [Wikipedia - Kalman filter]
In short kalman filter uses a dynamic model that is composed of
prediction step and update step to track an object in an unsertain
enviroment.

More advanced option you can consider, if you have couple of objects you want to track at the same time is particle filter


Convert measured spectrum to impulse response

Can I just write: X = 10.^(mag/20) .* exp(jphi/3602*pi) ?

Yes. If your software returns the phase in degrees, that's exactly right.
The tricky part is, however, the interpolation. The phase wraps round at 360 degrees. For example, if you have 1 degree and 355 degree, the correct interpolation would 0 degree (since 355 is the same as -1) but regular interpolation will give you 180 degrees which is about as wrong as you can get.
In this case unwrapping the phase often helps. You would start somewhere in the middle of the frequency spectrum where you have good signal to noise ratio and than unwrap up and down from there. Unwrapping means to look for phase jumps of more than 180 degrees and add/subtract multiple of 360 degrees to the phase until the jump is less than 180 degrees.
Another potential problem is bulk delay in the measurement which is fairly common for acoustic measurements. Bulk delay leads to very steep phase gradients which makes both unwrapping and interpolation difficult. This can be addressed by "eye balling" the bulk delay, take it out of the measurement. Then do the interpolation and inverse FFT and put the bulk delay back into the impulse response.


Auto Correlation vs Cross Correlation vs Convolution and their applications

I can tell you of at least three applications related to audio.

Auto-correlation can be used over a changing block (a collection of) many audio samples to find the pitch. Very useful for musical and speech related applications.

Cross-correlation is used all the time in hearing research as a model for what the left and ear and the right ear use to figure out a sound's location in space (this is called sound source localization). In the case of two microphones you would cross-correlate the left channel with the right channel.

Convolution is used in simulating reverberation. A room's impulse response can be determined from measurements and that impulse response can be convolved with any sound source to simulate the reverberant response (at the impulse response recording's exact location).

I know this answer isn't complete but maybe it can give you some idea that there is in fact a practical use for auto- and cross- correlation!

So in general, auto-correlation can be used to extract properties of a signal, cross-correlation can exploit the information between two related signals, and convolution can be used to modify the properties of an incoming signal based on some time, frequency, and phase response specified by the impulse response you're convolving the source with.


How to use Kalman filter for altitude prediction based on barometer data?

You might get good results when you consider the acceleration as an input $u$, so the model could then be written as

$$
x[k+1] =
\begin{bmatrix}
1 & \Delta t \\
0 & 1
\end{bmatrix} x[k] +
\begin{bmatrix}
\frac{\Delta t^2}{2} \\ \Delta t
\end{bmatrix} u[k] + w[k]
$$

$$
y[k] =
\begin{bmatrix}
1 & 0
\end{bmatrix} x[k] + v[k]
$$

where $x[k]$ a vector with the position and velocity and $w[k],v[k]$ zero mean Gaussian white noise. Here $v[k]$ has a covariance $R$ equal to the variance of the barometer and $w[k]$ probably has the following covariance $Q$

$$
Q = \sigma_u^2\,
\begin{bmatrix}
\frac{\Delta t^2}{2} \\ \Delta t
\end{bmatrix}
\begin{bmatrix}
\frac{\Delta t^2}{2} \\ \Delta t
\end{bmatrix}^\top
= \sigma_u^2\,
\begin{bmatrix}
\frac{\Delta t^4}{4} & \frac{\Delta t^3}{2} \\
\frac{\Delta t^3}{2} & \Delta t^2
\end{bmatrix}
$$

where $\sigma_u^2$ the variance of the measured acceleration (this assumes that there are not disturbance forces acting on the system).


The I/O relation
$$y[n] = \cos(x[n])$$

represents a stable system.

Because by definition of $\cos(x)$, irrespective of the input $x[n]$, its output will always be bounded by $$-1 \leq \cos(x[n]) \leq 1 $$ for all $x[n]$ for all $n$. Even if $x[n]$ is undefined ( such as would be the case with $x[n] = 1/n$ at $n=0$ as Marcus pointed out) we for sure know that the system output is bounded (the value of input is undefined, but the output is known to be within $\pm1$.)


Matlab: What should be the BER performance for BPSK using Constant Modulus Algorithm equalizer

Well, I took a look at your code, and spent a lot of time on it, and I have discovered some mistakes, some practicals and some theoretics. Here are my answers:

(1) You can use the function filter. Just remember that it returns a vector of the same length of the input, while conv returns the $ N_x + N_h - 1 $, where $ N_x $ is the length of your input and $ N_h $ is the length of your channel. This difference needs to be take in consideration when you create your noise vector. You won't be able to add the noise vector to your signal if you don't pay attention to this. Read the documentation to get more information.

(2) As I said, there is a lot of mistakes at your code:

First I changed the vector noise to only add a real number to the signal, and not a complex number.
Second, I changed the equation

w = w - mu * e(i) * yy * x;

to

w = w - mu * x * e(i) * yy.'

that I took from here.

There is some mistakes in your CMA algorithm. Here are some lines and the coments:

yy = w(:,i).'*x; Remeber to conjugate w in this. In BPSK we only have real numbers so this is not a problem, but if you had complex numbers, you must conjugate w.

e(i) = abs(yy)^2 - 1; In your original code, you didn't use abs(y) to calculate the error. Again, it is not a problem if you are only working with real numbers, but it will become a problem with complex numbers.

The way that you are getting the vector x inside the CMA loop is not getting all the samples since the beginning. So I iniciated the vector x before the loop as:

x = zeros(Le, 1);

and inside the loop I used:

x = [r(i); x(1:end-1)];

It will start CMA from the beginning of the receiving symbols.

This line doesn't make sense:

sb=w'*x;

This only equalizes ONE sample. The correct way to perform the equalization once you have the filter coefficients is:

sb = filter(w(:,N+1), 1, y);

Remember what I said about the function filter and the function conv. Read the documentation.

The next problem is in your number of errors computing. This is a very trick, but you need to remember (or learn) that every equalizer generates a delay in your signal. So you cant do ip- ipHat directly. You need to compensate this delay. I used:

nErr_CMA(ii) = size(find([ip(1:end-Le+1) - ipHat(Le:end)]),2);

If you don't know about this delay, you need to study a litle bit more.

In the plot part, I used the following:

figure
semilogy(Eb_N0_dB,theoryBer,'b'),grid;
hold on;
semilogy(Eb_N0_dB,simBer_CMA,'r-','Linewidth',2);
%axis([0 14 10^-5 0.5])
grid on
legend('Theoretical BER','sim-CMA');
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for BPSK in ISI with CMA equalizer');

I also increased N to N = 1e5, so as to perform a Monte Carlo simulation properly.
I changed Le to Le = 5;
I changed the w initializations for w = ones(Le, 1); to avoid that the update part of the CMA equation results in zero.
Also, change your channel. I didn't have time to verify that, but I strongly suspect that the channel [0.2 0.9 0.3] has some properties that makes it difficult for the CMA to equalize. I used ht = [1 0.5 0.2]. Try using this channel.

Ok, I gave you a lot of answers. I don't know if there is any more mistakes, but here I was able to get a BER curve that looks right. If it is right or not, it is up to you. I suggest you to look for articles about communications systems using CMA algorithm and compare your results.

Let me know if I need to let anything I said more clear. Good luck on your research.


Your start looks right, although I think your speed equation is wrong. To continue, the FT has no time dependance so you can take the integral out.

$$G(u, v) = \int_0^T \mathcal F\{f(x-x_0, y-y_0)\} dt $$

$$ = F(u,v) \int_0^T e^{-2\pi i(x_0u + y_0v)} dt$$

Divide by $F$ to get the transfer function.

If I understand the question correctly $x_0(t)$ and $y_0(t)$ are positions such that $x_0(t) = x_0(0) + v_x(t)t$ and similar for $y$. You can then put in an appropriate version of the speed given and do the integration.


I can answer the first part of your question: Where formula 2(3) comes from.

Since the noise model is Laplacian, you are looking at an estimate of the standard deviation of a Laplace distribution. Compare the moments equation on this Wikipedia page: http://en.wikipedia.org/wiki/Laplace_distribution#Moments

As for what N is, it looks like they define it as a matrix operator in 2(2). So the sum in 2(3) is taken over pixels of an image generated by convolving the original image I(x,y) with the Laplacian operator N (which is denoted by $I(x,y)*N$, where $*$ is the convolution operator).


Relation between the DTFT and CTFT in sampling- sample period isn't as the impulse train period

We are given two signals:

$$x_c(t)=\sum_{n=-\infty}^{\infty}x(nT_1)\delta(t-nT_2)\tag{1}$$

and

$$x_d[n]=x(nT_1)\tag{2}$$

where $x(t)$ is some known signal whose Fourier transform $X(j\Omega)$ exists, and $T_1>0$ and $T_2>0$, but generally $T_1\neq T_2$.

We want to express the DTFT $X_d(e^{j\omega})$ of $x_d[n]$ in terms of the CTFT $X_c(j\Omega)$ of $x_c(t)$. The DTFT $X_d(e^{j\omega})$ is an aliased version of $X(j\Omega)$:

$$X_d(e^{j\omega})=\frac{1}{T_1}\sum_{k=-\infty}^{\infty}X\left[j\left(\frac{\omega-2\pi k}{T_1}\right)\right]\tag{3}$$

In order to compute the CTFT $X_c(j\Omega)$ of $x_c(t)$ we define a compressed ($T_1>T_2$) or expanded ($T_1<T_2$) version of $x(t)$:

$$\tilde{x}(t)=x\left(\frac{T_1}{T_2}t\right)\tag{4}$$

The CTFT of $\tilde{x}(t)$ is

$$\tilde{X}(j\Omega)=\frac{T_2}{T_1}X\left(j\frac{T_2}{T_1}\Omega\right)\tag{5}$$

With $(4)$ we have $x(nT_1)=\tilde{x}(nT_2)$, and $x_c(t)$ can be written as

$$x_c(t)=\sum_{n=-\infty}^{\infty}\tilde{x}(nT_2)\delta(t-nT_2)\tag{6}$$

The CTFT of $(6)$ is

$$X_c(j\Omega)=\frac{1}{T_2}\sum_{k=-\infty}^{\infty}\tilde{X}\left[j\left(\Omega-\frac{2\pi k}{T_2}\right)\right]\tag{7}$$

and with $(5)$ this becomes

$$X_c(j\Omega)=\frac{1}{T_1}\sum_{k=-\infty}^{\infty}X\left[j\left(\frac{T_2}{T_1}\Omega-\frac{2\pi k}{T_1}\right)\right]\tag{8}$$

Comparing $(3)$ and $(8)$ we finally obtain

$$X_d(e^{j\omega})=X_c\left(j\frac{\omega}{T_2}\right)\tag{9}$$


How to find peak value of an analog signal efficiently after sampling in the digital domain?

You are unlikely to be able to evaluate correctly an ideal sinc interpolation filter in real-time (or even with a small time lag), since the future samples of your data are unknown. Plus, when you received your last sample $x[n]$, you generally don't know whether you are on  a peak.  So, it somehow boils down to

how real-time you want to be;
how fast you can estimated you are passed a peak.

There are many local possible interpolations. You can also think about  fast option, and close to real-time, only require to collect the next sample $x[n+1]$, and if knowing $x[n-1]$ you consider that a peak is located between  $x[n-1]$ and $x[n+1]$ (and there values are positive), you can use the standard frequency peak  estimators, which only use the above three samples. For instance:
$$n_{\textrm{peak}} = n+\delta$$
with
$$\delta = \frac{(x_{n+1}−x_{n-1})}{(4x_{n} − 2x_{n-1} − 2x_{n+1})}$$
A primer is provided by famous people here: Fast, Accurate Frequency Estimators, Eric Jacobsen and
Peter Kootsookos, IEEE signal processing magazine, May 2007. The above equation is only one of those this paper provides.


When negative frequencies cannot exist in practice, why do we calculate power in negative frequency?

Negative frequencies do exist. The vector base for Fourier transforms (and therefore power spectral density) is the set of complex exponentials $e^{j\omega t}$ where the frequency $\omega$ can be positive or negative.

A real sinusoid can be decomposed into complex exponentials:

$$\cos(\omega t) = \frac{e^{j\omega t} + e^{-j \omega t}}{2}$$

There you have your practical signal (the cosine wave) expressed in terms of the Fourier vector base (complex exponentials). Both sides of the spectrum must be considered in order to calculate power.

This is similar for any other real signal. In general, for real signals, the power spectrum will be symmetric and you must consider power in both sides of the spectrum.


How many moments are needed for a discrete distributions

You can do $E\{ x^n \}$ as much as you want. How much is necessary is up to you but the unique characterization of a distribution is not the same thing as a cutoff point of a moment expansion.  Some distributions can be very compact with respect to their moments but others are not.  My experience is that moment (taylor like) expansions can converge very slowly.

For something like a binomial distribution, all I need to know to represent the pdf is $p=$ success and know that the trials are independent.   The probability density for successes given a number of trials flollows from the mathematical formula.


There's simply no inherent relation between carrier frequency and data rate.

I assume a classical comms engineering education here, but since that most often introduces the concept of equivalent baseband early, a good explanation would be:

  This becomes obvious when you think about the complex baseband representation of a passband being the same, regardless of the passband's center frequency, and that the complex baseband carries identical information to the passband signal.

  Thus, carrier frequency cannot have direct influence on the information transport capabilities of any system, given that you mix the same signal to these carriers.

Indirectly, there might be influences: Bandwidth is proportional to capacity; it's impossible to get 1 GHz bandwidth centered around 433 MHz, but cheap around 60 GHz.

Noise and interference are the other limit to capacity; these are different for different bands and technologies.


What is a good FFT window function to reject DC?

If you're concerned with doing spectral analysis on a signal with a large DC component, and you want to suppress that DC peak, then a window function is not what you want. As some other answers noted, a highpass filter (or, viewed differently, a notch filter with the notch at zero frequency) is an appropriate solution.

To understand why, you need to think about what applying a window function does to the frequency response of each DFT output. The DFT is defined as:

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

One interpretation of how the DFT works is as a bank of filters at $N$ equally-spaced frequencies between $-\frac{f_s}{2}$ and $\frac{f_s}{2}$. Recast the sum above as follows:

$$
X[k] = \sum_{n=0}^{N-1} x_k[n]
$$

where:

$$
x_k[n] = x[n] e^{\frac{-j2 \pi n k}{N}}
$$

So, the $k$-th DFT output is generated by first taking the input signal $x[n]$ and multiplying it by a complex exponential at frequency $\frac{-2\pi k}{N}$ to yield a downconverted signal $x_k[n]$. The resulting signal is then summed over the $N$-sample window to yield the DFT output $X[k]$. This is effectively a moving average filter (sometimes called a boxcar filter), whose impulse response can be described as:

$$
b[n] =
\begin{cases}
1,\ x = 0, 1, \ldots , N-1 \\
0,\ \text{otherwise}
\end{cases}
$$

The magnitude response of the boxcar filter can be found by taking the discrete-time Fourier transform (DTFT) of that impulse response:

$$
|H(f)| = \left|\frac{\sin\left(N \pi\frac{f}{f_s}\right)}{\sin\left(\pi\frac{f}{f_s}\right)}\right|
$$

This is a Dirichlet kernel, and is sometimes referred to as a "periodic sinc" since it looks a bit like a sinc function but repeats periodically, which a sinc does not. This expression gives the magnitude response of each DFT output, where $f$ is measured as the frequency offset from the center frequency of the respective output bin. This illustrates the spectral leakage effect; each DFT output has a frequency response that covers some continuous swath of the input signal's spectrum, not just the discrete center frequency of each output.

Now consider how things change if you apply a window function to the input signal $x[n]$ before performing the DFT:

$$
\begin{align*}
X[k] &= \sum_{n=0}^{N-1} w[n] x[n] e^{\frac{-j2 \pi n k}{N}} \\
&= \sum_{n=0}^{N-1} w[n] x_k[n]
\end{align*}
$$

With the window function in place, the downconverted $x_k[n]$ is effectively passing through an FIR filter with an impulse response described by the window function. So, the per-output magnitude response of the DFT is:

$$
|H(f)| = |W(f)|
$$

where $W(f)$ the DTFT of the window function $w[n]$. Now note that if you chose a window function that had a zero at DC and used it to premultiply $x[n]$ before the DFT, you would actually cause the unintended effect of nulling out not only DC in the resulting spectrum, but the center frequencies of every one of the DFT outputs. This is probably not what you want.

So, if you truly just want to cancel the signal's DC component, removing it via some other type of pre-processing, not time-domain windowing, is the way to go. You could use a linear highpass filter with a very low cutoff frequency or subtract the estimated mean from the signal first, for example. Choosing between these methods should be based upon what other constraints your system has.


Question about convolution in time and frequency domain

Normally, the variables $t$ and $f$ are used for continuous time and frequency. But from your question I understand that you're talking about discrete time and frequency, and their relation via the Discrete Fourier Transform (DFT).

If you have discrete-time sequences $a[n]$ and $b[n]$ and their DFTs $A[k]$ and $B[k]$, then the DFT of the linear convolution $c[n]=a[n]\star b[n]$ is NOT equal to $A[k]\cdot B[k]$. The product of the DFTs corresponds to circular (or cyclic) convolution in the time domain.

The linear convolution of two $16$-point sequences has indeed $2\cdot 16-1=31$ points, whereas the circular convolution of two $16$-point sequences also has $16$ points.


Expectation and Auto-correlation of an Output from a Matched Filter.

First you should find $Y(t)$. In this case, it can be calculated as:

$$Y(t)=X(t)*h(t)=\int_{-\infty}^{\infty}h(t)X(t-\tau)\ \mathrm{d}\tau$$

But you have more information about $X(t)$ and $h(t)$:

$$Y(t)=\int_{-\infty}^{\infty}S(-t)[S(t-\tau)+N(t-\tau)]\ \mathrm{d}\tau$$

If we assume that the integral converges (which makes sense, because if not $Y(t)$ would not exist), then we can do the following:

$$\begin{align}
\mathbb{E}[Y(t)]&=
\mathbb{E}\left[\int_{-\infty}^{\infty}S(-t)[S(t-\tau)+N(t-\tau)]\ \mathrm{d}\tau\right]\\
&=\int_{-\infty}^{\infty}\mathbb{E}\left[S(-t)S(t-\tau)+S(-t)N(t-\tau)\right]\ \mathrm{d}\tau
\end{align}
$$

But $S(t)$ is deterministic, so:

$$\mathbb{E}[Y(t)]= \int_{-\infty}^{\infty}S(-t)S(t-\tau)+S(-t)\mathbb{E}[N(t-\tau)]\ \mathrm{d}\tau$$

The problem states that the mean of the noise process is zero, thus:

$$\mathbb{E}[Y(t)]= \int_{-\infty}^{\infty}S(-t)S(t-\tau)\ \mathrm{d}\tau$$

So, the expectation of $Y(t)$ is the convolution between $S(t)$ and its reversed version (also known as autocorrelation of $S(t)$):

$$\mathbb{E}[Y(t)]=S(t)*S(-t)$$

Because $S(t)$ is known, then $\mathbb{E}[Y(t)]$ is, too.

To find the autocorrelation, you have to do something similar to this procedure. I'll leave that to you.


Fourier transform of a periodic/aperiodic signal

Equidistant Dirac impulses in the spectrum imply a periodic time domain signal.

As pointed out in a comment, in continuous time, the signal $x(t)=\sin(\omega_0t)$ is always periodic, regardless of the value of $\omega_0$.

Your question about the spectrum being periodic is unclear to me.

A one-sided spectrum with equidistant Dirac impulses also implies a periodic time-domain signal, but that signal is necessarily complex-valued because its spectrum has no Hermitian symmetry, as required for real-valued signals.


Equation Derivation: Help with intermediate steps

This is not directly related to signal processing, but similar calculations are ubiquitous in DSP. So, here we go:

$$\begin{align}e^{-n/D}\frac{1-e^{(n+1)/D}}{1-e^{1/D}}&=\frac{e^{-n/D}-e^{1/D}}{1-e^{1/D}}\\&=\frac{e^{1/D}}{e^{1/D}}\cdot\frac{e^{-(n+1)/D}-1}{e^{-1/D}-1}\\&=\frac{1-e^{-(n+1)/D}}{1-e^{-1/D}}\end{align}$$


How can I identify Arabic letters in an image using Java?

You have to do following methods:

First step: Number plate extraction

You can use OpenCV for this purpose. It has Java api. Also it is ported to Android.

Second step : OCR

OCR is used to recognise the letters on notice board.

You can train your own OCR for arabic letters and digits using OpenCV.

Or you can use tesseract-OCR. I think it supports arabic also.

Following paper will be useful for you: Paper 1

Also, this site includes complete code for number plate extraction, but in C++. But it will be a good start.


Minimum Output Samples needed to crack a "Gold Code" Generator (Dual LFSR)

For a single polynomial generator, the minimum number of samples is 2N+1 where N is the order of the polynomial. If noise free this would result in 10 equations with 10 unknowns, sufficient to solve the linear equations involved. For the case of the Gold Code given above, the generator shown could also be done with a single polynomial which is determined by convolving the two polynomials shown, resuting in a 20th order polynomial (it just would no longer be irreducible and primitive, but would still generate each Gold Code by changing the initial state to 20 consecutive values that are in the sequence desired). Therefore the answer is 21 chips, and the Berleykamp-Massey algorithm can quickly decipher the generator polynomial that could recreate every Gold code.


Analytical Calculation of the Impulse Response of an Arbitrary Filter

Your expression is the Laplace transform of the exponential function $e^{-pt}$. The Laplace transform of the dirac unit is $1$. So the inverse Laplace transform is exactly $e^{-pt}$. Maybe the fact that this function is noncausal adds a little complexity here.


Clarifying matrix notation from an ICA-CMN paper

The matrices $\hat{\mathbf{M}}$ and $\tilde{\mathbf{M}}$ are constructed in such a way that the relation $\mathbf{M}\mathbf{x}=\mathbf{y}$ implies $\hat{\mathbf{M}}\hat{\mathbf{x}}=\hat{\mathbf{y}}$ and $\tilde{\mathbf{M}}\tilde{\mathbf{x}}=\tilde{\mathbf{y}}$.

Consequently, for constructing the matrix $\tilde{\mathbf{M}}$,  each element $m_{kl}$ of $\mathbf{M}$ must be replaced by a $2\times 2$ sub-matrix

$$\begin{bmatrix}m_{kl}&0\\0&m_{kl}^*\end{bmatrix}$$

In a similar way, the matrix $\hat{\mathbf{M}}$ is formed by replacing each element $m_{kl}$ of $\mathbf{M}$ by the $2\times 2$ sub-matrix

$$\begin{bmatrix}m_{kl}^R&-m_{kl}^I\\m_{kl}^I&m_{kl}^R\end{bmatrix}$$

where the superscripts $^R$ and $^I$ denote real and imaginary parts, respectively.


Controlling one signal with another one (Rocket engine effect - Subtractive synthesis of sound)

I have already solved this problem:
I have a white noise generated as an array of doubles:

public static double[] whiteNoise(int length) {
    double[] out = new double[length];
    for (int i = 0; i < length; i++) {
        out[i] = (Math.random() * 2 - 1.0)/2.0;
    }
    return out;
}

I obtain the red noise from the white noise in the following way:

public static double[] redNoise(int length, int m) {
    double out[] = whiteNoise(length);

    for (int i=0; i<length/m-1; i++) {
        int k = i*m;
        int l = ((i+1)*m <= length-1) ? (i+1)*m : length-1;
        double d = (out[l] - out[k])/((double)(m));
        for (int j=0; j<m; j++) {
            out[k+j] = out[k] + d*j;
        }
    }
    return out;
}

Where the m parameter is the number of interpolated samples (we take the first sample, and the m+1 sample, and interpolate values of all samples between these two).

And now the SHOW BEGINS

I generate the red noise with the interpolation number equals to let's say 10000. According to the value of consecutive samples, I alter this red noise with different values obtained from these samples taken as the interpolation number parameter. In this way a rocket sound effect is obtained.

int interpolatedframes = 10000;
    double[] rednoise1 = SigGen.redNoise(duration, interpolatedframes);
    while(i<rednoise1.length-1)
    {
        selectedNumber = (int)Math.floor(Math.abs((100*(rednoise1[i]+0.5)+100)));
        if(i+selectedNumber < rednoise1.length-1) {
            rednoise1 = SigGen.alterRedNoise(rednoise1, i, selectedNumber);
        }
        else
        {
            break;
        }
        i+=selectedNumber;
    }

With the following code the array with the noise is altered in a way, that it interpolates "interpolationNumber" of samples beginning from the startingIndex.

public static double[] alterRedNoise(double[] input, int startIndex, int interpolationNumber)
{
        int k = startIndex;
        int l = ((startIndex + interpolationNumber)<input.length-1?(startIndex+interpolationNumber):input.length-1);
        double d = (input[l] - input[k])/((double)(interpolationNumber));
        for (int j=0; j<interpolationNumber; j++) {
            input[k+j] = input[k] + d*j;
        }
    return input;
}

Sorry for a bit messy code, but I just wanted to publish it ASAP.

#EDIT 1

The sound which is supposed to be generated is the sound of a rocket being launched. Take a look at this movie (maybe from 1:45): https://www.youtube.com/watch?v=OnoNITE-CLc

To have the effect of the rocket flying over, an envelope generator might be used. It can be implemented as a simple multiplication of the samples obtained from the described method in such way, that it will be gradually suppressed.


Calculate aliasing of $x_a(t) = \cos{(2\pi300t)} + \cos(2\pi600t)$ when sampled with $F_s = 1000$

Peter said it all in his comment, but I'll try to explain it here in a bit more detail.

Your thinking is correct, but incomplete. Of course, the component at $f=300$ (unit depends on normalization of $t$) will not cause aliasing, so there will be no additional component below Nyquist, but what you forget is that sampling makes the spectrum periodic with period $f_s$ (sampling frequency). So the $f=\pm 300$ component results in components at $nf_s\pm 300$, $n\in\mathbb{Z}$. For $n=1$ you get the components at $f=1000\pm 300$, one of which (the one at $f=700$) is the one you see in the figure at normalized frequency $0.7$.


What happens when we Fast Forward the Video

What happens to the video:

In the most simple implementation, which is suitable for large integer fast forward rates, only 1 frame out of N, where N is the fast forward rate is played. For example, with a normal play rate at 30 FPS, frames 0 to 29 are played in sequence over a duration of one second. In fast forward with a 3x rate, frames 0, 3, 6 ... 87 are played in one second; and frames 1, 2, 4, 5... 88, 89 are discarded. This can be quite computationally expensive as the video decoder has to provide 3 times as many frames per second as in normal play (though the display continues receiving 30 frames per second).
When the hardware permits it and/or when the fast forward rate is small, the frame rate is simply increased. For example, in fast forward with a 1.1x rate, frames 0 .. 32 are played in one second.
There are applications when a video stream has to be sped up/down by a small factor (say 5%) while maintaining a target frame rate (for example when a broadcaster wants to speed up a show to make room for twenty extra seconds of ads). In this case, interpolation/resampling methods have to be used - this is the same kind of techniques as those used for telecine conversion.

What happens to the audio:

In some implementations, the audio is read at a higher sample rate, and decimated to the target sample rate. For example let us consider the example where the original audio is at 44.1kHz, and the fast forward rate is 3x. In 1s, you need to play 3s of the original audio. So you read 44100 x 3 samples from the original audio stream, low-pass filter it with a filter which has a 1/3 normalized cutoff frequency; and keep only every third sample. This gives you 44100 samples of sped-up audio, which are played in a time span of 1s at 44.1kHz. This has the unwanted effect of shifting up the pitch of all sounds, giving voices a "chipmunk effect". For small shift up/rates the effect is acceptable to the common public (though musically trained people might notice a well-known song has been shifted in pitch and is not in its original tonality).
In some better implementations (eg Apple Quicktime player), a time-stretching algorithm is used to multiply the speed of the input audio without altering the pitch. For playback / scrubbing, the sound quality rarely matters and the time-stretching algorithms are rather crude (you can have a look at the source code of soundtouch).


Literature on noise removal in OCR applications

Google research has some excellent papers, see for example:

An Overview of the Tesseract OCR Engine

Also, it seems that stackoverflow has a similar question:

https://stackoverflow.com/questions/4180629/ocr-and-image-preprocessing-techniques

The most powerful image filtering techniques that I'm aware are graph-cuts based,
which run in the following steps:

Calculate a (sparse) distance matrix between pixels (based on their intensity)
Spectral Clustering, keep only the lowest 3-8 eigenvectors
K-means clustering of the eigenvectors


There are two problems here that need to be solved. The first is the delay of the measured step response, i.e. how many samples from the beginning you need to discard. This is critical in your example because you chose a second-order system to model the data, and such a system is not very flexible in adding or subtracting delay. Delay could be added by adding more zeros to the system, i.e. by increasing the order of the numerator. If you don't want that then you will need to optimize the delay of the step response such that it can be optimally approximated by a second-order system.

The other problem is the fact that the designed filter doesn't reach the final value of the measured step response. The reason is that the impulse response required by Prony's method is derived by computing the difference of adjacent step response samples without taking into account that the step response was truncated. Since several samples of the measured step response are discarded, you also discard the same number of samples of the impulse response, and if you were to compute a step response from that truncated impulse response (by computing the cumulative sum), you would get a step response that is different from the measured one (and that would be close to the one implemented by the designed filter). The solution to this problem is very straightforward. Just truncate the step response and add a value of zero at the beginning. Then compute the corresponding impulse response by taking differences, and the first value of the impulse response will equal the first (non-zero) value of the step response, i.e., it will represent the accumulated past of the original impulse response (before truncation).

I did that and found as optimal starting value for the step response an index of $530$ for a filter order of $2$. The figure below shows that the step response of the designed filter closely approximates the measured step response:


Optimal receiver for a binary system employing waveform $s=A\cos^{2}(\pi f_c t \pm \alpha)$

Awesome, your transformation into a sum of a constant and a cosine is wonderful!
Since both signals are cosines shifted from the origin by the same $\frac A2$, that doesn't chang their euclidean distance.
Thus, subtracting $\frac A2$ doesn't change your decider's error rate. Also, it's technically super easy – just a high pass filter that removes the constant, no matter what it was.
The you're left with
\begin{align}
\tilde s_1 &= \frac A2 \cos(2\pi f_c t + 2\alpha)\\
\tilde s_2 &= \frac A2 \cos(2\pi f_c t - 2\alpha)& \beta := 2\alpha\\[1.5em]
\tilde s_1 &= \frac A2 \cos(2\pi f_c t + \beta)\\
\tilde s_2 &= \frac A2 \cos(2\pi f_c t - \beta),&0<\beta<\pi
\end{align}
Now, IQ-mixing this down with $f_c$ will indeed give us some phase-modulation-type thing, with the two constellation points with $\pm \beta$ as phase – but that's not a BPSK, since that is not always a phase difference of 180°!
However, it's always symmetrical to the I-axis in your constellation diagram. So, you can just remove the real part of your baseband signal and end up with something where the sign contains all the information.
Note that this, geometrically, is a Maximum Likelihood decision only if your noise distribution (which haven't told us) is is symmetrical after the transformation – which is the case for the usual complex Gaussian noise we assume, but I'd be careful whether that's the right assumption here - such non-symmetric PSK constellations are (surprisingly) actually used in a couple of cases, most of which involve interesting residual carrier systems with nonlinearities that shape the noise in a non-symmetric way.


What are the advantages of higher order Kalman Filters like EKF, UKF?

They should improve the convergence properties of the state estimate, i.e. if you have an initial value far from the trajectory, they could improve the region of attraction. They should also reduce the bias in the estimate.
In the EKF and UKF, there are no direct ways of including higher order moments of distributions, in both cases the state distributions are assumed to be Gaussian. Presumably, you can approximate non-Gaussian distributions with additional modeling. For (some) higher order approximations, it should be possible to include higher order moments. [However, as far as I understand, the additional computational complexity is usually not worth the performance gain (if any), since there are so many sources of uncertainty when modeling a system.] The go-to solution for determining the distributions of state-estimates is Monte-Carlo simulation.


The manual entry does not say that $\alpha = 1/\sigma$. It says $\alpha$ is proportional to the inverse of $\sigma$.

If you look at the example on that page you'll see

stdev = (N-1)/(2*alpha);

which, for your values says

stdev = (100-1)/(2*16.6667) = 825

you should be scaling the relationship by $N-1$.


Matlab produces two unknown spikes in custom FFT

The following is the corrected code. It seems you have the problem in the twiddle factor and the selection of even and odd samples of x[n]...

N=128;
f1=16;
num_cycles=32;
fs=f1*N/num_cycles;
x_time=0:1/fs:num_cycles/f1-1/fs;
x=sin(x_time*2*pi*f1);
plot(x_time,x);

% split inputs into even and odd samples and compute fft of each division
X_o=x(2:2:N);            % odd samples begin at x(2) --> x[1] in sequence
X_e=x(1:2:N);            % even samples begin at x(1)  --> x[0] in sequence

fft_x_o = fft(X_o, N/2);
fft_x_e = fft(X_e, N/2);

% Generate base twiddle factor
WN=exp(-1i*2*pi/N);

% Combine fft even and odd with twiddle factors to produce final output
for k=0:N-1
    if k<N/2
        X(k+1) = fft_x_e(k+1)+(WN^k)*fft_x_o(k+1);
    else
        X(k+1) = fft_x_e(k+1-N/2)+(W3N^k)*fft_x_o(k+1-N/2);
    end
end

% plot butterfly fft and matlab fft
FFT_xaxis=0:fs/N:fs-fs/N;
figure
plot(FFT_xaxis,abs(X))
title('Butterfly FFT')
xlabel('Frequency')
ylabel('Magnitude')
matlab_fft=fft(x);
figure
plot(FFT_xaxis,abs(matlab_fft))
title('Matlab FFT')
xlabel('Frequency')
ylabel('Magnitude')


Conditions for unique (not exact) recovery from noisy compressed measurements

"Recovery of Exact Sparse Representations in the Presence of Bounded Noise", by J. Fuchs, deals with the case you ask. From the abstract:

  The purpose of this contribution is to extend some recent results on sparse representations of signals in redundant bases developed in the noise-free case to the case of noisy observations. [..] We consider the case $b = Ax_0 + e$ [...] and seek conditions under which $x_0$ can be recovered from $b$ ...

I think theorem 2 and 3 are what you're looking for (at least in part).


What does the exponential term in the Fourier transform mean?

It's a complex exponential that rotates forever on the complex plane unit circle:
$$e^{-j\omega t} = \cos(\omega t) - j \sin(\omega t).$$
You can think of Fourier transform as calculating correlation between $f(t)$ and a complex exponential of each frequency, comparing how similar they are. Complex exponentials like that have the nice quality that they can be time-shifted by multiplying them with a complex number of unit magnitude (a constant complex exponential). If the Fourier transform result at a particular frequency is a non-real complex number, then the complex exponential of that frequency can be multiplied by that complex number to get it shifted in time so that the correlation to $f(t)$ is maximized.


How Could I extract the signal in the below Image

Normally you could do some valley detection. Scipy does this (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html) but this requires knowing a chunk of the signal in advance so you couldn't do it real time.

For real time maybe you could look at the derivative of the signal and notice when it goes from negative to positive, it would find valleys although it may not be as smooth as your black line. Maybe also check out this: https://stackoverflow.com/questions/22583391/peak-signal-detection-in-realtime-timeseries-data


So, first of all, we need to be clear that you mean FIR filtering not IIR.

By "high-order filtering", I presume you mean an FIR with a long length:  $L \gg 1$.  That is not unusual.  FIR filters have to be pretty long in order to get anything useful done.  The sharper the selectivity of the filter, the larger $L$ must be.

Most often, the best practice for using the "windowing method" of FIR design, is to use a Kaiser Window because it's optimized to give you, for a fixed $L$, the best stopband rejection or the sharpest transition-band region (what goes in between the passband and the stopband).  There is this $\beta$ parameter that let's you adjust that tradeoff.

So, if you're doing the windowing method (say, using MATLAB), first choose and FFT length, $N$, that is much much larger than your expected FIR length.  Like your FIR length might be $L=1024$, but then set $N=65536$ or something.  Then draw out your target frequency response and inverse FFT that to see what the impulse response might look like.  Remember the latter half is really negative frequency before the iFFT and negative time after the iFFT, so you may have to use fftshift() to swap the halves.  Then apply the window to that longer FIR to make it shorter, but for the purpose of analysis, keep the zeroed samples that became zero because of the windowing.  Then FFT back to see how bad the windowing affected your frequency response.  Adjust the Kaiser $\beta$ to see what your tradeoff looks like.

Rinse and repeat.  If you can't get the frequency response you want, you may have to increase the FIR length $L$.


In discrete time, distinct frequencies only have sense in a range of $2\pi$ (where the units are radians per sample).

A sinusoid (or complex exponential) of frequency $2\pi$ is indistinguishable from one with frequency 0. If in doubt, evaluate $x[n] = \cos(2\pi \times n)$ and $x[n] = \cos(0 \times n)$. Same thing with any two frequencies separated $2\pi$ or multiples of $2\pi$.

For example, you can easily show that $\cos( \theta_0 \times n ) = \cos((\theta_0 + 2k\pi) \times n)$ for any integer $k$.

So a range of size $2\pi$ is enough not to repeat frequencies.

In general, the range chosen is $-\pi \cdots \pi$, where $\pi$ is the maximum frequency ( $\cos(\pi \times n ) = (-1)^n$ ).

In the $z$ domain, the frequency response is the transfer function ($H(z)$) evaluated on the unit circle ($z = e^{j\theta}$), so here again you have the range $2\pi$ before repeating your path.


On generalizing Gabor-Heisenberg uncertainty in a noisy environment

This uncertainty thing is often misapplied.  Your statement for #2 is such a misapplication.  The bell curve, aka Gaussian, is an eigenfunction of the FT, (only an approximation in the DFT, but that is a different story).  The uncertainty principle is saying when a shape is taken back and forth from the FT and the narrowest cluster of information is desired, this is accomplished by the Gaussian.  It has nothing to do with parameter estimation, but noise indeed does.

You are conflating different concepts.

My latest blog article gives the actual code for calculating the frequency, amplitude, and phase of a pure tone exactly in the single noiseless pure tone case.  It is also robust in the presence of noise.

A Two Bin Solution

Early on, when trying to show "experts" that I had found exact formula frequencies for tones in the DFT, I was often quoted what you said and they would refuse to even look at my math as if I were proposing a perpetual motion machine.

So many misconceptions exist about the DFT, even among supposed learned folks, and that is scary.

Don't believe me?

Here is an actual quote from an expert who shall remain unnamed.

----------------------------------------------------------
An "exact" computation of frequency based on FTs is not possible for
some very simple reasons:  The Fourier Transform is itself an
estimator, and the time-frequency uncertainty principle means
frequency can't be determined exactly without an infinitely long
observation window.  So essentially all frequency computations of time
series or signals are estimates, it's just that some are more accurate
or more efficient than others.

So claiming that you have an exact solution just isn't going to get
the attention of everyone interested in the field.
----------------------------------------------------------

The only truth in there is the last sentence.  So, I wrote it up in a blog where the math is indisputable and anyone who wants to can verify it for themselves.

No, I've never gotten an apology.


Mean of input based on acf of filtered output

Yes. Mathematically, the limiting value of the autocorrelation function $R_X(\tau)$ of a wide-sense-stationary process is $$\lim_{\tau\to\infty}R_X(\tau) = \lim_{\tau\to-\infty}R_X(\tau) = \mu_X^2$$ where $\mu_X$ is the mean of the process. So it is not too surprising that res1 and res2 have identical shapes but res2 is shifted upwards in comparison to res1 and the above equation tells you how much the shift is.


How to relocate face points in opencv / face distortion?

Sure, that is one way of accomplishing your task so long as you do blending/interpolation to smooth the pixels in between. See these slides (esp. the last slide) or this survey for some literature on image warping and registration.

Thin-plate spline or "surface spline" problem
One reasonable approach would be establish point correspondences $r_1, \ldots, r_N$ and  $w_1, \ldots, w_N$ as you have done already then solve the following problem:

Find the 2D transformation model $\mathbf{p}=\{c_1,...,c_N,A,t\}$ s.t.
$$w_i=f(r_i)=Ar_i+t+\sum_{k=1}^{N}c_k\phi(d_i),$$
where $\phi(d_i)$ are radially symmetric basis functions (RBF). An example RBF choice is $\phi(d_i) = d_i^2\ln[d_i^2]$, where $d_i$ is a distance of your choice, e.g. Euclidean $d_i=\|r_i-r_k\|_2$.

Under certain thin-plate spline conditions not detailed here, this problem can be cast as a linear system $\mathbf{L}\mathbf{p}=\mathbf{w}$:
$$\begin{bmatrix}\mathbf{K} & \mathbf{P}\\\mathbf{P}^T & \mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{c}_{N\times2}\\\mathbf{A}_{3\times2}\end{bmatrix}=\begin{bmatrix}\mathbf{r}_{N\times2}\\\mathbf{0}_{3\times2}\end{bmatrix},$$
where $\mathbf{K}_{ij\,\,[N\times N]}=\phi\left(\left\|(r_i^x,r_j^y)-(r_j^x,r_j^y)\right\|_2\right)$ and $\mathbf{P}_{N\times3}=\begin{bmatrix}\mathbf{1}_{N\times1}\left|\mathbf{w}_{N\times2}\right.\end{bmatrix}$, and then solve for $\mathbf{p}$.

You can then apply the transformation model $\mathbf{p}$ to every pixel in the original/control image, and use 2D interpolation to get the final result.

And to answer your second question, assuming you are using OpenCV's CvPoint structure, you can just do the following:

cv::Point newPoint( oldPoint.x + xoffset, oldPoint.y + yoffset )


Dynamic range compression - Need pseudo algorithm for normalizing a signal

from your posted waveform, i am assuming that this is a unipolar signal.  that is

$$ x[n] \ge 0 $$

in audio, it would be the same, except that we would be working on $|x[n]|$ instead.

so first you want a sliding maximum of your signal, where the window length is $L$.

$$ x_1[n] = \max_{0 \le i < L} \Big| x[n-i] \Big|  $$

since your input signal $x[n]$ appears to be unipolar, you can leave off the absolute value operation.

then you want to apply a low-pass filter (LPF) to that sliding max.  you want the gain of the LPF at DC (0 Hz) to be 1 (or 0 dB gain).  a simple first-order LPF is

$$ \begin{align}
 x_2[n] &= (1-p) \cdot x_1[n] + p \cdot x_2[n-1] \\
\\
        &= x_1[n] + p \cdot (x_2[n-1] - x_1[n]) \\
\\
        &= x_2[n-1] + (1-p) \cdot (x_1[n] - x_2[n-1]) \\
\end{align} $$

$p$ is the pole value of the LPF and

$$ 0 < 1-p \ll 1 $$

so

$x_2[n]$ will be an envelope for your input signal $x[n]$ and will be delayed by about half of the max window length plus about 4 times the "time constant" of the LPF:

$$d = \frac{L}{2} - \frac{4}{\log(p)} \quad \quad \text{samples}$$.

so, to normalize, you want to invert (compute the reciprocal of) the envelope $x_2[n]$ and multiply your input signal by that inverted envelope.  the inverted envelope is

$$ x_3[n] = \frac{A}{x_2[n] + \epsilon} $$

$A$ is the normalized amplitude you want (it can be $A=1$ or $A=$ any other positive number that you like).  $\epsilon$ is a tiny number, much smaller than most of your non-zero $x[n]$ that you need to add to the denominator to keep from dividing by zero (which is a bad thing).

but you should line up the signal and the delayed inverted envelope, so your normalized output is also delayed by $d$ samples (defined above):

$$ y[n] = x_3[n] \cdot x[n-d] $$

the most computationally expensive operation is the sliding max.  we were just talking about this sliding max in the music-dsp mailing list.  i will dig up C code for it and post that in a following answer.


How to determine a suitable threshold when using amplitude thresholding method for spike detection?

For threshold-based spike detection, the threshold defines your performance in the detection vs false alarm performance. That is to say, a lower threshold will improve your probability of detecting (smaller) spikes but eventually, some of the noise may be considered a spike, increasing the probability of a false alarm, detecing a spike where there is none. A common term for showing this trade-off is ROC.
Therefore, it is highly application-dependent how to choose the threshold. If you have a specific false alarm rate that is tolerable, you can use the noise distribution to determine the threshold that gives you the desired false alarm rate (refered to as the constant false alarm rate (CFAR) approach). Same is true if you have a tolerable missed detection rate (àla I need to find at least 99.99% of the spikes).
If not, you should consider the threshold a variable and do your evaluations as a function of the threshold. This allows you to pick a threshold that optimizes your final figure of merit.


Can the inverse system of a stable system be unstable?

The ROC is determined by a condition on the real part of $s$, so for the given transfer function $H(s)$ you could define the ROC as $\Re\{s\}>1$, but then the system wouldn't be stable because the ROC does not include the $j\omega$ axis. You could as well have the ROC $-1<\Re\{s\}<1$. In this case the system would be stable but not causal.

And now for your main question: yes, the inverse of a stable system can be unstable. If we consider systems described by rational transfer functions (i.e. lumped elements) then the location of the system's poles determine its stability. The poles of the inverse system are the zeros of the current system, so there's no relation between the stability of a system and the stability of its inverse, because the zeros do not at all influence the system's stability, but they do determine the stability of the inverse system.

For your example, the inverse system has a pole at $s=2$, and there are two possible ROCs: $\Re\{s\}>2$, which corresponds to a causal but unstable system, and $\Re\{s\}<2$, corresponding to a stable but anti-causal system.


Selecting the Nyquist rate for a combination of signals

According to the lowpass sampling theorem, the required minimum sampling-rate (a.k.a Nyquist rate) is twice the highest frequency present in the frequeqncy spectrum of the bandlimited signal.
Note however that, the lowpass sampling theorem assumes there is no impulse (no pure sinusoidal component) at the highest frequency (a.k.a Nyquist frequency) of the signal being sampled. If there's an impulse, then you should choose a rate arbitrarily larger than the Nyquist rate.
Therefore, for your combined (linearly superposed) signal, the highest frequqency is 30 Hz; and thus the Nyquist rate is 60 Hz, however, since there's a pure sinusoidal at 30 Hz, then your sampling rate should be arbitrarily larger than 60 Hz.
In practice a 10% tolerance is already applied for various reasons. Hence you may use about 66 Hz for sampling your combined signal.


What your distortion box does is apply a non-linear transfer function to the signal: output = function(input) or y = f(x).  You're just applying the same function to every individual input sample to get the corresponding output sample.

When your input signal is a sine wave, a specific type of distortion is produced called harmonic distortion.  All of the new tones created by the distortion are perfect harmonics of the input signal:

If your transfer function has odd symmetry (can be rotated 180° about the origin), then it will produce only odd harmonics (1f, 3f, 5f, ...).  An example of a system with odd symmetry is a symmetrically-clipping amplifier.
If your transfer function has even symmetry (can be reflected across the Y axis), then the harmonics produced will only be even-order harmonics (0f, 2f, 4f, 6f, ...)  The fundamental 1f is an odd harmonic, and gets removed.  An example of a system with even symmetry is a full-wave rectifier.

So yes, if you want to add odd harmonics, put your signal through an odd-symmetric transfer function like y = tanh(x) or y = x^3.

If you want to add only even harmonics, put your signal through a transfer function that's even symmetric plus an identity function, to keep the original fundamental.  Something like y = x + x^4 or y = x + abs(x).  The x + keeps the fundamental that would otherwise be destroyed, while the x^4 is even-symmetric and produces only even harmonics (including DC, which you probably want to remove afterwards with a high-pass filter).

Even symmetry:

Transfer function with even symmetry:

Original signal in gray, with distorted signal in blue and spectrum of distorted signal showing only even harmonics and no fundamental:

Odd symmetry:

Transfer function with odd symmetry:

Original signal in gray, with distorted signal in blue and spectrum of distorted signal showing only odd harmonics, including fundamental:

Even symmetry + fundamental:

Transfer function with even symmetry plus identity function:

Original signal in gray, with distorted signal in blue and spectrum of distorted signal showing even harmonics plus fundamental:

This is what people are talking about when they say that a distortion box "adds odd harmonics", but it's not really accurate.  The problem is that harmonic distortion only exists for sine wave input.  Most people play instruments, not sine waves, so their input signal has multiple sine wave components.  In that case, you get intermodulation distortion, not harmonic distortion, and these rules about odd and even harmonics no longer apply.  For instance, applying a full-wave rectifier (even symmetry) to the following signals:

sine wave (fundamental odd harmonic only) → full-wave rectified sine (even harmonics only)
square wave (odd harmonics only) → DC (even 0th harmonic only)
sawtooth wave (odd and even harmonics) → triangle wave (odd harmonics only)
triangle wave (odd harmonics only) → 2× triangle wave (odd harmonics only)

So the output spectrum depends strongly on the input signal, not the distortion device, and whenever someone says "our amplifier/effect produces more-musical even-order harmonics", you should take it with a grain of salt.

(There is some truth to the claim that sounds with even harmonics are "more musical" than sounds with only odd harmonics, but these spectra aren't actually being produced here, as explained above, and this claim is only valid in the context of Western scales anyway.  Odd-harmonic sounds (square waves, clarinets, etc.) are more consonant on a Bohlen–Pierce musical scale based around the 3:1 ratio instead of the 2:1 octave.)

Another thing to remember is that digital non-linear processes can cause aliasing, which can be badly audible.  See Is there such a thing as band-limited non-linear distortion?


Seperation of wideband and narrowband - Adaptive Filter

A narrowband signal seems like (almost) periodic as indicated by

$$ x[n] = m[n] \sin( w_0 n) $$ where the message $m[n]$ has such a low bandwidth that the peak amplitude (the envelope) of the carrier sine wave changes very slowly compared to how fast the sine wave oscillates between those +/- envelope limits. This makes its autocorrelation sequence also to be long tailed and almost periodic. And it means that the narrowband signal is highly correlated for long lags (long delays).

The opposite is true for the wideband signal, its correlation sequence ceases very rapidly even for a small amount of lag.

The Wiener filter works by estimating the supplied desired response signal from its filter input signal, whose success depends on the amount of correlation between the desired signal and the filter input.

If there is large correlation between the desired signal and the filter input, then the filter will successfully estimate the desired signal and the error (which is desired signal minus the filter ouyput) signal will be small. On the other hand if there is not enough correlation between the two, then the estimation will be poor.

Based on this operational principle of the Wiener filter, and based on the fact that the narowband and the wideband signals have different correlation characteristics, from the given adaptive filter scheme, you can explain its operation as follows.

First, you see that the filter input signal and the desired response signals are the same (except the delay) and are a mixture of an auto-correlated narrowband component plus an uncorrelated wideband component.

The delay is essentially used to decorrelate the desired response and the filter input. However, by carefuly chosing the amount of delay (the lag), you can insure that the wideband component is fully decorraleted whilst the narrowband component is still corelated enough. If you chose the delay to too long, both components will lose their cross-correlation to the desired response. If you chose the delay too short, you would not be able to suppress the wideband component enough. Some intermediate delay will be sufficiently decorrelating the wideband component while still maintaining the cross-correlation between the narrowband components at the desired signal and the filter input.

When this is the case, the filter will successfully predict the correlated narrowband component in the desired response from the delayed version of it at its input, and it will therefore be subtracted from the supplied desired response signal which yields as the narrowband component as the error signal.


Object Recognition using bag of features using matlab

The, for now free available computer vision book, by Richard Szeliski contains an excellent introduction to this algorithm. Page 612.

EDIT:

I answered a similar question here on DSP with some example code.


I have no good explanation for why this happens, except that perhaps quadratic signals are nonstationary in that their mean is changing (or unbounded).

If I do something to ensure the mean is bounded (and not changing) such as:

a = sin(t.^2);
b = sin((t-15).^2);

or

a = rem(a,100);
b = rem(b,100);

then I get the right answer: 15 samples.

This is what happens when I use the first:

This is what happens when I use the second:

both modified signals generate the same (correct) estimate of delay.

Scilab code to generate it below.  I also modified your zero padding etc. because it seems overly complicated.

Code Only Below

// 25867
t = [1:1000];
a = (t.^2);
b = ((t-15).^2);
a = a - fix(a/200)*200;
b = b - fix(b/ 200)*200;

c = [a zeros(1,length(b) - 1)];
d = [b zeros(1,length(a) - 1)];

mu = mean([mean(c) mean(d)])

c_fft = fft(c - mu);
d_fft = fft(d - mu);

cross_spectral_density = c_fft.*conj(d_fft);
cross_correlations = fftshift(ifft(cross_spectral_density));
xc = xcorr(c- mu,d- mu,length(t)-1);

[max_xcorr shift_idx] = max(abs(cross_correlations));
shift = shift_idx - max(length(a),length(b));

clf;
subplot(211);
plot(a);
plot(b,'r');
subplot(212)
plot(cross_correlations);
plot(xc,'g:')
plot(shift_idx,cross_correlations(shift_idx),'r.')
xtitle('Delay is ' + string(shift))


This is not a standard 5-point derivative formula, the corresponding transfer function of which is

$$H(z)=\frac{1}{12}\left(-z^2+8z-8z^{-1}+z^{-2}\right)\tag{1}$$

The figure below shows the magnitude responses of an ideal differentiator (red), of the standard 5 point approximation given by (1) (green), and of the approximation in your question (blue):

As you can see, the non-standard approximation approximates the ideal differentiator only in a frequency region up to about $1/10$ of the Nyquist frequency (half the sampling frequency), whereas the standard approximation approximates the ideal differentiator quite well up to about a normalized frequency of $0.3$.

The differentiator in your question is a maximally flat lowpass differentiator. It avoids amplification of high-frequency noise. It achieves a maximally flat approximation of the ideal differentiator at $\omega=0$ and it has a given number of zero derivatives at $\omega=\pi$ (Nyquist), i.e. it is has a flat response at Nyquist, which causes high frequencies to be attenuated.

The design formula is derived in this paper. There is also a Matlab program available to design such maximally flat differentiators: link. Your specific filter is obtained by the command h=lowdiff(3,0);