Convolution and windowing using a buffer - how do I do overlap add?

normally i like to define windows to have even symmetry about $t=0$ or $n=0$, but this time i'll do it like MATLAB (except i know how to count from $0$).

Hann Window of width $B$ samples:

$$ w[n] \triangleq \begin{cases} \frac12 \left(1 - \cos\left(\frac{2 \pi}{B} (n+\frac12)\right) \right) \quad \text{for } 0 \le n < B \\
0  \quad \text{for } n < 0 \text{ or } B \le n \\
 \end{cases} $$

another way to write it is

$$ w[n] = \tfrac12 \left(1 - \cos\left(\tfrac{2 \pi}{B} (n+\tfrac12) \right) \right) \cdot \left(u[n] - u[n-B] \right) $$

where $u[n]$ is the discrete unit step function

$$ u[n] \triangleq \begin{cases}
1 \quad \text{for } n \ge 0 \\
0 \quad \text{for } n < 0 \\
\end{cases} $$

all $B$ samples within the window width are non-zero.  let's assume $B$ is even, so that $\frac{B}2$ is an integer.  turns out that the Hann window is "complementary" which means the downslope of one Hann window adds to $1$ with the upslope of the window adjacent to its right if they are staggered by $\frac{B}2$ samples.  in fact this is true:

$$ \sum\limits_{m=-\infty}^{\infty} w\left[n-m\tfrac{B}2 \right] = 1 $$

that means

$$ \begin{align}
x[n] & = x[n] \sum\limits_{m=-\infty}^{\infty} w\left[n-m\tfrac{B}2 \right] \\
     & = \sum\limits_{m=-\infty}^{\infty} x[n]  w\left[n-m\tfrac{B}2 \right] \\
     & = \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{B}2]   \\
\end{align} $$

where

$$ x_m[n] \triangleq x\left[n+m\tfrac{B}2 \right]  w[n] $$

so we broke up the input into overlapping little "wavelets" (my term, not in the sense of Daubechies) or "grains", $x_m[n]$, each delayed, in the summation, by $\frac{B}2$ samples in relation to their preceding neighbor $x_{m-1}[n]$.

now we know that for the $m$-th delayed grain, $x_m\left[n-m\tfrac{B}2 \right] \ne 0$ only for $m\tfrac{B}2 \le n < m\tfrac{B}2 + B$

FIR convolution

we have an FIR filter with causal finite impulse response, $h[n]$, with length $L$.  so $h[n]=0$ for $n<0$ or $L\le n$.  the output of this FIR filter is

$$\begin{align}
y[n] & = h[n] \ \circledast \ x[n] \\
     & \triangleq \sum\limits_{i=-\infty}^{\infty} h[i] x[n-i] \\
     & = \sum\limits_{i=0}^{L-1} h[i] x[n-i] \\
\end{align}$$

OLA convolution

one way to write this:

$$\begin{align}
y[n] & = h[n] \ \circledast \ x[n] \\
     & = h[n] \ \circledast \ \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{B}2] \\
     & = \sum\limits_{m=-\infty}^{\infty} h[n] \ \circledast \ x_m[n-m\tfrac{B}2] \\
     & = \sum\limits_{m=-\infty}^{\infty} y_m[n-m\tfrac{B}2] \\
\end{align}$$

where

$$ \begin{align}
 y_m[n] & \triangleq h[n] \ \circledast \ x_m[n] \\
        & = \sum\limits_{i=0}^{L-1} h[i] x_m[n-i] \\
\end{align} $$

another way to say it is

$$\begin{align}
y[n] & = h[n] \ \circledast \ x[n] \\
     & = \sum\limits_{i=0}^{L-1} h[i] x[n-i] \\
     & = \sum\limits_{i=0}^{L-1} h[i] \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{B}2 - i]  \\
     & = \sum\limits_{m=-\infty}^{\infty} \sum\limits_{i=0}^{L-1} h[i] \ x_m[n-i - m\tfrac{B}2]  \\
     & = \sum\limits_{m=-\infty}^{\infty} y_m[n-m\tfrac{B}2]  \\
\end{align}$$

$y_m[n]$ defined the same way as above.

the potentially non-zero length of the block of input samples $x_m[n]$ is $B$ samples, the length of the FIR is $L$ samples, and that makes the potentially non-zero length of the block of output samples $B+L-1$.  you can see that

$$ y_m[n] = 0 \quad \text{for }n<0 \text{ or }n \ge B+L-1 $$

$y_m[n]$ could be non-zero for $0 \le n < B+L-1$.

unlike the overlapping input segments $x_m[n]$, the length of the output grains $y_m[n]$ need not be the block length $B$, but can be (and are) longer.  so instead of 50% overlap (which is $\frac{B-\tfrac12 B}{B} = \frac12$ or an overlap of $B - \tfrac12 B$ samples) like you might have with the Hann-windowed input, your overlap with the output blocks will be $B+L-1 - \tfrac12 B$ samples.

FFT convolution

the convolution of each block

$$ y_m[n] = \sum\limits_{i=0}^{L-1} h[i] x_m[n-i] $$

can be done with an FFT as long as the FFT size $N$ is at least the length of the output block.

$$ N \ge B+L-1 $$

in this case

$$\begin{align}
X_m[k] & = \mathcal{DFT} \left\{ x_m[n] \right\} \\
       & = \sum\limits_{n=0}^{N-1} x_m[n] \ e^{-j 2 \pi \frac{nk}{N}}  \\
\end{align}$$

$$ Y_m[k] = H[k] \cdot X_m[k] \quad \text{for } 0 \le k < N $$

where

$$\begin{align}
H[k] & = \mathcal{DFT} \left\{ h[n] \right\} \\
       & = \sum\limits_{n=0}^{N-1} h[n] \ e^{-j 2 \pi \frac{nk}{N}}  \\
\end{align}$$

$$\begin{align}
Y_m[k] & = \mathcal{DFT} \left\{ y_m[n] \right\} \\
       & = \sum\limits_{n=0}^{N-1} y_m[n] \ e^{-j 2 \pi \frac{nk}{N}}  \\
\end{align}$$

you need only compute $H[k]$ once for use with all of the blocks. and you get the time-domain samples for each output block with the inverse FFT.

$$ \begin{align}
y_m[n] & = \mathcal{IDFT} \left\{ Y_m[k] \right\} \\
       & = \frac1N \sum\limits_{k=0}^{N-1} Y_m[k] \ e^{+j 2 \pi \frac{nk}{N}}  \\
\end{align} $$


The way that I would implement this would be to have a circular buffer with a slower output pointer than input pointer (or rather an output pointer that considers overlapping frames).

Depending on how quickly you need to perform the FFT and what other operations are happening, the size of the buffer may vary, but it will need to be at least 1.5 times the FFT size. The extra half period can be used to simplify the wrap around issue (assuming you are calling an external library for you FFT and not implementing it yourself). It may be more convenient to use a buffer size that is twice the FFT size or four times the FFT size to maintain a power of 2 for the buffer length.

You are using 2048 in your example, I am going to reduce this to use an FFT size of 8 (for simplicity). I will consider each number in the following sequence to be a single sample (not half a frame as you consider). The input sequence will be $x_1, x_2, x_3, \ldots$ and so on. When the first sample comes, we skip half an FFT size, and put it in there. We retain an input pointer and an output pointer. The input pointer starts at half the FFT size. The output pointer starts at half the FFT size.

So after 4 samples are input, your buffer looks like this.
$$ \left\{ 0, 0, 0, 0, x_1, x_2, x_3, x_4, 0, 0, 0, 0 \right\} $$

Now, 4 more samples will enter. These samples should be mirrored at the end and the beginning of the buffer. So, they will look like this.
$$ \left\{ x_5, x_6, x_7, x_8, x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8 \right\} $$

When the final sample enters, you call your FFT function on the part of the buffer beginning with $x_1$ and ending with $x_8$ (hopefully it isn't an in-place algorithm, if it is then copy to a temporary buffer first).

Then, the next 4 samples will be input. These replace the first four so your buffer will look like this.
$$ \left\{ x_5, x_6, x_7, x_8, x_9, x_{10}, x_{11}, x_{12}, x_5, x_6, x_7, x_8 \right\} $$

Now, call your FFT routine on the part of the buffer that starts with $x_5$ and ends with $x_{12}$. Then input $4$ more (mirroring them to the front and back), and call the routine again.

I hope this clarifies things.


Effect of sign change on ROC & z transform?

You need to understand that a function such as $F(z)=\frac{z}{z+1}$ doesn't have a region of convergence (ROC). Only an infinite series has a ROC. E.g., the $\mathcal{Z}$-transform of a sequence $x[n]$ can be expressed by the series

$$\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{1}$$

and this expression only makes sense if the series converges for certain values of $z$. The ROC is the region of the $z$-plane for which the series $(1)$ converges.

If you turn the problem around and if you're given a function of $z$ like $F(z)$ defined above, you can ask yourself if $F(z)$ is a valid $\mathcal{Z}$-transform of some sequence $f[n]$. It turns out that $F(z)$ is the $\mathcal{Z}$-transform of two different sequences, but with different ROCs. So a $\mathcal{Z}$-transform is only completely defined by the expression for the function $and$ the corresponding ROC.

Coming back to the function $F(z)$, it's straightforward to show that it is the $\mathcal{Z}$-transform of the sequence $f_1[n]=(-1)^nu[n]$ with the ROC $|z|>1$, and of the sequence $f_2[n]=-(-1)^nu[-n-1]$ with the ROC $|z|<1$.

Concerning your last question, as mentioned above, the expression for the $\mathcal{Z}$-transform including the ROC is unique, so the two sequences $u[n]$ and $u[-n]$ cannot have the same $\mathcal{Z}$-transform and ROC. It's actually pretty easy to come up with the $\mathcal{Z}$-transform of $u[-n]$:

$$\begin{align}\mathcal{Z}\big\{u[-n]\big\}&=\sum_{n=-\infty}^0z^{-n}\\&=\sum_{n=0}^{\infty}z^n\\&=\frac{1}{1-z},\qquad |z|<1\end{align}$$

which is of course different from the $\mathcal{Z}$-transform of $u[n]$. Note that the algebraic expression as well as the ROC are different.


advantage of using impulse response function for LTI systems?

A sum of shifted and scaled impulse responses would look more like

$ y[n] = \sum_{k=-\infty}^n x[k] \cdot h[n-k] $

The impulse response of an accumulator is

$ h[n] = u[n] $

where $u[n]$ is the unit step function,

$u[n] = \begin{cases}
1 & n \ge 0 \\
0 & \mathrm{otherwise}
\end{cases}$

so for your example case of

$ x[n] = \{1, 1, 2\} $

The result will be

$ \begin{align}
y[0] &= \sum_{k=-\infty}^0 x[k] \cdot h[-k] \\
     &= 1 \cdot 1 + 1 \cdot 0 + 2 \cdot 0 \\
     &= 1 \\
y[2] &= \sum_{k=-\infty}^2 x[k] \cdot h[2-k] \\
     &= 1 \cdot 1 + 1 \cdot 1 + 2 \cdot 1 \\
     &= 4
\end{align}$

and hopefuly you can see how $y[1]$ will be 2, and $y[n], n> 2$ will be 4, as expected.

Your equations are incorrect for the sum of scaled and shifted impulse responses, and for your first input signal $x_0$.


$X^*(j\omega)$ is the complex conjugate of $X(j\omega)$. So if

$$X(j\omega)=X_R(\omega)+jX_I(\omega)$$

then

$$X^*(j\omega)=X_R(\omega)-jX_I(\omega)$$

and

$$X^*(-j\omega)=X_R(-\omega)-jX_I(-\omega)$$

To gain some more understanding, consider the definition of the Fourier transform:

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

So we have

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

and

$$X^*(-j\omega)=\int_{-\infty}^{\infty}x^*(t)e^{-j\omega t}dt=\mathcal{F}\{x^*(t)\}$$

as claimed in the quote of your question.


What is the ideal FIR length, given a specific pulse length?

I'm converting some of my comments from above into this answer:

If you are trying to detect the presence of some pulse shape $p(t)$ on the AWGN channel, then the optimum detector uses a filter that is matched to the pulse shape $p(t)$ (aptly named the matched filter); this maximizes the signal-to-noise ratio (SNR) at the output of the filter and thus provides the best detection statistic. This approach is equivalent to a sliding cross-correlation of the pulse shape with the observed signal, expressed mathematically as follows:

$$
d(t) = x(t) * p(t)
$$

where $x(t)$ is the observed signal and $d(t)$ is the resulting detection statistic.

The main problem, therefore, consists of selection of an appropriate threshold that can be used to determine where the pulses of interest occur in $x(t)$. Specifically, one would indicate a detection when $d(t) > T$, where $T$ is a threshold used to balance between two opposing performance metrics: probability of detection $P_d$ and probability of false alarm $P_{fa}$. Here is a link to a previous answer where I talk about the tradeoff a bit more. Target values for these metrics would be chosen according to the requirements for your specific application.

For this case, we can come up with general expressions for $P_d$ and $P_{fa}$ pretty easily:

Probability of false alarm: A "false alarm" indicates a case where the detector reports the presence of the target pulse $p(t)$ when it is in fact not present. Since we've defined the channel to be AWGN, that means that for a false alarm, the input signal $x(t)$ is a white Gaussian noise (WGN) process. Without loss of generality, we will assume that noise to be zero-mean with variance $\sigma^2$.

To determine $P_{fa}$, we are concerned with what the signal at the output of the matched filter looks like. Recall that $d(t)$ is defined as:

$$
d(t) = x(t) * p(t) = \int_{-\infty}^{\infty} x(\tau) p(t-\tau) d\tau
$$

Assuming that $x(t)$ is a white Gaussian noise process with variance $\sigma^2$, it can be shown that $d(t)$ will also be Gaussian, with a variance equal to:

$$
\sigma_d^2 = \sigma^2 \int_{-\infty}^{\infty} |p(t)|^2 dt
$$

That is, the variance at the output of the matched filter is just scaled by the total energy of the pulse waveform $p(t)$. The detection statistic $d(t)$ is then just a Gaussian process with variance $\sigma_d^2$. The probability of a false alarm at any given time $t$ is equal to the probability that the detection statistic exceeds the threshold $T$. Using the properties of the Gaussian distribution, we can write this as:

$$
\begin{align}
P_{fa} &= P(d(t) > T\ | \text{ no signal present}) \\
&= 1 - F_d(T) \\
&= Q\left(\frac{T}{\sigma_d}\right)
\end{align}
$$

where $F_d(d)$ is the cumulative distribution function (CDF) of the Gaussian distribution and $Q(x)$ is the Q-function.
Probability of detection: This case differs from the false-alarm case in that the signal of interest is present. Specifically, we examine the situation where the matched filter is perfectly aligned with the pulse of interest. The same noise component that we analyzed previously is present, but the autocorrelation of the desired pulse shape gives it a non-zero mean. This mean is equal to the total energy of the pulse waveform:

$$
\mathbb{E}\left(d(t)\right) = \int_{-\infty}^{\infty} |p(t)|^2 dt
$$

The probability of detection, therefore, is the probability that the detection statistic exceeds threshold:

$$
\begin{align}
P_{d} &= P(d(t) > T\ | \text{ signal present}) \\
&= 1 - F_d(T) \\
&= Q\left(\frac{T - m_d}{\sigma_d}\right)
\end{align}
$$

The design process would then look like this:

Select an operating range for your detector, defining a minimum signal-to-noise ratio (or equivalently, a maximum noise variance $\sigma^2$) at which you will operate.
Assuming the worst-case conditions (i.e. the maximum noise level), select a threshold $T$ that meets either your required $P_d$ or $P_{fa}$ (whichever is more important to you).
Plug the resulting value for $T$ into the other equation to determine your predicted performance metrics.

This is a pretty high-level treatment of the problem, and if you get down to business trying to build something that works practically, you'll run into some other details of note:

One thing that may be relevant for your problem of detections sinusoids is that it's likely that the received pulses will be at some unknown starting phase $\phi$ relative to your template pulse shape $p(t)$. You would then observe a reduction in the correlation peak based on the amount of phase offset, which will wreak havoc with your detector's performance. If that is the case, a noncoherent detector is a better approach:

$$
d_1(t) = x(t) * p(t)
$$
$$
d_2(t) = x(t) * p_Q(t)
$$
$$
d(t) = \sqrt{d_1^2(t) + d_2^2(t)}
$$

where $p_Q(t)$ is a 90-degree shifted (or quadrature) version of the template sinusoidal pulse. The statistics of this case are slightly different and are left as an exercise for the interested reader.
The above treatment assumes that the pulse is received in $x(t)$ at the same power level as the template $p(t)$, which is almost certain to be untrue. Two ways to approach this complication come to mind: either use some sort of automatic gain control (AGC) process in order to steer the received power level to what you would expect, or you can make $T$ an adaptive threshold that adjusts itself relative to the observed signal (e.g. you could try to estimate the background noise variance $\sigma^2$ and then set $T$ appropriately).


Is there any relation between QAM and streaming?

QAM is a digital modulation scheme. As such it is one way of implementing a physical layer that allows to convey digital information over a given medium. QAM is frequently used in all kinds of systems, including wireless (cf. broadcast TV and yes, also WiFi) as well as wired (Ethernet uses some variations of QAM as well).

What kind of information you convey over this link does not matter to the physical layer. You could send documents, voice, videos, whatever you like, as long as you represent it in digital form. The physical layer does in general not need to know what kind of data is being sent by the application layer.

So yeah, if you want to stream video over some medium, chances are the physical layer uses QAM. If it happens to be an optical fiber, it might not, since optical modulation schemes work a bit differently. But you can view these two things more or less independently from each other.


STFT with varying window lengths (like the continuous wavelet transform)

This is sometimes done when creating a log frequency scaled spectrogram (which might better match human time versus pitch perception).
One issue with using multiple sizes of STFTs it that this results in multiple FFT outputs for the same graphical plot point (frequency vs. time).  So you have to decide how to select from, scale, mix, interpolate, and/or cross-fade the overlapping FFT information without adding (more) visual banding artifacts, and in the frequency domain as well as the time domain, where the resolution difference problems can be orthogonal.


Digital filtering based on frequency component amplitude... why not?

I can tell you why I avoid such methods. When you pick some frequencies above threshold and remove others, you are using a kind of hard limit. In DSP this is same as Rectangular Windowing in frequency domain. It's effect in time domain will cause ringing. Because when you multiple in frequency domain with rectangular window, in time domain the signal gets convolved with a sinc function. There is a whole wikipedia page dedicated to this (Gibbs Phenomenon) - https://en.wikipedia.org/wiki/Gibbs_phenomenon

Again it depends on our requirements though. For example, if you see your filtered signal, the rise and fall time are not steep compared to original signal because you have eliminated higher frequencies which contribute to this. If my requirement was to apply a filter without compromising on fast rise time, I would not choose this method.


Generating sound of particular frequency matlab

Try wave = amp*sin(2*pi*freq/Fs*t) instead. Since Fs maps to $2 \cdot \pi$, freq maps to $\frac{2\cdot\pi\cdot {\tt freq}}{\tt Fs}$.

Could be useful for reference on Wikipedia.


I can only answer your second question:

"How can the loop bandwidth in GNU Radio synchronization be configured as a percentage of the symbol rate?"

The tracking loop in the symbol synchronizer block operates at the symbol rate, estimating timing error and making a correction once per symbol.  So the sample rate of the error signal from the TED is at approximately 1 sample/symbol.  (I say approximately, because technically the block is continually estimating symbol clock error and adjusting its estimate of the symbol clock period, but its objective, nominal operating rate is 1 error sample/symbol.)

This means that in terms of normalized digital radian frequency, $\omega T_s$, for the loop filter, the symbol frequency corresponds to $2\pi$ radians/symbol.

The loop bandwidth parameter of the symbol synchronizer block is expected in units of the normalized digital radian frequency.  So when one specifies $\dfrac{2\pi}{200}$ for $\omega_n T_s$, one is specifying an approximate one-sided loop filter bandwidth of $2\pi \cdot 0.5\%$ radians/symbol, a one-sided filter bandwidth that is $0.5\%$ of the symbol frequency.

To express the loop bandwidth input of the symbol synchronizer block as a percentage of the symbol frequency, just use

$$\mbox{(one-sided) loop bandwidth} = \omega_n T_s = 2\pi \cdot n\%$$


Image processing based on wavelet transform in python

Three aspects should be taken into account. Suppose that the image is initially positive, coded on eight bits. Then:

standard wavelet oefficients become real, and are not converted to uint8 easily, because wavelet filter coefficients are often real (and not dyadic rationals)
since the filters are orthogonal, the coefficients risk to grow, and exceed the initial $[0,\ldots,255]$ range. So a mere rounding is dangerous
even if you use integer wavelets (by lifting), the same happen, and in 2D you may need $1$ more bit for each level
some detail coefficients will be negative, so uint8 will destroy the signs in a way you cannot recover afterward.

So, the inverse wavelet tranform is likely to be crappy. THe simplest way around:

divide all coeffciients by their maximal absolute value (and keep track of it),
code all the $[0,1]$ magnitudes on $7$ bit,
code the sign on the $8$th bit.


If the function is zero-valued throughout then so will be the output. But if you exclude this trivial homogeneous case, then your question can be reformulated as:

Given any function $f$ of compact support and any list $l$ of frequency domain zeros, does a function $g$ always exist that equals $f$ within the support of $f$ and has a Fourier transform that has the zeros given in $l$, with other zeros allowed?

Starting at $g_0 = f$, it would be enough to demonstrate a procedure $p$ such that $g_{i+1} = p(g_i, l_i)$ that adds a single given frequency domain zero $l_i$ to function $g_i$ without modifying it within or to the left of the support of $f$ and without altering the existing frequency domain zeros. We choose a causal filter as that procedure. The procedure being a filter guarantees that existing frequency domain zeros are preserved. To not change the function within or to the left of the support of $f$, the filter must have an impulse response that has a Dirac delta function followed by a zero-valued segment of the same length as the support of $f$. A comb filter consisting of two Dirac delta functions, the second one inverted and located after the zero-valued segment, matches this description. The comb filter can be made to have a zero at any wanted frequency by placing the inverted Dirac delta to an integer multiple of the period of that frequency. Choosing a large enough multiplier moves the inverted Dirac delta to outside the zero-valued segment.

If the functions are discretely sampled sequences rather than functions of a real argument, a Dirac delta function that lands between the sample points gives a sampled sinc function that has no compact support. To circumvent this problem, instead of a sinc impulse, a fractional, unity gain delay at any given frequency can be created by setting to certain values the two samples surrounding a multiple of the period of that frequency. The delay and gain at other frequencies will be off but this matters not.

So (unless I missed something), the answer to your question is: No, there are no limitations if the (allowed to be long enough) tail is thrown away.

Note that successive application of the above procedure does not generate the impulse response of your original question but an output (one of many possible) that enables a realizable impulse response.


Multiple-image dense point cloud reconstruction with camera extrinsics/intrinsics

First, a warm welcome to SE!
Basically, you have a calibrated 3D reconstruction problem. The typical approach follows a 5-stage pipeline:

Identify 2D features in each image along with the associated descriptors. Algorithms such as SURF, SIFT or AKAZE are heavily used and are available in many vision libraries such as OpenCV.
Match the extracted keypoints across the whole dataset by using their descriptors. Once again, this is a well studied topic. The keypoints and their multiview association will be called feature matches and will form the basis for 3D reconstruction.
Extend the pairwise matches to the multiview case by linking the cameras to one another, associating the ones that have enough common points. In MATLAB, the relevant function would be findTracks.
Use an N-view triangulation algorithm (typically a multi-ray intersection) to find an initial estimate of the 3D points. In MATLAB one could use the triangulateMultiview to solve for the multiview intersection. Normally, prior to this stage, one would solve for the camera poses, but since you have it already, I'm skipping this part.
Now use bundle adjustment to refine the 3D structure. I am not sure how you obtain the poses, but it is strongly recommended to also adjust the camera extrinsics additionally to the points. Libraries such as Ceres Solver are capable of solving for both structure and motion, simultaneously.

Optionally, a post-processing step consisting of densification can be run. PMVS could be of help here:

  Accurate, dense, and robust multi-view stereopsis. Yasutaka Furukawa
  and Jean Ponce. IEEE Trans. on Pattern Analysis and Machine
  Intelligence, 32(8):1362-1376, 2010.

It is another possibility to benefit from existing multi-stage SfM pipelines, such as OpenMVG. You would simply replace the OpenMVG poses by your camera poses and keep the rest of the workflow the same. This way, you get some outlier treatment and post-processing steps for free. MATLAB also has an SfM toolbox that you could easily test.


2D Frequency Domain Convolution Using FFT (Convolution Theorem)

Similar to your question Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB the issue is what happens when you multiply in 2D in frequency domain.
So few remarks about that:

Multiplying in frequency domain for discrete signals with finite support is equivalent to applying convolution in spatial domain under the assumption of cyclic / periodic boundary conditions.
In image processing we usually define per kernel the anchor pixel of the kernel. Usually it is marked as (0, 0) of the kernel. We also mostly set it as the center pixel (In Image Processing most kernels have odd length). When we pad the kernel to the size of the image we usually add zeros on its bottom and right. Which means its (0, 0) isn't aligned with the image.

The misalignment with the circular boundary extension yields the following for the naïve code:
clear();
close('all');

gaussianKernelStd       = 0.5;
gaussianKernelRadius    = ceil(5 * gaussianKernelStd);

mI = im2double(imread('cameraman.tif'));
mI = mI(:, :, 1);

numRows = size(mI, 1);
numCols = size(mI, 2);

vX = [-gaussianKernelRadius:gaussianKernelRadius].';
vK = exp(-(vX .* vX) ./ (2 * gaussianKernelStd * gaussianKernelStd));
mK = vK * vK.';
mK = mK ./ sum(mK(:)); %<! The Gaussian Kernel

mIFiltered = ifft2(fft2(mI) .* fft2(mK, numRows, numCols), 'symmetric');
figure();
imshow([mI, mIFiltered]);

As you can see at the top and left filtered image (The right) has artifacts which are the result of the circular extension and the misalignment. How to fix it?
Well, padding the image correctly and padding with circular extension the kernel.
I showed it in Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB.


Spectrogram, time vs frequency localization, and length of signal

The frequency resolution stays the same because increasing the sampling rate increases the maximum frequency that can be resolved. For example, doubling the sampling rate would double the maximum frequency.
Adjust the size of the spectrogram window, and set the number of FFT points based on that. The command is

[S,F,T] = spectrogram(X,WINDOW,NOVERLAP,NFFT,Fs)

X is the signal
WINDOW is the size of the window, try 1000 for starters
NOVERLAP is the amount of overlap as the window is slid along. signal length / (WINDOW-NOVERLAP) is approximately how many time points you will have.
NFFT is the number of FFT points. The maximum I would use is half the WINDOW value, say 500.
Fs is the sampling frequency. Input this and F will come out as the frequency points and T as the time points with no need for your own calculations.

You should not be running out of memory - I routinely calculate spectrograms on data with millions of points, you just have to tweak the input parameters. For example, putting WINDOW = 1000 and NOVERLAP = 999 will give you one time point for every sample and the output spectrogram will be, in your case, 74767 wide! Try the above settings and see how it goes.

If you make WINDOW large your frequency resolution will increase but you time resolution will decrease.


Terminologies - lags, order in time series model

Q1: should the model generate a time series of length 'N=16` i.e, would the output of the above model $\mathbf{y} = [y_1,y_2,\ldots,y_N]$ contain 16 elements where $n = 1,2,\ldots,16$?

If one thinks about your question in terms of FIR filtering, then an input signal, $x$ of length $N$ filtered with an FIR filter of length $p$ ($p$ coefficients) will have an output of length $p + N - 1$.

Note: there are ways to define the output so that only $N$ outputs are obtained, but you need to understand that is what you are doing.

  Since the order p=2 would the physical interpretation  be that the channel can take in only 2 data values as  input?

It means the channel does not change the input much at all. It's a two coefficient FIR filter.

  Also does p=2 mean that the channel has 2 paths through which the data input travels?

One way to interpret the FIR (MA) nature of the filter you are using is that it involves $p$ paths with a time delay of $nT$ (where $T$ is the sample period).

  Q2: Is there a relationship with the number of input elements and the number of lags?

In general, no, there will not be a relationship.

  (1) Is this related to "multipath" often used in Rayleigh fading? If there are $p$ paths then would the transmitted signal travel in two directions?

The usual method of modeling multi-path is to say that the receiver receives:
$$
y^c(t) = \sum_{i=1}^{p} h^c_i x^c(t - d_i)
$$
where this is the continuous-time version so the $d_i$ are not integers and $d_1 = 0$.

Moving to the discrete-time version means you need to sample everything using $t = nT$.  So if you stick with $p$ components, you can only really have $d_i = (i-1)T$.

  (2) What is the general way to set $nT$? Using Matlab, I did not do any thing in my code. Can you say what is $nT$ for my example?

$T$ is determined by your modeling. How did you relate the actual continuous-time problem to the discrete case?


Multiplication of two impulse function $\delta(t)\cdot \delta(t)=?$

The property

$$f(t)\delta(t-t_0)=f(t_0)\delta(t-t_0)\tag{1}$$

is only valid for a function $f(t)$ that is continuous at $t=t_0$. Since the Dirac delta impulse $\delta(t)$ is not a function (it is a distribution) and since $\delta(t)$ is not continuous, property $(1)$ does not hold (and does not make sense) for $f(t)=\delta(t)$.

The product $\delta(t)\cdot \delta(t)$ is undefined. Also note that the quantity $\delta(0)$ is meaningless. Since $\delta(t)$ is not a function you cannot determine its value for any value of its argument.


what is the starting point of a contour as given by cvConvexityDefect() function in OpenCV

If you can obtain coordinates for the points you show in the right hand image, and you know the points relating to the shoulders, could you presume that the palms were further from the shoulder than the elbow is?

If you assume this to be true, then you simply need to find the distance from the shoulder to the elbow & palm and call the point with the larger distance the palm.

Of course, you could bend your arm in such a way that the palm would be closer to the elbow, so perhaps this isn't an ideal method...


How to derive an expression for the optimal importance distribution?

Hint

Well, it looks like the question tells you what the optimal importance distribution is! It's $p(\mathbf x_k | \mathbf x_{k-1},\mathbf y_{1:k})$. So the problem boils down to computing a closed form expression for this conditional density using the given state space model.

More Hints

The transition model tells you $p(\mathbf x_k | \mathbf x_{k-1})$ is a multivariate Gaussian. The observation model tells you $p(y_k|\mathbf x_{k})$ is another multivariate Gaussian.

Can you write these multivariate Gaussians in closed form?
Can you express $p(\mathbf x_k | \mathbf x_{k-1},\mathbf y_{1:k})$ in terms of $p(\mathbf x_k | \mathbf x_{k-1})$ and $p(y_k|\mathbf x_{k})$?

Almost complete solution

\begin{eqnarray}
p(\mathbf x_k | \mathbf x_{k-1}, \mathbf y_{1:k}) &=& p(\mathbf x_k | \mathbf x_{k-1}, \mathbf y_{1:k-1}, y_k) \label{1}\tag{1} \\
&=& p(\mathbf x_k|\mathbf x_{k-1}, y_k) \label{2} \tag{2}\\
&=& \frac{p(y_k|\mathbf x_k,\mathbf x_{k-1}) p(\mathbf x_k | \mathbf x_{k-1})}{p(y_k | \mathbf x_{k-1})} \label{3} \tag{3} \\
&=& \frac{p(y_k|\mathbf x_k)p(\mathbf x_k | \mathbf x_{k-1})}{p(y_k | \mathbf x_{k-1})} \label{4}\tag{4}
\end{eqnarray}
where (\ref{1}) is just another way of writing $\mathbf y_{1:k}$, (\ref{2}) follows from Markov property (Eq. 4.2 of the book), (\ref{3}) follows from Bayes' theorem and (\ref{4}) follows from conditional independence of measurements (Eq. 4.3 of the book).

Using the model given in Eq. (7.49) in the original question we have
$$
p(y_k|\mathbf x_k) \sim  N((1 \;\; 0) \mathbf x_k,10^2) \label{5}\tag{5}
$$
and
$$
p(\mathbf x_{k}|\mathbf{x}_{k-1}) = N\left(\left(\begin{array}{}
1 & 1\\
0 & 1
\end{array} \right)\mathbf x_{k-1}, \left(\begin{array}{}
1/10^2 & 0\\
0 & 1
\end{array} \right)\right). \label{6}\tag{6}
$$
To compute $p(y_k | \mathbf x_{k-1})$ we express $y_k$ in terms of $\mathbf x_{k-1}$ as follows:
\begin{eqnarray}
y_k &=& (1\;\;0)\mathbf x_{k} + r_k \\
  &=& (1\;\;0)
\left(\begin{array}{}
1 & 1\\
0 & 1
\end{array} \right)\mathbf x_{k-1} + (1\;\;0)\mathbf q_{k-1} + r_k \\
&=& (1 \;\; 1) \mathbf x_{k-1} + (1\;\;0) \mathbf q_{k-1} + r_k
\end{eqnarray}
which implies that
$$
p(y_k | \mathbf x_{k-1}) \sim N\left((1 \;\; 1) \mathbf x_{k-1}, \frac{1}{10^2}+10^2\right).\label{7}\tag{7}
$$
Finally you should use Wikipedia or Eq. (A.1) from Appendix A.1 of the book for the expression for a multivariate Gaussian density function.

Time for grunge work. Plug (\ref{5}), (\ref{6}) and (\ref{7}) into (\ref{4}) and finish with some algebra.


How to reliably compute the group delay of a comb filter

Since this is an FIR, the group delay is D=(N-1)/2=20 samples.

No, since this is a linear phase (i.e. symmetric or anti-symmetric) filter, the group delay is half the length! (being a FIR isn't sufficient.)

  The issue is that I get too peaks in the cross correlation, one at zero lag and another at 20 lag.

Write down the formula for auto-correlation at zero lag. Compare that to the formula of "energy of a signal". They are identical!

This really shouldn't surprise you!

  This method work really well with other filters like averaging filters.

This method works with anything that has a non-zero zero-lag coefficient.

  Does anyone know while I get a the peak at zero which is overshadowing the real peak?

Yes, because autocorrelation at zero is simply the energy. And since correlation is a linear, and your system passes through the original signal, plus a delayed version of it, you get the sum of the auto-correlation of the input signal and the cross-correlation of your delayed signal and the input signal.

The 20-lag peak is no "realer" than the 0-lag peak.

  Is this normal for comb filters?

This is normal for any linear time-invariant system.

  Is there another method to compute delays using the filtered and unfiltered signal other than cross correlation?

The group delay is really defined as derivative of the phase of your signal over frequency. If in doubt, estimate the spectrum of your system, and derive its phase. You'll notice that only a few specific systems (linear-phase, see above) have constant group delay.

Hence, I'm not sure your cross-correlation had much to do with group delay to begin with.


Relationship between windowed fourier transform and wavelet transform

Wavelets commonly factor the frequency of the basis vectors into their temporal widths, some with more than one temporal location.  Windowed FFTs fix the temporal width of each basis vector to the length of the FFT, thus independent of any bin frequency, and with the basis location fixed to the center of the FFT.

I suppose one could consider the FFT to be similar to a wavelet transform with a bad set of wavelets (too many high frequency ones and with poor temporal/locality resolution, and too few low ones and with poor frequency resolution).


Fourier Transform minimum sample requirement

Please define "specific frequency".

Any signal with no frequency components above the Nyquist limit can be exactly recorded and re-constructed.

However, the number of samples determines the frequency resolution.

if you mean "i want to be able to tell the difference between silence and a 1000hz signal, then you would need 4 samples taken less 1/2000th of a second apart - using that you can solve for the 0 frequency signal (dc offset), and the strength of the other frequency.

if you mean "i want to be able to separate a 900hz signal from a 1000hz", then you need to resolve 11 frequency windows, ie, you'll need 22 samples taken at at least 2000 samples per second; then you'll have measured "frequencies around 0, 100,200,300...900,1000 hz", and you can answer "do i have a 900hz signal".

if you want to measure to the nearest hz the frequency of a N hz signal, then you would need to take at least 2N samples, taken at atleast 2N samples/second.


RMS is essentially a three step process on a signal:  Square it, mean it, root it.

Squaring is not so difficult in software.  Two things to keep in mind here.  First, squaring fixed point data will double the data width, so make sure you have enough precision to avoid overflow.  Second, squaring a signal will double its bandwidth.  Make sure you have enough margin in your sample rate to either eliminate aliasing, or that the aliasing doesn't affect your baseband signal.

A square root can be computationally expensive, but it isn't exactly difficult to do.  I does have an impact on bandwidth, but this could be treated as a data representation problem, for which you would not consider it to affect bandwidth.  I wouldn't worry about it.

A mean technically requires the full data set, so it doesn't apply to real time processing.  To get around this, we can instead apply low pass filtering.  Low pass filtering is a big enough topic that I don't think it warrants much discussion here, but you could any standard filtering techniques, such as IIR, FIR.  I would recommend trying some and maybe asking another question on here if you need more information.

It isn't explicitly required, but it is likely that you would need to decimate your signal at some point.  For example, your sample rate may be 1MHz, but the user interface might only refresh at 100Hz.  This would happen either after or during your filtering stage.  You could apply a pair or second order maximally flat IIR low pass filters with a cutoff of say 25Hz, then throw away 9,999 out of every 10,000 samples.  You could also save 10,000 samples in a buffer at a time, and calculate the average of this set.  This is (I believe) how an oscilloscope does its math, but keep in mind that the accuracy wouldn't be great, because the filter is not aggressive enough for the decimation.  You could 'smooth' the data by applying a window to the buffered data, but you would still get some aliasing if your buffers do not overlap.

Lastly, the bandwidth of the input signal.  At 1MHz, 500kHz is at Nyquist, so it is the absolute maximum bandwidth you can sample without aliasing.  However, this doesn't provide any margin, so your sampling is likely to be sub par because of real world sampling limitations, such as quantization error, and noise.  Additionally, squaring a 500kHz would create content at 0Hz as well as 1MHz.  The 1MHz content would alias back down to 0Hz, causing distortion of the 0Hz signal.  As such, you would want to build in some margin to avoid aliasing.

All this being said, I don't know how much proper signal processing technique goes into your off the shelf DMM.  In fact, it is possible to do a fair bit of this process in the analog domain to offload processing power.  It seems more likely that someone designing a DMM would be inclined to throw high precision, high bandwidth ADCs at the problem and brute force a solution that provides a number.  That number might be accurate for a steady state sine wave at a given frequency range, but inaccurate for real world signals.  However, it is unlikely that the user would either notice or care.


Carry the units and it becomes understandable.

T =  0.03125 seconds per sample

Thus your sampling rate is

1/T = 32 samples per second

Each bin corresponds to the cycles per frame.  Your frame has four samples.  Suppose the bin index is k.

k (cycles per frame) * 32 (samples per second) / 4 (samples per frame) = 8k cycles per second = 8k Hz

So the frequencies for your bins, in Hz, become

k   f_k
0   0
1   8
2   16
3   24

However, the DFT wraps around, and choice of the principal frequency among many aliases is a matter of convention.  It is mathematically sounder to consider the upper half of the DFT as negative frequencies.  For real valued signals, the upper half is just the complex conjugate of the lower half and thus provides no new information so it is usually ignored.  So, if you want to consider f_3 as 24 or -8, it doesn't really matter.  For a real valued signal, a negative frequency is the same as a positive frequency.


How would I go about detecting and following edges in time lapse data?

Deep learning might have something to offer in the long term but at the expense of training the network with very detailed information about what constitutes an "outer wall". A deep learning network could then iterate all frames and markup those blocks that are very likely to constitute the "outer wall".

In the short term, a technique based on normalised cross correlation might offer more satisfactory results. The idea in this case is simple: Define an area of the image (let's call it a Region of Interest (ROI)) that you want to track and then try to "find" that area within the next frame and the next frame and so on.

The primary assumption here, of course, is that the ROI doesn't change significantly from frame to frame. And in this particular case, this takes us back to what exactly does an "outer wall" mean?

Here is an example of tracking a point on the "outer wall":

And here is the data of the track (frame num, X, Y).

This is the best tracking I could do out of approximately 8 attempts (I will come back to this) and it was done using Blender.

Blender is using the normalised cross correlation method to track "markers" from a video feed. Depending on the type and number of markers tracked, you can do a number of things such as stabilise video or solve for the camera movement and position in order to synchronise a "virtual" and a real scene.

In this particular case, I am only tracking Location. I could track Location, Rotation, Scale or any combination of those but I figured that this phenomenon is happening more or less on a narrow depth of field. So, just tracking for location will not distort the results.

The tracked coordinates are the green and red timeseries, green for a "Y" coordinate and red for a "X" coordinate.

As you can see, it's not a simple linear movement. The tracked marker seems to be prescribing a tiny little ellipse around its position.

And this brings us to the main "problem" you are going to have to deal with: WHAT is an outer wall? Here are a few observations:

The "outer wall" displacement varies across it's length. The "outer wall" of this membrane (?) is roughly parallel to the width of the video. So, the recorded displacement, really depends on WHERE you are going to track. If I tried to track somewhere else along the length of the track I would have had a different profile.
The "WHERE are you going track?" question is made harder, because the "outer wall" is semi-transparent and furthermore, its transparency throughout it's "pumping" cycle seems to depend on the phase of it. So, there are points in the cycle that there are very clearly defined edges (which is good for accurate tracking) but gradually, these points vanish and introduce drift to the tracking. Gradually, the marker will "slip" from the target it is supposed to be tracking and start tracking a slightly adjacent image patch that happens to look similar.

Therefore, what I would suggest is to set-up a similar script that opens a video, establishes a marker patch and then tracks the patch frame-by-frame, using information that is readily available in many different scientific platforms (for example, for another example).

BUT!, having now established tracking, you will have to go ahead and manipulate the successive frames of the video in order to preserve the features. This could be a simple contrast enhancement with a steep window, or a histogram equalisation or a normalisation step or possibly the addition of a high pass filter to make the edges more pronounced and improve tracking.

Finally, if you would like to replicate this tracking, using your video, you will need:

Blender
Some demo on tracking. This one is very good.
Specific parameters

This tracking was produced with: Tracking Settings(Preset, Blurry Footage, Motion Model: Loc, Match: Previous Frame). You can match either on a frame-by-frame mode, or to a keyframe, in which case, you get displacements relative to that particular frame. Might be useful for you if you want to measure everything with respect to a particular point in the cycle.

A script to extract the tracking data to a CSV file. Like this one.

(Of course, you can track as many markers as you like and their data will be put in a separate file).

Hope this helps.


periodic noise detection in image's frequency domain

Even though the lines are periodic, but they are closer to a train of impulses than a harmonic, so they will leak to many harmonics.
The best feature to distinguish the lines from the rest is that the lines are either horizontal or vertical, vertical lines are represented at the horizontal borders of the FFT, horizontal lines are represented at the vertical borders of the FFT.
Here is the result if you filter out (set to zero) the first four and the last three lines and column of the FFT.

Answering a comment
Here you can see the frequency domain log(fftshift(abs(fft2(image)))), before and after filtering.

Bonus
Just because maybe you would be interested, you can apply morphological filters and get a result like this :D

Offtopic
Yesterday the the author said this

He accepted the answer, but after unaccepted


Shift in Time Domain After DFT Based Convolution

Note that multiplying the FFTs and taking an IFFT is equivalent to cyclic/circular convolution. This is not the same as linear convolution. However, your signals have sufficient zero-padding for the two types of convolution to be equivalent. The correct result of the linear (and cyclic) convolution is what you obtained via FFT and IFFT, a triangle between indices $20$ and $60$, as shown in the figure. This result is easy to predict using graphical convolution, as explained in this animation (from Wikipedia).

The strange result of np.convolve with the option same can be explained as follows: the result is the center part of length $100$ of the linear convolution of the two signals. Since the linear convolution is a signal of length $199$ ($=N_1+N_2-1$) with $99$ trailing zeros, the same option cuts out the center part between indices $50$ and $149$. So the result is the right half of the cyclic convolution with $49$ zeros added at the end, as shown in the bottom right figure.


A non-rectangular window will remove "noise" from distant bins at the cost of adding more "noise" to the immediately adjacent bins to a narrow-band spectrum peak.  The sum of both these spectral leakage effects is greater than zero for a non-rectangular window. So if you count the raising of the level total of all adjacent bins as noise, then the S/N ratio is lowered.

Some people don't care about the bins immediately adjacent to a spectrum peak (their spectral peaks are a priori assumed to be widely spaced; and/or they interpret, interpolate, or phase-vocoder adjust the energy out of those adjacent bins back into the central peak bin), so for those purposes, the reduced far-side-lobe energy means less noise.

Another reason for a lower S/N ratio is that windowing of quantized data is an informationally lossy process, and these (re)quantization losses can also be considered a form of noise.

"Best" is relative to some weighting of quality metrics, and different users may have very different weightings.  If you don't have a set of prioritized design goals for which to optimize a window, then you may not have a strong reason to not just use a Von Hann window.

Depending on your data source and your needs, using just some windowed FFTs may not even be a good form of spectrum analysis, much less the best possible.  Or the opposite.


Removing periodic spike noise from ECG signal

In the first place I would suggest you to look for some electrical methods of avoiding those spikes... If it's impossible to do so, then you must find a close enough mathematical model of your noise generation process from which you may decide to choose an adequate method of noise removal...Provided you can do so...


What is the general form of a transfer function

You may write the second equation on the first form by using $b_i = 0$ corresponding to those poles appearing at $s = 0$. Hence, the only difference is that in the second form you know that there are $m$ poles at $s=0$ while in the first form they may still be there, but one will have to check the values of $b_i$ to determine if they are there.

Note that a system with poles at $s=0$ is a marginally stable system and that the poles are typically introduced by integrators, which have a transfer function $1/s$.


This is a digital state variable filter. Note that your pseudo-code of the filtering process is not correct. You replaced all += in the original code by =, which results in a filter that wouldn't work as desired.

The actual recursive equations implemented by the code are

$$\begin{align}l[n]&=l[n-1]+fb[n-1]\\
h[n]&=q\big(x[n]-b[n-1]\big)-l[n]\\
b[n]&=fh[n]+b[n-1]\end{align}\tag{1}$$

where $l[n]$, $h[n]$, and $b[n]$ are the low-pass, high-pass, and band-pass outputs, respectively. The equations $(1)$ implement a second-order IIR filter (actually three).


Calculating frequency shift of short bursts

Try using the instantaneous frequency

$$f [n] =\dfrac {F_s/2} {\pi}Arg (x [n] \cdot x^*[n-1]) $$

and average it over the burst.


Audio Steganography: Inaudible Audio Watermarking for Source Identification

Petrus. First of all, thanks for attaching my github library as an update into your question for other readers. I have just seen your question.

According to my own experiences, echo hiding methods are fragile towards additional echoes, so are towards reverb. A simple spread-spectrum model using a PN Sequence, which introduces some noise, embeds data into high-frequency bands, so hidden data is usually audible and not robust against low-pass filters.

Transformations such as cosine transforms and wavelet transforms are often used in audio watermarking. Applying a spread-spectrum method without a PN Sequence onto Approximation Coefficients of the third level wavelet decomposition (A3) would be very robust against many common attacks except high-pass filter, and would be almost completely imperceptible. On the other hand, hiding data onto Detail Coefficients of different levels, gives different balance between robustness & transparency. These decompositions are done very easily in Matlab if you have Wavelet Toolbox. However, I have not uploaded my related works into github yet.

We have not worked with real-time audios, but I know that all of these techniques can be applied for them as well. I am not sure whether you already solved your problem, but if you have any further questions you can point out in a comment. I am not an expert in digital signals unlike many of you on this site, so different opinions are also very welcomed.

Best Regards.


Compensation for Irregular Time Step for DFT (FFT)

One option, in addition to the Non-uniform DFT that has been mentioned, is simply interpolating the data on a regularly-spaced grid by using a timestep approximately the same size as the smallest timestep, and then taking the DFT on the interpolated data.


Calculating the 3dB frequency of an nth order low pass filter and nth order high pass filter

In general there is no straightforward analytical solution. As you know, you need to solve
$$\left|H(e^{j\omega_c})\right|=\frac{1}{\sqrt{2}}\tag{1}$$
for $\omega_c$, where it is assumed that the maximum filter gain equals $1$.
For Butterworth filters, the specified cut-off frequency always equals the $3\textrm{ dB}$ frequency. This is not the case for other filter types.
In general, you have to either estimate $\omega_{c}$ from the magnitude plot, or you solve $(1)$ numerically.

EDIT: For Chebyshev filters it's relatively straightforward to compute an analytical expression for the $3\textrm{ dB}$ cut-off frequency.
As an example I'll use a Chebyshev II highpass filter, as mentioned in the question. The magnitude of the normalized frequency response of the analog prototype filter is given by
$$\left|H(\Omega)\right|=\frac{1}{\sqrt{1+\frac{1}{\epsilon^2T_n^2(\Omega)}}}\tag{2}$$
where $T_n(\Omega)$ is the $n^{th}$-order Chebyshev polynomial, and $\epsilon$ is a parameter determining the ripple size in the stop band. The normalized $3\textrm{ dB}$ cut-off frequency $\Omega_c$ can be determined from the equation $\left|H(\Omega_c)\right|=1/\sqrt{2}$, which is equivalent to
$$T_n(\Omega_c)=\frac{1}{\epsilon}\tag{3}$$
From the trigonometric definition of the Chebyshev polynomials
$$T_n(\Omega)=\cosh\big[n\; \textrm{arcosh}(\Omega)\big],\qquad\Omega>1\tag{4}$$
we can solve $(3)$ for $\Omega_c$:
$$\Omega_c=\cosh\left[\frac{1}{n}\; \textrm{arcosh}\left(\frac{1}{\epsilon}\right)\right]\tag{5}$$
Now we just need to transform $\Omega_c$ to the discrete-time domain according to the warping introduced by the bilinear transform:
$$\omega_{c}=2\arctan\left[\Omega_c\tan\left(\frac{\omega_s}{2}\right)\right]\tag{6}$$
where $\omega_s$ is the specified stop band edge of the discrete-time Chebyshev II high pass filter.
Note that we can also rewrite $(6)$ to determine the required value of $\omega_s$ given a specified $3\textrm{ dB}$ cut-off frequency $\omega_c$:
$$\omega_s=2\arctan\left[\frac{1}{\Omega_c}\tan\left(\frac{\omega_c}{2}\right)\right]\tag{7}$$
The following Matlab/Octave script illustrates the procedure of computing the $3\textrm{ dB}$ cut-off frequency of a Chebyshev II high pass filter:

% example filter specs of Chebyshev II high pass filter
n = 6;          % filter order
Rs = 60;        % stop band ripple in dB
ws = .5 * pi;   % stop band edge of Chebyshev II filter

% Cheby II: stopband ripple = ep / sqrt( 1 + ep^2 )
ep = 1 ./ sqrt( 10^( .1 * Rs ) - 1 );

% compute normalized 3dB cut-off frequency of analog prototype
Wc = cosh( acosh( 1 / ep ) / n );

% transform 3dB cut-off frequency to discrete-time domain
wc = 2 * atan( tan( .5 * ws ) * Wc );       % bilinear transform

I designed a discrete-time $6^{th}$-order Chebyshev II high pass filter according to above specifications. The resulting $3\textrm{ dB}$ cut-off frequency obtained from above script is $\omega_{c}=0.69\pi$. This is the correct value as shown in the plot below:


Marginally stable means that an otherwise stable system has one or more simple poles on the unit circle (in discrete time), or on the imaginary axis (in continuous time).

The consequence of that is that transients don't decay, but they also don't grow without bounds. Marginally stable systems are unstable in the bounded-input bounded-output (BIBO) sense.


Role of the rank of the filter mask matrix in image processing?

Separable just means you can do it in the x-direction and then in the y-direction and have it come out the same as if you did it in both dimensions simultaneously to begin with. It's not too hard to see that this will work for an average filter. If the filter is averaging over a 3x3 grid then in the 2-d case you take an average of nine values. In the separable filter case you first take three averages of three values. Then you average those three averages together. In both cases you get the same answer.

The separable case is much faster because you get to reuse some of the work you did in the x-dimension when you are doing the y-direction. In other words, each average of three values you computed in the x-direction will be used multiple times when filtering in the y-direction. The filters that can be made separable are precisely those whose matrix rank is one.


This's a too broad question including continuous time systems, discrete time systems, control systems, convolution operation, difference operation and differentiaiton operations.

However your main focus is on convolution. Therefore I will state in plain simple terms that the operation of convolution defines and calculates the output of a linear time invariant system to a given input from the impulse response of the system. This is performed by the convolution integral for continuous time systems and convolution sum for discrete time systems. Therefore convolution forms the core of signal processing. Just as integration or differentiation forms the core of calculus.

Control theory is an application of extensive signal processing, feedback theory, linear algebra, calculus, matrix theory, transform theory and probability theory. Therefore it will naturaly utilize convolutions in its many formulations. Nevertheless convolution is nothing more than a tool for the control field point of view.

A physical process requires an analog (continuous time) system to control it. But modern technology enables digital computers (sampled data systems) to be used instead. As a consequence, integrations and differentiations of the continuous time system are replaced by sums and differences in digital software implementation of analog systems. Whereas the link between the digital and continuous time systems can be defined via sampling theorem or some other mapping relations.


Does one ever store audio in mp3s with one channel being inverted?

Well, Audacity has an option to do that on digital files, so I assume that means it is done.

From that page:

  Invert flips the audio samples upside-down, reversing their polarity. The positive samples (above the horizontal zero line in the Audacity Waveform) are moved below the zero line (so becoming negative), and negative samples are made positive.
  Invert does not usually affect the sound of the audio at all, but it can be used for audio cancellation. If Invert is applied to one track and that track is mixed with another uninverted track that has identical audio, the identical audio is cancelled out (silenced).


How is data compression possible for audio?

there are at least two different meanings of "compressed audio", one is, as your title mentions, compression of data and the other is level compression.  they're two totally different functions.

about compression of audio data, there are two general classes: lossless compression and lossy compression.  in a crude sense of the word, signal quantization is a form of lossy compression, but i would not call word-width reduction the same topic as compression.  so let's assume the audio data has already had its word-width (or "bit depth") reduced to its production length (say, 16 or 24 bits or whatever).

as the names might suggest, losslessly compressed audio data will have the audio data restored exactly when uncompressed.  the original waveform will be restored.  for lossy compression, when the audio is uncompressed the original waveform is not restored, but the expectation is that the audio will sound the same.

to do lossless compression of audio, there are two main steps, the first is what is commonly called "Differential Coding".  "Delta Modulation" is the simplest example, but one can employ Linear Predictive Coding (LPC) to combine the past N samples in such a way as to predict the next incoming sample.  if the prediction is good, the difference between the actual sample and the predicted sample will be small and that reduces the number of values to encode significantly.  say, for 16-bit, instead of having 65536 values to encode, your difference signal might be not much bigger than 128 and you'll need only 8 bits to encode that difference.

what about the rare times you need more than 8 bits because the difference is greater than 128?  this is where entropy coding or Huffman coding comes in.  this is essentially what file compression programs like pkzip do.  words that appear often are assigned fewer bits than words that appear rarely.  this all comes out of Shannon's Information Theory.

that's, in a short form, what lossless compression is about.  to learn what lossy compression is about, there is a whole additional science about perception of audio, what audio features can be heard and what are not heard.  it's hard to explain, but there are ways of separating what can be heard from what cannot and for the separated data that is deemed not audible, no bits are allocated to it.  it's a crude explanation.


Filters that affect the real and imaginary parts equally, and thus
  have no effect on the phase, are appropriately called zero-phase-shift
  filters."

The filter itself is zero-phase, but the input (and therefore the output) waveform can be complex. In this case the real and imaginary components of the input would be effected equally since the filter is not adding any additional phase shift.


The formula you are looking for is

  (desired freq / Sampling Frequency) * Length of samples = sample number

you can see this works out even by verifying the units

  (Hz / Hz )* sample = sample

depending on your exact frequencies, the result may or may not be an integer number. If it is not an integer, that means the exact frequency wasn't captured (due to Fs). For example if your result was 23.8 you would have to look at either sample 23 or sample 24 to get the data you required, or do some sort of interpolation.

Good luck

edit

as your comment specified, this doesn't work exactly, the reason is becasue matlab does 1 indexing, not 0 indexing. So our DC component 0Hz is at element 1, not element 0 but the above formulas are correct. However implementing in matlab is a little different

desired_freq = 120; %Hz
desired_index = (desired_freq / Fs) * L + 1; %+1 is because matlab uses 1 indexing not 0 indexing

again desired index may not be an integer number so you can't simply say mag(desired_index) one workaround would be mag(round(desired_index))


First of all very similar questions have been asked here and even I have answered some of them, but now I will provide an answer using yet another approach.

Our aim is to show the effect of down-sampling on the DFT.
Let $x[n]$ be a sequence of length $N$. And let's define $y[n]$ by downsampling $x[n]$ by an integer factor $M$ as shown:

$$ x[n] \longrightarrow \boxed{ \downarrow M } \longrightarrow  y[n] = x[Mn]  \tag{1}$$

$y[n]$ is of length $K = N/M$. At the moment assume $K$ is integer.

I will make use of the relation between DTFT $X(\omega)$ and DFT $X[k]$ for finite length sequence which states:

$$ X[k] = X(\omega)|_{w = \frac{2\pi}{N}k } = X( \frac{2\pi}{N}k  ) ~~~,~~~ k = 0,1,2,...,N-1 \tag{2}$$

From the theoretical DTFT definition, it can be shown that the DTFT $Y(\omega)$ of $y[n]$ is given as:

$$ Y(\omega) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{ \omega + 2\pi k}{M} ) ~~~,~~~ -\pi \leq \omega < \pi \tag{3}$$

And we define the $K$-point DFT of $Y[k]$ as:

$$ Y[k] =  Y(\omega)|_{\omega = \frac{2\pi}{K} k } = Y(\frac{2\pi}{K} k) ~~~,~~~ k = 0,1,...,K-1  \tag{4}$$

Note the range of DFT index $k$ for $Y[k]$. Since $y[n]$ is a $K$-point sequence we have defined a $K$-point DFT of it.

Finally plug Eq(4) into Eq(3)

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{ \frac{2\pi}{K} k + 2\pi m}{M} ) $$

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{KM} k + \frac{2\pi}{M} m ) $$

Now $KM = N$ and we replace $\frac{2\pi}{M}$ with $\frac{2\pi}{N}(N/M) = \frac{2\pi}{N}K $ to get

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{N} k + \frac{2\pi}{N} Km ) $$

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{N} ( k + Km) ) $$

And finally we recognize that $X( \frac{2\pi}{N} ( k + Km) )$ is actually the $N$-point DFT $X[r]$ of $x[n]$, evaluated at $r = k + Km$ and insert it into the equation:

$$ \boxed{ Y[k] = \frac{1}{M} \sum_{m=0}^{M-1} X[ k + Km] ~~~,~~~k = 0,1,...,K-1 } \tag{5}$$

Note that DFT index k for $Y[k]$ is of modulo-K where as DFT index $r = k + Km$ for $X[r]$ is of modulo-N. Hence Eq(5) can alo be written like:

$$ \boxed{ Y[(k)_K] = \frac{1}{M} \sum_{m=0}^{M-1} X[ (k + Km)_N] } \tag{6}$$

The following is a MATLAB/ OCTAVE code to demonstrate the result:

N = 32;        % Length of x[n]
M = 4;         % Downsampling factor (integer)
K = N/M;       % length of y[n] (assuming N divisble by M)

x = randn(1,N);     % x[n]
y = x(1:M:end);     % y[n] = x[Mn]

X = fft(x,N);       % N-point DFT X[r] of x[n]

Y = zeros(1,length(y));  % K-point DFT Y[k] of y[n]
k = 0:length(Y)-1;       % DFT index range k

for m = 0:M-1            % implements Eq.6 above
    Y = Y +  X( mod(k+K*m,N)+1 );
end
Y = Y/M;

Y2 = fft(y,K);          % Comput Y[k] directly from y[n] for         checking.

figure,stem(abs(Y))     % display that they are the same..
hold on
stem(abs(Y2),'r+');


How is LMS/FXLMS noise cancelling different than simple polarity inversion?

Your polarity inversion method simply applies the same phase shift (delay) and gain to all frequencies, and is presumably not adaptive (or even closed loop). It may work reasonably well for some limited cases, where the phase and gain have been tailored for that particular environment, but any change to the environment will tend to reduce its effectiveness.

LMS however is an adaptive method which constantly "tunes" a filter to give you optimum cancellation using an appropriate gain and phase response which varies with frequency, and which adapts to a changing environment.


Why does OFDM use cyclic-prefix while QPSK uses root-raised cosine filters?

ISI, or intersymbol interference, means different things in the context of PSK and OFDM signals.  In PSK signals the symbols almost always have tails that extend, in the time-domain, into the times of other symbols.  This is what they mean by "intersymbol interference".

Unfortunately they have to do this to reduce the bandwidth of the signal.  They eliminate the negative effects of intersymbol interference by very carefully shaping the symbols so that when the receiver filters them through a matched filter, the other symbols have zero output at the sampling times of the other symbols.

With OFDM signals there is no overlap in the time-domain between symbols, so there is no intersymbol interference in the PSK sense of the term.  What they mean by ISI in the context of ISI is when time-delayed multi-path signals cause the time-delayed symbols to overlap with the following symbol in the non-delayed signal.  The cyclic prefix is added to the symbols to give the receiver enough un-corrupted data to recover the signal.

For instance, let's say that there are 64 sub-carriers and when we FFT them we get 64 samples of time-domain data.  Let's also say that we think the worst multi-path we will generally see will be a time delay of 7 samples.  If we make the cyclic-prefix 8 samples long then we should always have at least 65 "clean" samples, which gives us enough to recover the symbol.


Since what interests you is the "embedded system" part, and since you have a low budget (this excludes anything that requires proprietary compilers), I'd recommend building yourself a board with an ARM MCU and a codec, like this one. There's less than $50 of parts - the processor, the codec and the bare minimum to get them to work.

I'm recommending this because you can easily grab a FOSS toolchain for arm-none-eabi, because these chips have a serial bootloader which allows them to be programmed with a $2 FTDI dongle, and because they have enough flash and RAM onboard for many applications - so you don't have to deal with the complications of an external RAM or flash chip. For a beginner, it's good to have a small self-contained, system you can entirely master!

You'll learn how to set up a toolchain for an embedded processor, how to use a JTAG or serial interface to flash the chip, how to configure the peripherals and hardware for your application, and you'll write code in a relatively basic environment. It'll be a fun experience! The STM32F4 is good enough to implement audio effects and 6-10 voices of audio synthesis (Used in: the Owl stompbox, Audio Damage and Make Noise Eurorack modules, DIY synth projects like PreenFM2, Audiothingies P6, or Sonic Potions LXR). UAD needs a lot of computing power because they emulate analog hardware at the component level. You are not going to do that. To give you an idea of what's possible, the STM32F4 can run about 250 biquad filters or 100 wavetable oscillators or a phase vocoder with 75% overlap. This is a general purpose architecture, but there are DSP-ish things like saturation or MAC instructions that make the generated code relatively compact on audio DSP applications.

What this won't teach you are the specifics of VLIW architectures, and the art of writing assembly for such architectures. For that, you'd need a proper DSP board from AD or TI, but this costs more, and once you build something with their devboards it's more expensive and complicated to replicate it on your own board. You can partly capture this experience with a dsPIC.

Another cheap toy you could play with: ADAU1701 board. It has a built-in AD/DA. Downside: it can only be programmed through AD's drag and drop tool, with building blocks like filters, dynamic processors... Elektor also has a spin FV-1 board - it's a very peculiar minimalistic chip that can be mastered in a few hours. Might be worth exploring if you like vintage reverbs.


variance of the product of two samples with awgn

There is a paper, which is more general than Gaussian distributions.

  Bohrnstedt, George W., and Arthur S. Goldberger. "On the exact
  covariance of products of random variables." Journal of the American
  Statistical Association 64.328 (1969): 1439-1442.

https://www.jstor.org/stable/pdf/2286081.pdf?casa_token=m1_rSG4hiwIAAAAA:sGjisD6dekanCat5EBGl0yvdn1EJ5LY783fq4VKdxNhd0RA4z4izkcyDJ6yU7MylDke_qD335SWUNu1gJ5opTzXFWIVQWDrjoMvxcFgvMBEgG6r8En0W

Their equation (9)
$$
 V(xy) = E^2(x)V(y) + E^2(y)V(x) + V(x)V(y). \quad (9)
$$
reduces to your result


Time delay estimation of oscilloscope signals using cross correlation

@NickS

Since it is far from certain that the second signal in the plots is in fact a solely delayed version of the first, other methods besides the classical cross-correlation have to be attempted. This is because the cross-correlation (CC) is merely a maximum likelihood estimator if your signal(s) are delayed versions of each other. In this case, they are clearly not, to say nothing about the non-stationarity of them either.

In this case, I believe what may work is a time-estimation of the significant energy of the signals. Granted, 'significant' can or cant be somewhat subjective, but I believe that by looking at your signals from a statistical point of view, we will be able to quantify 'significant' and go from there.

To this end, I did the following:

STEP 1: Compute the signal envelopes:

This step is simple, as the absolute value of output of the Hilbert-Transform of each of your signals is computed. There are other methods to compute envelopes, but this is pretty straight forward. This method essentially computes the analytical form of your signal, in other words, the phasor representation. When you take the absolute value, you are destroying phase and only after energy.

Furthermore since we are pursuing a time-delay-estimate of the energy of your signals, this approach is warranted.

STEP 2: De-noise with edge-preserving non-linear Medial Filters:

This is an important step. The objective here is to smooth out your energy envelopes, but without destruction or smoothing out your edges and fast rise times. There is actually an entire field devoted to this, but for our purposes here, we can simply use an easy to implement non-linear Medial filter. (Median Filtering). This is a powerful technique because unlike mean filtering, medial filtering will not null out your edges, but at the same time 'smooth' out your signal without significant degradation of the important edges, since at no time is any arithmetic being performed on your signal (provided the window length is odd). For our case here, I selected a medial filter of window size 25 samples:

STEP 3: Remove Time: Construct Gaussian Kernel Density Estimation Functions:

What would occur if you looked at the above plot sideways instead of the normal way? Mathematically speaking, that means, what would you get if you projected every sample of our denoised signals onto the y-amplitude-axis? In doing this we will manage to remove time so to speak, and be able to study the signal statistics solely.

Intuitively what pops out of the figure above? While noise energy is low, it has the advantage that it is more 'popular'. In contrast, while the signal envelope that has energy is more energetic than the noise, it is fragmented across thresholds. What if we considered the 'popularity' as a measure of energy? This is what we will do with (my crude) implementation of a Kernel Density Function, (KDE), with a Gaussian Kernel.

To do this, every sample is taken and a gaussian function constructed using its value as a mean, and a pre-set bandwidth (variance) selected a-priori. Setting the variance of your gaussian is an important parameter, but you can set it based on noise statistics based on your application and typical signals. (I only have your 2 files to go off on). If we then construct the KDE Estimation, we get the following plot:

You can think of the KDE as a continuous form of a histogram so to speak, and the variance as your bin-width. However it has the advantage of guaranteeing a smooth PDF that we can then perform first and second derivitave calculus on. Now that we have the Gaussian KDEs, we can see where the noise samples peak in popularity. Remember that the x-axis here represents  the projections of our data onto the amplitude space. Thus, we can see which thresholds the noise is the most 'energetic' in, and those tell us which thresholds to avoid.

In the second plot, the first-derivative of the Gaussian KDEs is taken, and we pick the abscissa of the first sample after first-derivative after the peak of the mixture of Gaussians to attain a certain value close to zero. (Or first zero-crossing). We can use this method and be 'safe' because our KDE was constructed of smooth Gaussians of reasonable bandwidth, and the first derivative of this smooth and noise-less function was taken. (Typically first-derivatives can be problematic in anything but high SNR signals since they magnify noise).

The black lines show then at what thresholds we would be wise to 'segment' the image at, such that we avoid the entire noise floor. If we then apply to our original signals, we attain the following plots, with the black lines indicating the start of energy of our signals:

This thus yields a $\delta{t} = 241$ samples.

I hope this helped.


What base image size to choose to assist easier downsizing?

There is no objectively correct answer, there are too many variables to consider such as the nature of the image you are compressing (how much entropy/how 'compressible'?), how much you would need to reduce the jpeg quality to achieve the desired filesize etc... If it is possible to use a lossless compression format then I would consider that as a solution.

If you want to take an empirical approach you can take some representative images, go through the process of converting them using your two methods, then compare the results. You can use a tool like compare for a qualitative metric or visual comparison.


First, recognise that the fundamental period of the signal x(t) (sum of two sinusoidals with periods 2/3 s and 1/50 s) is 2s. And the fundamental frequency, $f_0$, is 0.5 hz.

Then considering the CT Fourier Series expansion of the periodic signal x(t) wrt its fundamental frequency $f_0$ given as $\sum_{-\infty}^{\infty}{a_k e^{jk2\pi f_0t}}$ you can  therefore find the coefficients corresponding to $2\cos(3\pi t)$ and $\sin(100\pi t + \pi/3)$ with k=3 and k=100 respectively.

Then, by expressing the sinusoidals in terms of complex exponentials using the Euler's identity this becomes $a_3=1 ,~ a_{-3}= 1 , ~~a_{100} = e^{j\pi/3}/2j, ~~a_{-100} = -e^{-j\pi/3}/2j $  which can be further manipulated if you want a rectengular representation of the complex coefficient instead of the above mixed-polar one.

Assuming no errors made during the arithmetic.


Why is distorted / quantized white noise still white noise?

White noise is defined as having a flat spectrum over the entire frequency range. Using uncorellated random numbers is just one of many methods to generate this, but by far not the only or the best one.
Also, whiteness implies no gaussian or whatever distribution, so your first sentence is wrong twice.


What factors should I consider in choosing an edge detection algorithm?

There are lots of edge detection possibilities, but the 3 examples you mention happen to fall in 3 distinct categories.

Sobel

This approximates a first order derivative. Gives extrema at the gradient positions, 0 where no gradient is present. In 1D, it is = $\left[ \begin{array}{ccc} -1 & 0 & 1 \end{array} \right]$

smooth edge => local minimum or maximum, depending on the signal going up or down.
1 pixel line => 0 at the line itself, with local extrema (of different sign) right next to it. In 1D, it is = $\left[ \begin{array}{ccc} 1 & -2 & 1 \end{array} \right]$

There are other alternatives to Sobel, which have +/- the same characteristics. On the Roberts Cross page on wikipedia you can find a comparison of a few of them.

Laplace

This approximates a second order derivative. Gives 0 at the gradient positions and also 0 where no gradient is present. It gives extrema where a (longer) gradient starts or stops.

smooth edge => 0 along the edge, local extrema at the start/stop of the edge.
1 pixel line => a "double" extremum at the line, with "normal" extrema with a different sign right next to it

The effect of these 2 on different types of edges can be best viewed visually:

Canny

This is not a simple operator, but is a multi-step approach, which uses Sobel as one of the steps. Where Sobel and Laplace give you a grayscale / floating point result, which you need to threshold yourself, the Canny algorithm has smart thresholding as one of its steps, so you just get a binary yes/no result. Also, on a smooth edge, you will likely find just 1 line somewhere in the middle of the gradient.


Extraction of the phase lag form cross-correlation function

The rotation of the constellation is due to carrier phase offset, which is not the same as a time delay. What the OP is seeing with the additional peaks is evidence that portions of the waveform itself repeat every 0.3x10^6 samples. Since the entire correlation is going down toward the edge of the plot, I assume the OP has done a linear cross-correlation and we are just seeing the effect from less overlap between the two waveforms that are being correlated.
Regardless of that, the carrier phase offset (constellation rotation) can be determined from the cross-correlation with zero time lag (the direct correlation function as a sum of complex conjugate products) as demonstrated below:
Consider such a correlation between the received $y[n]$ and known $x[n]$ signal where the only difference between the two signals is a phase rotation with
$x[n] = e^{j\theta[n] n}$
$y[n] = e^{j\theta[n] n}e^{j\phi}$
Where $\theta[n]$ is the known modulation of the carrier, while $\phi$ is an unknown phase offset (rotation).
The generalized correlation $r_{xy}[n]$ is a sum of complex conjugate products (the cross correlation function re-computes this at all delay offsets, but delay is NOT phase):
$$r_{xy}[n] = \sum_{n=0}^{N-1}x[n]y^*[n]$$
$$ = \sum_{n=0}^{N-1}e^{j\theta[n] n}e^{j\phi}Ae^{-j\theta[n] n}$$
$$ = Ne^{j\theta}$$
As we see, the phase of $r_{xy}[n]$ is the carrier phase offset we are trying to determine. Further, we have the advantage of processing gain in the correlation to optimally (in the presence of white noise) increase the SNR and therefore our estimate of the phase. The process is to use the cross-correlation function to align the received and transmitted waveforms, and once aligned (as evidenced by the largest peak) we can continue to use the generalized correlation function at that offset to determine the phase offset. Since the OP's cross-correlation function is maximum at offset = 0, the transmit and receive waveforms are already aligned, so the generalized correlation function can be used directly to determine phase (which would match the phase of that peak correlation at zero, since these would be the same).
I demonstrate this in this graphic pasted below showing correlation versus phase and frequency offset, with the added contributions of noise.


How to reconstruct RF signal using IQ data

The RF signal $r(t)$ is obtained from the complex baseband (IQ) signal $s(t)=x(t)+jy(t)$ in the following way:

$$r(t)=\text{Re}\{s(t)e^{j\omega_0t}\}\tag{1}$$

where $\omega_0$ is the carrier frequency (in rad/s). There are two other equivalent representations of $(1)$, where I use $s(t)=x(t)+jy(t)=a(t)e^{j\phi(t)}$:

$$\begin{align}r(t)&=x(t)\cos(\omega_0t)-y(t)\sin(\omega_0t)\tag{2}\\
r(t)&=a(t)\cos(\omega_0t+\phi(t))\tag{3}\end{align}$$

Eq. $(2)$ shows that the I- and Q-components are both modulated by orthogonal carriers, whereas $(3)$ shows that $r(t)$ generally exhibits amplitude modulation as well as phase modulation.

The quote in your question suggests that in the respective text the signals are all represented as discrete-time signals, hence the decimation and interpolation stages.


How can I smoothly interpolate between 2 position?

When you talk about bezier curves, it sounds like you think about them from an "illustrator" point of view, which is not totally right when it comes to spline interpolation (most probably what you are looking for).

Splines are piecewise curves that pass thrhough points. Bezier curves are third degrees splines between two points, and a series of them can form what "looks like" a single curve in a drawing aplication, but are actually lots of individual curves tied together (it counts as if the intermediate knots are doubled).

Illustration concepts aside, I think you could choose any 1-D spline interpolator from some popular library. I use Python, so I'd go with Scipy, take a look:

http://docs.scipy.org/doc/scipy/reference/interpolate.html#d-splines


How moving part pixel intensity values of video frames becomes dominant compared to stationary part intensities in reconstructed frames?

Q1. As shown by Oppenheim's experiment, the phase spectrum contains most of the structural information about the image. In 2D this are things like lines and edges. In 3D it is things like lines and edges but also movement. Instead of 2D frames and time, imagine the video as a 3D solid where the z-axis is the frame number. If you took a slice along the z-axis (3rd dimension), movement appears like an edge in the signal.

When you reconstruct using just the phase spectrum, you give all the component sinusoids the same magnitude. This basically normalises the brightness everywhere in the image. Try it in MATLAB:

I = double(imread('YOURIMAGE.TIF'));
f = fft(I);
f = f./abs(f);
I = real(ifft(f));
imagesc(I);

In a 3D image, such as a video, this normalisation is also in the time axis. Hence some smaller amplitude moving parts becoming greater in amplitude (and the opposite for high amplitude parts).

Q2. One video frame in your case is the same as 1 pixel in the x or y directions. The phase change depends on the movement.

Q3. Again that depends on the movement. Slow movement will fall into the lower frequency of the spectrum, high movement into the higher frequency. The spectrum is a global response though. It can only tell you about the entire signal. If you want to do analysis locally (i.e. restricted in location and time) you have to use quadrature filters if you want phase information.


This is a tougher problem than it appears at the first glance. This may be confusing at first, but read carefully, and I hope that it will give a good idea what to do. Here is what is my experience and how I resolved a similar problem:

apply a very good band-pass filter to get the frequency range of interest.
apply a very tight DFT to get a good resolution in frequency domain. You may wish to pad the original signal with zeros if it is not long enough.
determine the bandwidth of the notch filter that you want to use. Be very careful with notch filters, since they have some inconvenient characteristics.
make the running window sum of the DFT results using the width of your notch filter.
set the threshold to decide what is a peak in the results.
the running windows sum should give you values where the peaks are.
apply the notch filter designed for the frequency in mind. If the frequencies are changing, and the filtering has to be done in a real time, adaptive filter has to be applied, but that is yet another question.


For symplicity, consider the SISO linear system
\begin{align*}
\dot x &= ax +bu\\
y &= cx
\end{align*}
with $x$,$u$ and $y$ taking values in $\mathbb R$. Assume that you want to stabilize the system by making $y$ to converge to some fixed set point $\bar y$. Then necessarely there must be a point $(\bar x,\bar u)$ such that
\begin{align*}
0 &= a\bar x+b\bar u,& \bar y=c\bar x
\end{align*}
since $\bar x$ must be an equilibrium and such that $y=\bar y$. Trivially one has $$\bar u = -\dfrac{a}{bc}\bar x$$ One can define the error variables as $\tilde x=x-\bar x$, $\tilde y=y-\bar y$ and $\tilde u=u-\bar u$ that yield the error system
\begin{align*} \dot{\tilde x} &= a\tilde x + b\tilde u,& \tilde y&=c \tilde x
\end{align*}
One can stabilize it by taking $\tilde u = k\tilde x$ such that $a+bk<0$. That in fact ensures that asymptotically $\tilde y=0$ which implies $y=\bar y$. The overall control law is
$$
u = \tilde u + \bar u = k(x-\bar x) -\dfrac{a}{bc} \bar x = kx + \dfrac{1}{c}\left(-k-\dfrac{a}{b}\right)\bar y
$$
Thus as you can se, simple proportional is not enough, you also need a feed-forward control action given by
$$
u_{ss} = \dfrac{1}{c}\left(-k-\dfrac{a}{b}\right)\bar y
$$
If you define an integrator system like $\dot\sigma = y-\bar y = C(x-\bar x)$ then also $\sigma$ will have a set point $\bar \sigma$, and with $\tilde \sigma=\sigma-\bar\sigma$, the extended error system becomes
\begin{align*}
\begin{bmatrix}\dot {\tilde x}\\\dot{\tilde\sigma}\end{bmatrix} &= \begin{bmatrix}a & 0\\c & 0\end{bmatrix}\begin{bmatrix}\tilde x\\\tilde{\sigma}\end{bmatrix} + \begin{bmatrix}b\\0\end{bmatrix}\tilde u, & \tilde y&=C\tilde x\end{align*}
Take $\tilde u=k_p \tilde x + k_i \tilde\sigma$ (proportional+integral), then the closed loop system reads
\begin{align*}
\begin{bmatrix}\dot {\tilde x}\\\dot{\tilde\sigma}\end{bmatrix} &= \begin{bmatrix}a+bk_p & bk_i\\c & 0\end{bmatrix}\begin{bmatrix}\tilde x\\\tilde \sigma\end{bmatrix}, & \tilde y&=C\tilde x\end{align*}
You can chose $k_p$ and $k_i$ to make the closed loop matrix to be Hurwitz, and thus to stabilize the error system and you don't need any feedforward action, in fact the control law is
$$
u = \tilde u + \bar u = k_p(x - \bar x) + k_i (\sigma -\bar \sigma)
$$
Then with $$ \bar \sigma = \dfrac{k_p}{k_i}\bar x$$ the control law to becomes
$$
u = k_p x + k_i \sigma \qquad\qquad (1)
$$
and that's all. Moreover without integral control you need a feed-forward action which depends on the system parameters. If you have uncertainties on the parameters for sure you wil have an error on your set-point. with the integral action your final control law (1) does not depend on the parameters! and as long as the closed-loop system is stable you have null error at the steady state. Oscillations depend on how you chose $k_p$ and $k_i$. In fact if you chose them such that the closed loop matrix has complex eigenvalues then you can have socillations. For MIMO system the reasoning extends trivially.


Capturing images on fast moving conveyor belt

Try reducing exposure time if your camera allows that. With the short exposure time the conveyor belt should be well-lit.

If that is not enough, a common solution is to use a flash which is synchronized to the camera. If a sync signal is not readily available from the camera and you can't trigger the camera with your own sync signal, you probably could find a pin on your prototype camera's circuit board that carries a signal you can sync the strobe light to.

Without synchronization, setting the strobe rate to about 2/3 of the frame rate may be a reliable solution that gives good frames separated by at most one bad frame. Whether that works depends on whether there is a rolling or a progressive shutter and on the exposure time. Other ratios closer to 1 might work better in case of a rolling shutter and/or a short exposure time.

About DSLR cameras, you should use a camera that uses an electronic shutter in the mode (still/video) you would use. Cameras do not normally overheat but a mechanical shutter is only specified for that many cycles.

Also line scanners are a possibility, but I have no idea if those are commonly used. Think of a flatbed scanner where instead of moving the scanner head the target is moving. Those never blur the image but require bright illumination along the line being imaged. That may be easier to arrange than illumination over a large area, but it's a more exotic solution than a normal area camera.

More of an anecdote: A synchronously rotating prism with parallel opposite sides, located in front of the camera, can approximately compensate for linear motion.


RoC of Z transform of signal consisting of 3 values

When we speak about convergenge in all the $Z$ plane, the points $z=0$ and $z=\infty$ are not considered.

Note that by time-shifting the signal, you include/exclude zeros and poles at these 2 special points.


Is it accurate to say in digital audio that, when a fader is down, then its value is "-$\infty$"?

If by "its" you mean the value of the fader, then yes - it's absolutely correct.
Fader defines the attenuation of the signal with respect to the reference level. The units are in logarithmic scale.

For example, no amplification (or in other words, multiplication of the signal by 1) corresponds to $0\; \mathrm{dB}$ since $10\log_{10} 1 = 0$.

If you want to completely attenuate your signal, then you have to multiply everything by $0$, hence $10\log_{10} 0 = -\infty$.


PSD - Correct Formula for removing Noise

What you need is a digital filter. Digital filters are means to shape or "filter" the frequency response of a signal. Specifically, you need something called a band pass filter which would leave the 60-80 frequency mostly alone, and would mostly get rid of everything else. Snoop around this site or Google them. Digital filter are DSP 101, so there are many examples and tutorials out there.


The skeleton is a connected set of ridges of the distance transform. It is not a local maximum. The image should be analyzed through derivatives. See Haralick's pioneering work on that:
http://haralick.org/journals/ridges_and_valleys.pdf
http://haralick.org/journals/topographic_primal_sketch.pdf

And some MATLAB source code is here:
http://www.cs.utah.edu/~tch/notes/matlab/tchlib/tpslab.m
(and here of course: http://www.cs.utah.edu/~tch/notes/matlab/analysis/)
The source code obtains a Topographic Primal Sketch as described by Haralick. Ridges are only a part of that.

Optionally, you could invert the depth map as in $D^{-1}=1/D$, where $D$ is the distance map and $D^{-1}$ denotes the inverse distance function. Then start from a point with the highest distance and compute the shortest distance map using $D^{-1}$ as weights of the undirected graph. You should then select the points on the boundary to compute the shortest path from/to.


IQ Sampling and Oversampling Differences? *Confused*

I think you're confusing IQ sampling with oversampling.  Your first scenario is oversampling, or taking twice the number of samples on a real signal.  Other than stating the sampling takes place 90 degrees apart (it would actually be 180 degrees), your analysis is sound.

In the second section, you describe IQ sampling.  This is usually done by mixing a signal with a nonzero center frequency to baseband using two mixers, which multiply the same signal by a sine and a cosine, respectively.  The two mixing signals, separated by 90 degrees, are generally output from the same oscillator.  The outputs of the two mixers are separately sampled to get the I and Q components of the sampled stream.  This is typically represented as a stream at the original sample frequency but with two values (I and Q) for each sample.

IQ sampling essentially allows one to differentiate between positive and negative frequencies in the baseband signal, so it does in fact double the useful bandwidth of the signal just as oversampling does, but in a different way.  Instead of doubling the Nyquist rate, it revitalizes the band from -Fn to 0.


the sub-range of circular and linear convolution

Given a finite length discrete-time sequence $x[n]$ nonzero for the range $ 0 \leq n < L_x$, and another finite length sequence $y[n]$, nonzero for the range $0 \leq n < L_y$. Their linear convolution will be
$$ z[n] = x[n] \star y[n] $$

which is nonzero for the range $0 \leq n < L_x + L_y - 1$. In other words, the convolution will have a length of $$L_z = L_x + L_y - 1$$

The $N$-point circular convolution of $x[n]$ and $y[n]$ (usually implemented through $N$-point DFTs) will be $w[n]$

$$w[n]   = x[n] \circledast_N y[n] $$

Since most typically the circular convolution is used to implement the linear convolution, then we are interested in the relation between $z[n]$ and $w[n]$; i.e., what's the range of $n$ for which they are the same ?

The answer depends on a relation between $L_z$ and $N$ and stated as follows:

if $~~ L_z \leq N ~~~ $ then $w[n] = \begin{cases} { z[n]  ~~~,~~~ 0 \leq n < L_z \\ ~0~ ~~~~~,~~~~ L_z \leq n < N }\end{cases} $

$$\\$$

if $~~ N < L_z ~~~ $ then $w[n] = \begin{cases} { \text{aliased}  ~~~,~~~ 0 \leq n < L_z-N \\ z[n]~ ~~~,~~~~ L_z-N \leq n < N }\end{cases}$

Note that as a consequence of the second case, if $L_z - N > N$ or $ 2N \leq L_z$ there will be no clean samples.

The following Matlab stem-plot shows the overlap between a linear and a circular convolution of two sequences of lengths $L_x = 32$ , $L_y=10$ and DFT length $N=25$ and the linear convolution of length $L_z = 41$.

Since $N < L_z$, the first $L_z-N = 16$ samples of the circular convolution in the range $0 \leq n < 16$ will be aliased and only the remaining $N-(L_z-N) = 2N-L_z = 9$ samples in the range $ 16 \leq n < 25$ will be clean; i.e. equal to the linear convolution samples.

Note that the initial $16$ aliased samples of the circular convolution for $0 \leq n < 16$ is set to zero for clarification, and plotted them on together with the linear convolution which has a length of $41$ samples. The circular convolution just has $N=25$ samples only. Hence the last plotted $16$ zeros of the circular convolution actually belong to the first $16$ samples of the next period of the periodic result of the circular convolution.


Hough Transform not working to recognize a line

What you want is just to find the function

$$y(x) = mx + q$$

that fits through the centroids you've found.

You wouldn't do a Hough transform on that (since Hough transforms generally perform badly in presence of only few representatives of the full line, and are pretty computationally intense), but simply look for something like the least squared error fit of a $y(x)$ through the $(x_i, y_i)$ tuples given by the coordinates of your centroids.


How to detect whether the object has triangular shape?

I think you can find the orientation of the arrow using principal component analysis (or fitting an ellipse and taking major/minor axis directions).

To locate the arrow head, use distance transform. You'll have to preprocess the image, like flood-filling the region. Because the arrow head base is usually the widest part of the image, you should get a peak there.

Now you can use this point along with the two orthogonal vectors calculated earlier. In the image below I've translated those vectors to go through the distance-transform peak(I didn't do any calculations, the image is just for illustration).

The blue axis (second component in PCA or the minor axis of ellipse) now divides the image into two halves, and if you take projection of each half on to the blue axis, the half that contains the arrow head should give you a triangular profile (blue and orange profiles in the image below).

EDIT

Another approach would be to use the distance transform peak to isolate the neighborhood of the arrow head. If we threshold the distance transformed image by some fraction of the peak value, for the arrow head we should still get a triangular shape. Now we can prepare a mask centered on the peak to extract this region. Below is a simple Matlab code that does this for a specific image:

im = imread('arrow.jpg');
gray = rgb2gray(im);
bw = im2bw(gray, graythresh(gray));
di = bwdist(bw);
[mx, mxidx] = max(di(:));
[ypeak, xpeak] = ind2sub(size(di), mxidx);
bw2 = di > mx*.5;
se = strel('disk', double(ceil(mx*1.5)), 0);
mask = zeros(size(bw2));
mask(ypeak, xpeak) = 1;
mask = imdilate(mask, se);
figure,
subplot(1, 3, 1), imshow(bw), title('binary')
subplot(1, 3, 2), imshow(bw2.*di, []), title('thresholded dist')
subplot(1, 3, 3), imshow(mask.*di.*bw2, []), title('masked')

And the result is:


Determining Causality and Time-Invariance of a system

Hint
Substitute $s = t-1$. That gets you the equations in a more standard form $y(s) = ...$ Then go through the same excercise.


Does the Quality of the Pictures Taken by a USB 3.0 Camera Depend on Cable Length?

No, the data in the USB cable is digital.
Namely there is either an error in the data or the data is the same at any place along the cable.