Elegant way of calculating minimal sampling frequency without aliasing for annoying signals

Let me start by pointing out that for a band pass signal with a single band, there is no need for an iterative procedure. Taking your first figure as an example, in order to avoid aliasing, the sampling frequency $f_s$ must satisfy

$$\frac{2f_2}{n+1}\le f_s\le \frac{2f_1}{n}\tag{1}$$

for some integer $n$. The largest integer $n$ for which the range defined by $(1)$ is non-empty results in the smallest possible sampling frequency. This largest value of $n$ is given by

$$n_{max}=\left\lfloor\frac{f_1}{f_2-f_1}\right\rfloor\tag{2}$$

where $\lfloor\cdot\rfloor$ denotes the floor function. So you just need to determine $n_{max}$ from $(2)$ and then evaluate $(1)$ with $n=n_{max}$ to determine the possible range for $f_s$.

For multi-band signals, you usually need an iterative procedure. Let's take your second figure as an example (assuming that the spectrum is only shown for positive frequencies, and that it is symmetrical around $f=0$). Let's define $f_0=1$, $f_1=5$, and $f_2=6$. First consider the band pass component in the range $[f_1,f_2]$. From $(2)$ we get $n_{max}=5$. From $(1)$ we get the following "range" for $f_s$:

$$2\le f_s\le 2$$

Let's also evaluate the other ranges for values $n<n_{max}$:

$$\begin{align}n=4:&\quad 2.4\le f_s\le 2.5\\
n=3:&\quad 3.0\le f_s\le 3.3\\
n=2:&\quad 4.0\le f_s\le 5.0\\
n=1:&\quad 6.0\le f_s\le 10.0\end{align}\tag{3}$$

Let's now consider the low pass component. In any case, due to $f_0=1$, the minimum sampling frequency must satisfy $f_s\ge 2$. However, the aliases of the low pass component must not interfere with the band pass component. In a manner completely analogous to the single band band pass case, this leads to the following inequalities:

$$kf_s+f_0\le f_1,\quad (k+1)f_s-f_0\ge f_2$$

resulting in

$$\frac{f_0+f_2}{k+1}\le f_s\le \frac{f_1-f_0}{k}\tag{4}$$

The maximum value of $k$ for which the range defined by $(4)$ is non-empty is

$$k_{max}=\left\lfloor\frac{f_1-f_0}{2f_0+f_2-f_1}\right\rfloor\tag{5}$$

For the given values this gives $k_{max}=1$, resulting in the range

$$3.5\le f_s\le 4.0\tag{6}$$

Comparing $(3)$ and $(6)$, the only overlap is $f_s=4.0$, as you've already figured out.

Since this "overlap" is only a single frequency at the edge of two ranges, this solution is a purely theoretical one. Any imperfection in the implementation will cause aliasing. In practice, you would choose a larger sampling frequency. Considering $k=0$ in $(4)$ (meaning that the upper limit becomes infinite), we get $f_s\ge 7$. Comparing this with the ranges in $(3)$, we can conclude that a practical sampling frequency should satisfy $7+\delta<f_s<10-\delta$ with some reasonably chosen $\delta>0$.

In sum, for a signal composed of a low pass and a band pass signal (as in your example), the procedure is as follows:

determine the maximum values $n_{max}$ and $k_{max}$ from $(2)$ and $(5)$
determine the possible ranges for $f_s$ from $(1)$ and $(4)$ for decreasing values of $n$ and $k$, until you find overlapping ranges.

A general iterative procedure for finding the minimum sampling frequency which guarantees zero aliasing for arbitrary multi-band signals can be found in this paper.


Detect non rectangular bounding box for objects

Let's label your corners like this:

1-------2
|       |
3-------4

For your destination.  Let, $P_1, P_2, P_3, \text{ and } P_4 $ be vectors representing the corners in pixel coordinates.
For your source, Let $H$ and $W$ represent the height and width, and $(x,y)$ be a point in that rectangle you want to project.
You can calculate your coefficients like this:
$$
\begin{aligned}
C_1 &= \left(1 - \frac{x}{W} \right)\left( 1 - \frac{y}{H} \right) \\
C_2 &= \frac{x}{W} \left( 1 - \frac{y}{H} \right) \\
C_3 &= \left(1 - \frac{x}{W} \right)\frac{y}{H} \\
C_4 &= \frac{x}{W} \cdot \frac{y}{H}
\end{aligned}
$$
Note: $C_1 + C_2 + C_3 + C_4 = 1 $
Then the projection on the screen is at:
$$ P = C_1 \cdot P_1 + C_2 \cdot P_2 + C_3 \cdot P_3 + C_4 \cdot P_4 $$
This is a first order approximation that should be adequate for your needs.  Especially if the TV image is a relatively small portion of the screen.
To locate the corners, a rough approach is to scan horizontally looking for black streaks.  Match them row by row.  The endpoints will define your TV edges.  Use linear regression to fit a line.  Next scan vertically in the same manner and do the same.  Calculate the intersections as your corners.  They need not be integer, and will be more accurate if you leave them fractional.
This is the easy way.  Depending upon what you are projecting will likely work fine.  The more proper way is to find the inverse mapping, then for every pixel on the screen, find the location in your source rectangle and do a weighted average there using the same splitting formula there to find the weights of the corners.


I'm guessing the other frequencies it gets are harmonics of the fundamental?  Like you're playing 100 Hz and it picks out 200 Hz or 300 Hz instead?  First, you should limit your search space to the frequencies that a guitar is likely to be.  Find the highest fundamental you're likely to need and limit to that.

Autocorrelation will work better than FFT at finding the fundamental, if the fundamental is lower in amplitude than the harmonics (or missing altogether, but that's not an issue with guitar):

You can also try weighting the lower frequencies to emphasize the fundamental and minimize harmonics, or use a peak-picking algorithm like this and then just choose the lowest in frequency.

Also, you should be windowing your signal before applying the FFT.  You just multiply it by a window function, which tapers off the beginning and end of the waveform to make the frequency spectrum cleaner.  Then you get tall narrow spikes for frequency components instead of broad ones.

You can also use interpolation to get a more accurate peak.  Take the log of the spectrum, then fit a parabola to the peak and the two neighboring points, and find the parabola's true peak.  You might not need this much accuracy, though.

Here is my example Python code for all of this.


Power Line interference with different Sampling frequency

That depends if the bandpass filter was applied before or after sampling. If it was before sampling (i.e. in analog domain) it worked as an anti-aliasing filter and you dont need to filter out the PL. If it was applied in digital domain, you'd need to suppress it, because then the aliasing would occur. Have a look at the FFT of your signal to see if there is a peak or not at 14Hz, which could come from the PL.

(I guess the filter was applied in analog domain, because it does not make much sense to do it after sampling, since the high cutoff is already the Nyquist frequency.)


Suppose you have signals $x(t)$ and $y(t)$ whose cross-correlation function $R_{x,y}(t)$ is not something you like; you want $R_{x,y}$ to be impulse-like.  Note that in the frequency domain,
$$\mathcal{F}[R_{x,y}] = S_{x,y}(f) = X(f)Y^*(f).$$
So you filter the signals through linear filters $g$ and $h$ respectively to get
$\hat{x}(t) = x*g$, $\hat{X}(f) = X(f)G(f)$, and $\hat{y} = y*h$, $\hat{Y}(f) = Y(f)H(f)$, and now their cross-correlation function is
$R_{\hat{x},\hat{y}}$ whose Fourier transform is
$$\begin{align*}
\mathcal{F}[R_{\hat{x},\hat{y}}] = S_{\hat{x},\hat{y}}(f) &= [X(f)G(f)][Y(f)H(f)]^*\\
&= [X(f)Y^*(f)][G(f)H^*(f)]\\ &= [X(f)Y^*(f)][G^*(f)H(f)]^*,
\end{align*}$$
that is, $R_{\hat{x},\hat{y}}$ is the cross-correlation of $R_{x,y}$ with
$R_{h,g}$.  More importantly, you want to choose $g$ and $h$ so that
the cross-spectral density $G(f)H^*(f)$ of $g$ and $h$ is the multiplicative
inverse of the  cross-spectral density $X(f)Y^*(f)$ of $x$ and $y$, or
something close to it.  If you have only one signal and one filter, then
you get the result given by Hilmar (with amendment as given by my comment there).
In either case, the issue of compensating for spectral nulls, or generally,
frequency bands where the signals have little energy still remains.


CIE-Lab changes on Light intensity changes

I hope I get it right. You want to reduce effect of the light sources and perform the image segmentation. So basically you want to perform color based image segmentation.
If you do want to do that on a variety of data the first step should be performing the white-balancing technique on your images, and adjust the color temperature on the same values. This ensures the same type of your illumination source, even if the sources change over time (e.g. color temperature of daylight changes). The type of the light source definitely affects your color values. Unfortunately this can be performed only on RAW image datatypes, but if you expect change of light sources and you are keen on accurate results, I would strongly recommend the white-balance process.

For image color segmentation I do recommend mentioned CIELAB color space, because of it's independence. It does not depend on the input devices as well. The big advantage here is the possibility of separation the Lightness channel and perform segmentation only on the color ones (a*, b*). That is why this color space is more useful than for example RGB, which includes information of lightness in it's R, G, B components and the results are significantly biased due to that. However the L channel of the CIELAB corresponds more to the psychometric lightness rather than the real luminance.  The problem might be the intensity of light source, EV as well. If you get saturated and also really dark parts of the image you will not be able to restore the right color values. Other settings of the camera might have impact on the colour values as well. So if you want seriously accurate measurements and results from image color segmentation it might be really tricky. On the other hand, with a "normal" conditions and removing L channel in CIELAB color space, you can reduce the effects of the light sources on colour values as much as possible and there are not many other ways to do that. So ideally take your images with the same settings, white-balance them, transform them to CIELAB, remove the L channel and than perform color based segmentation.


Averaging Different Length Speech Signal

Dynamic Time Warping is pretty well explained on this site. I'll use some of the diagrams from the PPT on that site to explain.

The idea is to divide the signals into segments (frames) and then compare frames sequentially through each signal.  As illustrated below, motion from a segment in one signal to the next segment depends on the similarity to the current segment in the other signal.

For your example of audio files, you usually have features that you extract from them (LPCs or MFCCs for example).  Extract these from each frame.  Use these features as the comparison features used in the DTW algorithm.

  1) Suppose I have 10-10 frames of my 2 signals and each frame has 13 MFCC coefficients so the 2 DTW signals will be all coefficients arranged sequentially?

Well, before DTW you'll have $MFCC_s^f$ for signal $s$ for each frame, $f$:

Signal 1: $MFCC_1^1$,$MFCC_1^2$, $\ldots MFCC_1^{N_1}$

Signal 2: $MFCC_2^1,MFCC_2^2, \ldots MFCC_2^{N_2} $

where $N_1 \not=N_2$ and $N_s$ is the number of frames for each signal, $s = 1,2$.

After DTW, you'll have the best matching path over $N_\mbox{min} = \min(N_1,N_2)$ frames:

Signal 1: $MFCC_1^1$,$MFCC_1^1$, $\ldots MFCC_1^{N_?}$

Signal 2: $MFCC_2^1,MFCC_2^2, \ldots MFCC_2^{N_?} $

Then the comparison will be done by running all $those$ MFCCs sequentially.

  2) How this shortest path can be used to decrease or increase the size of one signal so that they are of same size?

See above explanation: After DTW you'll only have the same number of MFCCs for each and every signal.  All other (unused) frames with MFCCs are discarded for comparison purposes.


Your FFT code initializes a complex array because that is what the Fourier Transform produces. You should understand the complex exponential, and the nature of the transform as a change of basis (i.e. a projection onto a basis of sines & cosines).

It is a big topic that is well travelled, start here:
https://stackoverflow.com/questions/10304532/why-does-fft-produce-complex-numbers-instead-of-real-numbers


non-uniform sampling on the unit circle - applications

Uniform sampling is the definite standard as it simplifies not only mathematical derivations but also practical applications. On the other hand, there are some applications (such as array of sensors located in a nonuniform grid) in which taking uniformly distributed samples from a signal of interest is impossible. In such cases nonuniform sampling in either time or space is performed and the usual approach is to convert those nonuniform samples to uniform ones. Then follows the standard operations on those precomputed uniform samples.


Frequency components much less prominent at the receiver for multiple tones

The problem was from the fact that I was using a rectangular window for the ifft at the transmitter.

Using any other window at the ifft ensures there are no harmonics outside the 200kHz (which are interrupting the signal).

Here is an example of the signal using the rectangular window

And using a hamming window

Four tones

Received

Transmitted


Getting introduced to signal processing for remote sensing

@AndrewP. You have a long road ahead you. But the road is well lit and the rewards will are greater and greater as you travel along. You need to learn digital signal processing (DSP). May I humbly suggest you buy the two books "The Scientist & Engineer's Guide to Digital Signal Processing"

http://www.amazon.com/Scientist-Engineers-Digital-Signal-Processing/dp/0966017633/ref=sr_1_2?ie=UTF8&qid=1438652461&sr=8-2&keywords=scientist%27s+guide+signal+smith

and "Understanding Digital Signal Processing."

http://www.amazon.com/Understanding-Digital-Signal-Processing-3rd/dp/0137027419/ref=ntt_at_ep_dpi_1

They're widely considered to be the best ways to learn DSP for beginners. Do NOT sign up for any free Internet online DSP courses!  They will only discourage you at this point in your learning. You'll need to obtain either (1) a 'Student' version of Matlab software or (2) Octave software so you can experiment with DSP on your computer. Do not take any math classes at this point--you can learn selected mathematical topics as necessary as you travel your long road. Andrew if you have enough energy and enthusiasm, to quote Susan B. Anthony, "Failure is impossible."


The issue is resolved: the butter() function automatically does bilinear transform, you just need to specify the frequencies normalised to the Nyqyist freq.


It looks like using the polartocartesian transformation is destroying the characters (look at the "M", "E" before and after...). I think it will be a better idea to think about some rotations mechanism using imrotate.

Any way, I'd try the following additional pre-proc:
1. you might want to perform some morphological operations on the image after the adaptive threshold.
2. I would also try labeling the image (bwlabel) and removing too large objects (most of the background noise in your image is connected so it will be a huge object - in terms of pixel counting)

Have fun!


Why harmonic components appear only after a certain level when a signal is clipped?

Is this a well-known phenomenon?

Yes, of course. You will see harmonics as soon as your clip point is lower than the maximum amplitude in the time domain. The latter is a function of the relative phases between the harmonic components. In your case the max amplitude is indeed 2.5 (plus whatever the noise adds).

If you change the phases you will get a different clip point. For example, if you use cosine instead of sine, you will get a clip point of 3.0


How to Determine Covariance Matrix $Q$ and $R$ in Kalman Filter

R depends on the sensor sensitivity. If this is a real world problem this can be obtained from the manufacturer. If not use the identity matrix multiplied by a scalar that is less than 1.
Q is the covariance of the process noise. Again if this is a real world problem this can be obtained in the noise level in the states of the system at steady state. if not you can assume Q is zero matrix. A non zero Q helps achieve good convergency characteristics as explained here.


The three main reasons appear to be these:

1) In AFSK, it's not just jumping back and forth between frequencies. The tones must also be continuous phase. In other words, when the tone changes, there can't be any jump in phase. For example, if you're sending a 1200 hz tone and the waveform is at its peak when you switch to 2200 hz, the waveform is still at its peak, so it can't be just any two frequencies. It has to be at the peak. It can't start back at zero, or halfway through; it can't start anywhere but at the peak. Additionally, the length of each tone must be the same. Because of these two points, the amount of wavelengths that share compatible periods is finite. For instance, if you use 1000 hz, your next compatible is something like 1500 hz (my apologies if my numbers are off; I'm worn out).

2) As a follow-on from 1, since AFSK is audio, the range of compatible frequencies is not that high.

3) The amount of bandwidth required for the transmission increases for each additional tone, making data efficiency greater, but power and bandwidth efficiency less.

That doesn't meant that multiple frequencies can't be used; in fact MFSK has been used since the mid-60s. However, MFSK has several drawbacks that make it untenable in some situations. Delay spread, when broadcasting to multiple locations, and fading make it less desirable for some applications.


That text is confusing on its own as the spectrum spreading as transmitted is NOT affected as described. This is suggesting that if you acquire on any arbitrary 20 ms sequence of received data, your acquisition results will vary since a data transition can occur equally likely anywhere in that 20 ms given you select a random starting point for your acquisition. You avoid this whole issue by effectively performing acquisition using a 20 ms sliding correlation on 40 ms of received signal (assuming you want to maximize processing gain over the full 20 ms available).

The BPSK modulation includes the C/A Code specifically making the carrier change back and forth 180°. The C/A code is 1023 chips long; when a C/A code chip is = 1, the carrier prior to data modulation (such as L1 at 1575.42 MHz) is transmitted with 0° phase reference. When a C/A code chip = 0, the carrier is transmitted with a 180° phase reference. The C/A code is at a 1.023 MSps rate, so there are 1023 phase transitions in each 1ms C/A code sequence.

Then there are 20 C/A code sequences in one data bit:  Each data bit additionally toggles the carrier between 0 and 180°: So the C/A code is sent 20 times to send a "1" data bit, and the negative of the C/A code (0110... instead of 1001) is sent 20 times to send a "0" data bit. Thus the data rate is 50 bps.

Below is an example of a correlation to SV24 for a single 1023 length sequence associated with that satellite from actual data captured from a GPS antenna with a digital sampling scope. The correlation was done by sliding the sample-matched 1023 length sequence through the received signal that in every other way appeared as white noise. With frequency offsets removed, the result appears as in the plot below showing the correlation peaks every 1 ms at every point in time that the local C/A code aligned with the transmitted C/A code. What we see in the plot is three data symbols (either 1 1 0 or 0 0 1 since phase ambiguity hasn't been removed).

This is the same correlation result plotted on a complex plane (the above plot is the real part of this versus time). This shows how the carrier offset was removed and suggests carrier tracking approaches.


That's not true, it's not better.
The thing is: the matched filter just implements the projection in the signal vector space, onto the signal vector itself (or a multiple thereof). (You'll find correlation is just an inner product in that space.)
The line through that vector is the signal subspace, the plane to which that vector is normal is the noise space. Matched filtering is a method of linear algebra, and not different than using the SVD to find the subspaces. It's just that for matched filtering, the signal space is known.


Sharpening Using Laplacian Kernel - How to Handle the Output (Specifically Negative Values)?

@SOFUser, Welcome to our community.

Since the High Pass Filter, such as the Laplacian Kernel, output has zero mean (As the DC Component is removed) it indeed can (Should) have negative values as the output oscillates around zero.

What to do with it depends on the objective:

Display Purpose
In case on wants to display the output of the Sharpening / HPF Filter process it is very common to add 50% gray to it. In case of 8 [Bit] image it means adding 127 / 128. This is what's don in Adobe Photoshop.
Processing Purpose
Usually we keep the negative values as they might be important in the next steps. For instance in Laplacian Pyramid if one won't keep the negative values the reconstruction won't be accurate. Another example would be sharpening where we add the output of the Laplacian Filter scaled by the sharpening factor to the image to have sharpened image.
Generating Image
If one wants to generate image and archive it in one of the UINT based formats (JPEG, PNG, etc...) usually one would clip data into the compatible range.


What are the known use cases for different wavelet families?

Answer 0: ask yourself if you really need wavelets

Say yes. Let us concentrate on 2-band real discrete wavelet first. JPEG 2000 is a special case, where CDF9/7 and 5/3 biorthogonal wavelets are used. For the first one, the analysis wavelet is moderately regular, but its moments nicely concentrate piecewise regular parts of the data. The synthesis wavelet is smoother, and visually compensates quantization artifacts. The 5/3 works similarly, with the addition of exact computations on integers. There are almost symmetric, which is somehow better for edge location. The paper Mathematical Properties of the JPEG2000 Wavelet Filters gathers some of the most notable.

Answer 0.5. apart from that example, there is very few "fit-for-all signals" optimality results.

The quality of wavelets depends a lot on the nature of the data: signals, images, on the task, on the resource constraints. But the most important is the mastery of the wavelets you use. Because most often, you are trying to enhance or restore data, but the target is unknown, and quality metrics are often imperfect,  so you can only get informed guesses on simulated datasets. Theoretical properties are guides, but you need models for signals or noise at least (eg spline wavelets with a fractional Brownian noise).

Orthogonal wavelets, says Daubechies', are widely used because orthogonality is useful to prove stuff, and their "shortest" support interval for given moments is nice. On natural signals, a rule of thumb is that it is rare to observe polynomial parts above degree 3. Hence short DB should suffice. However, longer ones are often used to compensates for the stronger asymmetry of the shortest.

The (overall) symmetry of the wavelet is also observed: anti-symmetric to act like derivatives, symmetric to mimic Laplacians.

There is a no-go theorem of discrete 1D 2-band wavelets: except for Haar, they cannot be real, orthogonal with symmetry and finite support. Haar seems crude, but it is very fast, and since you can almost compute everything (delays, etc.) it can do a marvelous job.

But if you go deeper, you can find that $M$-band wavelets $M=4$ for instance, can have all the required properties, with more degrees of freedom. And they have better aliasing cancelation between channels, which makes them better at textures, for denoising.

Answer 1: a good knowledge of  wavelets properties restricts the candidate list, and a little trial-and-error is often unavoidable

But you can also use a set of different wavelets to enhance the result of a single one: perform the processing with each wavelet, and (wisely) combine the result. One of the most efficient technique is based on a pair of Hilbert-related wavelets: take any wavelet basis you like (for its properties), and build a second basis with an (approximate) Hilbert transform relationship. The resulting primal and dual wavelets share the same moment property,  have an opposite symmetry, and form together a dual-tree wavelet transform, that provides almost shift invariance at little cost. But you can also use different shifted versions (time-invariant wavelets), symmetrized combinations, etc.

Answer 2: a union of wavelet bases is often a powerful combination.


Why do butterworth filter sometimes make signals diverge?

Basically, this is bound to happen with higher order filters and finite precision. The filters' poles are getting closer and closer to the unit circle for sharper transition bands. But poles very close to the unit circle lead to unstable filters for some level of precision.

Things you could do:

Increase the precision of your numeric type (e.g. double instead of single precision)
Factor the filter into a series of 2nd order sections

Matlab and Octave allow you to directly compute the poles (p), zeros (z) and the overall gain (g) instead of the filter coefficients b and a:

[z,p,g] = butter(...);

and the result can be converted to a SOS matrix (second order sections):

[sos,g] = zp2sos(z,p,g);

Each row in the sos matrix represents a biquad, a 2nd order IIR filter. Using this second order section factorization -- meaning: a couple of 2nd order IIR filter in series -- is less subspectible to rounding errors as far as I know. Actual filtering can be done with sosfilt instead of filter:

y = sosfilt(sos,x)*g;

And if you still have trouble getting it stable, you should probably replace butter with ellip and allow a bit of passband ripple and/or do multiple passes of the same lower order filter.


The function gives you the path loss in dB. This number is positive, i.e. a path loss of $3$ dB corresponds to a gain of $-3$ dB. So for a given path loss $x$ (in dB), the signal is attenuated by a factor $10^{-x/20}$. Note that multiplication with a dB value hardly ever makes any sense.


Texture detection and region properties on an image

Here's a simple process:

Assign texture measurements to each region in the image.
Use a region labeling / connectivity algorithm (or region growing algorithm) to join adjacent regions having the same texture measure.
Implement a simple edge-following algorithm to trace the contour of each region.

Law's Texture Measures are an older but still useful technique for determining the texture in an image, and they may be sufficient for you to distinguish texture A from texture B throughout the image. See the section "Laws Texture Energy Measures" in the Wikipedia article:

http://en.wikipedia.org/wiki/Image_texture

As a first step, calculate all texture measures and determine which particular measure (e.g. Edge or Spot) allows you to distinguish one texture from the other most easily. (If you post some pictures I could help you identify a texture measure.)

If you have only two textures, A and B, then you can treat them as foreground and background, and a standard region labeling algorithm will work. To make it easier to see what's happening in the processing, you could generate a new image by assigning A texels (texture elements, little chunks of texture) to the color white, and B texels to the color black. Region labeling and/or contour-following algorithms would then find the connected white and black regions. The findContours() function in OpenCV will work well.

http://en.wikipedia.org/wiki/Connected-component_labeling

The same Wikipedia article includes both the traditional multipass algorithm as well as a single-pass algorithm. I haven't implemented the single-pass algorithm described there, but I have worked with the single-pass algorithm described in the paper "A Component-Labeling Algorithm Using Contour Tracing Technique" by Chen and Chang. The paper by Chen and Chang also describes a standard contour-following algorithm that can be implemented quickly.

If you have more than two textures, then you might use a watershed or mean shift algorithm to cluster regions together after you've remapped textures to colors. Although this remapping from texture to color isn't necessary, it does make the process easier to debug and understand.


How to blend images on a spherical surface / cube map

I have eventually found a simple idea for the solution on the paper called "Blending Images for Texturing 3D Models".

The idea is not doing the complete blending for each new image added to canvas, but rather build separate lowpass and highpass canvas (for two-band blending). These are then merged as a final step of the rendering.

Finding seam is still a bit tricky and not used in the paper. It requires reprojecting data from canvas to image, computing the seam and then projecting the blending mask to the canvas.

The method described in the paper is simple and does not leave any blending artifacts.

The only drawback is that more data have to be kept in memory.

Update

This solution does not work if the rendered canvas is smaller than original image. Furthermore gradient-domain blending is not easily applicable as well.


When and how can multiple EQ filters be merged into one?

This concept generalizes to any series cascade of linear filters. If I have a collection of linear filters with frequency responses $H_1(f), H_2(f), H_3(f), \ldots$ that are applied in series (i.e. an input signal is applied to $H_1$, its output is applied as the input to $H_2$, and so on), that series cascade could be replaced by a single filter whose overall frequency response is equal to $H(f) = H_1(f)H_2(f)H_3(f) \ldots$. This follows directly from the properties of LTI systems.

With that said, I'm not too familiar with some of the audio EQ terminology that you're using, so I'm not sure. If you have two filters with the same frequency response apart from a constant gain factor, like:

$$
H_1(f) = k_1 S(f)
$$

$$
H_2(f) = k_2 S(f)
$$

where $S(f)$ defines the shape of the frequency response, then the frequency response of the cascade of the two filters is:

$$
H(f) = H_1(f) H_2(f) = k_1 k_2 S^2(f)
$$

Note the squaring of $S(f)$; this is an important distinction. Applying two filters in series with the same frequency response shape (apart from a constant) will have a different effect than if you just created one similar filter with a different gain. Specifically, due to the squaring of $S(f)$, the passband will experience twice the ripple and the stopband will experience twice the attenuation than you would observe had you just used one of the two filters.


Why is the second half of the FFT negative frequencies

If you say that plugging in a negative value for $k$ would make sense to you then have a look at this:

$$\begin{align}X[N-1]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi (N-1)n/N}\\&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi n}e^{j2\pi n/N}\\&=\sum_{n=0}^{N-1}x[n]e^{j2\pi n/N}=X[-1]\end{align}$$

The DFT $X[k]$ is periodic with period $N$, so it doesn't matter if you consider the interval $k\in[0,N-1]$, or $k\in[-N/2+1,N/2]$ (assuming $N$ is even).


How to Generate the Following Image Effect (Modulation)?

The image seems to be modulated, namely multiplied by Harmonic Signal along its columns.

This is the MATLAB Code:

mInputImage = im2double(imread('2vE0G.png'));
numRows = size(mInputImage, 1);
numCols = size(mInputImage, 2);

signalFreq = floor(numRows / 1.25);

vModulationSignal = 1 + (0.5 * cos(2 * pi * (signalFreq / numRows) * [0:(numRows - 1)].'));

mOutputImage = bsxfun(@times, mInputImage, vModulationSignal);

figure();
imshow(mOutputImage);

This is the result:

You can play with the parameters.

Enjoy...


Does an asymmetric Bernoulli matrix satisfy the RIP?

As others have stated in the comments, the answer is "No".  The non-zero mean of the matrix dictates that a nonzero mean vector (say, all ones), will have substantially higher gain than a random vector with zero mean (say uniformly random +1,-1).

Consider the squared norm of A times a constant vector y is expected to be n*(p*N)^2.  (iteration of expectations)

The squared norm of A times a vector x drawn uniformly from (-1,+1) is expected to be n*(p*N).  (calculable by sum of variances of Binomial distribution)

The norms of x and y are the same, but the expectation of transformed norms differ by a factor of p*N -- diverging as the dimensions grow large.

Here's matlab code to help demonstrate.

n=2000;
N=1000;
p=.9;
A=double(rand(n,N)<p);
x=sign(randn(N,1));
y=ones(N,1);
Ex_normSqAx = n*(N*p);  % E[ squared norm of A times random signs ]
Ex_normSqAy = n*(N*p)^2; % E[ squared norm of A times constant vector ]
normSqAx = norm(A*x)^2;
normSqAy = norm(A*y)^2;


Yes, you could. When you're using spread spectrum for multiple access, what's important is the cross-correlation between any two users' spreading signals. Any cross-correlation between two users' spreading sequences will cause "leakage" from one user's channel into another's; this causes an artficial noise floor that will limit your system's performance.

You'll find that an up-chirp and down-chirp have very little correlation with one another, so they would be usable for this purpose.


Is there a rule of thumb for zero padding in image processing?

See Press et al. "Numerical Recipes in C++". Chapter 13 on "Fourier and Spectral Applications", Section 13.1 with a subsection entitled "Treatment of end effects by zero padding".  This is perhaps the best summary of zero-padding anywhere in the literature (but it's not referenced in either the contents or index).  My edition of this book is the 2002 second edition and the page number there is 545.


What is impulse response if i know input and output?

If $x[n]$ is non-zero for $n=0\ldots2$, and $y[n]$ is non-zero for $n=0\ldots3$, then $h[n]$ must have non-zero components for $n=0$ and $n=1$. This is a property of the convolution.

So now you have to solve $y[n] = h[0]x[n] + h[1]x[n-1]$ for the relevant values of $n$, and there you go.

Start with the edges ($n=0$ and $n=3$).


How are the constellation point coordinates determined in digital modulation

Constellation diagrams exist in what is called signal space which is an abstraction used to describe finite-energy signals. The coordinate axes, even if they are marked $x$ and $y$ as in Marcus Muller's answer, really represent unit-energy signals such as $$s_I(t) = \sqrt{\frac 2T}\cos(2\pi f_c t)\mathbf 1_{t \in [0,T)}(t),$$ $$s_Q(t) = -\sqrt{\frac 2T}\sin(2\pi f_c t)\mathbf 1_{t \in [0,T)}(t)$$ where $T = \frac{n}{f_c}$ and $\mathbf 1_{t \in [0,T)}(t)$ means a rectangular pulse of duration $T$. Note that both sinusoidal pulses consist of and integer number $n$ periods of a sinusoid of frequency $f_c$. Verify for yourself that $s_I(t)$ and $s_Q(t)$ are indeed unit energy signals and that they are orthogonal, that is, verify that
$$\int_0^T [s_I(t)]^2 \mathrm dt = \int_0^T [s_Q(t)]^2 \mathrm dt = 1;\quad
\int_0^T s_I(t)s_Q(t) \mathrm dt =0.$$
If we use $\pm As_I(t)$ as our PSK signals, then the signal energy is $A^2$(why?) and so if we use $E_b$ to denote the energy per bit, then
with respect to the signal space with axes $s_I(t)$ and $s_Q(t)$, the two signals used in PSK can be represented by the constellation points $(\sqrt{E_b},0)$ and $(-\sqrt{E_b},0)$. Left-handed folks might prefer to use PSK signals $\pm As_Q(t)$ in which case the signal constellation would have points $(0,\sqrt{E_b})$ and $(0,-\sqrt{E_b})$.
Mealy-mouthed ambidextrous people might want to use
$$\pm A[s_I(t)\cos(\theta)+s_Q(t)\sin(\theta)] = \pm A\sqrt{\frac 2T}\cos(2\pi f_c t+\theta)\mathbf 1_{t \in [0,T)}(t)$$
as PSK signals in which case in the signal space with axes $s_I(t)$ and $s_Q(t)$, the constellation points would be $\pm (\sqrt{E_b}\cos\theta,\sqrt{E_b}\sin\theta)$. But if we had chosen axes $\sqrt{\frac 2T}\cos(2\pi f_c t+\theta)$ and
$-\sqrt{\frac 2T}\sin(2\pi f_c t+\theta)$ instead (this just corresponds to a rotation of axes by $\theta$ from the previous signal space), then the constellation points would once again be $(\sqrt{E_b},0)$ and $(-\sqrt{E_b},0)$.

Thus, in signal space in which the coordinate axes represent unit-energy orthogonal signals, the energy of the signal corresponding to a constellation point equals the squared distance from the origin. Note that in (antipodal) PSK considered above, this squared distance from the origin is $E_b$ for each of the two constellation points. For QPSK (which I have ignored for the most part but you can read about it  in this answer of mine), each of the four constellation points are at squared distance $2E_b$ from the origin which makes sense in that $E_b$ is the energy per bit and a QPSK signal carries two bits.

Finally, with regard to the importance of distance as a concept, given two constellation points in signal space, the probability of error (in an AWGN channel with two-sided power spectral density $N_0/2$) is $Q\left(\frac{d}{\sqrt{2N_0}}\right)$ where $d$ is the distance between the two points. For the PSK examples considered above, $d = 2\sqrt{E_b}$ which gives the familiar result
$$P_e = Q\left(\frac{d}{\sqrt{2N_0}}\right) = Q\left(\frac{2\sqrt{E_b}}{\sqrt{2N_0}}\right) = Q\left(\sqrt{\frac{2E_b}{N_0}}\right)$$


CNN

Let's say you have $P$ phones $p_0,p_1,\ldots,p_{P-1}$ and the CNN generates the posterior probabilities $x_i = x(p_i)$ for $i=0,\ldots,P-1$.

HMM

A Hidden Markov Model is a system that has $S$ states $s_0,\ldots,s_{S-1}$ that can produce $O$ observations $o_0,\ldots,o_{O-1}$.

The probability of jumping from state $s_i$ to state $s_j$ is $a_{ij}$, and this is captured by the state transition (probability) matrix $\mathbf{A} = [a_{ij}]$.

The output $o_i$ is generated by the state $s_j$ with probability $b_{ij}$ to form the output matrix $\mathbf{B} = [b_{ij}]$.

How to connect them?

So let's say that your phones $p_i$ are the states of the HMM $s_i$ (so that $S=P$).

But after that, I can't really say. The $x_i$ aren't really the $a_{ij}$. Nor are they the $b_{ij}$.

This paper contains the image below... but their explanation isn't really clear to me.

I'll read the paper again in the morning and see if that sheds any light.


With every polyphase filter bank I have worked with, the first block in the analysis phase is an IFFT, and the block in the synthesis phase is a DFT. These operations essentially cancel one other, so it should be fairly intuitive. Think of the DFTs acting as sinusoids modulators in this case rather than an operation to convert to the frequency domain. By having complementary DFTs, the modulation is essentially performed in analysis and then removed in synthesis. From looking at the figure in the text you linked, it appears that they’re may be a misprint, at least in comparison to my industry experience implementing this types of filter banks.

The following is a great, compact and concise resource on polyphase filter banks, multi-rate identities, and DFT based filter banks:
https://apps.dtic.mil/dtic/tr/fulltext/u2/a457390.pdf


From what I found online, this paper is a pretty good intro to curvelets:
Curvelets A Surprisingly Effective Nonadaptive Representation For Objects with Edges.

Another decent one I found was Curvelets and Ridgelets.

A quick Google search can lead you to a lot of good resources, and provided sufficient background, you can learn most of these things on your own. I'm sure there's more literature freely available out there.


My guess is that the currently popular and new and likely robust way to solve this detection problem is to feed a sequence of audio fingerprints (such as MFCCs) to an RNN machine learning algorithm that was trained on a large wide range of rhythm tracks mixed with increasing levels of realistic background noise.

Feeding audio stream samples directly to a deeper CNN+RNN ML inference engine might work even better by detecting humanly resolvable beat details that audio fingerprinting misses, but that would be computationally far less efficient and would take a lot longer to train.


applying single pole filter to simple signal

If your pole is at $z=p$ (real or complex), this will do:

Y = filter( [1], [1, -p], S );

Free bonus: if you also want a zero at $z = c$:

Y = filter( [1, -c], [1, -p], S);

This last expression corresponds to the filter:

$$H(z) = \frac{1 - c z^{-1}}{1 - p z^{-1}}$$


Advice needed for separating 3-channel signals

I understand that you receive A, B, and C and you want to separate it into A1, A2; B1, B2; C1, C2.

From this, you can do the following: Model your expected signal A1 as some pulse shape (here, I assume its a tanh function, which is close to your example). Then, normally you can perform matched filtering on the received signal to find peaks at the position of your two pulses. However, since your pulses are non-finite (i.e. they remain constant non-zero), I decided to use their derivative as the signal of interest. Then, you can find something like this:

n = -10:0.01:40;

noise = 0.5;

ntest = -5:0.01:5;
testpulse = tanh(ntest);
testpulse_diff = diff(testpulse);

% Amplitude of first pulse
A1 = -4;
start1 = 0;
% Amplitude of second pulse
A2 = 11;
start2 = 3;
% The input signal (sum + some noise)
sumOfBoth = A1*tanh(n-start1) + A2*tanh(n-start2) + sqrt(noise)*randn(1, length(n));

subplot(2,2,1);
plot(n, sumOfBoth)
title('sum signal');

subplot(2,2,2);
plot(testpulse_diff);
title('Diff of required pulse');

subplot(2,2,3);
plot(diff(sumOfBoth));
title('1st derivative of input');

subplot(2,2,4);
metric = xcorr(diff(sumOfBoth), testpulse_diff);
plot(metric);
hold on;
[peaks,locs] = findpeaks(abs(metric),'MinPeakHeight',0.01,'MinPeakDistance',100);

plot(locs, metric(locs), 'or');
hold off;

title('Metric with the detected peaks that correspond to the two signals');

Which gives you the following output:

From the two detected peaks you can derive both the amplitude and the time position of the both peaks. Note that, even though the first derivative of your input signal looks like pure noise, it contains all the information that is needed. This information is extracted by means of the cross-correlation between the expected and the received peak.


Binary classifier for a small range of the audio spectrum

It could have been Goertzel's algorithm, though it looks at a single frequency rather than a specific band.

Another approach is to apply modulation techniques to shift the central frequency of your range of interest into the baseband, aka "zoom FFT".

Your intuition about the max to average ratio is good. Another "peakedness" metric is the ratio of geometric mean to arithmetic mean.


why type 2 FIR cannot be used to design high pass filter

Type 2 filters are even and symmetric.  A signal with a frequency of $f_s/2$ is just a scaled version of [+1,-1,+1,-1,+1,...].

That means, at $f_s/2$, the response of a type 2 filter will always be zero as the coefficients of the same magnitude are an even number of samples apart so they cancel each other out.

A filter with a magnitude frequency response of zero at $f_s/2$ makes a very poor high pass filter.


What part of complex number of inverse discrete Fourier transform?

To apply IFFT you need back the signal do complex numbers, you need use magnitude and phase information to rebuild correctly.

The real part is = magnitude * cos(phase)

The imaginary part is = magnitude * sin(phase)

You can use square roots of −1 (sqrt(-1)) to get Imaginary unit.

Now multiply imaginary unit with imaginary part and sum with real part, OK now are you done to apply IFFT !

At the end I apply a mat2gray function to convert the matrix to the intensity!

here how it is really done in matlab:

x=imread('C:\Users\Eder\Pictures\download.jpg');
figure(1);imshow(x);
%Make FFT
y=fft(x);
%Amplitude of the FFT
mx=abs(y);
%get Phase Information
ma=angle(y);
%back the signal to complex
y2= mx .* ( cos(ma) + sqrt(-1) * sin(ma) );
%Apply Inverse FFT
x2=real(ifft(y2));
result=mat2gray(x2);
figure(2);imshow(result);


Geometric representation of a signal using basis functions

Feels like a bit of a trick question. Your answer is indeed correct. $\Psi_1$ ...$\Psi_3$ are an orthonormal basis but it's an incomplete basis. Any signal you can construct with this basis has 4 degrees of freedom, but you only have three basis functions. In order to have a complete basis, you need 4 functions. Turns $x(t)$ (scaled properly) would make the 4th basis function, since it's orthogonal to the other three.


Simple fft to Gaussian pulse with MATLAB

You simply don't plot what you want to see. Note that your time domain signal is complex-valued and you modulate by a negative frequency. So the range of frequencies where things are happening are the negative frequencies (or - by periodicity - the frequencies in the range $[f_s/2,f_s]$), but those you don't plot. If you change the definition of your frequency vector and the corresponding plot command you'll see what you expect to see (i.e. a Gaussian in the frequency domain centered at $f_s-f=3.37e15$ Hz):

f = fs*linspace(0,1,NFFT); % (full range) Frequency Vector
...
plot(f,2*abs(X))

OR, simply change the negative modulation frequency to a positive one

x=(exp(-(t-50e-15).^2/(2*sigma)^2)).*exp(1i*2*pi*f*t);

and leave everything else unchanged. Then the Gaussian in the frequency domain is centered at $f=374.7$ THz. This is probably what you expected to see in the first place.


Valuable results can be obtained using Mel Frequency Cepstral Coefficients (MFCC) and Dynamic Time Warping (DTW).

Divide the signal in multiple segments. I recommend to use an overlap of 50%. Compute the MFCCs for each frame, then you will obtain a vector of MFCCs for each frame. Finally you will have a matrix of [No of MFFCs, Number of frames], for each signal. To compare the similarity between sequences with different time length, compute the Euclidean distance between the two matrices. Also you can use the standardized Euclidean distance (check scipy.spatial.distance.cdist for Python).

cost_matrix = cdist(x, y, metric='euclidean')

Thus, you will obtain a cost matrix. Applying the DTW algorithm on this cost matrix, you will find the score which will tell you the similarity between the input signals. I recommend you to see how a different number of MFCCs affects the results and also take in consideration the length of the signals.

I recommend this link to read more details about MFCCs and there is a useful video tutorial about DTW. Also you can find many articles about DTW, it is a well known measure for the similarity between two time series.


The transcription engine has its own background noise removal however it isn't trained on telephone audio.

There are many engines which are trained on bad quality audio and which perform just right.

It is pretty hopeless to repair the quality when it is already lost, it is better to find an accurate engine.


Weird results out of fourier transform of unknown signal, what am I doing wrong?

My guess: nothing to see here. Move along.

You mention that "the signal appears to be a helix in [the] frequency domain". I assert that what this helix really is, is a complex exponential. And there is nothing wrong with that, because complex exponentials are what you get when you time-delay a signal: $f(x-a)\Leftrightarrow \hat{f}(\xi)e^{-2\pi i a \xi}$. The helix pattern simply reflects that most of the power in the signal you are acquiring isn't located at $T=0$. (And it shouldn't be there anyway!)

Eyeballing, the period of oscillation in your frequency plot appears to be ~28 samples. At a sampling rate of 20hz, if what I'm saying is true, then most of the energy in your time-domain signal should be centered roughly around $t\approx\frac{28}{20}=1.4$ seconds. Was I close?

The real problem here is likely to be your plots. It seems like you're plotting raw real/imag DFT outputs. Don't do that. Instead, compute magnitude/phase from real/imag values, and plot that.


differences in fundamental frequency/pitch between music and speech

I don't think pitch information is relevant for what you want to do.

The variation of pitch during speech is known as intonation, and can convey emotions, indicate if a sentence is a question etc. However, there is no universal rule as to how pitch variation patterns are mapped to meaning - this is quite language dependent ; and some languages sound "flatter" in terms of pitch contour than others.

While there are some exceptions (say a solo flute performance), music is polyphonic in that several instruments are playing at the same time - and to make things worse some of those instruments might be playing several notes at the same time. As a result, at a given point in time, music is unlikely to contain just a single fundamental frequency ; and when applied to music, Praat's f0 estimator is very likely to yield a noisy and unreliable pitch contour, sometimes jumping from one instrument to the other.

Regarding automatically annotating movies from the audio, here are a few ideas:

Looking for sharp increases in loudness to detect important events as you suggested.
Use speech/music discrimination to help detect scene boundaries.
Train classifiers for common sound effects (explosions, gunshots, punches), and use this to detect action/violence scenes.
Use speaker diarization techniques (along with some face recognition) to detect all lines spoken by a particular actor.


How to reconstruct an image from s-parameters in Matlab?

Each complex value can be described as $$A \cdot e^{i\cdot \phi}  $$.
Where A is the magnitude and $\phi$ the angle.
With this knowlege you can use the ifft as you allready mentioned.
If you consider that the image is a real signal you can take the real part by using real().

So overall it the part after reading the magnitude and angle would look like

IMG = uint8(real(ifft2(mag .* exp(i * angle))));

Where IMG is the resulting image, mag the array of magnitudes (in linear scaling), angle the array of angles.

The only problem with that is, that you have to know the dimensions of the image. The spectrum of an image has (after FFT) the same dimensions as the picture (113 x 200 pixel picture -> 113 x 200 matrix after FFT).
So your mag and angle must have the dimension of the resulting picture (e.g. 113 x 200).


How do I show that the total system has the impulse response $h[n] = aδ[n] + bδ [n - 1] + cδ [n - 2]$ and determine the constants a,b and c?

Since this is a homework type problem I will not present the complete solution, but I'll help you to solve it yourself.

First of all, don't write out the convolution as an infinite sum, but simply use the fact that for any sequence $x[n]$ the following holds:

$$x[n]\star\delta[n-k]=x[n-k]\tag{1}$$

The lower path of the given system has the impulse response

$$h_2[n]\star(h_3[n]+h_4[n])=h_2[n]\star h_3[n]+h_2[n]\star h_4[n]=h_a[n]+h_b[n]\tag{2}$$

and if you use $(1)$ to evaluate $h_a[n]$ and $h_b[n]$ in $(2)$ you'll quickly see that both are (very simple) finite impulse response.


GNU Radio - PSK Mod block - unexpected constellation diagram

What you are seeing is the transitions from one constellation point to another.  In order to reduce the signal bandwidth, the baseband signal is low-pass filtered.  This causes the transitions to not be instantaneous (i.e. the I and Q are not square waves), so they take some time.  You are simply seeing those transitions.  The low-pass filtering also causes the overshoot that you see at the constellation points.

Just to see, I tried putting this into R. This is not going to get quite the same result because I am not sure of the parameters GNURadio uses, but it illustrates Jim's point.

As you increase the samples per symbol, the number of "transition" samples increases, making the edges of the constellation clearer.


Most efficient way to extract frequency components of audio signal with near-human accuracy in real time

First, with the classic short-term Fourier transform approach, there are alternative to interpolation - in particular techniques making use of phase information to recover the instantaneous frequency (See this question) which can give you very accurately the position of a spectral peak without an increase of FFT size. The drawback, as you correctly said, is that you are not increasing the ability of the system to discriminate adjacent peaks - but this is already a great improvement compared to using the central frequency of the FFT bin index.

Regarding your brute force approach... FFT has this resolution limitation (with bins at $\frac{sr}{FFT\_size}$) because it is "probing" the signal with complex exponentials that have a null integral over the analysis window - and these must have a period which is an integer division of the length of the analysis window. If you try to naively use the same approach to "probe" frequencies at other ratios, it won't work because you'll be "probing" your signal with functions that do not contain a complete cycle. If you wanted to make it work, you would have to increase your analysis window to the least common multiple of the periods of all the signals in your brute-force search, and this is precisely what you want to avoid!

There is another brute-force approach that works: "probe" your signals with windowed complex exponential (Gabor wavelets). These are characterized by a center frequency, a center time, and a bandwidth (which measures how the wavelet is spread over time or over frequency). You'll have to evaluate many, many, many correlations between your signal and these wavelets at as many time offsets, frequencies, and bandwidths you want. The result will be the same as a very flexible "tiled" STFT in which an optimal window size is selected for each time-range and each frequency-band. Besides the computational cost, the downside is that there is no efficient algorithm, and no causal algorithm (you will need to know as many samples in advance as the longest wavelet in your dictionary). If you want to experiment with these techniques, have a look at MPTK.

A last class of approaches working quite well for music signals are subspace tracking / high resolution methods. They work the problem from a model estimation angle - "assuming this signal is a sum of $k$ exponentially damped sinusoids, what are the best amplitudes/frequencies giving the least square error?", and as such the number of samples required to get an accurate estimate is only dependent on the noise level. There are some hurdles to make them work, though:

The number of components $k$ must ideally be known in advance. Order estimation methods like ESTER have been proposed to optimally pick the model order.
They perform well in the presence of white noise - this requires the signal to be whitened before analysis; performing the analysis in individual channels of a filter bank helps too.

These are computationally expensive, but they can work online, with short windows if model orders and/or noise are low.


Compressing a time-domain signal in the frequency domain Matlab

Everything depends on what you call compressing. What you are proposing is Fourier or spectrum thresholding. You will generate a $50$-length vector with some zeroes, hopefully. Let us ignore "complexness" that requires two floats for now,  you can deal with it using Hermitian symmetry.

Your main problem is actual compression. Forget about quantization first. The main issue is index location. Assume I have a transformed domain vector $T=[5, 1, 3 ,4, 0.5]$. Suppose I use a threshold on $2$. Now I give you the compressed version  $T_c=[5,  3 ,4]$, a ratio of $3/5$.

Can you retrieve $T_d=[5, 0, 3 ,4, 0]$? Probably not, you do not know where the zeroes were. This problem somehow prevented matching pursuit to be used for image compression, because coding the index penalized the final compressed file too much.


Running Kalman smoother is not same thing with running Kalman filter backwards in time. They are different. In particular, RTS smoother does not take care of observations, it uses filtered estimates. Kalman filter is an online estimation algorithm and estimates the filtering distribution $p(x_t | y_{1:t})$, i.e., the posterior distribution over the hidden states given the observations up to time $t$. However, smoothing recursions give the posterior of the states given the whole observations: $p(x_t | y_{1:T})$.

Even if your observations are reversible in time (if that is what you imply by saying 'symmetric'), in that case whenever you process the filtering algorithm, the smoother will work backwards and since filtered estimates are incomplete and does not know future observations, they will typically be different.


It is not clear from your question why you are even taking an FFT; knowing that may also lead to additional insights- but the main point to understand is frequency resolution in the FFT; assuming you are not doing any windowing (shaping the envelope of the time domain signal before taking the FFT), the frequency resolution is 1/T where T is the length of your FFT block in seconds. From this you will see that the frequency resolution is one FFT bin, and is equal to $f_s/N$ where $f_s$ is the sampling rate used and N is the number of samples in the FFT.

For additional info on this see What happens when N increases in N-point DFT and the additional referenced posts there.


Here are some graphics I have to add to the other good answers, as well as typical correction schemes for DC Offset and IQ Imbalance.
DC OFFSET
Below is an example DC Offset on a 16 QAM signal (with some added AWGN). In this example the offset if shifted by an amount such that the power in the fixed vector representing the shift is 10 dB below the total power in the signal. Notice that the shift can be at any magnitude and phase and the fact that it does not change with time makes it "DC". Also notice that even though it is 10 dB BELOW the total power in the signal, it appears significantly higher in the spectrum plot which may confuse some observers: this is because the waveform's power is spread over the modulation bandwidth, so the power in each FFT bin is lower than the power in the DC offset.

Below is a simple approach for eliminating DC offset. This filter is also introduced in this post: Transfer function of second order notch filter, where the gain constant $\alpha$ controls the bandwidth, getting tighter as $\alpha$ approaches 1.

IQ IMBALANCE
IQ Imbalance is demonstrated on the same 16 QAM signal in the next figure below, showing how it can result from either amplitude imbalance or phase imbalance or both. With amplitude imbalance, the scale of one axis is changed from the other, with phase imbalance, the axis are no longer in perfect quadrature.

The correction for IQ Imbalance is done as in the figure below with an $\alpha$ gain factor shown for amplitude imbalance and a $\beta$ gain factor shown for quadrature imbalance.


Signal strength is only important in determining signal quality insofar as it affects the signal-to-noise ratio (SNR).  As the name implies, SNR is determined by dividing the signal power by the noise power in the same bandwidth ($\frac SN$).

Claude Shannon, the pioneer of information theory, came up with a formula that gives an upper limit to how much information (you can think of this as bits) can be transmitted for a given SNR and bandwidth $B$ (I will note, though, that Wikipedia calls it the Shannon-Hartley theorem).  The formula is as follows-

$$
C = B\log_2\left(1 + \frac{S}{N}\right)
$$
where $C$ is the channel capacity in bits/s, $B$ is the bandwidth, and $\frac{S}{N}$ is the SNR.  Modern transmitters/receivers try to get as close to this limit as they can through a number of means.  One way is to increase the bits per symbol as the SNR improves (e.g. going from QPSK to 16-QAM).  Another is to reduce the error correction coding as the SNR improves, and increase it when it worsens.


General form of dilations and translations of wavelet function

Yes, they are equivalent, since $j\in\mathbb{Z}$. Note that your first definition becomes equal to the second when $j$ is substituded by $-j$. It is just a convention, if you consider a positive $j$ to be a up- or downscaling.

So, when writing a paper or reading a book, it is always important to understand the definition of the transform used.


How to find the lower eyelid in an image?

There is a section in this paper¹ that deals with eyelid contour approximation. It should help you out.

¹ Pernek, Ákos and Hajder, Levente (2014) Precise 3D Pose Estimation of Human Faces. In: VISAPP 2014. Proceedings of the 9th International Conference on Computer Vision Theory and Applications. Vol. 3. SciTePress, Lisszabon, pp. 618-625. ISBN 978-989-758-005-5


The science of finding values in between those points where you actually have samples is called Interpolation.

It's a rich field, and every signal processing toolkit (be it Matlab or OpenCV, or …) has functions for it.

Picking the right interpolator has a lot to do with the nature of your input, and the purpose of you doing the interpolation  – so, no general advice can be given on what specific interpolator to use.


Convolution with an impulse response does not give the same result as the frequency response

Explicitly making the imaginary part of your FFT coefficients 0 will change your magnitude response.

It is very clear if you consider this:
$$
X_k = a_k + jb_k\\
|X_k| = \sqrt{a_k^2 + b_k^2}
$$
Now, you're explicitly setting new coefficients to be $X'_k = a_k$, and hence $|X'_k| = |a_k|$.
Taking the IFFT of your new frequency response will obviously not be the same as the impulse response you expect.

If instead you made $X'_k = |X_k|$ and then performed an IFFT, you would get your expected filter taps. This is because your magnitude response is enforced to be the same and your phase is enforced to be zero.


Using Total Variation Denoising to Clean Accelerometer Data

Apart from Total Variation Denoising you could try a first much simpler approach: a median-filter. You just move a window along your data and replace the current input value by the median of all data in the window. You just have to optimize the window length (by experimenting).

By the way, the equations you copied into your question are for 2-dimensional data, but your data are 1-D.


How do I find the fundamental period of the given signal?

If your top equation is really
$$
x(t) = 2\cos\left(\frac 45 \pi t\right)\sin^2\left(\frac{16}{3} t\right)\tag{1}
$$
You gonna have a hard time getting the fundamental period/frequency as the there isn't an exact integer relating the two periods/frequencies. You can make approximations of the $\pi$ multiplier/divisor but the errors accumulates and this doesn't cut it. However if your equation is:
$$
x(t) = 2\cos\left(\frac 45 \pi t\right)\sin^2\left(\frac{16}{3} \color{magenta}\pi t\right)\tag{2}
$$
You could then easily extract the fundamental frequency $f_0$, and simply get the inverse $\displaystyle\frac{1}{f_0}$  to get the fundamental period $T_0$. For a signal $x(t)$

$$
x(t) = \sum_{k=1}^N a_k\cos\left(2\pi f_k t+ \phi_k\right)
$$

The fundamental frequency is the greatest common divisor (GCD) of all the frequencies; or
$$\mathrm{GCD}\left\{f_k\right\}\quad\big\vert \quad k = 1, \ldots, N$$

With non-integer values of $f_k$, you get the frequencies to the same common denominator or their least common multiple (LCM) and $f_0$ is the GCD of the numerators over this LCM, their common denominator. In your case, you have the form:
$$
x(t) = \prod_{k=1}^N a_k\cos\left(2\pi f_k t+ \phi_k\right)
$$
But the same principle applies. So, looking at the second-last equation you have two distinct frequencies: $\displaystyle f_1 = \frac 25$ and $\displaystyle f_2 =\frac{16}{3}$.

So in the end, you get $\displaystyle f_0 = \frac{2}{15}\implies T_0 = \frac 1f_0 = 7.5$.

Below is a plot illustrating this with the four waveforms below:

\begin{cases}
\color{cyan}{y_1} &= 2\cos\left(\frac 45 \pi t\right)\\
\color{green}{y_2} &= \sin^2\left(\frac{16}{3}\pi t\right)\\
\color{blue}{y_3} &= 2\cos\left(\frac 45\pi t\right)\cdot \sin^2\left(\frac{16}{3}\pi t\right)\\
\color{red}{y_4} &= 2\cos\left(2\pi\frac{2}{15}t\right)\\
\end{cases}

You see that the fundamental period of your signal $y_3$ ($\color{blue}{\rm blue}$) is at exactly $7.5 \ \rm seconds$ in the plot, matching exactly that of a pure sinusoidal wave $y_4$ ($\color{red}{\rm red}$) at $\displaystyle f_0 = \frac{2}{15}$.


In a mathematical sense it doesn't but (almost) any nonlinear interaction, say modulation or demodulation, will cause additional signals, leakage/out of band, signals.
The point is that when you do the evaluation you must include the modulator or demodulator; as well as things like transformers which will be biased off center and thus generate extra distortion.  It's the non-linearities in the overall signall path that induce extra signals.  By definition of the Fourier/Laplace transform superposition and linearity properties a linear system will not generate extra signals.


Trouble with time-syncing two radio signals using cross-correlation

It will work when you take the 2nd gradient of the signals:

import numpy as np
from scipy import signal

s0 = np.gradient(np.gradient(s0))
s1 = np.gradient(np.gradient(s1))
np.argmax(signal.correlate(s0, s1)) -> 525358

That corresponds to a shift of 1071 which is close to your expected 1069
Interestingly the minimum (most negative correlation) is close by at 1069

Plotting the cross-correlation of 2nd gradient gives a reasonably sharp peak:

Zomming in to the accordingly shifted data shows that this worked:

Using the most negative correlation results in the almost same picture:

What might be surprising is that the min and max are just two samples apart and also that the signal seems to be significantly phase shifted. On the other hand the start is matched perfectly. My interpretation of this is that the sync is not so much done by what we see as signal but by some high frequency parts of noise received in both channels. Since this 2 sample shift between min and max is the typical result when matching two channels that received the same random noise.

I can not explain easily why the correlation of the original data works in some cases. But here is an simple example that gives an idea:

a = np.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0], dtype=float)
b = np.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0], dtype=float)
np.argmax(signal.correlate(a, b)) -> 16

which is correct

but when adding an offset it fails

a += 5
b += 5
np.argmax(signal.correlate(a, b)) -> 20

the correlation of the first gradient still gives the right result but it will fail also if a significant ramp is added. Though for the 2nd gradient the correct result is achieved

a += np.array(range(0,len(a)))*5
b += np.array(range(0,len(b)))*5
np.argmax(signal.correlate(a, b)) -> 20
a = np.gradient(a)
b = np.gradient(b)
np.argmax(signal.correlate(a, b)) -> 20
a = np.gradient(a)
b = np.gradient(b)
np.argmax(signal.correlate(a, b)) -> 16


Proof of First Difference Property for Fourier Series

You have the result already written down in your question. If $a_k$ are the Fourier coefficients of $x[n]$, then your third formula is the Fourier series of $y[n]=x[n]-x[n-1]$. So the Fourier series coefficients of $y[n]$ are

$$b_k=a_k-a_ke^{-jk\omega_0}=a_k(1-e^{-jk\omega_0})$$

meaning that

$$y[n]=\sum_{k=0}^{N-1}b_ke^{j\omega_0 nk},\qquad\omega_0=2\pi/N$$


Finding audio signal time delay: cross-correlation vs. impulse response?

Your description of the cross-correlation method is not entirely correct. What you describe is actually a convolution. For a cross-correlation you would perform an FFT on both signals and then multiply one with the complex conjugate of the other, and finally perform an IFFT.


How can I calculate the total power? It is "the variance plus the mean squared", so I should just calculate $(\sum_{i}^N x[i]) / N$ for the signal of length $N$?

$(\sum_{i}^N x[i]\color{red}{^2}) / N$ is the average power, not the total energy. Just the sum $\sum_{i}^N x[i]\color{red}{^2}$ is what you're looking for.

Or just sum the PSD (like in this link)?

and

Should I perform this calculation in the time domain (just time series like above), or should I calculate FFT and PSD (Power Spectral Density) first and calculate on PSD (which is also just an array of numbers, if I understand correctly)?

The sum of the power in time domain (as above) and in frequency domain are the same, so that's just another way to compute the same (probably just more labor-intense).

I have seen quite a few formulas for the total power. For finite number of samples and discrete measurements I have seen (∑Nix[i]) divided by N, N+1, 2N or 2N+1. Which one should I use? Does it have something to do with bias/unbiased sampling, i.e. in statistics we add +1 to the divisor to get unbiased calculations?

Since you're really not looking for the power but the energy, neither. You don't have to divide at all.


Image distortion and alignment on 3D mesh / pointcloud

If you're only left with matching 3D points to their corresponding 2D parts, this is a fairly task and called PnP pose estimation. OpenCV has a lot of options to solve this problem. Here is one tutorial, which you could maybe benefit.


How to decompose a mixed signal into its original signals?

An FFT will only decompose a mix of sinusoids to the original constituents if all the constituent sinusoids are exactly integer periodic in the FFT's length.  One of your signals is not integer periodic in the number of samples in your FFT's length, so (1) its "leakage" (really windowing artifacts) will contaminate the magnitudes of other signals, and (2) you will need to interpolate (with a Sinc interpolation kernel) between FFT result bins to get a more accurate estimate of its magnitude.


What is the Asymptotic Equipartition Property (AEP)?

When you are sampling a stochastic process $n$ times, the larger you make $n$, the higher the probability that the series of samples is contained in the so called strongly typical set of outcomes of length $n$, where all members have roughly the same probability to be realized. If you let $n\rightarrow \infty$ it holds that
$$-\frac{1}{n}\log p(X_1,X_2,...,X_n)\rightarrow H(X)$$
with $H(X)$ being the entropy rate of the process.

This should be intuitively clear: the longer the sequence, the more possible outcomes, each with possiblity near zero but positive. The number of typical outcomes with $p\approx\epsilon$ is growing faster than the number of non typical outcomes as $n$ is growing. As $n\rightarrow \infty$, the "leftover probability" for non typical outcomes gets smaller and smaller, thus making the above equation valid.


How i can a 1 Hz resolution fft, from a limited data (between 0 to 100 msec)?

Your DFT's resolution is always $\frac{f_\text{sample}}{L_\text{FFT}}$, by definition of what the DFT is.

So, for any given frequency resolution, you'd need to transform the inverse of that frequency worth in samples. So, for 1Hz resolution, you need 1/(1Hz) = 1s worth of samples; you're giving it less than 100ms, i.e. 1/10 of that – you can zero-pad, but for an FFT-based frequency estimate, that few samples simply won't be enough to give you a low-variance estimate.

So, garbage in, garbage out.


impulse responses, $h(t)$, of continuous-time systems (those with output that is the same species of animal of the input) have dimensions 1/time.

consider a wire which has impulse response

$$ h(t) = \delta(t) $$

the area of the $\delta(t)$ is equal to the dimensionless $1$.  the width of the nascent delta is in dimension of time, so the height must be in dimension of 1/time.

the impulse response for a simple RC low-pass filter is

$$ h(t) = \tfrac{1}{RC} e^{-^{\, t}/_{RC}} \ u(t) $$

the dimension of $h(t)$ is in 1/time because $RC$ is in units of time.


How can the order of a transfer function be derived from its equivalent state space representation?

One trick will be whether there are any pole-zero cancellations in the overall transfer function.  For that, you will probably need to check for the minimality of the system (i.e. check whether the system is both controllable and observable).

The question of "determine the order of the transfer function for each input" seems odd, to me.  If the system is minimal (it can't be represented by a smaller order system), then you need it as-is.  Pulling out one of the $m\times 1$ systems on its own may be of interest, but you can't reduce it without destroying the overall system.

Can you give more information about why you want to examine just one of the $m$ transfer functions?


Rectangular 8 QAM bit error rate calculation

Neglecting square terms of the Q function, your expression is correct, but it is the symbol error probability, not the bit error probability:

$$P_s=2.5\,Q\left(\frac{d}{2\sigma}\right)\tag{1}$$

With Gray encoding, the bit error probability is closely approximated by

$$P_b\approx\frac13Ps\tag{2}$$

because there are $3$ bits per symbol, and the most likely symbol errors (i.e. between neighboring symbols) cause only one bit error.

Furthermore, if $N_0$ is the one-sided noise density, then $\sigma=\sqrt{N_0/2}$, and the energy per bit is $E_b=E_s/3$, where $E_s$ is the symbol energy. It remains to express the distance $d$ between the symbols in terms of $E_s$. The energy of the four inner symbols is $d^2/2$. The energy of the corner symbols is $5d^2/2$. So for equally probable symbols, the average symbol energy is

$$E_s=\frac12\cdot\frac{d^2}{2}+\frac12\cdot\frac{5d^2}{2}=\frac{3d^2}{2}\tag{3}$$

which gives

$$d=\sqrt{\frac{2E_s}{3}}=\sqrt{2E_b}\tag{4}$$

Plugging (4) into (1) gives for the symbol error probablity

$$P_s=2.5\,Q\left(\sqrt{\frac{E_b}{N_0}}\right)\tag{5}$$

which, combined with (2), gives for the average bit error probability

$$E_b\approx\frac{5}{6}\,Q\left(\sqrt{\frac{E_b}{N_0}}\right)\tag{6}$$

Note that Eq. (6) for the bit error probability was derived under the assumption that the SNR is relatively high. This assumption is implicit in Eq. (1), where quadratic terms of the Q function were neglected. The exact expression for the symbol error probability is

$$P_s=2.5\,Q\left(\frac{d}{2\sigma}\right)-1.5\,Q^2\left(\frac{d}{2\sigma}\right)\tag{a}$$

For $d/2\sigma\rightarrow 0$ the expression (a) approaches $7/8$. This is a consequence of all symbols being detected with equal probability, and since only one of them is the correct symbol, the symbol error probability must be $7/8$.

The expression (2) for the bit error probability assumes that each symbol error causes only one bit error, which is true if Gray encoding is used and if it is assumed that symbol errors only occur between neighboring symbols, the latter being a reasonable assumption at relatively high SNRs. However, for $d/2\sigma\rightarrow 0$ any symbol will be confused with any other with equal probability. There is only one symbol (the correct one) where there are zero bit errors. There are $3$ symbols causing one bit error, $3$ symbols causing $2$ bit errors, and one symbol causing $3$ bit errors. Since for very bad SNRs all symbols are detected with equal probability, the average number of wrong bits per detected symbol is

$$\frac18\cdot 0+\frac{3}{8}\cdot 1+\frac{3}{8}\cdot 2+\frac18\cdot 3=\frac{3}{2}$$

And since there are $3$ bits per symbol, the average bit error rate for $d/2\sigma\rightarrow 0$ is

$$E_b=\frac{3}{2}\cdot\frac13=\frac12$$

as expected.


Hey I've been meaning to write this program, too.  :D

The way you use the initial conditions is to just pass them from one stage to the next.

Also lfilter returns (y, zf), not (zf, y).

So it should look more like:

zi = signal.lfilter_zi(b, a) * xn[0]
z, zi2 = signal.lfilter(b, a, xn, zi=zi)
z2, zi3 = signal.lfilter(b, a, z, zi=zi2)
z3, zi4 = signal.lfilter(b, a, z2, zi=zi3)

Though you don't need to rename each one, of course.  Initialize it with

zi = signal.lfilter_zi(b, a) * xn[0]

and then in each loop do something like

for x in chunks(signal):
    y, zi = signal.lfilter(b, a, x, zi=zi)
    look_for_glitches_in(y)


Unit step function present in convolution result

A step by step derivation of this convolution would start with the following:

$$
\begin{align}
f(t)
  &= \int_{-\infty}^\infty x(\tau) h(t-\tau) \,d\tau \\
  &= \int_{-\infty}^\infty e^{-\alpha\tau} u(\tau) e^{-\beta(t-\tau)} u(t-\tau) \,d\tau
\end{align}
$$
And this is probably where you probably went on and tried to simplified the bounds on the integral with:
$$
\begin{align}
  f(t) &= \int_0^t e^{-\alpha\tau} e^{-\beta(t-\tau)} \, d\tau \\
\end{align}
$$
forgetting that this is only true if $t \geq 0$.
Indeed for $t < 0$ the product $u(\tau) u(t-\tau)$ is 0 for all values of $\tau$, and the result of the integral is thus 0.

So, the correct expression would be:
$$
\begin{align}
f(t)
  &= u(t) \int_0^t e^{-\alpha\tau} e^{-\beta(t-\tau)} \, d\tau \\
  &= \cdots \\
  &= u(t) \frac{e^{-\alpha t} - e^{-\beta t}}{\beta-\alpha}
\end{align}
$$


What is the spectrogram of an impulse response used for?

I think you are thinking of waterfall plots and/or cumulative spectral decay plots.  They are spectograms in the widest sense of the word since their windowing function tends to be exponential.  They are popular for plotting the response of speaker systems as they show group delay for various frequencies a lot more, well, graphically than the impulse response itself or a magnitude/phase plot of the transfer function would.  Resonators can be used for straightening out an otherwise wobbly (for whatever unrelated reason) frequency response magnitude at the cost of introducing delays.

The resulting lack of "impulse accuracy" which leads to problems with sound clarity and spatial perception is not easily apparent with other representations.


2D Convolution as a Doubly Block Circulant Matrix Operating on a Vector

The point is that circular convolution of two 1-D discrete signals can be expressed as the product of a circulant matrix and the vector representation of the other signal.

The circulant matrix is a toeplitz matrix which is constructed by different circular shifts of a vector in different rows. For example, consider two signls $h[n]$ and $g[n]$, each of length $4$. If we assume $\mathbf{h}=\begin{bmatrix}
h_0 & h_1 & h_2 & h_3
\end{bmatrix}$, then the circulant matrix denoted by $\mathrm{circ}(h)$ is $$\mathrm{circ}(\mathbf{h})=\mathbf{H}=\begin{bmatrix}
h_0 & h_3 & h_2 & h_1\\
h_1 & h_0 & h_3 & h_2\\
h_2 & h_1 & h_0 & h_3\\
h_3 & h_2 & h_1 & h_0
\end{bmatrix}.$$
So we can calculate $h[n]*g[n]$ by evaluating the product $\mathbf{H}\mathbf{g}$, where $\mathbf{g}=\begin{bmatrix}
g_0 & g_1 & g_2 & g_3
\end{bmatrix}^T$.

This can be extended to 2-D signals. I explain it through the given example.

The size of signals are $2\times2$ and $3\times3$. The size of convolution is $(2+3-1)\times(2+3-1)=4\times4$. So we need to pad zeros to adjust the size of 2-D signal $x_2(m,n)$:
$$\mathbf{A}=\begin{bmatrix}
 0 & 0 & 0 & 0\\
0 & -1 & 0 & 0\\
-1 & 4 & -1 & 0\\
0 & -1 & 0 & 0
\end{bmatrix}$$
Now construct four circulant matrices corresponding to the four rows of $\mathbf{A}$ denoted by $\mathbf{a}_i$ as follows:

$$\begin{align}
\mathbf{X}_0&=\mathrm{circ}(\mathbf{a}_0)=\mathrm{circ}(\begin{bmatrix}
0 & -1 & 0 & 0
\end{bmatrix})\\
\mathbf{X}_1&=\mathrm{circ}(\mathbf{a}_1)=\mathrm{circ}(\begin{bmatrix}
-1 & 4 & -1 & 0
\end{bmatrix})\\
\mathbf{X}_2&=\mathrm{circ}(\mathbf{a}_2)=\mathrm{circ}(\begin{bmatrix}
0 & -1 & 0 & 0
\end{bmatrix})\\
\mathbf{X}_3&=\mathrm{circ}(\mathbf{a}_3)=\mathrm{circ}(\begin{bmatrix}
0 & 0 & 0 & 0
\end{bmatrix})
\end{align}$$
and construct the $16\times16$ doubly circuland matrix $$\mathbf{X}=\begin{bmatrix}
\mathbf{X}_0 & \mathbf{X}_3 & \mathbf{X}_2 & \mathbf{X}_1\\
\mathbf{X}_1 & \mathbf{X}_0 & \mathbf{X}_3 & \mathbf{X}_2\\
\mathbf{X}_2 & \mathbf{X}_1 & \mathbf{X}_0 & \mathbf{X}_3\\
\mathbf{X}_3 & \mathbf{X}_2 & \mathbf{X}_1 & \mathbf{X}_0
\end{bmatrix}.$$
Now we need to construct the vector corresponding to $x_1(m,n)$. Pad zeros from left and bottom to adjust indices considering $m,n$, and pad zeros from top and right to make the size $4\times4$:
$$\mathbf{B}=\begin{bmatrix}
 0 & 0 & 0 & 0\\
0 & 3 & 4 & 0\\
0 & 1 & 2 & 0\\
0 & 0 & 0 & 0
\end{bmatrix},$$
then construct a $16\times 1$ vector $\mathbf{b}=\begin{bmatrix}
 \mathbf{b}_0 \ |& \mathbf{b}_1 \ |& \mathbf{b}_2 \ |&\mathbf{b}_3
\end{bmatrix}^T$ by concatenating the rows of matrix $\mathbf{B}$.

Finally, the 2-D circular convolution is $$\boxed{\mathbf{G}=\mathbf{X}\mathbf{b}}$$


What are some easy to detect reference mark shapes and corresponding detection methods?

These marks (in the electronics manufacturing industry) are usually called fiducials or fiducial markers.

This paper from Accutron suggests:

Though I have no reference (or even measure) of the "optimality" of the disk / circle.


How does Correlation for 1D Signals work?

For the simple $3$-tap filter given in the example, each output value $r(j)$ is computed as
$$\begin{align}r(j)&=f(-1)I(j-1)+f(0)I(j)+f(1)I(j+1)\\&=-\frac12 I(j-1)+\frac12 I(j+1)\tag{1}\end{align}$$
So each input value is replaced by the (scaled) difference between its neighboring values.


Weighting function for output error IIR filter design

The reason for the worse fit at low frequencies (without additional weighting) is the low density of frequency points at low frequencies due your choice of a logarithmic frequency vector. I'm pretty sure that you will get a better fit at low frequencies if you choose a (more) linear frequency grid.


Intuition behind the scaling property of Fourier Transforms

Let $g(x) = f(ax)$.  Then, using a change of variables $y = ax$,
$$\begin{align*}
G(u) &= \int_{-\infty}^{\infty} g(x)\exp(-j u x)\,\mathrm dx\\
&= \int_{-\infty}^{\infty} f(ax)\exp(-j u x)\,\mathrm dx\\
&= \frac{1}{|a|}\int_{-\infty}^{\infty} f(y)\exp(- j (u/a) y)\,\mathrm dy\\
&= \frac{1}{|a|}F\left(\frac{u}{a}\right)
\end{align*}$$
not $\displaystyle \frac{1}{|a|}F\left(\frac{u}{|a|}\right)$ as you say (cf. Henry Gomersall's
comment on your question).
Now, if $f(x)$ is a periodic function with period $2\pi$, then its
Fourier transform exists only in a generalized sense (involves impulses).
The Fourier series
representation $\displaystyle \sum_n c_n \exp(j nx)$ of $f(x)$
has Fourier coefficients
$$\displaystyle c_n = \frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)\exp(-j n x)\,\mathrm dx,
 ~ -\infty < n < \infty$$
and these also show up in its Fourier transform
which is given by
$\displaystyle F(u) = \sum_n c_n \delta(u -n)$. Note that
the inverse Fourier transform is
$$\int_{-\infty}^\infty F(u)\exp(jux)\,\mathrm du
= \sum_n \int_{-\infty}^\infty c_n \delta(u-n)\exp(jux)\,\mathrm du
= \sum_n c_n \exp(j nx) = f(x).$$
The function $g(x) = f(ax)$ of period $2\pi/a$ thus has Fourier transform
$$G(u) = \frac{1}{|a|}\sum_{n=-\infty}^{\infty} c_n \delta(u/a -n)
= \frac{1}{|a|}\sum_{n=-\infty}^{\infty} c_n \delta\left(\frac{u-an}{a}\right).$$
But, $\delta(\alpha x) = \frac{1}{|\alpha|}\delta(x)$, and so with
$\alpha = a^{-1}$, we have
$$G(u) = \sum_n c_n \delta(u - an)$$
whose inverse Fourier transform is
$$\int_{-\infty}^\infty G(u)\exp(jux)\,\mathrm du
= \sum_n \int_{-\infty}^\infty c_n \delta(u-an)\exp(jux)\,\mathrm du
= \sum_n c_n \exp(j nax) = g(x) = f(ax).$$

So, you are not missing anything. Scaling an ordinary
function scales both
the magnitude and the frequency axis of the Fourier transform,
but when impulses are involved, scaling the argument of the
impulse leads to a magnitude scaling that exactly compensates
for the magnitude scaling imposed on the Fourier transform,
and leaves the Fourier coefficients unchanged.
frequency scaling also scales the amplitude


What's the difference between "Discrete Fourier Series" and "Discrete Fourier Transform"?

They are the same thing. Just the same.

But for some reasons, we want to call the first period of the DFS sequences $\tilde{x}[n]$ and $\tilde{X}[k]$, as the new DFT sequences $x[n]$, $X[k]$.

They are both periodic. The periodicity of DFS is natural extension of periodicity of complex exponential and periodicity of time sequence, but the periodicity of DFT is also forced by the modulus operator.

Be careful! DFT is also linked with finite length aperiodic sequences $x[n]$ and their DTFT via the frequency sampling relation; i.e., $X[k] = X(e^{j \frac{2\pi}{N}k})$.

Because of this latter relation of the DFT $X[k]$ to the DTFT $X(e^{j \omega})$ of aperiodic (and finite length) sequences $x[n]$, in order to remember and imply the inherent periodicity in DFT (due to its definition through DFS) we shall enforce the modulus operator into the arguments of DFT sequences.


How to Estimated the Noise Standard Deviation (STD - $ \sigma $) as a Function of Brightness from an Image?

I didn't read the article you referred to but I can try giving you some idea.

Run along the image for each pixel consider its 9 x 9 neighborhood.
For each pixel's neighbor hood calculate the STD and the Mean.

Use the Mean as the brightness value.
Plot the STD as a function of the mean value.

Some ideas for improvement:

Set the neighborhood size in an adaptive manner to only include data from the same brightness level.
For patches with the same value, aggregate the STD and average it.


Why is White Noise so important in System Identification or Adaptive Filters

I believe the point of feeding white noise into the system is for the filter to adapt its coefficients before actually generating the signal $x[k]$. This would mean there are two "operating modes" for the system: coefficient adapting mode (in which white noise, a broadband signal, is used to adapt the filter to the feedback path), and performing mode (where the signal $x[k]$ is input into the system with the filter already adapted).

So the answer to your first question would be that the noise is output by the loudspeaker and follows the whole signal path. The system adapts it's coefficients to be similar to $F(q)$, and once that process has ended the white noise can be turned off. Afterwards, when generating $x[k]$ the system will already reduce feedback using the $\hat{F}_0(q)$ transfer function that was adapted previously, and the coefficients will only change when the forward path $G(q)$ is modified.

The purpose of using white noise (or a pulse train) is to have a broadband signal that has frequency components throughout the whole audible spectrum so that the adapted filter works for any kind of input signal $x[k]$. If the system was used without first feeding the broadband signal into it, every time $x[k]$ presents new frequency components, the coefficients may need to adapt during the performance/speech, which may have audible effects depending on the convergence speed of the algorithm.


Prefactors computing PSD of a signal with numpy.fft VS. scipy.signal.welch

I would like to post here the link to the accepted answer that I obtained to the question in StackOverflow, as I think it is relevant to this forum.

The original answer by francis is here: post in StackOverflow


Difference between "scaled white noise " vs. "unscaled white noise"

Wavelets, as bases or frames, are linear like-decompositions. As such, signals and processes can be decomposed into a linear or weighed sum of coefficients multipling (wavelet) vectors:

$$s[n] = \sum a_n e[n]$$

The indexing above is not wavelet-specific. Traditionally, there was some natural order: low $n$ denoting low-scales, or high frequencies, and vice-versa. So traditional linear filtering often consisted in keeping indices above or below some cut-off indices (with weighting), assuming that data above and below could be noise (typically low-pass filtering).

This linear (order-related) scheme remains practically useful. However, and some cases, a non-linear selection can prove more efficient. Especially when indices, across time and scale, don't exhibit clear ordering anymore.

$$s[n] = \sum a_{j,k} \psi\left[\frac{n -k2^j}{2^j}\right]$$

The main idea is to select "cleaned signal" coefficient not within indices, but with respect to magnitude: the highest the magnitude, the most likely it belongs to data instead of noise.

Selecting the highest coefficients amounts to setting a threshold below which you discard coefficients. This threshold should be selected as cleverly of possible. nd here is the catch. In the early works on wavelet thresholding, for proofs and practice, it was often considering that the noise was Gaussian with (known) unit variance. This would yield nice theories and practical algorithms.

Yet, this requires a knowledge on the noise, which is often unknown (in shape and amplitude). If we forget the shape (Gaussian or not), if the noise level is $\sigma$ (instead on one), this should be taken into account, esp. by changing the threshold accordingly. Caveat, estimating this $\sigma$  remains tricky.


Origin of brightness artifacts in intensity resolution compressed images

I've taken the original tiff and just set a low threshold for 1bit quantization (conversion to black and white):

so we can clearly make out that these artifacts look very circular, and placed in specific patterns around the image borders.
I think these are artifacts could stem from a resizing of the image with imperfect anti-aliasing, though I'm certainly confused by the shape and the lack of "peaks" along the x and y axes. That would actually suggest lens artifacts, as these tend to be circularly distributed. But I think we can rule this out, since the "dot array" at the edges of the picture clearly indicate this happened after cropping the picture digitally.
Another interesting find: quantizing at the other end yields this:

Looks like "cathode ray burn in" used to look – I would, however, not know how this could happen with a digital camera or a scanner. Maybe it's actually a feature of the face of that clock?


Within the context of optical image/signal processing, the Fresnel transform mathematically describes the manifestation of a physical phenomenon, known as diffraction, on the image formation process of a practical camera system using a finite aperture size.

Diffraction causes blurring on the image, and defines the physical limit of achievable image resolution, given an ideal lens system (which has no distortions of any kind on the image formation process). No better resolution is possible using that aperture size and optical imaging system.

Diffraction can be analysed through the Fourier transformation of the aperture function.


The comments actually clearly put forward the complexity of a serious answer. But I'ld like to go simple and put mine for your example filters.

As you've stated, if two IIR filters of $70$th order and $4$th order (of whatever type) can satisfactorily do the job for you, then you should select the $4$th order filter because it will take much less processing power and processing time to compute its output.

That being said, those two filters (the elliptic and the Butterworth) have different frequency domain characteristics. They both have their own magnitude and nonlinear phase responses.

The Butterworth filter will have a monotonic frequency response in both of its pass and stop bands therefore in the design process it will be an overkill to specify a maximum allowed error ripple (from the ideal response) which will be satisfied at the cutoff freqencies but otherwise will be oversatisfied on the rest of the band due to monotonicity of the magnitude curve. That's why the Butterworth filter requires a much larger order for the same given filter specs.

The elliptic filter is known to have equiripple oscillatory characteristics on both the pass and stop bands which distributes the error energy accross all of the bands, therefore reducing the maximum (peak) of the allowed error ripple for a given filter order or minimizes the order of the filter for a given allowed peak ripple and transition width.

Since you don't want to consider the effects such as coefficient quantization or implementation architectures for high order recursive filters, then there isn't much remains to be considered for those two specific filters.


Slope of 2nd order low pass vs band pass Butterworth filters

By convention, the pole count refers to the number of poles in the design polynomial (Butterworth, Chebyshev, etc) and this defines the amount of roll off in the stop band.

If we maintain this convention, then all filters roll off at the rate of 6 dB / octave / pole in the stop band, whether it is a low pass, high pass, band pass, or notch.

Also, a nth order Elliptic has the same ultimate roll off rate in the stop band as a nth order Butterworth. Of course, the Elliptic's transition band will be much narrower.


How can a zero-padded length n signal be truncated to n coefficients from which a signal can be reconstructed?

From a theory standpoint, that ought to be possible. It's not bound to be simple, however. Zero padding will scale and interpolate the FD representation of your signal. If you remove coefficients the right way, you should be able to reconstruct what's missing.

Let's say you doubled the signal length by padding. You can then remove every other bin of the fourier transformed (starting with the second) and get exactly the result you would have without padding.


I have found the bug !
n2 is indeed a integer in my case, but its not an int()

print(type(n2)) gives:   <class 'numpy.int64'>   !!!

In the previous version of scipy, it was already a numpy.int64 but it was working correctly.

Maybe a test slightly too stringent...

Thanks for helping everyone.


Setting Parity Check Matrix for LDPC channel coding

As specified in documentation, using 'ParityCheckMatrix' you can configure the Parity Check Matrix (PCM) during the constructions of the encoder/decoder objects.
The syntax is encoder = comm.LDPCEncoder('ParityCheckMatrix',pcm) or simply
encoder = comm.LDPCEncoder(pcm);

where pcm is the desired PCM which must be sparse type.
An example for (probably poor performance) 1024x2048 PCM:
pcm = sparse([eye(1024) eye(1024)]);
encoder = comm.LDPCEncoder(pcm);
decoder = comm.LDPCDecoder(pcm);

%data
m = randi([0 1],1024,1);
%encode then decode. Note the convention of the output of comm.LDPCDecoder
y = 1 - decoder(encoder(m));
%error which equals to 0
biterr(y,m)


What's the best criterion for determining a frequency peak?

First convolve the FFT sequence in the given range by a fairly wide gaussian window. In case there is only one peak, this will just widen it. In case there are many peaks (as in the second example), this will merge them into a uniform and less spiky mass.

Then you can use one of the following peakedness metrics:

Kurtosis.
Ratio of geometric to arithmetic mean


How to get a Standard Deviation weighted average in image?

The $k$ and $l$ indices are to be understood with indexes under the $\sum$ sign. They can be called "free variables". Suppose that the block $k$ has only two pixels, then $p=1$ and $q=2$. So $\overline{L\_red(k)} $ is just a weighted average (a kind of center of mass) of the $2$ red components from the two pixels  of the block:
$$ \frac{L\_weight(1,1)red_{1,1}(k)+L\_weight(1,2)red_{1,2}(k)}{L\_weight(1,1)+L\_weight(1,2)} $$
which gives you a single value for each block.
Formula (5) has the same interpretation,  you perform the weighted average of those averaged values over blocks, with weights $SD\_red(k)$.

it is quite easy to use Matlab's block processing tools.


Ok, question solved. I found out that indeed L1 pre/post constellations are not rotated.

The constellation above comprises:

L1 pre (BPSK)
L1 post (64-QAM)
Data PLP from P2 symbol (rotated 16-QAM)

Also P2 cells are interleaved (ETSI EN 302 755 @ 8.5).


Contradictory results from power spectral densities and autocorrelation in shot noise

Your Zener diode will have an intrinsic/paracetic capacitance (you might also have connected a separate capacitor in your circuit as well, across zener' output or nearby) and assocaited to it will be a RC time constant with the resistance across it. For the noise to show it's random nature accross this circuit's tapping point, we would have to wait for a time that is  certain order of the time greater than the time constant. That is why a slower sampling rate is able to better capture that.

That is why you observe that when the sampling rate was 1MHz, the autocorrelation was high all accross and when the sampling interval is of the order of Milliseconds, you have very much acheived a good correlation value because you are nearing the typical time constant of an RC circuit, for the half second sampling the noise fluctuations have indeed taken effect across the tappping point as now the time is much greater than the time constant and you observe true randomness, that is short noise that is white.


It is a million-dollar question.
It depends on the modality and the task the radiologist has to accomplish.

The gold standard is to give radiologists your images, assign them a task and measure how effective they are. However, this requires time, money and clinical images. Typically, this is used as a last resort to convince the regulatory agencies (the FDA for instance) that your medical imaging system is effective.

There are other ways of doing this without radiologists, with "model observers". Basically, you define an algorithm (often a variation on linear template-matching) with is supposed to mimick the performance of a human observer, and you measure the performance of the algorithm. However, not all tasks and modalities are suitable to such an approach.

You can also use metrics measured on the image such as the Signal Difference to Noise Ratio, or metrics measure on the whole imaging system, such as the Detective Quantum Efficiency for X-ray systems.

For a good overview, read this: http://www.ncbi.nlm.nih.gov/pubmed/9419600
Image quality, the ideal observer, and human performance of radiologic decision tasks. Burgess A.


How to Improve Signal Processing Skills?

One of the best ways is to hang out on this community.
Read questions and answers.
Try to replicate results on answers.

Once you have more knowledge, try answering questions of others.
If you encounter a question on a field you don't know, find resources to read about it until you can answer the question (Assuming it is not deep and requires deep knowledge).

Other classic approaches would be books.
Yet, I recommend Hand On approach.
Namely if you read a book, choose those which has code inside them to replicate results and have better understanding of real world models.

You may even look for some courses on Coursera:

Digital Signal Processing.
Audio Signal Processing for Music Applications.
Fundamentals of Digital Image and Video Processing.

There are few options on edX as well:

Discrete Time Signal Processing.
Signals and Systems, Part 1.
Signals and Systems, Part 2.
Discrete Time Signals and Systems, Part 1: Time Domain.
Discrete Time Signals and Systems, Part 2: Frequency Domain.

For more classic recommendation you may refer to What Resources Are Recommended for an Introduction to Signal Processing (DSP)?

Enjoy the ride!