Sequential Form of the Least Squares Estimator for Linear Least Squares Model I'm currently working on a project in which I need to find the tilt of a surface. Let's assume I'm only concerned with a single dimension tilt (i.e. slope) to begin.

I currently have the ability to calculate parameters of Affine Functions as described in Line of Best Fit (Least Square Method). However, this requires me to batch all the data prior to performing the calculation.

Is there instead a sequential / iterative form for computing the least squares linear fit such that I'm not required to batch the data prior to performing the calculation? Instead I would like to continuously update the least squares slope for each new data point that is received. The motivation to this is that I could avoid saturation / overflows due to batching data.

If you still don't understand what I'm asking, refer to this web page here which outlines a recursive form for calculating the mean - Heiko Hoffmann - Unsupervised Learning of Visuomotor Associations - PhD Thesis - Iterative Mean. The advantage of using this method is that you don't need to keep any large sums going thus avoiding potential overflows, and instead, you simply weight the new data against the previously calculated mean.

Slope from all samples obtained

To summarize the question's problem, you want to calculate the slope based on all samples obtained thus far, and as new samples are obtained, update the slope without going through all the samples again.

On the page you cite is the equation for calculation of the slope $m_n$ that together with $b_n$ minimizes the sum of squares of errors of the first $n$ of the $y$-values approximated by the linear function $f_n(x) = b_n + m_n x$. Converting the equation to zero-based array indexing for $x$ and $y$, we get:

$$m_n = \frac{\displaystyle\sum_{i=0}^{n-1}(x_i-\overline X)(y_i-\overline Y)}{\displaystyle\sum_{i=0}^{n-1}(x_i-\overline X)^2},\quad n \ge 2,\tag{1}$$

where $\overline X$ is the mean of the $x$ values and $\overline Y$ is the mean of the $y$ values. The problem with Eq. 1 is that the sums are actually nested sums:

$$m_n = \frac{\displaystyle\sum_{i=0}^{n-1}\left(x_i-\frac{\sum_{j=0}^{n-1} x_j}{n}\right)\left(y_i-\frac{\sum_{j=0}^{n-1} y_j}{n}\right)}{\displaystyle\sum_{i=0}^{n-1}\left(x_i-\frac{\sum_{j=0}^{n-1} x_j}{n}\right)^2}, \quad n \ge 2.\tag{2}$$

As $n$ is incremented, the inner sums can be updated recursively (iteratively), but by direct inspection the outer sum would need to be recalculated from the beginning. Least squares is an old and well-studied problem, so we don't try to bang our heads but look elsewhere. Converted to zero-based array indexing, Eqs. 16–21 of the Wolfram page on least squares fitting state:

$$\sum_{i=0}^{n-1}(x_i-\overline X)(y_i-\overline Y) = \left(\sum_{i=0}^{n-1} x_iy_i\right) - n\overline X\overline Y,\tag{3}$$

$$\sum_{i=0}^{n-1}(x_i-\overline X)^2 = \left(\sum_{i=0}^{n-1} x_i^2\right) - n\overline X^2.\tag{4}$$

These allow to rewrite Eq. 1 as (also given as Eq. 15 on the Wolfram page):

$$m_n = \frac{\left(\displaystyle\sum_{i=0}^{n-1} x_iy_i\right) - n\overline X\overline Y}{\left(\displaystyle\sum_{i=0}^{n-1} x_i^2\right) - n\overline X^2}, \quad n \ge 2.\tag{5}$$

It is easy to update the sums as $n$ is increased, and also updating $\overline X$ and $\overline Y$ is easy, by keeping and updating an accumulator containing the sum of all $x_i$ and another containing the sum of all $y_i$. The accumulated sum is then divided by the current $n$. The factor $n$ in the numerator and the denominator of Eq. 5 cancels together with one such division each, giving a formula for you to use:

$$m_n = \frac{\left(\displaystyle\sum_{i=0}^{n-1} x_iy_i\right) - \left(\displaystyle\sum_{i=0}^{n-1} x_i\right)\left(\displaystyle\sum_{i=0}^{n-1} y_i\right)/n}{\left(\displaystyle\sum_{i=0}^{n-1} x_i^2\right) - \left(\displaystyle\sum_{i=0}^{n-1} x_i\right)^2/n}, \quad n \ge 2.\tag{6}$$

Uniform sampling

If the samples are uniformly distributed, then further optimizations are possible. If the $x$ values are successive ascending integers, for example:

$$x_i = i,\tag{7}$$

then we get:

$$m_n = \frac{12\left(\displaystyle\sum_{i=0}^{n-1} i y_i\right) - 6(n - 1)\displaystyle\sum_{i=0}^{n-1} y_i}{n^3 - n}, \quad n \ge 2.\tag{8}$$

Avoiding large sums

Based on your comment to the question and to @M529's answer, perhaps this is not enough and you want to only ever store normalized sums that do not grow indefinitely with increasing $n$, similarly to what was done here. We can rewrite Eq. 8 as:

$$\begin{gather}m_n = 12 A_n - 6 B_n,\\
A_n = \frac{\displaystyle\sum_{i=0}^{n-1} i y_i}{n^3 - n},\quad
B_n = \frac{\displaystyle\sum_{i=0}^{n-1} y_i}{n^2 - n},\quad
n \ge 2,\end{gather}\tag{9}$$

which has recursions for $n \ge 3$:

$$\begin{gather}A_n = A_{n-1} + \frac{y_i}{n^2 + n} - \frac{3 A_{n-1}}{n + 1},\\
B_n = B_{n-1} + \frac{y_i}{n^2 - n} - \frac{2 B_{n-1}}{n},\end{gather}\tag{10}$$

that were found by simplifying the boxed parts of:

$$\begin{gather}A_n = A_{n-1} + \boxed{\frac{A_{n-1}\left((n - 1)^3 - (n - 1)\right) + (n - 1)y_{n-1}}{n^3 - n} - A_{n-1}},\\
B_n = B_{n-1} + \boxed{\frac{B_{n-1}\left((n - 1)^2 - (n - 1)\right) + y_{n-1}}{n^2 - n} - B_{n-1}},\end{gather}\tag{11}$$

where:

$$\begin{gather}A_{n-1}\left((n - 1)^3 - (n - 1)\right) = \sum_{i=0}^{n-2} i y_i,\\
B_{n-1}\left((n - 1)^2 - (n - 1)\right) = \sum_{i=0}^{n-2} y_i.\end{gather}\tag{12}$$

The recursive (iterative) algorithm can be tested in Python against the direct algorithm by:

from random import random
N = 10000  # total number of samples
n = 1
y = 2*random()-1  # Input y_0
sum1 = 0  # Direct (also needed to init recursive)
sum2 = y  # Direct (also needed to init recursive)
print("n="+str(n))
n = 2
y = 2*random()-1  # Input y_1
sum1 += y  # Direct (also needed to init recursive)
sum2 += y  # Direct (also needed to init recursive)
A = sum1/(n**3 - n)  # Recursive A_2
B = sum2/(n**2 - n)  # Recursive B_2
m = m_ref = 12*A - 6*B  # Recursive and direct m_2
print("n="+str(n)+", m=m_ref="+str(m_ref))
for n in range(3, N + 1):
    y = 2*random()-1  # Input y_2 ... y_(N-1)
    A += y/(n**2 + n) - 3*A/(n + 1)  # Recursive A_3 ... A_N
    B += y/(n**2 - n) - 2*B/n        # Recursive B_3 ... B_N
    m = 12*A - 6*B;                  # Recursive m_3 ... m_N
    sum1 += (n-1)*y  # Direct (not needed for recursive)
    sum2 += y        # Direct (not needed for recursive)
    m_ref = 12*sum1/(n**3 - n) - 6*sum2/(n**2 - n)  # Direct m_3 ... m_N
    print("n="+str(n)+", A="+str(A)+", B="+str(B)+", sum1="+str(sum1)+", sum2="+str(sum2)+", m="+str(m)+", m_ref="+str(m_ref))

The recursive algorithm seems stable, with the state variables A and B remaining small with a total of N=10000 independent random samples from a uniform distribution $-1 \le y < 1$. At the last sample the state variables and the output are:

n     = 10000
A     = 6.246471489186524e-07
B     = 5.20978725548091e-07
sum1  = 624647.1426721811
sum2  = 52.09266276755386
m     = 4.369893433735283e-06
m_ref = 4.369893433735271e-06

where m is $m_n$ calculated by Eq. 9 using the recursions of Eq. 10 for $n \ge 3$ and m_ref is the direct calculation of $m_n$ by Eq. 9. The two ways of calculating $m_n$ agree to 12 decimal digits. The direct sums sum1 and sum2 seem to be growing without bound. Of course, n is also growing without bound, and I don't think there's anything that can be done about that if this is what you want to calculate.

I did some testing against an arbitrary precision implementation using Python's mpmath (data not shown), and found that the root mean square error of the direct and the recursive methods did not differ significantly for $N = 10^5$ random samples. Personally, I would prefer the direct summation method using a large (in number of bits) fixed-point accumulator, if necessary. With direct summation, especially if it is not necessary to compute $m_n$ for all $n$, many slow division instructions are avoided, compared to the recursive method.

Recursion without auxiliary state variables

For what it's worth, I also found a recursion that doesn't require auxiliary state variables (such as $A$ and $B$), but requires remembering the previous two $m$ and the previous $y$:

$$m_n = \frac{2m_{n-1}(n - 2)}{n + 1} - \frac{m_{n-2}(n - 2)(n - 3)}{n^2 + n} + \frac{6y_{n-1}(n - 3)}{n^3 - n} - \frac{6y_{n-2}(n - 3)}{n^3-n}.\tag{13}$$

Slope from the most recent $N$ samples

Earlier I had thought that you want to calculate the slope based on $N$ most recent samples. I will keep the following stub of an analysis about that, in case it might be useful for someone.

For a least squares fit polynomial $f(x) = c_0 + c_1 x$ with $x=0$ at the center of the length-$N$ window, you can calculate the "intercept" $f(0) = c_0$ using a Savitzky–Golay filter, and the "slope" $f'(0) = c_1$ using a Savitzky–Golay differentiation filter, with both filters parameterized to have a polynomial degree of 1. These are linear time-invariant filters each described by an impulse response, easy to obtain for example using the function sgolay in Octave:

pkg load signal;
N=9;  # Number of samples in the window
sg0 = sgolay(1, N, 0)((N+1)/2,:);  # Savitzky-Golay filter
sg1 = sgolay(1, N, 1)((N+1)/2,:);  # Savitzky-Golay differentiation filter

The filters obtained have the following coefficients (which when reversed give the impulse response):

sg0 = 0.11111   0.11111   0.11111   0.11111   0.11111   0.11111   0.11111   0.11111   0.11111
sg1 = -6.6667e-02  -5.0000e-02  -3.3333e-02  -1.6667e-02   3.5832e-18   1.6667e-02   3.3333e-02   5.0000e-02   6.6667e-02

Figure 1. Savitzky–Golay (blue crosses) and Savitzky–Golay differentiation filter (red o's) coefficients, with $N=9$ and polynomial order 1. Swap the sign of the horizontal axis for impulse responses.

The first filter is a moving average. The second filter has a truncated polynomial impulse response. I think such filters can be implemented recursively (see my question about that), perhaps elaborated in this pay-walled paper: T.G. Campbell,  Y. Neuvo: "Predictive FIR filters with low computational complexity" IEEE Transactions on Circuits and Systems (Volume: 38, Issue: 9, Sep 1991).


Doppler effect on Baud Rate I'm having a hard time wrapping my head around doppler and how it affects a digital communications signal's baud rate.  The examples that I've seen explain that the center frequency of a signal is shifted when either the receiver and/or transmitter are moving relative to each other.  What I don't understand is the why doppler also affects baud rate?  Can someone point me to some tutorials on this topic or explain?  What I'd like to understand clearly is:

Intuitively understand why this happens?
Mathematically understand why this happens?
I'd like to be able to calculate when this is something that a communications receiver has to take into account?  My guess is that based on the doppler shift that will be experienced by the receiver, in combination with the signal's baud rate, but I guess I need some hints or help to get me along with this.

This is a follow-up to this question already posted: How to estimate and compensate for doppler shift in wireless signals?

I'm sure that you are familiar with the idea of how Doppler affects a sine wave.  When the transmitter and/or receiver move towards the other the wave arrives more quickly, which makes the wave appear to be a higher frequency.  In regards to the baud rate, it also arrives more quickly.  If you think about it you will see that this must be so.  After all, the data stream cannot be separated from the carrier wave that modulates it!

So, you have a signal that is exactly identical to the "non-Doppler" signal except it arrives faster or slower than the non-Doppler signal.  This can be modeled exactly, in the discrete time domain, through resampling.

So why is Doppler usually modeled as a simple frequency shift when it should really be modeled as resampling?  Two reasons- first, the frequency shift effect is almost always vastly larger than the baud rate change.  The reason for this is that the Doppler effect is proportional to the frequency of the original wave ($\Delta f = \frac{\Delta v}{c}f_0$).  The frequency of the carrier wave is vastly greater than the baud rate, so the effect is much larger on the carrier wave.

The second reason most people don't model it as resampling is that it is difficult to do very small arbitrary sample rate changes, like 3013/3014.  If either the numerator or denominator is a prime number or has really large factors then it can become impractical to do traditional fractional resampling.  In that case polynomial interpolation is probably the best way to go.


How to extend hough lines over the image space I need to extend lines over the image space. I used Hough transform to finde vanishing lines, but i need to get lines longer and see where is vanishing points.
How can i do it? Please help.
This is my code:

clear; clc;

I = imread('22pp.jpg');

%rotI = imrotate(I,33,'crop');
bw_I = rgb2gray(I);
BW = edge(bw_I,'canny');

figure, imshow(BW); %pause;

BW=edge(BW,'canny',(graythresh(I)*(0.3)),'horizontal');

figure, imshow(BW); %pause;

[H,T,R] = hough(BW);
imshow(H,[],'XData',T,'YData',R,'InitialMagnification','fit');
xlabel('\theta'), ylabel('\rho');
axis on, axis normal, hold on;
P  = houghpeaks(H,15,'threshold',ceil(0.3*max(H(:))));
x = T(P(:,2)); y = R(P(:,1));
plot(x,y,'s','color','white');
% Find lines and plot them
lines = houghlines(BW,T,R,P,'FillGap',200,'MinLength',200);
figure, imshow(I), hold on

for k = 1:length(lines)
   xy = [lines(k).point1; lines(k).point2];
   plot( xy(:,1) ,  xy(:,2),  'LineWidth',2,'Color','green');

   % Plot beginnings and ends of lines
   plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
   plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
end

Fit a line through the end points using standard line model $ax+by+c=0$. Then draw them on the image from -10000, 10000. The points remaining on the image will give you what you are asking for.


What is the reason behind complex conjugate pairs in Linear Phase FIR filter analysis from the Pole Zero plot When we say that a Generalized Linear Phase System must satisfy the pole zero plot with the condition that a complex zero not on the unit circle exist's in a pair of 4. Then I understand that I need the conjugate pairs to have real coefficients. Why do we need the complex conjugate reciprocal pairs and what role do they play for the FIR filter? I am aware that we use complex conjugate reciprocal pairs in case of all pass systems but how does that relate to the Linear Phase System if it even does.

The impulse response of a linear phase FIR filter is odd or even symmetrical. This has as a consequence that if $z_0$ is a zero, $1/z_0^*$ must also be a zero. If in addition to the linear phase property you also have real-valued coefficients, complex zeros always occur in complex conjugate pairs. So if a real-valued linear phase filter has a complex root $z_0$, you automatically get $z_0^*$, $1/z_0$, and $1/z_0^*$ as roots.
Complex roots on the unit circle just occur in complex conjugate pairs (because $z_0=1/z_0^*$, and real-valued roots just occur in pairs $z_0$ and $1/z_0$, unless they're on the circle (i.e., $z_0=1$ or $z_0=-1$).
In sum, for real-valued linear phase FIR filters, there are four types of "clusterings" of roots:

complex roots with $|z_0|\neq 1$ occur in groups of four: $z_0$, $z_0^*$, $1/z_0$, $1/z_0^*$
complex roots on the unit circle occur as complex conjugate pairs $z_0$ and $z_0^*$
real roots with $|z_0|\neq 1$ occur in pairs as $z_0$ and $1/z_0$
real roots on the unit circle ($z_0=1$ or $z_0=-1$ ) don't imply another root

Mathematically, for a real-valued length $N$ FIR filter the linear phase condition in the frequency domain is
$$H(z)=\pm z^{-(N-1)}H\left(\frac{1}{z}\right)\tag{1}$$
where $N$ is the filter length (number of taps). From $(1)$ it is clear that if $z_0$ is a zero of $H(z)$ then also $1/z_0$ must be a zero.


Calibration (intrinsic) parameters of a virtual camera with some known information I need the calibration matrix of my simulated camera, I am trying to fully comprehend and follow the steps from Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers as well as Estimating Plane Pose Without Knowledge of Intrinsic Camera Parameters, but I still don't "get it". Ultimately, I want to use K to find the essential matrix and compute the transformation between two images.

The information available to me are:

The images.
Images are planar scenes, camera is perpendicular to the ground. Transformation between them includes rotation (around ground-perpendicular axes), scaling and translation.
Known distance from image (altitude), FOV, pixel width.

Assuming square pixels and no skew, this is what I am after:

$$
K = \left[\begin{matrix} \alpha_u&0&u_0\\0&\alpha_v&v_0\\0&0&1 \end{matrix}\right]
$$

$\alpha_u$ and $\alpha_v$  are the scale factor in the $u$ and $v$ coordinate directions, and are proportional to the focal length $f$ of the camera: $\alpha_u = k_u f$ and $\alpha_v = k_v f$. $k_u$ and $k_v$ are the number of pixels per unit distance in $u$ and $v$ directions.
$c=[u_0,v_0]^T$ is called the principal point, usually the coordinates of the image center.

Is there an "easy" approach to computing these variables, given the assumptions and the available info about the images?

Going through dsp, I am overwhelmed by the amount of links and information on the matter, usually very technical papers that are hard to follow for beginners in the field, and I am unable to understand what corresponds to my case.

Since your camera is simulated and the scenes are synthetic, you already have this information: Check out my updated answer to this question.


Noise cancellation by averaging multiple images together As the title suggests , I want to know how noise cancellation is achieved by adding multiple images together . I know the theory behind it , but simple adding the pixel values of 100 noisy images of the same sample would simply produce a white image . Can anybody elaborate , how it is normally done.

This works simply by virtue of the central limit theorem. For a given pixel, you are given 100 observations of its true value + random noise; and if you assume that the noise has zero mean, then the average of these will have a normal distribution centered at the true pixel value, with a variance decreasing as more observations are added. The central limit theorem also tells you that this is a rather terrible denoising method: for a gain in SNR of $6dB$ (reducing noise by a factor of 2) you have to increase the number of averaged samples by a factor of 4.

Regarding your problem, are you sure that you are averaging (dividing by the number of images) rather than just summing? And are you sure that you are using the correct data type for your computations? These are implementation/programming questions and are out of topic here.


Identifying impulse based on experimental data I have conducted an experiment in which the impulse was exerted on sound bowl (it was hit by a hammer). As a result I have obtained acceleration from accelerometer on 3 axis, calculated net acceleration ($\sqrt{x^2+y^2+z^2}$) and tried to obtain frequency components of the response.

What I need to do is to identify input impulse but I have no idea on how to do that.

Compute the position function from the accelerometer data using a scaled aggregator of some kind presumably, with initial conditions about resting position and velocity of the sound bowl. If you don't want to model the hammer's impact exactly, and want to think of it as an impulse, then the position function is the impulse response. Otherwise, go to 2.
Model the impact of a hammer as a non-ideal impulse function in time. Deconvolve the position function and the hammer model to get the impulse response. You could try to use the 3d fft on both these functions and find their ratio, and the compute its inverse 3d fft.


Centre channel derivation using correlation I would like to know the algorithm for  how centre channel is derived from rejecting left and right channel in audio upmixing process.
Can someone explain intuitively rather than hard hitting mathematical equations to describe the derivation process ...

Audio up-mixers are fairly complicated and the "better" ones are either protected IP or trade secrets.
To extract a center channel you would look at the correlation between the left and right input channels and dynamically steer the correlated content towards the center and the uncorrelated towards left and right. The devil is in the details though: it's tricky to do this without creating undesirable steering artifacts. There are lot of implementation issues tow work out: time constants, energy detectors details, correlation metrics, steering curves, frequency band splicing, etc.
Here is on old patent you can plow through if you are interested: https://www.freepatentsonline.com/7016501.pdf


Difference of Gaussian filter but using Gaussian CDFs not PDFs A very common filter in signal / image processing is the Laplacian of Gaussian filter, and its approximation, the difference of gaussians filter.

The DoG filter is the difference of 2 Gaussian PDFs that have slightly different standard deviations.

My question is: what would a filter that uses the difference of 2 Gaussian CDFs (instead of PDFs) look like and what properties would it have?

You need to ask yourself why do we use the difference of Gaussians from the first place? The reason is because the difference will give us a measurement for the change in value around the point we apply it to as a function of the variances of the Gaussians. If we have a big change the difference between the Gaussians it means that we have some frequencies that were low passed by the Gaussain with the bigger variance. If the change is low it means that there i no information lost. So we can think of it as a band-pass filter.

A simplified reason why it acts like that is because of the shape of the Gaussian PDF as can be seen in the following figure.

If we would use the CDF of the Gaussains and not the PDF it will give us a positive difference on one side and a negative difference on the other side.


Synchrosqueezed STFT phase transform Synchrosqueezed STFT paper defines STFT phase transform as:
$$
w(u, \xi) = \xi - \frac{\partial_tS_{g}f(u, \xi)}{j2\pi S_gf(u, \xi)}
$$
which differs from synchrosqueezing original paper's (per CWT) two-fold:

The $\xi\ - W$, where original's is simply $W$
The $2\pi$, absent in original's

Why the differences; what motivates STFT's definition? And why use 'modified' STFT?

TL;DR:

CWT's $\partial_t$ is taken with respect to wavelet timeshifts, STFT's with respect to window time samples, with CWT's FT taken along rows, and STFT's along columns.
The $2\pi$ is per assuming a different 'oscillatory factor' of Fourier Transform's kernel, i.e. $e^{-j2\pi f}$ vs. $e^{-j\omega}$. Numeric implementations can work with either.
"Why modified STFT" is a separate question, addressed in second answer.

Motivating phase transform
To derive CWT's phase transform, original paper uses $f(t) = A \cos(\Omega t)$ as a motivating example, or for our purposes equivalently $A e^{j\Omega t}$, then works backwards from the result of $\partial_t W_f / W_f$ to define $w$ such that $\Omega$ is extracted. The idea is, this $f(t)$ is the basis function (building block, atomic unit) of the Fourier Transform, so results here will transfer to  any other $h(t)$ the more $h(t)$'s (D)FT is "well-behaved" (more specifically, $h(t)$'s "components" obey synchrosqueezing assumptions detailed in followup paper), which spans a large variety of real-world signals.

Deriving phase transform
The SSQ STFT paper considers $e^{j\omega t}$ directly, but omits a derivation; I will derive it here for both the 'standard' and 'modified' STFT, along their derivatives, from which the phase transform will arise naturally.
Let $f(t) = e^{j\xi_0 t}$, and assume a real, even-symmetric ($g=g^*,\ \hat{g} = \hat{g^*}$), window $g(t)$. We have:
$$
\begin{align}
S_gf(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g(t-u) e^{-j\xi t} dt \\
& = \int_{-\infty}^{\infty} e^{j\xi_0t} g(t-u) e^{-j\xi t} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) \hat{g}(\omega) e^{-ju\omega} d\omega \\
& = \hat{g}(\xi - \xi_0) e^{-ju(\xi - \xi_0)} \\
\end{align}
$$
where Plancherel-Parseval's Theorem was applied, and similarly
$$
\begin{align}
S^{\text{mod}}_g f(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g(t-u) e^{-j\xi (t-u)} dt \\
& = \int_{-\infty}^{\infty} e^{j\xi_0t} g(t-u) e^{-j\xi (t-u)} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) \hat{g}(\omega) e^{-ju(\omega - \xi)}d\omega \\
& = \hat{g}(\xi - \xi_0) e^{ju\xi_0} \\
\end{align}
$$
The derivatives are:
$$
\begin{align}
S_{g'} f(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g'(t-u) e^{-j\xi t} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) j\omega\hat{g}(\omega) e^{-ju\omega} d\omega \\
& = j(\xi - \xi_0)\hat{g}(\xi - \xi_0) e^{-ju(\xi - \xi_0)} \\
\end{align}
$$
$$
\begin{align}
S_{g'}^{\text{mod}} f(u, \xi)
& = \int_{-\infty}^{\infty} x(t) g'(t-u) e^{-j\xi (t-u)} dt \\
& = \frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi \delta(\omega -(\xi - \xi_0)) j\omega\hat{g}(\omega) e^{-ju(\omega - \xi)}d\omega \\
& = j(\xi - \xi_0)\hat{g}(\xi - \xi_0) e^{ju\xi_0} \\
\end{align}
$$
Then, for both cases, to retrieve $\xi_0$, we consider the ratio
$$
\begin{align}
\frac{\partial/\partial t( S_{g}f(u, \xi) )}{S_gf(u, \xi)} = \frac{S_{g'}f(u, \xi)}{S_gf(u, \xi)} = j(\xi - \xi_0)
\end{align}
$$
which leads us to define
$$
\begin{align}
& \boxed{ w(u, \xi) = \xi - \frac{S_{g'}f(u, \xi)}{jS_gf(u, \xi)} } \\
&= \xi - \frac{j(\xi - \xi_0)}{j} = \xi_0 \\
\end{align}
$$
Note the phase transform being same for both is one consequence of the division "reducing influence of windowing".

Verifying derivation
$S_g f(u, \xi)$ expression can be verified with an example, comparing against Wolfram's output (set normalization=1, oscillatory factor=-1):
$$
\text{Let}\ f(t) = e^{jt},\ g(t) = e^{-\pi t^2}, \text{ and thus}\ \hat{g}(\xi) = e^{-\xi^2 / (4\pi)}.\ \text{Then}, \\
$$
$$
\begin{align}
S_gf(u, \xi)
&= \hat{g}(\xi - 1) e^{-ju(\xi - 1)} \\
&= \exp{\left( -\frac{1}{4\pi}(\xi - 1)^2 - ju(\xi - 1) \right)} \\
&= \exp{\left(\frac{1}{4\pi}( (\xi - 1)(-j4\pi u - \xi + 1 )\right)} \\
\end{align}
$$
More completely, ssqueezepy's ssq_stft on a variety of signals yields expected behavior.

Why the $2\pi$?
The TL;DR pretty much says it, which can be verified in the original MATLAB Toolbox's paper. In code, with frequency-domain differentiation, the $2\pi$ makes it so we define $j \omega$ to span $[0, ..., \pi]$ and $(-\pi, ..., 0)$, and without uses $0.5$ instead of $\pi$, reinterpreting the derivative's dimensionality as linear cyclic fraction of sampling frequency, rather than radian. See code.

Why the $\xi\ -$?
CWT's $\partial W_f$ takes $\partial / \partial b (\psi((t-b)/a))$, where $b = $ timeshift, whereas STFT's takes $\partial / \partial t (g(t))$. This is per how each is computed:

CWT takes DFT along rows; for every scale (row), we use frequency-domain convolution to correlate the wavelet with the signal at every timeshift $b$ (column).
STFT takes DFT along columns; for every timeshift (column), we apply DFT to correlate the windowed signal with cisoid bases at every frequency (row).

Following the derivations from motivating examples reveals the $\xi \ -$ stemming from this difference.
The definitions can be interchanged or made the same for both; why exactly they differ I'm unsure. My bet is, CWT is best computed for every timeshift (hop_length = 1), and the scales are already sparse, whereas in STFT we often apply large hop_length (small overlap), so convolution theorem is very compute-wasteful.


error variance of frequency and phase estimation using DFT What is the error variance of estimating the (frequency) and (initial phase) of a sinusoid using DFT.
please provide references if possible.

Since I don't have sufficient reputation to comment i will quote and add to Dsp_user's answer.

Summary: The error variance will depend how well you can estimate the frequency of your signal. This in turn will be limited by the amplitude of the DFT peak vs noise level. Spectral leakage will likely also skew your estimates.

  the latter can't be assigned a physical meaning (please someone correct me if I'm wrong on this point).

this is incorrect, calculated phase from the DFT is linearly correlated with phase shift of a sinusoid. When generating a cosinus curve and shift it between 0:2*pi i get errors in the range of (10^-15) between signal phase shift and phase estimate. Phase estimates from DFT can thus be said to compare a signals shift to a non-shifted signal. Simulated curve is shown bellow:

In my opinion, phase estimation error should only be dependent on how well the "sinusoidal-fit" from the DFT fits the actual signal. One of those limitations (frequency resolution/spectral leakage) was explained by Dsp_user. If your signal is between frequency bins you will estimate the wrong signal amplitude (correlated with how close the frequency is to the middle between two bins). Since phase estimation by DFT is linerly dependent on signal frequency i would guess that you will introduce a higher/lower derivative when trying to estimate real time shifts between signals but I admit here that I am not 100%.

Added: In case of spectral leakage, if you have other sinusoids close to your peak of interest, their leakage could pour onto that and introduce errors.

Your noise floor (if it is white) will in real cases likely be your largest source of error. It will cause additive errors to your phase estimates whose size will depend on DFT-Peak Amplitude vs noise floor amplitude (noise power spread out over the signal length). Bellow you can find a figure of this as well as the code to generate both in my answer

t=0:1000^-1:.999;
phi=0:1000^-1:2*pi-1000^-1;
signal = cos(2*pi*10*t+phi');

FFTSig= fft(signal,[],2);
figure; plot(phi,squeeze(atan2(imag(FFTSig(:,11)),real(FFTSig(:,11)))),'linewidth',2)
xlim([0 2*pi])
ylim([-pi pi])
xlabel('Signal Phase Shift')
ylabel('Phase shift estimate from DFT')
title('Phase Shift Estimate - No Noise')

signal2 = cos(2*pi*10*t+phi');
for iter=1:6283
   signal2(iter,:)= squeeze(signal2(iter,:))+2*randn(length(signal(1,:)),1)';
end

FFTSig2= fft(signal2,[],2);
figure; plot(phi,squeeze(atan2(imag(FFTSig2(:,11)),real(FFTSig2(:,11)))))
xlim([0 2*pi])
ylim([-pi pi])
xlabel('Signal Phase Shift')
ylabel('Phase shift estimate from DFT')
title('Phase Shift Estimate - Additive White Noise')


In Search of Voice Activity Detection Algorithms - What happens to the frequency or pitch or speech features as speaker moves away from microphone Actually i am in  search of a Voice Activity Detection Algorithm which could distinguish between voice and non-voice

Roughly speaking it must not detect even a bullet sound,even a foot stepping and other non speech activity should only detect people conversation  or any one shouting

in that search this question arises in my mind and i want know the effect of noise and distance between speaker and microphone on the speech features like pitch,frequency,cepstrum,zerocrossing rate,power spectral density,entropy  etc

if some component wont get distorted i would like to extract that feature and do the activity decision on that

Can any one help me in extracting dominant parameter of speech which would differentiate it from other common sounds even in Lower SNR conditions <0dB

Note:my algorithm  expects voice activity happens at a distance of at least 10m away from microphone and continuous generator hum as background noise

The intrinsic fundamental frequency (pitch) of a speech signal do not change when the speaker moves away from the microphone (unless it moves away quite fast in which case it's a case of Doppler effect!) - it would be strange if we heard people's voice transposed up or down as they move away or as ambient noise increases!

However, as the distance from the microphone increases, the signal to noise ratio will increase, and this might affect the reliability of a naive pitch detector, so you might indeed get different result when estimating the f0 of a speech signal and the same signal with more noise. It doesn't mean that the f0 is different, just that the noise makes the estimate unreliable.

I am not aware of any naive feature which allows noise-robust speech/non-speech discrimination with a simple threshold rule. I think supervised machine learning would be the way to go for your problem.


Compensate delay of non-symmetric (non-linear phase) filter Lets assume I want to apply a matched filter h which is non-symmetric to my signal x and the output is y in matlab:

y=filter(h,1,x);

Now I am interested in comparing the information between my filtered signal yand my input x. For this I want to compensate for the delay of the filter h.

When googling this topic I found the following tutorial which handles linear-phase fir filters. Obviously in such a case the group delay is $\frac{N-1}{2}$ given that the filter has the length $N$.

%this is the example from matlab
delay = mean(grpdelay(h,1,128));
x = x(1:end-delay);
y(1:delay) = [];

However, I am aware that the group delay of a non symmetric filter is frequency depending and I am not sure how to proceed to align my signals in time to compare them. As I am interested in real time application the function filtfilter is not applicable.

So my question is, how do I estimate in the best way the time delay to align (sync?) my signals in time given the coefficients of a non-symmetric filter h?

delay = .... h ....?
x = x(1:end-delay);
y(1:delay) = [];

The desired output should look like the red curve in the following figure

Here the signal in blue is the input signal x and in green the filtered signal y. However, I would like to to compensate for the delay of the green signal so that the input and output align in time as in the case of the red curve.

If your filters have a low pass characteristic (which at least the one in your example appears to have), then a possibly useful measure for the delay is the center of gravity of the impulse response:

$$\hat{\tau}=\frac{\displaystyle\sum_nnh[n]}{\displaystyle\sum_nh[n]}\tag{1}$$

For more general impulse responses, other than low pass, you can modify $(1)$ by either taking magnitudes or squares of $h[n]$. Note that $(1)$ is the filter's group delay evaluated at DC.


How to make the impulse response sparse? How does one know that the channel is sparse? I am new to sparse channel estimation algorithms and reading research articles. One such paper is blind sparse channel estimation using a modification of the BOMP technique titled, "Blind Acoustic Source Separation Via
System Identification for Leak Detection in
Pipelines."

How are the impulse response/coefficients made sparse?

Dankers, Jalilian, and Westwick have explained the way they have done parameterization under their algorithm on page 4 and (10). How does one do parameterization of impulse response to decide for sparse/non-sparse?

In general, how is parameterization applied to make the impulse response to zero and in general how does one know which coefficients are sparse?

How you parameterize your sparsity will depend on your application.  The authors of that paper, in a paragraph on page 231 say:

which is why they clump the coefficients together in $P$ blocks of start time $t_{B_k}$ of duration $n_{B_k}$.

For one impulse response this is shown in their figure 1(b).

The overall sparsity will depend on the reflections and delays that they have when this impulse response is used by the direct signal and all its indirect reflections.

  How are the impulse response/coefficients made sparse?

You need to decide which coefficients in the overall impulse response are zero. To do that, you need to decide why they are zero (what is the cause?). Perhaps there is a long dead-time between the onset and the first secondary response.

  How does one do parameterization of impulse response to decide for sparse/non-sparse?

Once you've decided on the form, you decide on the structure of the non-sparse parts of the response and how they are separated from each other.

Suppose your non-sparse parts are all defined by $\alpha_k h[n]$, where $\alpha_k$ is the gain of the $k^\mbox{th}$ non-sparse segment and $h[n]$ is the FIR response of all segments (modulo a gain term). $h[n]$ is of duration $H$.  Then the sparsity will come about by how you space out the $\alpha_k h[n]$:
$$
g[n] = \sum_{k=1}^K \alpha_k h[n-t_{B_k}]
$$
where $t_{B_k} \gg H$ for sparsity to be true (so there is lots of space between the $\alpha_k h[n]$.

  In general, how is parameterization applied to make the impulse response to zero and in general how does one know which coefficients are sparse?

This comes down to how do you choose $t_{B_k}$ and $H$. That will depend on what system you're trying to model.


Identify the correct spot to place a label In the image shown below, I have a 2D data set where I have identified four clusters labeled [0,1,2,3]. I'm looking for an algorithm to place the labels in a natural* way for each shape. My first guess was to place them in the "center-of-mass" of the cluster, which is shown below. For contiguous clusters this works fine. For clusters that take a shape as in cluster 0 however, the approach fails. Without using a legend, what would be a better method to place the labels on this image?

* natural here is subjective to some extent, but the point of the labels is to help the viewer associate a certain region in the x-y plane with a number.

What about putting the label at the innermost point of the segment? Let's define innermost by the maximum of the distance transform of the segment's mask.

With software systems like Mathematica and the sort, it is straightforward to achieve.

The mask for one segment, and its distance transform:

After repeating for each segment and positioning labels where the individual distance transforms are maximum:


Decoupling two signals where one is a delayed version of the other I have a signal $y(t)$ that consists of the sum of two signals where one is a time delayed version of the other.

$ y(t) = x(t) + x(t-t_{d}) $

$x(t)$ is known but the time delay $t_d$ is not.

I would be interested to know what (if any) techniques could be employed to recover the time delayed signal from $y(t)$.

I believe for techniques such as ICA it is required to have two receivers to decouple the signals fully however neither signal would be known in advance for this case.

You can employ cross correlation of the received signal with the signal you know that is contained.

Fs = 1000;
t = np.arange(0, 2, 1./Fs)
s = lambda t: np.cos(2*np.pi*10*t*t)

t0 = 0.5

single = s(t)
signal = s(t) + s(t-t0)

plt.subplot(1,2,1)
plt.plot(t, signal)
plt.plot(t, single)

plt.subplot(1,2,2)
plt.plot(t, np.correlate(signal, single, 'same'));

The delay between the peaks is the time distance between both signals.


Frequency response of a filter - fft vs. freqz I have designed my filter using firhalfband(N,Fo,A), which gives me the filter coefficients. Then I try to see the frequency response of it. I use these two function: fft, and freqz, but they are giving me different result in my figure, why is that? I wonder how these two MATLAB functions operate when taking the frequency response of a signal. I am expecting them to be the same.

Thank you

If you have the filter coefficients, freqz is the way to go - the various optional input arguments and the expected outputs are documented in MATLAB which you can see using help freqz. Snippets given below:

" freqz Digital filter frequency response.
    [H,W] = freqz(B,A,N) returns the N-point complex frequency response
    vector H and the N-point frequency vector W in radians/sample of
    the filter:

[H,F] = freqz(B,A,N,Fs) and [H,F] = freqz(B,A,N,'whole',Fs) return
    frequency vector F (in Hz), where Fs is the sampling frequency (in Hz).  "

If you have filter's impulse response - which will be different from the filter coefficients - you can use it to find the frequency response of the filter using the FFT.


Designing Butterworth filter I'm looking to implement the designing of a Butterworth filter based on the implementation in MATLAB/Scipy into C# but I've ran out of knowledge and can't quite seem how to compute the coefficients b and a  based on the order and the normalised cutoff frequency.

I think that I am calculating z, p and k correctly using:

var z = new double[0];
var m = new List<double>();

var start = -order + 1;
var step = 2;
while (start <= order)
{
    m.Add(start);
    start += step;
}

var complexP = new List<Complex>();

for (var i = 0; i < m.Count; i++)
{
    complexP.Add(-Complex.Exp(Math.PI*m[i]/(2*order)));
}

var k = 1;

I believe that to calculate b and a using z, p and k I have to find the coefficients of a polynomial based on the root provided by p. But I'm struggling to understand how you calculate a coefficient of a polynomial based on a root value.

I understand the result of this should always return 1 as the first value of the a vector but I can't find anything that describes how to calculate the coefficient, perhaps I've found the answer but I've brushed past it.

For the b value I believe that this is k * polynomialCoefficients(z).

Also to note is that I haven't taken into consideration the cutoff values in to how they relate to z, p and k but this doesn't look too difficult and I don't think it will affect the question.

I'm not looking for this is how you do it in C#, but a general run down of how to calculate the coefficient would help.

Figured out how to correctly change from the zero-poles coefficients to the numerator/denominator coefficients.

Essentially this simply requires a convolution of the negative of each value with 1 being added as the first value in the results vector.

The implementation of this for the application of convolution in the instance of filter design is:

var result = new List<Complex>
{
    1
};

for (var i = 0; i < roots.Length; i++)
{
    var value = -roots[i];
    result.Add(result[result.Count - 1] * value);
    for (var j = result.Count - 2; j >= 1; j--)
    {
        result[j] = result[j] + value * result[j - 1];
    }
 }

To visualise this we start with the roots as being:
[0.954697396423354, 0.8855894564017569, 0.86046511627906985, 0.8855894564017569, 0.954697396423354]

Before iteration:

result = [1]

First iteration:

result = [1, -0.9546974]
result[length-1] = 1 * -0.954697396423354

Second iteration:

result = [1, -1.84028685  0.84546995]
result[length-1] = -0.9546974  * -0.8855894564017569
result[length-2] = -0.9546974 + -0.8855894564017569 * 1

Third iteration:

result = [1, -2.70075197  2.42897259 -0.7274974]
result[length-1] = 0.84546995 * 0.86046511627906985
result[length-2] = 0.84546995 + 0.86046511627906985 * -1.84028685
result[length-3] = -1.84028685 + 0.86046511627906985 * 1

And so continues until you reach the end of the input vector.

The main thing to note is that the initial value (1) must never change in the result vector so it is required in the inner loop to make sure you skip the last value, as you've just calculated and added it, and also the first value.


Amplitude response of Lynn's low pass filter at the 3dB cutoff frequency I am referring to the calculations in https://courses.cs.washington.edu/courses/cse474/18wi/labs/l8/QRSdetection.pdf for the second order Lynn's low pass filter.
On the third page, the author mentions the transfer function used as:
$$
H(z)=\frac{(1−z^{-6})^2}{ (1−z^{-1})^2 }
$$
The amplitude frequency response for this is given as:
$$
|H(e^{j\omega T})|=\frac{\sin^2(3\omega T)}{\sin^2(\omega T/2)}
$$
Now, he has calculated the 3dB cutoff frequency to be around 11 Hz for the sampling frequency of 200 Hz (i.e. sampling time = 0.005 seconds). Substituting these values in the above equation for amplitude response, I get:
$$
|H(e^{j\omega_c T})|=\frac{\sin^2(3 \cdot 11 \cdot 0.005)}{\sin^2(11 \cdot 0.005/2)} = \frac{0.026979}{0.000756} = 35.69
$$
At the 3dB cutoff frequency indicated, shouldn't the amplitude response work out to around 0.707 and not 35.69? Where is my understanding incorrect?

Here's your mistake, you forgot that $ \omega = 2\pi f$


DAC implementation in OFDM In an OFDM transmitter, after IFFT & CP block, a DAC block is used.Since the output of IFFT is complex, and the input to DAC should be an integer, how can we integrate the two blocks? I implemented OFDM (upto CP block) in MATLAB but I could not integrate DAC with it.

Regards
Ehsan

In a real system you actually need two DACs: one for the real part and one for the imaginary part. The implementation of the DAC model depends on what should be included in your simulation. There are two main aspects: quantization and pulse shaping.

Quantization - if quantization effects of the DAC should be modelled, the floating point output of IFFT/CP insertion has to be mapped to a discrete number of values. Quantization can, for example, be implemented by scaling and rounding. Say your DAC has a vertical resolution of 8 bit. Then the IFFT output signal should be scaled such, that its maximum value is 127. Afterwards the signal is rounded to integer numbers. Note that OFDM signals are often clipped before quantization to reduce quantization noise, especially for lower resolutions (1).
Pulse shaping - if a physical channel should be modeled (for example a wireless channel), the analog output of the DAC must be simulated. This is implemented by converting every input sample into an impulse that is weighted with the sample's value. In the simulation the impulse itself consists of $L$ samples, of course, that mimic the analog signal. $L$ is the oversampling factor (not to be confused with the OFDM oversampling that is achieved by zero padded subcarriers). The simplest impulse shape is a rectangle, i.e. every sample in the simulation is just repeated $L$ times.

If the channel model is very simple, for example an AWGN channel, DAC and ADC can be omitted. It really depends on how realistic the simulation model should be.

(1) Michael Bernhard et al.: Analytical and numerical studies of quantization effects in coherent optical OFDM transmission with 100 Gbit/s and beyond


How to generalise the Fourier transform? The Fourier transform takes a signal and splits it into a series of sine and cosine waves.

I am told that it's supposed to be possible to split a signal into some other set of functions. My question is: How do you do this?

I'm presuming that the set of functions you use would have to have certain properties for this to work. (E.g., you have to have "enough" different functions to capture all the information of the original signal.) How do you figure out whether your set of functions is suitable? And then how do you do the actual splitting?

The Fourier Transform is just one of so many different transforms that alter the representation of (usually) a time-series from the time domain, to another domain (usually a frequency domain but other representations exist for other transforms such as time / frequency, time / scale and others).

You can find much more information about transforms in general from this Wikipedia list of articles that lists some popular and often used transforms. (You might want to focus on the Discrete and Integral Transforms at first)

Alternatively, you can check out this recent discussion on how the Wavelet transform, achieves a decomposition similar to that of the Fourier transform.

Finally, when you have the luxury to have acquired simultaneously many different time series from the same phenomenon you can even employ techniques such as Principal Component Analysis (PCA) and Independent Component Analysis (ICA) which go to the point of transforming a signal to a sum of elementary waveforms that are actually extracted from the signal itself (rather than being pre-set as it is done in the Fourier (and related transforms) or Wavelets).


Image analysis - fiber recognition I am completely new to image analysis. Do you know how to binarize this image in such a way to get the fibers only?

A step by step procedure would be awesome, but any hint is appreciated as well. I can use Matlab, .NET and I am open to any other tool.

PS: At this address you can find the same question

Hough transform is your friend. Basicaly it deals with straight lines (and almost all fibers in your example are straight). You can find open source implementation in OpenCV (here is doc + result example).

Moreover there is some generalizations that able to detect custom shapes.
P.S. Since you're using C#, OpenCV wrapper could be helpfull for you.


Phase synchronization between 2 PLLs If 2 PLLs as the one shown are driven from the same reference/crystal and the ratio between divider N and R is an integer (note that divider R is outside the loop), does one need a sync pulse for divider R or not and why? what if the division ratio N/R is not an integer?

Yes, in order to create multiple PLLs that are phase-coherent with one another, then you need to specifically design in some functionality to implement that synchronization.

The issue is that each of the dividers in the block diagram that you showed can have multiple stable operating phase offsets (phase offsets where the input/output frequency relationship is met). If you don't specifically design the system such that those operating points are deterministically chosen in some way, then there will be some residual phase offset in the output based on the randomly-occurring phase offsets in the divider outputs.


3D effects like HTC One Dimension plus Its true that HTC one M8 uses dual camera to achieve the 3D Dimension plus effect, but can anybody tell what is it called, if its a series of image manipulation?

Demonstrated here,
Video

The technique demonstrated in the video appears to be "image morphing" which is not specific to 3D.  They have "morphed" a left-eye view to a right-eye view and back.  The resulting animation is a form of "wiggle" stereoscopic presentation.

See http://en.wikipedia.org/wiki/Morphing

Here's an example I just made from an old stereo pair
of a Civil War cannon, from the Library of Congress.  I used
ImageMagick's "-morph" option to rock back and forth between
the two images:

http://www.simplesystems.org:/users/glennrp/stereo/loc_cannon.gif

See also the "wiggle method" of stereoscopy, with some more examples,
at http://en.wikipedia.org/wiki/Wiggle_stereoscopy
This method typically just alternates between the left and right images, but sometimes "morphing" or filling in with intermediate photos are used.


Polyphase decomposition of IIR filter Is there a way to decompose IIR filter into set of poly phase components, that are All pass. E.g. a an IIR filter  of order M, can be decomposed into a cascade of All zero filter of order 2*M( which can be implemented as polyphase) and all pole of same order( as Direct form), by rationalizing the Denominator. But, I have not been able to find a way to do it directly as sum of all pass polyphase components.

any help?

Using matlabs tf2cl, it may require the DSP toolbox, you can turn the transsfer function (numerator and denominator coefficients) in to a cascade allpass lattice coefficients, which is a two phase filter.

% Generate some IIR filter coefficients
[b,a]=cheby1(9,.5,.4);

% Coupled allpass decomposition
[k1,k2]=tf2cl(b,a);


Using large DFT for OFDM to avoid the guard interval In OFDM, we usually insert the guard interval to avoid the inter-symbol interference at the expense of achieved data rate. that guard interval can be either zeros or CP guard interval, which it must be longer than the channel delay.
I am wondering, what's about using long symbol to avoid that guard interval ? for example instead of having number of subcarriers $N= 64$ or $N = 1024$, we can use $N = 5120$, what's the impact of using a high number of $N$ ??

Sure,  you can !
Practically, some systems are using OFDM system with number of subcarriers 4096, but the PAPR will be an issue in that case, and if your channel is fast-fading, the time-domain equalizer will be high-complex since one-tap frequency-domain equalizer will be useless.
Usually, choosing the length of the symbols depends on the channel delay, e.g., if your channel delay is 128,  you must use at least length of symbol 1024 and so on.


$\mathcal Z$-transform ROC Let's say I have a $\mathcal Z$-transform that represents some transfer function and its has some ROC.

My question is how do I know if this system is causal?

I know that if the ROC contains the unit circle in the complex plane the system is stable and it has DTFT but I don't know how I can evaluate the causality. This also applies to the continuous domain, I believe the system is causal if the ROC is right sided, for example, $\Re(\sigma)>1$ but I also don't know why.
(All of this for LTI systems)

If the ROC is outside a circle in the complex $z$-plane ($|z|>a$), then the corresponding system is causal. If it is inside a circle ($|z|<a$), the system is anti-causal. If the ROC is a ring ($a<|z|<b$), the impulse response of the corresponding system is two-sided, i.e., the system is non-causal. In all cases, if the ROC contains the unit circle, the system is stable.

Consider the causal sequence $x[n]=a^nu[n]$ (where $u[n]$ is the discrete-time unit step sequence). The bilateral $\mathcal{Z}$-transform is defined by

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

Since $x[n]$ is causal we get

$$X(z)=\sum_{n=0}^{\infty}a^nz^{-n}=\sum_{n=0}^{\infty}\left(\frac{a}{z}\right)^n=\frac{1}{1-\left(\frac{a}{z}\right)}=\frac{z}{z-a},\qquad \left|\frac{a}{z}\right|<1\tag{2}$$

where I've used the formula for the geometric series. Note that the series only converges for $|a|<|z|$, i.e., outside a circle in the complex plane.

In a completely analogous manner we get for an anti-causal sequence such as $x[n]=a^nu[-n]$

$$X(z)=\sum_{n=-\infty}^{0}a^nz^{-n}=\sum_{n=0}^{\infty}a^{-n}z^n=\sum_{n=0}^{\infty}\left(\frac{z}{a}\right)^n=\frac{1}{1-\left(\frac{z}{a}\right)}=\frac{a}{a-z},\qquad \left|\frac{z}{a}\right|<1\tag{3}$$

which converges for $|z|<|a|$, i.e., inside a circle in the complex plane.

Finally, a two-sided sequence can be split up into a causal and an anti-causal sequence, with an ROC that is the overlap of the two individual ROCs, which, if it is not empty, is a ring in the complex plane.


How to estimate the frequency at certain time of linear sine sweep? experts.
I have a question concern "linear sine sweep".
Here is the equation of linear sine sweep(from Wikipedia Chirp).
$$ x(t) = sin\biggl[ 2\pi( f_0 +\frac{f_1-f_0}{2T}t)t\biggr] $$
where,

$f_0$ is start frequency
$f_1$ is end frequency
$T$ is total length of excitation time.

In my case, I would like to calculate the excited frequency at certain time frame.
Let me calculate it through an example.
$f_0 = 20 \text{ Hz}$ , $f_1=200 \text{ Hz}$, $T = 60\text{ s}$.
In this case, above equation can be written as below,
$$ x(t) = sin\biggl[ 2\pi( 20 +1.5t)t\biggr]. $$
If we estimate the frequency at time 30 s, I think we could calculate as 65 Hz.
For checking this is right or not, I checked the power spectrum by selecting time through 29.5 s  - 30.5 s for frequency resolution 1 Hz.
However, from the power spectrum, the frequency was 112 Hz.
There is difference between simulation value and calculated value almost twice.
Also, I just checked 1 second(29.5 s - 30.5 s) Fourier transform absolute value, it also displayed the peak frequency as 112 Hz.
What is the correct value?
If the calculated value, 65 Hz, is wrong, is there any comment for correcting this?
Thank you in advance.
If there is any comment or question, please let me know.

Frequency is the first derivative of the phase so in this case that's f0+(f0-f1)*t/T. If you plug in the numbers you get 110Hz, which is what your measurement shows. Not sure how you got 65Hz for that ??

It's intuitively clear as well. The sweep starts at 20 Hz and ends at 200 Hz after 60s. After 30s you should be at the half way point between 20 and 200, which is 110 Hz.

$$\varphi (t) = 2 \pi \cdot (f_{0}+\frac{f_{1}-f_{0}}{2T}\cdot t)\cdot t = 2 \pi \cdot (f_{0}\cdot t+\frac{f_{1} -f_{0}}{2T}\cdot t^{2})$$

$$f(t)=\frac{1}{2 \pi}\cdot \frac{\partial \phi}{\partial t}=f_{0}+\frac{f_{1}-f_{0}}{T}t$$


Does 46 dB gain of a filter for a frequency imply 200 times more amplitude? Consider we have a IIR filter or order 2 and it has gain of say 46.02 dB for frequency $f_0$. Now if we apply a time domain signal of frequency $f_0$ to the filter for a long time, can we expect to get 200 times higher amplitude in the output of the filter comapre to the input of the filter? The 46.02 db and 200 comes from the fact that $20\log_{10}(200)= 46.0206$.

can we expect to get 200 times higher amplitude in the output of the filter comapre to the input of the filter

Yes.

For a sine wave input the amplitude of the steady state output is simply the amplitude of the input multiplied with the magnitude of the transfer function at that frequency.


Complex Morlet coefficient A complex Morlet wavelet looks like this:

$$\psi(t) = C \cdot e^{i \omega t} \cdot e^{-t^2/2}$$

Here $\omega$ is the frequency and $C$ is some normalization constant. The first exponential represents the oscillation, and the second exponent is the Gaussian-like envelope. However, what should $C$ be?

First, why is $C$ even necessary? Admissibility of a wavelet is based on finite support (which the Morlet wavelet has practically, if not really) and its integral over its whole range summing to $0$, which is does regardless of $C$'s value.

Torrence and Compo's A Practical Guide to Wavelet Analysis, on whose work mlpy's wavelets are based, argue for the constant $\pi^{-\frac{1}{4}}$ to "ensure a total energy of unity". Wikipedia's Morlet wavelet page gives the more complicated coefficient $\pi^{-\frac{1}{4}} \left(1 + e^{-\sigma^2} - 2e^{-\frac{3}{4}\sigma^2}\right)^{-\frac{1}{2}}$, where $\sigma$ is a parameter allowing some kind of tradeoff between time and frequency resolution. I'm having a hard time rectifying these two different approaches.

Finally, does $C$'s value have any relation to whether you want the resulting frequencies to be on a logarithmic scale, or a linear scale?

Your $C$ is here for the wavelet energy normalization: it is chosen such that:

$$ \int \psi(t)^2dt = 1\,,$$

the most  common framework for wavelets. It ensures that wavelet coefficients can more or less directly be interpreted in energy. The problem is that, as is, the wavelet is not admissible. If you take $\omega \simeq 0$, it reduces to a Gaussian, and does not  oscillate enough to have zero-mean.

However, when $\omega$ becomes "big enough", some  say $\omega>5$ (a common choice is $\omega = \pi\sqrt{2/\log2} $), one usually considers that it is "numerically admissible" (with high precision).

A more proper version (with $\omega_0$ instead of $\omega$ to avoid strange interpretations) is:

$$ \psi_{\omega_0} = \frac{C_{\omega_0}}{\pi^{1/4}} e^{-t^2/2}(e^{i\omega_0 t} -e^{-\omega_0^2/2})$$

and then you get:

$$ 1/C_{\omega_0}^2 = 1+e^{-\omega_0^2} -2e^{-3 \omega_0^2/4} \,,$$

which is very close to $1$ when $\omega_0$ is big enough.


correcting delay between two samples I'm sampling two signals with fundamental frequency of 50 Hz. I want to sample them at the same instant, but there's a delay between the first and the second conversion. I don't keep the previous values of the signals in an array, at most I can store 10/20 samples.

I want to correct this delay. I thought of a FIR fractional delay, but i'm not sure is the right way to do it.

Assuming you are sampling above the Nyquist rate, then an interpolation kernel or equivalent fractional FIR filter on one or both channels (both to match up any non-zero processing delays) might be an appropriate solution.


Simple crosstalk cancellation (or signal separation)? To start with, I'm pretty new to DSP, so forgive me failing to understand any basic concepts.

I'm working with an iOS app that depends on the input from the analog audio jack whilst monitoring the signal through output (headphones, etc). Unfortunately, iOS devices suffer from crosstalk given a high enough volume due to the proximity of the output and input leads.

Garageband on iOS provides something called "Crosstalk prevention" for monitoring guitar audio playback through its amp modelling. I've no idea how they do it but it's something I'd like to do in my app.

I've access to the output signal, and I perform a FFT on the input signal (to perform spectral noise gating) anyways, so anything that uses fourier transforms might be good since I'm doing that operation. It's difficult for me to find any solutions as DSP is pretty new to me.

If I understand correctly what you are trying to achieve is usually called Acoustic Echo Cancellation. In many cases, so-called adaptive filtering algorithms, such as the (normalized) least-mean-squares ([N]LMS) algorithm, can be used effectively. In this context, I can really recommend the book Adaptive Signal Processing by Widrow and Stearns. I'm sure you'll find a number of ready-to-use implementations of the (N)LMS algorithm to experiment with on the internet. (Have an eye on the filter order - it needs to be high enough to capture the length of your crosstalk path.)

Furthermore, as you mentioned that you are already working in the frequency domain, you can have a look on algorithms that are called (you guessed it) frequency-domain LMS algorithms. They originate, via so-called block LMS algorithms, from the classical LMS algorithms. In this context, I really liked the paper by Shynk: Frequency-domain and multirate adaptive filtering.


Recursive filter with repeated poles This post follows this previous resolved post where I was trying to find the inverse Z-transform of a more simple filter output (that I used as an example to get the methodology). The present filter implements an additional integrator and I tried to use the same method. The problem is, I am not sure about how to deal with repeated poles.

$H(z)$ and $X(z)$ are given by:

$ H(z)= \frac{(1-z^{-K})(1-z^{-L})}{{(1-z^{-1})}^{2}}$
$ X(z) =\mathcal{Z}\{x[nT]\}= \mathcal{Z}\left\{Ae^{-\frac{nT}{\tau}}u[n]\right\} = A\dfrac{z}{z-e^{-\frac{T}{\tau}}}= A\dfrac{1}{1-e^{-\frac{T}{\tau}}z^{-1}} $

I want to express $y[n]$ as a function of the pulse amplitude $A$, $L$, $K$, $n$ (sample number), $T$(sampling period), and the time constant $\tau$; my goal being to retrieve the pulse amplitude $A$.

I used polynomial long division to develop $H(z)$ and i came up with something like this, assuming $K>L$:

$(1)\space H(z)=1+2z^{-1}+3z^{-2}+...+(L+1)z^{-L+1}+Lz^{-L}+Lz^{-L-1}+..+Lz^{-K+1}+(L-1)z^{-K}+(L-2)z^{-K-1}+...+z^{-K-L+2}$

Now I describe $Y(z)$ using partial fraction expansion: $Y(z) = H(z).X(z) = A.\left[\frac{(1-z^{-K})(1-z^{-L})}{{(1-z^{-1})}^{2}(1-e^{-\frac{T}{\tau}}z^{-1})}\right] = A.\left[\frac{B_0}{(1-e^{-\frac{T}{\tau}}z^{-1})} + \frac{B_1}{(1-z^{-1})} + \frac{B_2}{{(1-z^{-1})}^2}\right] $
$$\begin{align*}
& (1-z^{-K})(1-z^{-L})=B_0.{(1-z^{-1})}^{2}+ B_1.(1-z^{-1}).(1-e^{-\frac{T}{\tau}}z^{-1})+B_2.(1-e^{-\frac{T}{\tau}}z^{-1}) \\
& H(z)(1-z^{-1})^{2}= B_0.{(1-z^{-1})}^{2}+ B_1.(1-z^{-1}).(1-e^{-\frac{T}{\tau}}z^{-1})+B_2.(1-e^{-\frac{T}{\tau}}z^{-1})\\
\end{align*}$$
We obtain:
$$\begin{align*}
&B_0 = H(z) = 1+2z^{-1}+3z^{-2}+...+(L+1)z^{-L+1}+Lz^{-L}+Lz^{-L-1}+..+Lz^{-K+1}+(L-1)z^{-K}+(L-2)z^{-K-1}+...+z^{-K-L+2}\\
& B_2 = 0 \\
\end{align*}$$
Not sure about $B_1$:
$$\begin{align*}
&(2)\space B_1 = \frac{d}{dz^{-1}}\left[H(z).(1-z^{-1})^{2}\right]_{z^{-1}=1} = \frac{d}{dz^{-1}}\left[\frac{1-z^{-K}-z^{-L}+z^{-K-L}}{1-e^{-\frac{T}{\tau}}z^{-1}}\right]_{z^{-1}=1} = \left[\frac{(Kz^{-K}+Lz^{-L}-(L+K)z^{-K-L})(1-e^{-\frac{T}{\tau}}z^{-1})-(1-z^{-K}-z^{-L}+z^{-K-L})(e^{-\frac{T}{\tau}}z^{-2})}{(1-e^{-\frac{T}{\tau}}z^{-1})^{2}}\right]_{z^{-1}=1} \\
\end{align*}$$
If $B_0$ is correct (I am not sure), the inverse z-transform of $\frac{B_0}{(1-e^{-\frac{T}{\tau}}z^{-1})}$ can easily be retrieved. As for  $B_1$, I get:
$$\begin{align*}
&(3)\space \frac{B_1}{1-z^{-1}}=-\frac{(1-z^{-K})(1-z^{-L})e^{-\frac{T}{\tau}}z^{-2}}{(1-e^{-\frac{T}{\tau}}z^{-1})^{2}(1-z^{-1})}=-\frac{(1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1})e^{-\frac{T}{\tau}}z^{-2}}{(1-e^{-\frac{T}{\tau}}z^{-1})^{2}}{}
\end{align*}$$

Now according to this table http://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html, my expression looks similar to the follow one.

However I tried to apply it on a concrete example (using Excel), and it looks like I have missed something or got something wrong.

Any help would be greatly appreciated

I assume that by "inverse Z-transform" you simply mean the impulse response of the filter.

That's much easier done in the time domain directly. The first section has an impulse response of
$$
h_1[n] =\begin{cases}
1 & \text{ if } n= 0\\
-1 & \text{ if } x= K\\
0 & otherwise
\end{cases}
$$
So basically a positive spike at n=0 and a negative spike at $n = K$

The integrator picks up the first spike and keeps it until you hit the negative spike, which cancels the positive one and everything is zero again. So the resulting impulse response is a rectangular window
$$h_2[n] =\begin{cases}
1 & \text{ if } 0 \leq n < K \\
0 & otherwise
\end{cases}$$

The second half of the filter has exactly the same structure as the first half and so it produces a rectangular window of length $L$.

The overall impulse response is then simply the convolution of both rectangular windows. This will be a triangle for $K=L$ and a trapezoid otherwise. The ramps of the trapezoid have the length of the smaller of K and L and the length of the flat part is the difference between K and L. Total length is $L+K-1$

Assuming K is the smaller one we could write this as

$$h[n] =\begin{cases}
n+1 & \text{ if } 0 \leq n < K \\
K & \text{ if } K \leq n < L \\
K-1-(n-L) & \text{ if } L \leq n < L+K-1 \\
0 & otherwise
\end{cases}$$


Convolution Theorem using DCT I am trying to digest this paper.
[2002] Anna Usakova. Using Of Discrete Orthogonal Transforms For Convolution
So far, I've been successful with DFT and DHT.
Could someone give a hand with DCT?
I am trying to reproduce the convolution theorem for DCT-I and DCT-II. But my results never match the same as the result I get with DFT and DHT.

I am reading reference 9 of the paper

RAO, K. R.—YIP, P.: Discrete Cosine Transform — Algo- rithms, Advantages, Applications, Academic Press, Inc., Lon- don, 1990.

But it puts the conv. theorem of DCT-II in other terms.
I have also read other works, such as

[1994] Martucci. Symmetric Convolution and the Discrete Sine and Cosine Transforms
[2011] Roma. A tutorial overview on the properties of the discrete cosine transform for encoded image and video processing

and others. But I am interested on the Conv. Theorem on a similar fashion as the first two references.

I did it for the DCT-I and DCT-II.
At first I thought it is about circular convolution, but it is not. After desperate attempts to do it by myself I found the article: Convolution Using Discrete Sine and Cosine Transforms (Alternative Link).
In this paper there is a derivation of circular convolution using DTT's, but it is constructed from sum of $C_1^{-1}$ and $S_1^{-1}$, so, it's not what we are interested in.
What I suggest is to look at Fourier transform first, in order to get the idea:
Let we have
$\displaystyle F[t] = \sum_{k=0}^{N-1} f[k]e^{-\frac{i2π}{N}tk}$,
$\displaystyle G[t] = \sum_{k=0}^{N-1} g[k]e^{-\frac{i2π}{N}tk}$.
So then
$$\displaystyle ℱ\{f*g\}[t] = \sum_{x=0}^{N-1} f[x]e^{-\frac{i2π}{N}tx}\cdot\sum_{y=0}^{N-1} g[y]e^{-\frac{i2π}{N}ty} = \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f[x]g[y]e^{-\frac{i2π}{N}t(x+y)}$$
Now, let $x+y = m$, $m = 0\ldots2N-2$ and then group all $(x,y)$ which share the same $m$. So you will get:
$$\displaystyle \sum_{m=0}^{2N-2}\biggl(\sum_{0\leq x,y\leq N-1\\x+y=m}f[x]g[y] \biggl)   e^{-\frac{i2π}{N}tm} \tag{1}$$
Because $x, y$ are both positive and less than $N$, then sum inside is just linear convolution at $m$. Let me denote it as $A[m]$, then everything is going to be rewritten in that form
$$\sum_{m=0}^{2N-2} A[m] e^{-\frac{i2π}{N}tm} \tag{2}$$
DFT kernel is N-periodic that is why we can group terms another time:
$$\displaystyle \sum_{m=0}^{N-1} (A[m] + A[m+N])e^{-\frac{i2π}{N}tm}$$
$A[m]$ captures $A[0]\ldots A[N-1]$, $A[m+N]$ captures $A[N]\ldots A[2N-1]$. But remember that $A[m]$ is linear convolution of two sequences of size $N$ and therefore size is $N + N - 1$ then $A[2N-1] = 0$. So these two sums are equivalent. Also it is a DFT of $A[m] + A[m+N]$. What is it? Right, circular convolution, as expected.
Now plug in DHT kernel instead of DFT into $(1)$ and you will get the same result ─ circular convolution. That is because DHT and DFT kernels are quite similar.
Well, if you have 2 signals $f$, $g$ with lengths $Q$ and $P$, and $A = f*g$, with $Q+P-1 = L$. Now, add zeros to $f$ up to $L$, do the same for $g$. Now $f$, $g$, and $A$ are $L$ periodic signals, so we can take the DFT of them. DFT of $A$ is going to be:
$$\sum_{m=0}^{L-1} A[m] e^{-\frac{i2π}{L}tm} = \sum_{m=0}^{L-1}\biggl(\sum_{0\leq x,y\leq L-1\\x+y=m}f[x]g[y] \biggl)   e^{-\frac{i2π}{L}tm} = \sum_{x=0}^{L-1}\sum_{y=0}^{L-1}f[x]g[y]e^{-\frac{i2π}{L}t(x+y)} \tag{3}$$
So, via $(2)$ we can get the rule and expression for the convolution and $(3)$ shows that zero-padding up to convolution length is required to get linear convolution using the same rule. We're going to apply it to the DCT's and DST's.
Because there are many varieties of DCT's and DST's, I will mention what I used:
DCT-I (N+1) points:
$\displaystyle C_1\{f\}[n] = \sum_{k=0}^{N} c_k f[k]\\cos\Bigl[\frac{π}{N}nk\Bigl], \text{where}\
c_k =
\left\{
\begin{array}{l}
\tfrac12, \qquad &k (\mod N)=0, \\
1,\qquad &\text{otherwise} \\
\end{array}
\right. \tag{4}$
DST-I (N+1) points:
$$S_1\{f\}[n] = \sum_{k=0}^{N} f[k]\sin\Bigl[\frac{π}{N}nk\Bigl]=\sum_{k=0}^{N} c_k f[k] \sin\Bigl[\frac{π}{N}nk\Bigl] \tag{5}$$
(N-1) for data: $f[1] \ldots f[N-1]$ it doesn't matter what is placed in $f[0]$ and $f[N]$.
DCT-II (N+1) points:
$$C_2\{f\}[n] = \sum_{k=0}^{N} f[k]\cos\Bigl[\frac{π}{N}\Bigl(k+\tfrac12\Bigl)n\Bigl] \tag{6}$$
N for data: $f[0]\ldots f[N-1]$ and $f[N] = 0$
DST-II (N+1) points, N for data:
$f[0]\ldots f[N-1]$ and $f[N] = 0$:
$$S_2\{f\}[n] = \sum_{k=0}^{N} f[k]\sin\Bigl[\frac{π}{N}\Bigl(k+\tfrac12\Bigl)n\Bigl] \tag{7}$$
Inverse DCT-I:
$$C_1^{-1} = \frac{2}{N}C_1 \tag{8}$$
Inverse DST-I:
$$S_1^{-1} = \frac{2}{N}S_1 \tag{9}$$
Inverse DCT-II (N+1) points:
$$C_2^{-1}\{F_{C_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{C_2}[n]\cos\Bigl[\frac{π}{N}n\Bigl(k+\tfrac12\Bigl)\Bigl] \tag{10}$$
$C_2^{-1}\{C_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $f[N-1]$.
Inverse DST-II (N+1) points:
$$C_2^{-1}\{F_{S_2}\}[k] = \frac{2}{N}\sum_{n=0}^{N} c_n F_{S_2}[n]\sin\Bigl[\frac{π}{N}n\Bigl(k+\tfrac12\Bigl)\Bigl] \tag{11}$$
$S_2^{-1}\{S_2\{f\}\}[k] = f[k]$ except for $k=N$ there it is going to be $-f[N-1]$. Also, they are not interchangeable.
Everything mentioned above is just Wikipedia ones, but modified to achieve more similar look for easier coding and derivation.
So, now again, applying the rule $(2)$ for the DCT-I we get:
$$\sum_{m=0}^{2N} A[m] c_m \cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{m=0}^{N} p_m(A[m] + A[2N-m]) c_m \cos\Bigl[\frac{π}{N}nm\Bigl]$$
$p_m = \tfrac12$ if $m = N$ else $p_m = 1$
from the other hand we have
$$\sum_{m=0}^{2N} \frac{A[m]}{c_m} c_m \cos\Bigl[\frac{π}{N}nm\Bigl] = \sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\cos\Bigl[\frac{π}{N}n(x+y)\Bigl]$$
$$\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\cos\Bigl[\frac{π}{N}n(x+y)\Bigl] = \\\sum_{x=0}^{N}\sum_{y=0}^{N} f[x]g[y]\biggl(\cos\Bigl[\frac{π}{N}nx\Bigl]\cos\Bigl[\frac{π}{N}ny\Bigl] - \sin\Bigl[\frac{π}{N}nx\Bigl]\sin\Bigl[\frac{π}{N}ny\Bigl]\biggl) = \\ \biggl(\sum_{x=0}^{N} f[x]\cos\Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]\cos\Bigl[\frac{π}{N}ny\Bigl]\biggl) - \biggl(\sum_{x=0}^{N} f[x]\sin \Bigl[\frac{π}{N}nx\Bigl]\biggl)\cdot\biggl(\sum_{y=0}^{N} g[y]\sin\Bigl[\frac{π}{N}ny\Bigl]\biggl)$$
It is now easy to see that it can be rewriten in a such way:
$$C_1\{\frac{f[m]}{c_m}\}[n]\cdot C_1\{\frac{g[m]}{c_m}\}[n] - S_1\{f[m]\}[n]\cdot S_1\{g[m]\}[n]$$
Let $\frac{f[m]}{c_m} = f^{/}$ then we can write everything in this form:
$$\frac{p_m}{c_m}(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
alternative version
$$δ_{0m}(A[0] + A[2N]) + (A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
or even
$$(A[m] + A[2N-m]) ⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1} - \frac{f[0]g[0] + f[N]g[N]}{2}$$
Now, remember the $(3)$? If you add zeroes, so the size of f and g is equal to the convolution's length, then you have:
$$(f*g)^{/}⟷ F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}$$
It looks a little bit different to what is presented in the paper, but it is because I used different version of DCT.
For DCT-II we have the rule (remember that $f[N] = g[N] = 0$):
$$\sum_{m=0}^{2N} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] = \sum_{m=0}^{2N-2} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl]  = \sum_{m=0}^{N-1} (A[m] + A[2N-1-m]) \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] $$
And the expression is going to be:
$$\sum_{m=0}^{2N} A[m] \cos\Bigl[\frac{π}{N}\Bigl(m+\tfrac12\Bigl)n\Bigl] =
\sum_{m=0}^{2N} A[m] \Bigl(\cos\Bigl[\frac{πn}{2N}\Bigl]\cos\Bigl[\frac{π}{N}kn\Bigl] - \sin\Bigl[\frac{πn}{2N}\Bigl]\sin\Bigl[\frac{π}{N}kn\Bigl]\Bigl) = \\
\cos\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G^{/}_{C_1} - F_{S_1}G_{S_1}) - \sin\Bigl[\frac{πn}{2N}\Bigl](F^{/}_{C_1}G_{S_1} + F_{S_1}G^{/}_{C_1}) = \\
=F^{/}_{C_1}(\cos\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1} - \sin\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}) - F_{S_1}(\cos\Bigl[\frac{πn}{2N}\Bigl]G_{S_1}+ \sin\Bigl[\frac{πn}{2N}\Bigl]G^{/}_{C_1}) = \\
=F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2} = G^{/}_{C_1}F_{C_2} - G_{S_1}F_{S_2} $$
Now we finally can write it as:
$$A[m] + A[2N-1-m] ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$
And again, if you add enough zeroes, then it will be:
$$f*g ⟷ F^{/}_{C_1}G_{C_2} - F_{S_1}G_{S_2}$$
You can play with it here: https://repl.it/@dogtyping/dct#main.py — to make sure that it does work.
Thanks for the question, it was a lot of fun.


QPSK, from Square Root Raised Cosine to Rectangular Impulse Shaping I have implemented a QPSK modulator and demodulator using classic text-book examples: Square Root Raised Cosine Impulse shaping and Matched Filtering.

The problem is that I've been told that I could only used square waves for the modulation instead of sine waves on the emitter side (only +/-Vmod) and this really degrades the performance of the communication.

After further investigation it seems that the root cause of the lack of performance is not really the use of square wave but the change of impulse shaping/matched filtering implied by it.

Basically my question is the following: Using only +/-Vmod (and 0V if necessary) on the emitter side, am I wrong in assuming that the best case scenario is using a Rectangular Impulse Shaping with the according Matched Filter or is there some other combination that could provide better results?

If you have any ideas or links I could learn from don't hesitate! I'm struggling because I don't know if my reasoning is wrong or if its my implementation.

You can use rectangular pulses, but it is not necessary to do so.  Rectangular pulses have advantages and disadvantages.  The advantage is that they are very simple and easy to use, both for the modulator and the demodulator.  It makes life easier for the demodulator because it makes symbol timing very forgiving- you can be off on where you sample the symbol and still get a good measurement.

There are two disadvantages to using rectangular pulses.  The first is that it uses a lot of bandwidth.  The Fourier transform of a rectangular pulse is a sinc function, which drops off very slowly.  The second disadvantage is related to the first- very few real world solutions use rectangular pulses because they use a lot of bandwidth.  If you want something realistic, in terms of what most designers would do professionally, you should not use rectangular pulses.

You could, perhaps, get it working with rectangular pulses, and then when it is solid reintroduce the root raised cosine pulses.


Division of complex signal by frequency response I have a measured signal (time domain) which I take the FFT of. I have a measured frequency dependent sensitivity (magnitude and phase (radians) as a function of frequency) which I need to apply to the measured signal. In effect perform a frequency dependent weighting to the measured signal's FFT.
I calculate the complex signal for the sensitivity magnitude and phase before dividing this into the FFT of the original signal.
I then take the IFFT of this which should give me something like the original signal shape but it does not.

Is my method correct?

There are several things that could be going wrong here.

Check your FFT lengths

If you do multiplication or division in the frequency domain, your frequency domain data needs to be at least as long as $N+M-$ where $N$ is your signal length, $M$ is your weighting length.

In matlab pseudo code this means you need:

N = length(x);
M = length(w);
freq_length = N + M - 1;

X = fft(x,freq_length);
W = fft(x,freq_length);

y = ifft(X./W);

So your resulting y will be longer than x.

If you don't do this, you will get time-domain aliasing because of the circular nature of convolution using the FFT.

Is your weighting real-valued in the time domain?

You may need to use conj(W) in the above, but you also need to check that w (the time domain version of the frequency weighting) is real-valued. If it's not, you will get a complex-valued signal.

Also note that you may get a small imaginary component.  Often times all that is required to remove this is:

yr = real(y);


how sinc pulses can be used to distinguish different symbols Let's say we want to transmit a sequence of numbers Using PAM, this sequence is transmitted with the signal $s(t)$
The pulses in $s(t)$ are transmitted at a rate $R_p=\frac{1}{T_p}$. This is called the "pulse rate" or the "symbol rate". We want $T_p$ to be short, to transmit information faster. However, this increases the bandwidth. The maximum spectral efficiency can be achieved using sinc pulses. I have the following questions please:

How can sinc pulses be used to distinguish the different numbers? If we are using 16-QAM, then we have 16 different symbols. Is it the amplitude of the sinc pulse that will change to distinguish the 16 symbols??
The bandwidth of a single sinc pulse is $\frac{R_p}{2}$, what is the bandwidth of the whole signal $s(t)$ consisting of several pulses, each of width $\frac{R_p}{2}$?

How can sinc pulses be used to distinguish the different numbers? If we are using 16-QAM, then we have 16 different symbols. Is it the amplitude of the sinc pulse that will change to distinguish the 16 symbols?

You send one sinc after the other, and each sinc is multiplied with one complex number (which is your 16QAM constellation point)

The bandwidth of a single sinc pulse is Rp/2, what is the bandwidth of the whole signal s(t) consisting of several pulses, each of width Rp/2?

The bandwidth of something is something that can be seen from the power spectral density, which in itself is an average.
Hence, the spectrum of your transmission, assuming uncorrelated data symbols, is exactly the spectrum of your pulse.

Your two questions are the most fundamental statements underlying all modern digital communications:

you can transport data by sending a symbol impulse train through your pulse shaping filter
the spectrum of the signal is defined only by pulse shape

It feels a bit like you learned about 16QAM before learning the basics of complex baseband and modulation – maybe go a chapter back in your book until you find a block diagram showing a data source, a symbol mapper, a pulse shaper, and a mixer, and then refresh your knowledge of complex baseband, digital comms, if you feel like you're missing basics.


Producing masks from mostly-static-background subtraction when given the background My data set contains a number of high resolution greyscale images containing a person's face. A single image containing the background without the face is provided. The background is mostly static in the image, except for some minor changes in lightning from frame to frame. The background pixel intensities can be similar to those of the face. I need an algorithm that would produce a pixel-perfect mask of the face, given that we know what the background looks beforehand.

(the face outline will not necessarily be as prominent as in this example - edge detection does not find it correctly all the time)

A naive approach I've tried is: run a median filter for noise removal, do a diff of image-background, threshold on low values to produce a binary mask, find biggest contour and fill it in.

However, this fails in cases where:

The background is skin-coloured as image subtraction gives near-0 difference in the face region
The background contains a skin-coloured object right beside the side of the face
The image has a very low contrast. The lower the contrast, the more sensitive it becomes to threshold values, requiring fine tuning for every image.

How would you approach this problem?

You can try working off a different representation of the image, other than color. You can diff the image gradient for example. It is very unlikely that the gradient is the same on the face as the background. A gradient based descriptor such as HOG or SIFT captured at a suitable scale over the image has worked well for me in the past. You can even combine different representations of the image, such as RGB value as well as HOG. IF the image is different in either mark as a change. You will still get the odd few pixels that don't come up as a change but you can smooth over the result with a smoothing kernel with a similar scale tot he size of the object you want to detect.


What is Fourier transform in terms of area under the curve? We know that integration gives area under the curve.  And that FT is also an integration.  How do we interpret FT in terms of the area under the curve especially because e^(jwt) is a complex term?  In general, while dealing with complex terms in integration can we relate to some area?

Since

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

you get for the real part of $X(\omega)$ (assuming that $x(t)$ is real)

$$X_R(\omega)=\int_{-\infty}^{\infty}x(t)\cos(\omega t)dt\tag{1}$$

and for the imaginary part

$$X_I(\omega)=-\int_{-\infty}^{\infty}x(t)\sin(\omega t)dt\tag{2}$$

So if you like, from (1) and (2) you can interpret the real part of the Fourier transform as the area under the curve $x(t)\cos(\omega t)$, and the imaginary part would then be the area under the curve $-x(t)\sin(\omega t)$.

Note that interpreting an integral as the area under a curve does not always help intuition. In the case of the Fourier transform it is more natural to interpret the integrals (1) and (2) as inner products.


Failed to Get Silence in Audio Streaming I want to detect silence in audio streaming. I've been following lot of answers from any website and i feel like know how to do that. but i doubt i'm on the right path coz the result doesn't seems right.

The stream is in mp3 format 16 bit 22050 hz, the data stream divided into small packets. So every packet received I'll check the data if its a header or not to build complete chunk of mp3.
Every complete packets(mp3 chunk) will be decoded into PCM. the PCM data consist of signed value between 0-255. i will split every PCM data into 1.16 seconds of audio to detect silence in that duration. I've read some answer to scale the signals to make it between -1 and 1 for further processing.

My first try is calculating RMS. the PCM data will be splited into 20 miliseconds and calculated to RMS for each number of samples (samples=PCM.length*0.02). Then changes the RMS result to dB with 20*log(rms)/log(10). The result is not as i expected. The more silence in the audio, i get less dB rate (-20dB to -30dB) and get (-50dB to -70dB) for more noise in audio. Then i try to compare it in chart but it doesn't really help at all. I think the result is somehow wrong.

My second try is using FFT. The scaled signals from previous data added with padding to make it 2^n length. The padding used to make signals can be calculated with FFT. After i got the FFT result, then i do A-weighting. I can see some similarity in noise audio. The silence audio gives random chart that can't be compared with noise audio. So i decided to use noise to compare with every a-weighting result.
First i try to compare it with calculating DTW of those signals and make the threshold is not more than 5k in differences, but sometimes i getting audio error(the accuracy is around 90%).
Then i try compare it with another method using correlation ifft(fft(noise) * fft(test)), but i don't know how to use it. I've been comparing same 4 signals and the result doesn't showing any number 1 in that list (the highest is 0.625)

could you tell me what my mistake is? i really confused where my mistake is. actually i want the result like ffmpeg RMS level but i didn't get the same result with ffmpeg in my first try. i can't use ffmpeg library coz my bos ask me to calculate it from memory which is pipe in and pipe out. i've try that and it seems ffmpeg can't be use for that scenario. and tell me if you need some sample of error. I'll provide you that.
thanks

the PCM data consist of signed value between 0-255.

That seems to be your problem. A 16 bit MP3 file decodes into signed integer values of -32768 to +32767. Perhaps you are looking at bytes and not at samples ?


Old style signal processing? I am wondering if "an array of overlapping band-pass ﬁlters" it is 'mostly' equivalent to a Short-Time Fourier transform. See this (figure down), where $s(t)$ is the signal;
LPF a second-order low-pass ﬁlter with a resonance frequency of 4 kHz; BPF an array of band-pass ﬁlters, where frequencies that are spaced equidistantly on
a critical band scale; HCM a hair cell model which is assumed identical in all channels, and finally a low-pass ﬁlter at 1250 Hz does an envelope extraction of the patterns in each channel. Processing a signal with this method produces a matrix of intesities that are vertically distributed like in STFT.

What are the differences between  them? Could be improved with wavelets as basis?

Thank you in advance. All comments will be rewarded!

there is a relationship between filterbank analysis and the STFT analysis, but it is not an identity.

perhaps it can be identical, where for filterbank, you are slicing the frequency domain into band slices and for STFT you are slicing the time domain into windowed segments (in audio we might call them "grains" or maybe sometimes "wavelets") in time.  but to make the two identical will require some thinking.  the BPF can't be simple regular-old BPFs.


How to find minimum length of a FIR symmetric filter if I am given DTFT I am practising for upcoming exams and came across this question.
Let $h[n]$ be an FIR filter such that $h[n] = 0$ when $|n| > M$ and $h[n] = h[−n]$. A plot of $H(e^{j\omega})$ (DTFT of $h[n]$) is shown below.
The eventual goal is to implement this filter digitally. Then, with the given information, is it possible to determine the minimum
value of $M$? Yes or no, with a suitable technical explanation.
What I though was that if I knew the order of this filter I could have found the length of the FIR filter as well. However order cannot be extracted from this graph.

HINT: How many zeros of $H(e^{j\omega})$ do you count? Does that say anything about the (minimum) order of the corresponding polynomial?


Ergodicity of joint process If we have two processes and both of them are ergodic. Does this mean that the joint proces is ergodic? Or other way around? If we have the dynamics for both components of the joint process what are the possible ways to show ergodicityof the joint process?

Joint behavior cannot always be deduced from individual behavior. For example, if $X$ and $Y$ are (nondegenerate) random variables with
finite means, then $P\{X<E[X]\}$ and $P\{Y<E[Y]\}$ both are nonzero. This
is because there must be probability mass to the left of the mean (and to the right of the mean too) except in the degenerate case when $X$ has value $E[X]$ with probability $1$. However, $P\{X<E[X],Y<E[Y]\}$ can be zero. Don't believe me? Consider $X\sim$ Bernoulli$(\frac 12)$ and
$Y=1−X$ which is also Bernoulli$(\frac 12)$. Both have means $\frac 12$
but $P\{X < \frac 12, Y < \frac 12\} = 0$.

Turning to ergodicity, one characteristic of an ergodic process is that
over long periods of time (say duration $T$), a typical sample path has value $\alpha$ or less for approximately $F_X(\alpha)\cdot T$ time units where $F_X(\cdot)$ is the first-order CDF of all the random variables
comprising the process. That is, we can get good approximations to the
value of the first-order CDF at $\alpha$
by observing what fraction of time a sample path spends
below level $\alpha$. So, going back to your question, if $\{X(t)\}$
and $\{Y(t)\}$ are ergodic processes, we cannot deduce that the
joint process is also ergodic, for reasons similar to the example
I gave in the first paragraph above. Now, if you impose conditions
such as $\{X(t)\}$ and $\{Y(t)\}$ are independent processes
in addition to being ergodic processes, then Yes, the joint
process is also ergodic.


Under what condition is the available bandwidth $f_s$ instead of $\frac{f_s}{2}$? Can someone help me with this interview question? It is about the sampling frequency rate.
Under what condition the available bandwidth is $f_s$ instead of $\frac{f_s}{2}$?

If your signal is real-valued, then it's spectrum is conjugate symmetric. That means, that negative frequencies (or frequencies from $\frac{f_s}{2}$ up to $f_s$) are mirrored. Thus we can always neglect frequencies above Nyquist range.

Although, if your signal is complex valued, then such symmetry won't exist, and frequencies above $\frac{f_s}{2}$ contain extra information. So you can indeed say that the bandwidth of your signal is $f_s$. This somewhat makes sense, because you might represent now both real and imaginary part of your signal on your frequency axis (twice as much space is necessary).

Here is some small example of sawtooth signal ($f_s=105Hz$), with fundamental frequency of $f_0=10 Hz$ and same signal modulated by the complex exponential with $f_{mod}=5Hz$. In such way we get time domain signal with complex values:
 $$x_{mod}(t)=x(t)e^{2\pi if_{mod}t}$$

As you can see, amplitude (phase also) spectrum for real valued signal is symmetrical around the Nyquist frequency (aliasing is also visible whilst aliased harmonics are not covered by genuine). On the contrary same sawtooth modulated by complex sinusoid (and cosinusoid) is no longer symmetrical. Probably you will notice another "feature" of this signal, while taking a closer look at peak frequencies and relating it to $f_{mod}$...

And traditionally, code:

%% Time domain parameters
fs = 105;
dt = 1/fs;
T = 10*pi;
t = 0:dt:T-dt;
N = length(t);
%% Sawtooth signal generation
f0 = 10;
x = sawtooth(2*pi*f0*t);
X = abs(fft(x))/N;
freqs = (0:N-1).*(fs/N);
subplot(211)
plot(freqs, X)
grid on
hold on
plot([fs/2 fs/2], [0 max(X)], '--r')
xlim([0 max(freqs)])
xlabel('Frequency [Hz]')
ylabel('Amplitude')
title('FT of real-valued signal')
l=legend({'sawtooth','$\frac{f_s}{2}$'});
set(l,'Interpreter','Latex');
%% Complex modulated sawtooth signal
f0_cplx = 15;
y = x.*exp(2*1i*pi*f0_cplx*t);
Y = abs(fft(y))/N;
subplot(212)
plot(freqs, Y)
grid on
hold on
plot([fs/2 fs/2], [0 max(Y)], '--r')
xlim([0 max(freqs)])
xlabel('Frequency [Hz]')
ylabel('Amplitude')
title('FT of complex-valued signal')
l=legend({'complex mod-sawtooth','$\frac{f_s}{2}$'});
set(l,'Interpreter','Latex');


A kernel density estimation based on your data:

Depend on your programming language, you may find the corresponding fitting functions with different distributions.

EDIT

I tried to use the allfitdist function with statistic toolbox in Matlab in your data fitting.

x = [24.988 25.224 25.212 25.066 25.310 24.963 24.826 24.944 25.490 25.176 24.740 24.988   ...
24.994 24.722 24.510 24.863 25.249 24.931 25.218 24.907 25.459 25.292 25.163 25.447 ...
25.072 25.360 24.976 25.103 25.237 24.845 25.237 25.206 25.109 24.969 24.976 25.194 ...
25.218 24.975 24.777 24.832 24.564 24.637 24.919 24.931 24.944 24.994 25.378 25.042...
25.182 25.151 25.097 25.103 25.030 25.018 24.783 24.722 24.667 24.963    ];

[D PD] = allfitdist(x,'PDF');

Unfortunately the fitted probability density is not good enough:

The curves for each possible distribution are all overlapped. Like I mentioned yesterday, you may need more sample number to obtain a better estimation.


Find Preamble (pattern) in signal, using Corelation I have a signal which always has the same preamble (4x 0xAA - 10101010...).

I also have a signal with additive noise, and I want to find my preamble in that a noisy signal using corelation.

U used xcorr, but I was a little confused by that. Can you help me? It looks like that signal was too noisy so that no preamble was found?

I have just used

corel=xcorr(preamble(:,2),u_noisy(:,2));

You should be able to see an peak correlation as from your plot and the length of your pattern it does not appear to be noisy enough such that you would not see a peak (unless your noise is completely in the time domain as in excessive jitter that I cannot visibly see).

I suspect that the pattern you are using is extending far beyond the length of the preamble such that your data payload is corrupting the final cross correlation result (or you are not using just the preamble pattern to correlate with your data). It goes without saying that the preamble and data need to be running at a similar rate (although you can have a significant clock offset and still correlate; correlation amplitude versus frequency offset is a Sinc function with the first null at a clock offset of 1/T where T the length of your correlation sequence; the longer the sequence the more accurate your clocks need to be). From your result if I read the numbers on the plots correctly your correlation result is close to 60 on average, and from your description you are correlating over 32 symbols (4x 0xAA).  If the sequences were not correlated, and the average for your data and pattern are 1/2, then this suggests your sampling rate is approximately 8 samples per symbol if you are indeed correlating over the full 32 symbols.

You will also benefit by removing the mean from your data (such that it is bipolar in the correlation computation) as you will see in my example below. Also as shown this does not appear to be the reason you are not getting any result (as shown by my result even with the mean present).

Trying to duplicate your conditions, here is the result I got when I correlated the 4x 0xAA pattern to a random sequence of 0's and 1's for 5000 samples: (although I reversed the correlation from what you did such that the preamble correlation would appear first):

And then with the preamble at the beginning and repeated in the middle of the random sequence (note the noise is very big in this case; as I added the random data to the preamble, just to show what you should expect to see in extremely noisy conditions:

And by removing the mean of the signals prior to correlation the detection becomes quite easy even in the presence of the very significant noise (the noise is a random pattern of 0's and 1's at sample rate):

Also for completion, below shows the result in the case of correlating over just a portion of your sequence (1x 0xAA) with the same preamble in the sequence (4x 0xAA); 20Log(4) less processing gain but showing that a clear detection is still possible with the amount of noise I had added (and such an approach may be done if the rest of the preamble is reserved for timing recovery; meaning by design we need to detect early in the preamble):

Not directly related to why you are having issue, but a challenge also to be aware of with this pattern is that while great for timing recovery (as there is a transition at your data rate on every symbol), your preamble does not have a unique autocorrelation as you can see by the cyclical pattern in the plots above.

The same processing gain in SNR is achieved as in other patterns at the position in cross-correlation where the delay between the patterns is zero , $R_{xy}(0)$, (the gain in SNR is 20Log(N) where N is the number of samples where the noise is still uncorrelated; usually due to symbol rate filtering this is the number of symbols) and if you use a circular correlation you can also achieve that SNR advantage at every offset where the sequence aligns as described above as well. For timing critical applications where it is important to know the actual position within the sequence, Barker codes or other similar pseudo-random sequences are used that have the property that correlation is maximum at shift = 0 and close to 0 everywhere else (such as maximum length sequences).  An example 11 bit Barker code is 10110111000.  If you treat the 0's as -1 you will see that this sequence will add to +11 when they are added and in a circular correlation will add to +1 everywhere else. Note however the challenge for clock recovery in such a sequence given the longer runs of ones and zeros compared to 1010101010.


Laplace transform of averaging operator I am studying dc-dc converter now. I got a problem with Laplace transform of the averaging operator as in the image below.

Can anyone help me derive the Laplace transform result $G_{av}(s)$ as in the image?

Here's the outline of the argument, feel free to fill in the details.

The averaging operator is like a convolution with a "square" pulse of height $1/T_s$ supported on the interval $[-T_s/2, T_s/2]$.

You can express the square pulse as a sum of two heaviside step functions.

Finally, recall the Laplace transform of a step function $\mathcal{L} \{H(t-T_s/2)\}(s) = \frac{e^{-sT_s/2}}{s}.$