How to measure/estimate symbol rate of a signal with unknown ( but probably QAM ) modulation? Is it possible to estimate the symbol rate of a signal whose modulation is not know exactly but it can only be guessed that it is probably some kind of QAM (probably 64 or 256).

What is the best way to get an estimate for the  symbol rate of a signal of unknown modulation?

For example by feeding the signal through some non-linearity and then looking at it's spectra, should that reveal something about the symbol rate ?

More precisely, I have a 10 MHz wide signal recorded with an SDR at 2.45GHz. It is emitted by a wireless camera.

The nonlinear trick is in fact often used to recover the symbol clock of $\frac \pi {2N}$ phase-shifted signals, e.g. BPSK, QPSK: Squaring samples $x[n]$ (or even taking $x^4$) leads to a removal of phase information, and you can then optimize timing with something that is pseudo-derivate based like Müller&Muller.

Now, the problem with QAM is usually that there's no two constellation points with the same phase; so an arbitrary amount of squaring can not reduce the number of phases to 1.

For such signals synchronization is complicated even with exact knowledge of what to expect; it's pretty common to see preamble-correlation based synchronization in receivers. Upside of that is that you can not only get timing, but also a fair amount of CSI that way, which leads to the possibility of equalization.

As a first step, however, some kind of autocorrelation-based detector surely sounds promising: Assuming that with a finite set of constellation points, there's also only finitely many sequences of symbols within the influence of a single pulse shaping filter impulse response, which also means that one might expect the signal to repeat after that amount of time. Hence, a peak on the autocorrelation¹ function of your signal probably signifies the length of $N,\, N\in\mathbb N$ symbols.

With that known, you'd still be on the lookout for $N$; remember that the spectrum of a data-modulated signal is the pulse shape convolved with a diraccomb with an inter-dirac distance of $\frac{1}{f_{symbol}}$. Hoping that symbol timing coincides with the pulse shape filter's time domain roots, you could just look at a sufficiently oversampled signal and count the spectral "peaks"; in reality, this will typically not be as clear.

A few other hints: If the spectral shape is wideband and rectangular, possibly even with a bit of "sincy" sidelobe on either side, it's probably OFDM. That'll make your life a little harder, but: OFDM typically uses cyclic prefixes, which means that you can find the OFDM symbol duration using autocorrelation; since the cyclic prefix needs to be a whole number of subsymbols (IDFT bins), you can narrow down the possible OFDM configurations quite a bit.

¹In fact, any periodicity-based estimator might work; have a look at thegr-specest toolbox. It has some parametric spectrum estimators, of which MTM might be very interesting, here.


Accelerometer BMI160 changing bias Goal: Obtain the 3-axis accelerations from the Bosch BMI160 IMU when measuring the accelerations of the bogie of a commercial train. The accelerometer used is installed in a PCB and used by an mbed microcontroller. The mbed manages to receive and save the data from the accelerometer at 50 Hz. The BMI160 has some options such as measurement range that can be changed: to avoid sensor saturation the limit is set to ±4g.

Issue: The Z-acceleration (gravity-like oriented) shows an important bias due to unknown reasons. Phenomenon in Figure 1 and 2:

The odd phenomenon can be more easily seen in this picture of the zoom Z-acceleration:

time stamps and $ Acc_z $ data can be found here: https://pastebin.com/raw/4PyAtVU8

It can be seen on it, that the average of that signal is not the gravity value. If it where due to the rotation of the object, a similar effect shall be seen in the other axis. Anyway, this does not seem to be the case, as the vehicle is a railway going forward.

If we calculate the low-pass-filtered signal ($ f_cut = 0.1\ Hz $), the signal is clearly non-oscillating. The expected result should be something more similar to what is shown in purple in the next figure:

Moreover, if we zoom in the last part of the signal, we can see that there are not any particular trends:

Question: Why the bias of the acceleration, which should be a slow-changing as a random variable, causes such a distorsion on the signal? It seems to affect only to big accelerations caused by impacts.

In next figure, you can see more information about the device layout:

The accelerometer (green component) is well fixed to the PCB with strong adhesives. The PCB itself is screwed to some wooden pieces which are also strongly sticked to the plastic box.

The whole device is fixed to the train through duck tape using metal plates for the contact. All componentes are strongly fixed one to each other and we have not observed any relative bouncing between components. The next figure shows the set up:

The blue axis are the vehicle reference system and, as it can be stated with the sign of gravity value, it does not coincide with the accelerometer axis.

I think the phenomenon may be due to an aliasing issue or just a problem with this low-cost accelerometer. Any insight from anybody?

the problem is asymmetric railing. It only occurs on the Z-axis because it is affected by gravity, when the vibration exceeds the accelerometer's range the asymmetric range causes a bias. When the vehicle is in motion vibration occurs. Solution would be with improved mechanical mounting or artificially clipping the data. See link below:

https://www.analog.com/en/technical-articles/vibration-rectification-in-mems-accelerometers.html


Is this piecewise-defined system time-invariant? \begin{equation}
  G_c \qquad y[n] =
    \begin{cases}
      u[n] & \text{if} & n<0\\
      2u[n] & \text{if} & n\ge 0
    \end{cases}
\end{equation}
I found this question in web: Whether or not is this time-invariant. The short answer is no. But I want to make sure I know how to prove this.
Let's take a certain $n<0$ then $y[n] = 2u[n]$. Now, let's put instead of $n \rightarrow (n + n_0)>0$. Since, it is time-invariant then $y[n+n_0] = u[n+n_0]$.
But $(n+n_0)>0$, so $y[n+n_0] = 2u[n+n_0]$. So the only chance is: $u[n+n_0] = 0$, but since $n$ and $n_0$ where randomly chosen then, this system can be time-varying as long as always $u[n+n_0]=0 \rightarrow u[n] =0$ and $y[n] =0$ which doesn't hold as a system. So it is not time-varying. Does this makes sense?

Let us suppose that the system with input $x$ is defined by:
\begin{equation}
  y[n] =
    \begin{cases}
      x[n] & \text{if} & n<0\\
      2x[n] & \text{if} & n\ge 0
    \end{cases}
\end{equation}
As answered by @TimWescott, the two-part  definition suggests a lack of time invariance. Finding a counter-example could be a good start. For instance a two-part signal whose behavior when multiplied by 2 is visible. Indeed a modification of the Heavide jump is interesting:
\begin{equation}
  x_0[n] =
    \begin{cases}
      -1& \text{if} & n<0\\
      \phantom{-}1 & \text{if} & n\ge 0
    \end{cases}
\end{equation}
along with the right-shifted version:
\begin{equation}
  x_1[n] =
    \begin{cases}
      -1 & \text{if} & n<1\\
      \phantom{-}1 & \text{if} & n\ge 1\,.
    \end{cases}
\end{equation}
You get:
\begin{equation}
  y_0[n] =
    \begin{cases}
      -1& \text{if} & n<0\\
      \phantom{-}2 & \text{if} & n\ge 0
    \end{cases}
\end{equation}
along with the output for the right-shifted version:
\begin{equation}
  y_1[n] =
    \begin{cases}
      -1 & \text{if} & n<0\\
      -2 & \text{if} & n=0\\
      \phantom{-}2 & \text{if} & n\ge 1\,.
    \end{cases}
\end{equation}
which are not time-shifted versions of each other.


FFTs of a complex signal - separating the real and imaginary parts I have a complex time varying signal at a single frequency x = a + jb where a represents the contribution from the cosine basis function and b represents the contribution from the sinusoid basis function.

I am trying to understand how the output differs if I was to take a complex FFT of x and inspect the real and imaginary components compared to taking two FFTs of the real and imaginary parts separately. I naively expected these to be equal when a=b. What am I missing?

% signal parameters
fs    = 1e3;
f0    = 50;
t     = (0:999)/fs;
wt    = 0.5; %weighting between cosine and sin basis functions
NFFT  = 1024;
freq  = linspace(-fs/2,fs/2, NFFT);
freq2 = linspace(0,fs/2, NFFT/2);
wnd   = hanning(length(t)).';

% Complex signal
x = wt*cos(2*pi*f0.*t) + (1-wt)*1j*sin(2*pi*f0.*t);

% 1 - Take FFT of complex signal and split into re/im components
X   = fft(x.*wnd, NFFT);
Xre = abs(X(1:NFFT/2));
Xim = abs(fliplr(X(NFFT/2+1:end)));

% 2 - Split the signal into Re/Im and compute FFTs
Yre = abs(fft(real(x).*wnd, NFFT));
Yim = abs(fft(imag(x).*wnd, NFFT));

close all;
figure;
ax(1)=subplot(221);
plot(freq2, Xre, 'b');
legend('Cmplx Re');
ax(2)=subplot(222);
plot(freq2, Yre(1:NFFT/2), 'r');
legend('Split Re');
ax(3)=subplot(223);
plot(freq2, Xim, 'b');
legend('Cmplx Im');
ax(4)=subplot(224);
plot(freq2, Yim(1:NFFT/2), 'r');
legend('Split Im');
linkaxes(ax,'xy');

figure;
plot(freq, fftshift(abs(X)));

Below is a chart I had of "Universal Fourier Transform Properties", that apply in either direction (going from time to frequency or going from frequency to time).  For example, a signal that is periodic in one domain, will be discrete in the other: a digitized time domain signal becomes periodic in frequency (such that we can concern ourselves with the spectrum from $-F_s/2$ to $+F_s/2$ only. Similarly a signal that is periodic in time will only have discrete frequencies (at the fundamental repetition rate and its harmonics in frequency).

Specific to your case I draw attention to the properties when a signal is only real, and when a signal in only imaginary:

This may help you to see what is going on, specifically the FT of a signal that is only real in the time domain, and the FT of a signal that is only imaginary in the time domain. A signal that is only real will have a spectrum that is conjugate symmetric: the positive half of the spectrum is equal in magnitude to the negative half but opposite in phase (Meaning the real portion of the FFT is "even" and the imaginary portion of the FFT is "odd", in the same fashion that a cosine in time is an even function and a sine in time is an odd function. Similarly a signal that is only imaginary will be have the FFT real components to be odd and the FFT odd components to be even.

So hopefully you see the difference between taking the fft of b, and taking the fft of jb, where to be clear b is real, so jb is imaginary (we could complicate this by allowing a and b to be complex numbers but my guess in your case they are indeed real).

This is also interesting in showing how causal time domain signals MUST be complex in frequency, and the imaginary and real components are related by the Hilbert transform: if you consider the real and imaginary components of the spectrum seperately, we see via the odd/even relationship that everything for t<0 cancels while everything for t>0 adds.

This last point is further elaborated by Sir Robert Bristow-Johnson in this post: What is the easiest, most straight-forward way to prove this about minimum-phase filters?


Bandwidth and BER of FM0 encoding What is the bandwidth (i.e., the lowest possible bandwidth) and BER for FM0 encoding?

To be precise, this is my example of an FM0 encoded waveform:

There is always a transition between each bit. For a zero there is an additional transition in the middle. For a one there is no transition (that's similar to differential Manchester code).

From the plot above I would assume that the minimum required bandwidth is 80kHz (if the waveform would be interpreted as NRZ waveform the bandwidth would be 80kHz). In my experiments I am surprised that 40kHz seems sufficient. I had expected that FM0 would just be an inefficient encoding.

Regarding the BER, I found the article "Some Interesting Observations for Certain Line Codes With Application to RFID" which derives the BER. However, from this article it is not clear if E/N0 relates to the symbol or to the bit. Actually, for this encoding it is unclear to me what exactly is the bit and what the symbol (for me, the bit is equivalent with the symbol).

I simulated the BER in MATLAB and found a deterministic discrepancy between simulation and the formula reported in the paper. The way I decode it in my simulations: I look at each halfsymbol. If two consecutive samples are the same, the output bit is a one; otherwise zero.

Due to this observation, I wonder if I can only get the proper BER if I employ a more sophisticated decoding scheme.

40kHz would be the two sided bandwidth, in which case sampling requirement is 40khz minimum


5.1 Rear To 5.1 Side mixing matrix I am searching for a matrix definition for converting a audio signal (5.1 with rear) to a audio signal (5.1 with side).
Currently I am using the definitions from: https://msdn.microsoft.com/en-us/library/windows/desktop/ff819070(v=vs.85).aspx.
But I am missing the yellow marked matrices:

Does anybody has got an idea where to get it from?
EDIT: Would be perfect if someone would have some other matrices else like quadrophonic,...

The commonly called 5.1 format uses only surround channels, which are defined as rear/side channels in ITU-R BS 775. The case you want to deal with (turning rear surround channels to side surround channels) is therefore not explicitly defined.

Notice that in the referenced ITU document, the case of changing the number of rear/sides loudspeakers reproducing a common rear/side signal is mentioned (see Note 5) :

  If more than two rear/side loudspeakers are used, the LS signal should
  be fed to each of the side/rear loudspeakers on the left side of the
  room and the RS signal should be fed to each of the side/rear
  loudspeakers on the right side of the room. In doing so, it will be
  necessary to reduce the signal gain such that the total power emitted
  by the loudspeakers carrying the LS (or RS) signal is the same as if
  that signal had been reproduced over a single loudspeaker.

In your case, as you are keeping the same number of loudspeakers, I would simply map Rear Left (resp. Rear Right) to Side Left (resp. Side Right) with a gain of 1.0.

The case is different if you want to downmix a 7.1 format (of which several mappings are mentioned in ITU-R BS 2159, see figure 20) to a 5.1 format, where you will need to take into account which mapping is used in the 7.1 format to define downmix coefficients.

Also notice that although there are commonly used coefficients when downmixing multi-channels formats, most codecs use metadata sets to allow the content creator to specify it's own coefficients in order to preserve artistic intent (see Dolby Metadata Guide, part 4 for example).


If you would write the frequency-response in terms of projections through inproducts. How would you do it? While studying for signals and systems I realised that the most intuitive way of understanding Fourier Series for me, was to see it as a projection (through inproducts) of a signal onto the orthogonal complex exponentials.
Is it possible to view the impulse-response and frequency-response of a signal in this same manner? Looking at the definition it seems clear that there should be such an 'intuitive' understanding, but I am having difficulties seeing it..

the impulse-response and frequency-response of a signal

Only systems have impulse responses and frequency responses, signals don't. I assume that what you mean here.
A System  describes the relation ship between its input signal and its output signal. For an LTI system, that relationship can be captured through either the transfer function or the impulse response. The idea is similar

Project the signal onto the some orthogonal basis
Run the basis function through the system
Assemble output basis function into the output signal. Since it's an LTI system you can simply sum or integrate.

Of course that only makes sense if the bases functions are "easier" to run through the system than some arbitrary signal. For the transfer functions the basis functions are complex exponentials and for the impulse response they are time shifted dirac deltas.
Update
Elaborating on step 2: The basis functions have an index or parameter. For example for a complex exponential, the parameter is frequency. Now you calculate the output of the system with the basis function as a an input using the parameter as a variable. Specifically you calculate the output for a complex exponential with frequency $\omega$  as a variable. The result is a function of frequency: That's the transfer function.
So if you have a signal $x(t)$ that can be expressed as the sum of some basis functions b_k(t) like
$$x(t) = \sum a_k \cdot b_k(t) $$
You can calculate the output as a weighted sum of the outputs of the basis functions.
$$y(t) = T\begin{Bmatrix}
x(t)
\end{Bmatrix} = T\begin{Bmatrix}
\sum a_k \cdot b_k(t)
\end{Bmatrix} = \sum a_k \cdot T\begin{Bmatrix}
b_k(t)
\end{Bmatrix} $$


Sobel Vs Gaussian Derivative What would be the most significant difference when using either a sobel filter or a gaussian-derivative filter, apart from the obvious difference of the size of the 2 filters?

The Sobel kernels can also be thought of as $3 × 3$ approximations to first-derivative-of-Gaussian kernels.

The Sobel filter is nothing but jointly applying the following:

Smoothing in the perpendicular direction of the derivative such as $h_s=[1,2,1]^T$
Simple central difference in the derivative direction such as $h_d=[1,0,-1]$

Then $$h_{sobel} = h_s  h_d$$

The smoothing factor is an approximate triangle shaped filter. A Gaussian is naturally a better replacement. In fact, if larger sizes of Sobel is desired, people first smooth the image with a Gaussian filter, then apply the Sobel kernel as is. That's because:

$$
\frac{\partial}{\partial x} {(I*G)} = I*\frac{\partial}{\partial x} {(G)}
$$

This approximation enables the Sobel operator to be computed quite quickly, while Gaussian derivatives might be more accurate.


How to find the impulse response of this system with complex roots? I've been working trying to analyze a causal stable system. Hopefully a numeric example such as the one I am working on, and the problems I face during it could be useful to others.

I have this system function:

$$H(z) = \frac{1-\frac{3}{10}z^{-1}+\frac{1}{3}z^{-2}}{1-\frac{3}{5}z^{-1}+\frac{38}{75}z^{-2}+\frac{2}{5}z^{-3}}$$

So I originally thought that since the numerator has the power $z^{-2}$ and the denominator has the power $z^{-3}$ I could use partial fractions, but now I believe that I am wrong to think that. I believe this case is the numerator has a higher degree than the denominator?

I multiplied the fraction through with $\frac{z^3}{z^3}$

$$H(z)=\frac{z^3-\frac{3}{10}z^2+\frac{1}{3}z}{z^3-\frac{3}{5}z^{2}+\frac{38}{75}z+\frac{2}{5}}$$

Well now the numerator has a degree that is the same as the denominator. I think I need to use long division.

I get,

$$H(z) = 1+ \frac{\frac{9}{10}z^2-\frac{13}{75}z-\frac{2}{5}}{z^3-\frac{3}{5}z^2+\frac{38}{75}z+\frac{2}{5}}$$

Factoring the denominator,

$$H(z) = 1+ \frac{\frac{9}{10}z^2-\frac{13}{75}z-\frac{2}{5}}{(z-(\frac{2}{5}+j.711)(z-(\frac{2}{5}-j.711)(z+\frac{1}{5})}$$

After much later I find the coefficients,

$$H(z) = 1+ \frac{1.12+j.273}{z-(\frac{2}{5}+j.711)}+\frac{1.12-j.273}{z-(\frac{2}{5}-j.711)}-\frac{.3805}{z+\frac{1}{5}}$$

Writing it in terms of delays,

$$H(z) = 1+ \frac{(1.12+j.273)z^{-1}}{1-(\frac{2}{5}+j.711)z^{-1}}+\frac{(1.12-j.273)z^{-1}}{1-(\frac{2}{5}-j.711)z^{-1}}-\frac{.3805z^{-1}}{1+\frac{1}{5}z^{-1}}$$

Now taking the inverse Z-Transform

$$h[n]= \delta[n]+ (1.12+j.273)(\tfrac{2}{5}+j.711)^{(n-1)}u[n-1]+(1.12-j.273)(\tfrac{2}{5}-j.711)^{(n-1)}u[n-1]-.3805(\tfrac{-1}{5})^{(n-1)}u[n-1]$$

That was pretty difficult and time consuming for me too. I ended up find some sort of solution but I am not sure if it is correct. Is there an easier way to solve more difficult IIR systems such as these with complex roots? What do the complex coefficients mean? Is there a standard to simplifying these even more? Further more, what would the input output equation of this type be?

UPDATE:
Using part of the example kindly provided below, noting that complex conjugate pairs can be simplified further to obtain an impulse function I was able to come up with the following: (Focusing on the complex conjugate part in particular)

$$H(z) = \frac{1.12+j.273}{z-(\frac{2}{5}+j.711)}+\frac{1.12-j.273}{z-(\frac{2}{5}-j.711)}$$

Which simplifies to,

$$H(z)= \frac{(2.24)-(.9381e^{-j.811}+.9381e^{j.811})z^{-1}}{1-.8z^{-1}+.665z^{-2}}z^{-1} $$

I am going to ignore the extra factor of the delay tacked on at the end. I found a Fourier Transform Pair that particularly looks like this! Not sure if this is what should be done but it looks hopeful!

$$r^ncos(w_0)u[n]<---DTFT--->\frac{1-r \cos(w_0)z^{-1}}{1-2r\cos(w_0)z^{-1}+r^2z^{-2}}$$

rearranging some things to make it fit better with the DTFT pair,

$$H(z)= 2.24\frac{1-(.8375)\cos(.811)z^{-1}}{1-2(.58)\cos(.811)z^{-1}+(.815)^2z^{-2}} $$

Well one thing that is alarming is that my r values (according to the DTFT pair) are slightly off.

The three different values:

$r=.58$ is worrisome.

$r=.8375$

$r=.815$

The other two r values seem like it could be due to averaging errors. But $r=.58$ seems a little too off.

What I think I would do is take the average $r$, though I don't know if this is the right approach, and use that to find the new impulse function for this complex conjugate pair.

It turns out that those terms with complex conjugate scalers and exponentials can be combined.  Rather than do it to your example, I would rather do it in the Z-domain.  Let

$$\begin{align}
H(z) &= \frac{A_1}{1-p_1\,z^{-1}} + \frac{A_2}{1-p_2\,z^{-1}} \\
\\
&= \frac{A_1(1-p_2\,z^{-1})+A_2(1-p_1\,z^{-1})}{(1-p_1\,z^{-1})(1-p_2\,z^{-1})}  \\
\\
&= \frac{(A_1-A_1 p_2\,z^{-1})+(A_2-A_2 p_1\,z^{-1})}{1-(p_1+p_2)\,z^{-1} + p_1 p_2\, z^{-2}}  \\
\\
&= \frac{(A_1+A_2) - (A_1 p_2 + A_2 p_1)z^{-1}}{1-(p_1+p_2)\,z^{-1} + p_1 p_2\, z^{-2}} \\
\\
\end{align}$$

Now it turns out that there are two ways for the impulse response, $h[n] = \mathcal{Z}^{-1} \{H(z)\}$, to be total real (with no imaginary component).  One way is for both poles $p_1,p_2$ to be real and that will also require the two coefficients $A_1,A_2$ to also be real.

$$\begin{align}
h[n] &= \mathcal{Z}^{-1}\{H(z)\} \\
\\
 &= \mathcal{Z}^{-1}\left\{ \frac{A_1}{1-p_1\,z^{-1}} \right\} + \mathcal{Z}^{-1}\left\{ \frac{A_2}{1-p_2\,z^{-1}} \right\} \\
\\
&= A_1 p_1^n u[n] + A_2 p_2^n u[n] \\
\end{align}$$

The other way for $h[n]$ to be real is for both poles $p_1,p_2$ to be complex-conjugate of each other and that will also require the two coefficients $A_1,A_2$ to also be complex-conjugate.

$$\begin{align}
h[n] &= \mathcal{Z}^{-1}\{H(z)\} \\
\\
 &= \mathcal{Z}^{-1}\left\{ \frac{(A_1+A_2) - (A_1 p_2 + A_2 p_1)z^{-1}}{1-(p_1+p_2)\,z^{-1} + p_1 p_2\, z^{-2}} \right\}  \\
\end{align}$$

(i need to return to this later.)


Help for baseline drift removal/correction I have signal as seen in the picture. It has high frequency noise and drifting baseline. I am able to remove the noise in MATLAB but I am not sure what can I do about the drift.

I would like the signal to have a stable baseline. Any help or hints on how I can do that would be most helpful.

Thanks

Like I already mentioned in the comment, you might find MATLAB's detrend function useful. Generally it is useful for removal of linear (or piecewise linear) trend:

y = detrend(x);

Another solution, that should work in your case, is to use the polynomial fitting. This will remove the non-linear trend (providing you want that). A minimal code example:

% vectors t,x are given (time and the signal)
degree = 5; % to be adjusted
[p, S, mu] = polyfit(t, x, degree);
xp = polyval(p, t, [], mu);
y = x - xp;

Personally I do not recommend any Moving Average filtering, since it is introducing delay and can distort your signal.


How to efficiently calculate only the low coefficients of a zero-padded FFT I've got an algorithm that zero pads a sequence to 4N, does an FFT, and only uses the lowest frequency N points out of the generated 4N.

This seems like a lot of wasted work, any ideas how this can be done faster?

If you only a few bins then the following may be very efficient for you:
1. Simply do the DFT of at each frequency you need.
2. Use the Goertzel algorithm for each frequency in question.


Direct and Transpose Polyphase Multirate Processing Polyphase implementations of upsampling/ interpolation and downsampling/ decimation, after having invoked the Noble identities, are presented as follows (taken from Proakis):

(Three-Channel Polyphase Decimation)

(Three-Channel Polyphase Interpolation)
The decimation always seems to use the direct form (delay before filtering), while the interpolation always seems to use the transposed form (delay after filtering). I understand that these two forms are equivalent (or believe that's the case). But I'm unable to find an explicit discussion on why this is canonical; why, for instance, doesn't the interpolation use the direct form? I've seen a couple of suggestions that it has something to do with the critical path, but I'm not clear on precisely what this means.
Could anyone clear this up? I suspect I'm being a bit obtuse...

The antialiasing filters here are FIR so there's no difference in the order of operations. If they were IIR there would be implications, see NUMERICAL ROBUSTNESS OF TDF-II. I think the classical representation of polyphase decimators and interpolators here is an elegant way to say that filtering has to occur at the lowest rate (before interpolation or after decimation).


A snippet from the MATLAB documentation of the plot command:

  plot(Y) creates a 2-D line plot of the data in Y versus the index of each value...If Y is complex, then the plot function plots the imaginary part of Y versus the real part of Y, such that plot(Y) is equivalent to plot(real(Y), imag(Y)).

The circular plot you're seeing is corresponding to the real part on the x-axis and the imaginary part of the y-axis.

If you're doing plot(s), then MATLAB is plotting each column in s.


Alias-free digital nonlinear filter design @Jazzmaniac has a good answer to the question of how to design an alias-free digital nonlinear time-invariant filter here: https://dsp.stackexchange.com/a/28787/18276

Basically, according to that answer, a digital nonlinear time-invariant filter is alias free if and only if it commutes with subsample translations. Meaning that it doesn't matter whether you filter and then translate, or translate and then filter. Sinc interpolation is required for perfect subsample translations, but of course you can always use a finite interpolator that is good enough.

This question is to elaborate:

How can we see the link between subsample translation invariance and aliasing?
Is there any easy way to see what these filters look like?
Do the filters have some standard form they can be put in?
Do we know what the alias-free version of the monomials look like? (i.e. the alias-free version of $y[t]=x[t]^n$ for some positive natural number $n$)
Are there any good references or published works on the topic of alias-free nonlinear filter design?

First, allow me to address the subsample shift property in relation to non-linear signal mappings. It is fairly straightforward to see time shifting is not as simple a property as for linear systems. Consider a discrete time signal given by the sequence $$\dots ,1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, \dots$$
This sequence is invariant under memoryless nonlinearities of the form $x\mapsto x^n$ for natural $n$. If the above sequence is shifted by a fraction of the sampling interval, the sample values will not be $0$ and $1$ anymore, and the resulting sequence will lose its invariance under the nonlinear map given.

Enforcing this invariance and allowing only nonlinear maps that preserve it is equivalent to removing aliasing, as I will be showing here when I find a little time to expand upon my answer.

Edit: Some more details.

For simplicity, we will look at a discrete time signal centred around $t=0$ with a sampling interval of $T=1$. The discrete signal can be expanded using a power series $$x(t) = \sum_{n=0}^{2N} D_{N}^{(n)} \frac{t^n}{n!}$$
which describes $x[t]$ in the interval $t\in [-N,N]$. The coefficients $D_{N}^{(n)}$ are chosen so that the polynomial is the minimal interpolating polynomial on this interval. They coincide with the $n$-th discrete derivatives of $x[t]$ at $t=0$, also on this interval. The interpolating nature of this expansion makes it a natural linear homomorphism to the space of continuous time signals. It is worth noting, that in the limiting case $N\to\infty$, the interpolation approaches the bandlimited interpolation of the $\mathrm{sinc}$ kernel.

Applying a differential time shift $dt$ to a smooth continuous time signal $x(t)$ can be written as $$x(t-dt)=x(t) - dt \frac{\partial}{\partial t}x(t)=(1-dt \frac{\partial}{\partial t})x(t)$$
For finite shifts $\delta t$, we can concatenate many smaller shifts and take the limit $$x(t-\delta t)=\lim_{n\to\infty}\left(1-\frac{\delta t}{n}\frac{\partial}{\partial t}\right)^n \,x(t)=\exp\left(-\delta t \frac{\partial}{\partial t}\right)\,x(t)$$
Therefore the linear operator $S(\delta t)=\exp\left(-\delta t\frac{\partial}{\partial t}\right)$ shifts smooth continuous time signals.

The interpolation polynomial above is also an expansion of the $\exp$ function into a power series, with the discrete derivative taking the place of the partial derivative of the shift operator. For finite orders of the interpolation polynomial, the continuous time shift operator therefore approximates shifts of the interpolating function. This approximation becomes exact in the limit of infinite order.

With this understanding, we can calculate an (approximately) alias free memoryless nonlinearity $f$ acting on $x[t]$. We only need to evaluate $f(x[0])$, all other times follow by integer sample shifts.

With the interpolated continuous-time signal $x(t)$ approximating a band-limited reconstruction of $x[t]$, we can apply the nonlinearity in continuous time and then band-limit the result using an approximation of a band-limiting kernel $b(t)$ to avoid aliasing. A sufficient condition for a feasible symbolic evaluation of the band-limited result is that $f$ and $b$ be polynomials. Then $f(x(t))$ is a polynomial and we can directly calculate

$$y(0) = \int_{-N}^{N} f\left(x(t)\right) b(t) dt$$

which, in the limit of large orders $N$, achieves both full translation invariance and perfect alias rejection. This is not a proof that both are equivalent, but a good starting point to understand how these two properties are linked.

Suitable choices for $b$ include polynomial expansions of the $\mathrm{sinc}$ function. For example

$$b_N(t)=\frac{(N^2-t^2) \prod_{n=1}^{2N}(t^2-n^2) }{N^2 \prod_{n=1}^{2N}n^2}$$

for an approximately bandlimited kernel on the interval $t \in [-N,N]$.

This much must suffice as theoretical motivation for now.

Example

The simplest, non-trivial example is that with the minimal neighbourhood involvement and a crude approximation to a band-limited kernel. Do not expect good anti-aliasing properties from it. It's only here to demonstrate the general procedure of creating an anti-aliased memoryless nonlinearity.

We use the lowest possible order, $N=1$, and arrive at the expansion

$$x(t) = D_1^{(0)} t^0 + D_1^{(1)} t^1 + D_1^{(2)} \frac{t^2}{2}$$

where

$$
D_1^{(0)} = x[0]\\
D_1^{(1)}=\frac{1}{2}(x[1] - x[-1])\\
D_1^{(2)}=x[1]-2x[0]+x[-1]$$

and as a single expression

$$x(t)= x[0] + \frac{1}{2}(x[1]-x[-1])\, t + \frac{1}{2}(x[1]-2x[0]+x[-1])\,t^2$$

The nonlinearity is assumed to be a monomial $x\mapsto x^k$ and we take the $b_1$ from above and get

$$y_k[0] = \int_{-1}^{+1} \left(x[0] + \frac{1}{2}(x[1]-x[-1])\, t + \frac{1}{2}(x[1]-2x[0]+x[-1])\,t^2\right)^k \, \cdot \frac{1}{4}(4-t^2)(t^2-1)^2 dt$$

We can evaluate this expression for $k=2$ and also remove the $t=0$ simplification made earlier to arrive at a nonlinear filter

$$y_2[n] = \frac{4}{3465}(688 x[n]^2 + 40 x[n-1]^2 + 40 x[n+1]^2 - 41 x[n-1] x[n+1] + 82 x[n](x[n+1]+x[n-1]))$$

Generalising for non-smooth nonlinearities

The argument so far has required that the non-linearity can be well approximated by a polynomial. If that is not the case, then the integral will generally be harder to evaluate and the existence of a closed form solution is not even guaranteed. This is where the equivalence of commutativity of the filter with sub-sample shifts comes in.

The most general form of a nonlinear filter on a neighborhood $[-N,N]$ around $t=0$ is that of a nonlinear map

$$ y[0] = F( x[-N],x[-N+1],\dots,x[N-1],x[N] )$$

For a memoryless non-linearity, we want to map constant input signals to constant output signals, according to the non-linear transfer function $f$. If the constant input is $X$, then we have the condition

$$F(X,X,\dots,X) = f(X)$$

and from the shift-invariance we have the condition

$$ F( S(\delta t) x[-N], S(\delta t) x[-N+1],\dots, S(\delta t)x[N-1],S(\delta t) x [N]) = S(\delta t)F( x[-N],x[-N+1],\dots,x[N-1],x[N] )$$ for all $\delta t \in \mathbb{R}$

In general, there is no symbolic solution for this problem. It is however approachable with numerical optimization methods. In many cases, the free parameters can be reduced further by restricting the form of $F$.

I believe I should have answered all your questions with exception of the request for literature pointers. I am not aware of any. I don't know if this theory has ever been presented. But I believe the idea is not too difficult to come up with, so it has probably been done before. If you find something, please do let me know.


IQ demodulation from specific rate We have an acquisition hardware in our lab that acquires and digitizes data from a source generating band-limited signal in the MHz range. It samples the signal at a rate of 4 times the center frequency $f_0$ of the signal. In the documentation, they say that this specific sample rate makes it easier to convert the RF samples to demodulated IQ samples.
How would $f_s=4f_0$ make the demodulation easier? I have been using Hilbert transform to achieve the demodulation (the old-fashioned way). But, if there is any shortcut offered by this specific sample rate, that would be really helpful.
Thanks.

To get IQ data from a bandlimited signal centered at f0, with a bandpass well inside f0-f0/2 to f0+f0/2, and using a sample rate of fs=4*f0, you can simply keep 2 successive samples out of every 4, and ignore the rest.  Take every 4th sample, put them in the I channel, take the samples right after those and put them in the Q channel, and throw away the 3rd and 4th sample out of every 4 samples (I Q ignore ignore I Q ignore ignore I Q ignore ignore ... etc.).  The result is an IQ signal at a sample rate of f0.


How can I filter a noise with pattern type in audio? I have a set of audio files rescued from a source of phonograph record like a gramophone. Very old family stuff records from my elders back in year 1946. A relative has digitalized these with a turntable and gave the records to me in lossless WAV Format. The audios are just conversations of my relatives in cotidiane situations. So it is just vocals materials but mixed with a heavy quantity of noise. Noise is mostly concentrated in the mid-high frequencies, like 1Khz-10Khz.

So the records are VERY noisy. SNR is like 0dB or worst. My mission is to filter out the noise.

One advantage is that noise exhibits a cyclical pattern. I tried some classic vinyl restoration stuff like decracking o declicking with Sound Forge, DART 24 Pro. Results were unsuccessful. The best result I had is the "Noise reduction" effect in Audacity, which did output relatively good. But this last one just filters the frequencies where the noise is with narrow FFT. So the vocals content loses the high-mid frequencies. Leaving only mid or mid-bass content. The noise got down but voices are unintelligible.

So I think perhaps I could turn to a Artificial Intelligence solution?. but I have no clue from where can I start. I have some understading of theory of AI, the Perceptron, neural networks, Backpropagation, but I never applied them to audio. Perhaps a Convolution Neural network could be a solution, but how could I train it?. If helps, I have programming skills in C# and knowledge about digital linear filters like IIR and FIR. But equalization based solutions *are not, what I need *.

Below I offer an example track for your consideration.

Audio track for example

What can I do for filter this kind of cyclical pattern noise?

Simple idea: let some pros do the work. Send the audio through a cellphone to another one via a WhatsApp, Duo or whatever voice call app, even classic network call works, dealer's choice. The noise reduction algorithms in cellphones are designed to do exactly what you want to do. Reduce noise while conserving human speech. The outcome will depend to some degree on the market, the phone's firmware was designed for. US customers like it, when the background noise is cut out completely, accepting some degradation of speech. European customers prefer more naturally sounding outcomes with some noise left.

I think this would be worth a shot.


Why does cel-shading sometimes look similar to rotoscoping? Why does cel-shading sometimes look similar to rotoscoping?

Rotoscoping: A Scanner Darkly:

Cel-shading: Wikipedia: Cel-shading, Teapot Example:

I checked out the Wikipedia page for rotoscoping and it says one way it was done is that a live action image was shown behind frosted glass and an animator would draw it.  The frosted glass detail obfuscation coupled with the fact that it was a cartoonist means to me that you are going to get solid colored areas as shading, like you'd see in a cartoon.  Cell shading is basically trying to make the same effect but using 3d rendered graphics as a source instead of an actual photographed image.  Seems like they are the same effect, just with different names!


Design a Discrete-Time System Could you help me with this theoretical question.

The purpose is to design a discrete-time system. What it does is, it finds the maximum value in a one-dimensional sequence whose elements consist only of "non-negative" real numbers and contains only a single global maximum with no local maxima.

Thanks in advance.

Since you have defined your signal $x[n]$ to be a set of non-negative numbers with a single global maxima, the algorithm to find the maximum is quite simple.

compute the first (causal , backward) difference
$$y[n] = x[n] - x[n-1]$$

Then look for the sign change of $y[n]$; i.e., at the extreme point you will have the sign of $y[n]$ changing from positive to negative,  assuming you had a single global maxima.

$$z[n] = \text{sign_change_detector} \{ y[n] \} $$

Don't forget to check the end points as well. As they are also the candidates to be maximum in addition to those critical point(s).


Determining the mean and standard deviation in real time What would be the ideal way to find the mean and standard deviation of a signal for a real time application.  I'd like to be able to trigger a controller when a signal was more than 3 standard deviation off of the mean for a certain amount of time.

I'm assuming a dedicated DSP would do this pretty readily, but is there any "shortcut" that may not require something so complicated?

There's a flaw in Jason R's answer, which is discussed in Knuth's "Art of Computer Programming" vol. 2. The problem comes if you have a standard deviation which is a small fraction of the mean: the calculation of E(x^2) - (E(x)^2) suffers from severe sensitivity to floating point rounding errors.

You can even try this yourself in a Python script:

ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]]
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2

I get -128.0 as an answer, which clearly isn't computationally valid, since the math predicts that the result should be nonnegative.

Knuth cites an approach (I don't remember the name of the inventor) for calculating running mean and standard deviation which goes something like this:

 initialize:
    m = 0;
    S = 0;
    n = 0;

 for each incoming sample x:
    prev_mean = m;
    n = n + 1;
    m = m + (x-m)/n;
    S = S + (x-m)*(x-prev_mean);

and then after each step, the value of m is the mean, and the standard deviation can be calculated as sqrt(S/n) or sqrt(S/n-1) depending on which is your favorite definition of standard deviation.

The equation I write above is slightly different than the one in Knuth, but it's computationally equivalent.

When I have a few more minutes, I'll code up the above formula in Python and show that you'll get a nonnegative answer (that hopefully is close to the correct value).

update: here it is.

test1.py:

import math

def stats(x):
  n = 0
  S = 0.0
  m = 0.0
  for x_i in x:
    n = n + 1
    m_prev = m
    m = m + (x_i - m) / n
    S = S + (x_i - m) * (x_i - m_prev)
  return {'mean': m, 'variance': S/n}

def naive_stats(x):
  S1 = sum(x)
  n = len(x)
  S2 = sum([x_i**2 for x_i in x])
  return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }

x1 = [1,-1,2,3,0,4.02,5]
x2 = [x+1e9 for x in x1]

print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)

print "stats:"
print stats(x1)
print stats(x2)

result:

naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}

You'll note that there's still some rounding error, but it's not bad, whereas naive_stats just pukes.

edit: Just noticed Belisarius's comment citing Wikipedia which does mention the Knuth algorithm.


Wavelet decomposition For the following code,

X=[2 2 3 100 4 0 98 100 90 2 3 67 98 0 6 6 89 9 21 78]
[C,L]=wavedec(X,N,'db1');

where N is the level , I get n(approximation and detail coefficients)>n(X)

What does this signify?
Why does this increase happen?

This may be caused by "border effects" related to subsampling and signal length with perfect reconstruction filter banks.

For 2-band wavelets, there are standard  non-expansive solutions (for instance with periodic border extension) when $2^N$ divides the length $L$ of your signal. In your case, one or two levels work without length increase, becauses $20=20/2+20/2 = (10/2+10/2)+10$. Then, starting at level $3$, you cannot divide $5$ into equal parts, since $2^3 = 8$ does not divide 20.

So:

It signifies that you have decomposed too deep to keep the number of samples constant.
The increase happens because of your signal length, your decomposition  level and your border extention.

For some applications (outside compression) this is not really an issue.


the sub-range of circular and linear convolution circular convolution  $x_{_3p}[n]$ = $x_1[n]~\circledast_N ~x_2[n]$

is a period version of the linear convolution  $x_{_3p}[n]=x_1[n] * x_2[n]$

The length of $x_1[n]$ and $x_2[n]$ are $L$ ($n\in[0,\ldots,L-1]$) and $P$ ($n\in[0,\ldots,P-1]$) points, respectively.

The minimum value of $N$ that makes $x_{_3p}[n] = x_{_3}[n]$ for $n \in [0,\ldots N-1]$
is when $N\geq L+P-1$, right?

My question is: If $N=L$, what is the sub-range of $[0,\ldots,N-1]$ that $x_{_3p}[n]=x_{_3}[n]$?

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

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

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

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

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

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

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

$$\\$$

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

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

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

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

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


Do integration/differentiation processes work as simple filters? How these processes are different from simple IIR 1 order filtering, FIR filters in terms of amplitude and phase characteristics?

Yes, integration and differentiation can be linear filters.
You can start from laplace properties that say:

$ \int_{0}^{t} {x(t)dt} \longrightarrow \frac{X(s)}{s} \\ \frac{d}{dt}x(t) \longrightarrow sX(s) $

So you can find transfer function of integration and differentiation:

$ H_{INT}(s) = \frac{1}{s} \\ H_{DIFF}(s)=s $

You can convert these transfer function in digital IIR filters for example by bilinear transform or others digitization techniques. However you should notice that $ H_{DIFF}(s) $ is not causal, then you must add a pole to the transfer faction far from useful signal frequency, and it begin:

$ H_{DIFF_{causal}}(s)= \frac{s}{\alpha s + 1}$ where $\alpha$ is the time constant of the derivative filter a small real $>0$.

Using bilinear transform $H_{INT}(z)$ begin a trapezoidal integrator, using Euler transofrm $H_{INT}(z)$ begin a rectangular integrator, you can see this difference and the digitization in the MATLAB PID page.

I wrote a simple PID program in C that computes these operations using Euler transform, a PID in closed loop doesn't need to be accurate so Euler works well.

In literature there are a lot of ways to implement derivative filters (IIR and FIR) that workaround in elegant way anti-causal problem, however in a lot of situations you can simply digitize the analog transfer functions to make IIRs and FIRs.


What Properties Make Certain Wavelets "Better" Than Others In Image Compression? I am trying to teach myself more about image compression using the wavelet transform method.  My question is: What is it about certain wavelets that make them preferable when compressing images?  Are they easier to compute? Do they produce smoother images? Etc...

Example: JPEG 2000 uses the Cohen-Daubechies-Feauveau 9/7 Wavelet...why this one?

Overview

The short answer is that they have the maximum number of vanishing moments for a given support (i.e number of filter coefficients). That's the "extremal" property which distinguishes Daubechies wavelets in general. Loosely speaking, more vanishing moments implies better compression, and smaller support implies less computation. In fact, the tradeoff between vanishing moments and filter size is so important that it dominates the way that wavelets are named.  For example, you'll often see the D4 wavelet referred to either as D4 or db2.  The 4 refers to the number of coefficients, and the 2 refers to the number of vanishing moments. Both refer to the same mathematical object. Below, I'll explain more about what moments are (and why we want to make them disappear), but for now, just understand that it relates to how well we can "fold up" most of the information in the signal into a smaller number of values. Lossy compression is achieved by keeping those values, and throwing away the others.

Now, you may have noticed that CDF 9/7, which is used in JPEG 2000, has two numbers in the name, rather than one. In fact, it's also referred to as bior 4.4. That's because it's not a "standard" discrete wavelet at all. In fact, it doesn't even technically preserve the energy in the signal, and that property is the entire reason people got so excited about the DWT in the first place! The numbers, 9/7 and 4.4, still refer to the supports and vanishing moments respectively, but now there are two sets of coefficients that define the wavelet. The technical term is that rather than being orthogonal, they are biorthogonal. Rather than getting too deep into what that means mathematically, I'll just review the factors which led to using non-energy-preserving biorthogonal wavelets in the first place.

JPEG 2000

A much more detailed discussion of the design decisions surrounding the CDF 9/7 wavelet can be found in the following paper:

  Usevitch, Bryan E. A Tutorial on Modern Lossy Wavelet Image
  Compression: Foundations of JPEG 2000.

I'll just review the main points here.

Quite often, the orthogonal Daubechies wavelets can actually result in increasing the number of values required to represent the signal. The effect is called coefficient expansion. If we're doing lossy compression that may or may not matter (since we're throwing away values at the end anyway), but it definitely seems counterproductive in the context of compression. One way to solve the problem is to treat the input signal as periodic.
Just treating the input as periodic results in discontinuities at the edges, which are harder to compress, and are just artifacts of the transform. For example, consider the jumps from 3 to 0 in the following periodic extension: $[0,1,2,3] \rightarrow  [...0,1,2,3,0,1,2,3,...]$.  To solve that problem, we can use a symmetric periodic extension of the signal, as follows: $[0,1,2,3] \rightarrow [...,0,1,2,3,3,2,1,0,0,1...]$. Eliminating jumps at the edges is one of the reasons the Discrete Cosine Transform (DCT) is used instead of the DFT in JPEG. Representing a signal with cosines implicitly assumes "front to back looping" of the input signal, so we want wavelets which have the same symmetry property.
Unfortunately, the only orthogonal wavelet which has the required characteristics is the Haar (or D2, db1) wavelet, which only as one vanishing moment. Ugh. That leads us to biorthogonal wavelets, which are actually redundant representations, and therefore don't preserve energy. The reason CDF 9/7 wavelets are used in practice is because they were designed to come very close to being energy preserving. They have also tested well in practice.

There are other ways to solve the various problems (mentioned briefly in the paper), but these are the broad strokes of the factors involved.

Vanishing Moments

So what are moments, and why do we care about them? Smooth signals can be well approximated by polynomials, i.e. functions of the form:

$$a + bx + cx^2 + dx^3 + ...$$

The moments of a function (i.e. signal) are a measure of how similar it is to a given power of x.  Mathematically, this is expressed as an inner product between the function and the power of x.  A vanishing moment means the inner product is zero, and therefore the function doesn't "resemble" that power of x, as follows (for the continuous case):

$$\int{x^n f(x) dx = 0 }$$

Now each discrete, orthogonal wavelet has two FIR filters associated with it, which are used in the DWT. One is a lowpass (or scaling) filter $\phi$, and the other is a highpass (or wavelet) filter $\psi$. That terminology seems to vary somewhat, but it's what I'll use here. At each stage of the DWT, the highpass filter is used to "peel off" a layer of detail, and the lowpass filter yields a smoothed version of the signal without that detail. If the highpass filter has vanishing moments, those moments (i.e. low order polynomial features) will get stuffed into the complementary smoothed signal, rather than the detail signal. In the case of lossy compression, hopefully the detail signal won't have much information in it, and therefore we can throw most of it away.

Here's a simple example using the Haar (D2) wavelet. There's typically a scaling factor of $1/\sqrt{2}$ involved, but I'm omitting it here to illustrate the concept. The two filters are as follows:
$$ \phi = [1,1] \\  \psi = [1,-1]$$

The highpass filter vanishes for the zero'th moment, i.e. $x^0 = 1$, therefore it has one vanishing moment. To see this, consider this constant signal: $[2,2,2,2]$.  Now intuitively, it should be obvious that there's not much information there (or in any constant signal).  We could describe the same thing by saying "four twos". The DWT gives us a way to describe that intuition explicitly. Here's what happens during a single pass of the DWT using the Haar wavelet:

$$
[2,2,2,2] \rightarrow_{\psi}^{\phi} \left\{ \begin{array}{rr}
\left[2 + 2, 2 + 2\right] = \left[4,4\right] \\
\left[2-2,2-2\right] = \left[0,0\right]
\end{array}\right.
$$

And what happens on the second pass, which operates on just the smoothed signal:

$$
[4,4] \rightarrow_{\psi}^{\phi} \left\{ \begin{array}{rr}
\left[4 + 4\right] = \left[8\right] \\
\left[4-4\right] = \left[0\right]
\end{array}\right.
$$

Notice how the constant signal is completely invisible to the detail passes (which all come out to be 0).  Also notice how four values of $2$ have been reduced to a single value of $8$.  Now if we wanted to transmit the original signal, we could just send the $8$, and the Inverse DWT could reconstruct the original signal by assuming that all the detail coefficients are zero. Wavelets with higher-order vanishing moments allow similar results with signals that are well approximated by lines, parabolas, cubics, etc.

Further Reading

I'm glossing over a LOT of detail to keep the above treatment accessible. The following paper has a much deeper analysis:

  M. Unser, and T. Blu, Mathematical properties of the JPEG2000 wavelet
  filters, IEEE Trans. Image Proc., vol. 12, no. 9, Sept. 2003,
  pg.1080-1090.

Footnote

The above paper seems to suggest that the JPEG2000 wavelet is called Daubechies 9/7, and is different from the CDF 9/7 wavelet.

  We have derived the exact form of the JPEG2000 Daubechies 9/7 scaling
  filters... These filters result from the factorization of the same
  polynomial as $Daubechies_{8}$ [10]. The main difference is that the
  9/7 filters are symmetric. Moreover, unlike the biorthogonal splines
  of Cohen-Daubechies-Feauveau [11], the nonregular part of the
  polynomial has been divided among both sides, and as evenly as
  possible.

  [11] A. Cohen, I. Daubechies, and J. C. Feauveau, “Biorthogonal bases
  of compactly supported wavelets,” Comm. Pure Appl. Math., vol. 45, no.
  5, pp. 485–560, 1992.

The draft of the JPEG2000 standard (pdf link) that I've browsed also calls the official filter Daubechies 9/7. It references this paper:

  M. Antonini, M. Barlaud, P. Mathieu, and I. Daubechies, “Image coding
  using the wavelet transform,” IEEE Trans. Image Proc. 1, pp. 205-220,
  April 1992.

I haven't read either of those sources, so I can't say for sure why Wikipedia calls the JPEG2000 wavelet CDF 9/7. It seems like there may be a difference between the two, but people call the official JPEG2000 wavelet CDF 9/7 anyway (because it's based on the same foundation?).  Regardless of the name, the paper by Usevitch describes the one that's used in the standard.


Is gradient vector of pole zero carry usful information? Consider a transfer function (TF) with $n$ number of poles $(p_1,..p_n)$ and $m$ number of zeros $(z_1,..z_n)$.  One can write the magnitude of the frequency responce of the TF in terms of poles and zeros.

This means one can take partial derivative of the magnitude of the TF with respect to each pole and zero for a paticular frequency. Now consider $f_1$ is the frequency at which TF is deviated most from the desired responce. Thus one can change the poles and zeros based on the gradient of the magnitude of TF at $f_1$ and optimize the poles and zeros for $f_1$. Can one can keep doing this procedure for each frequency until a desired result is obtained? Is this a valid process of filter design?

Any comment would be appreciated.

Does gradient vector of pole zero carry useful information?

Yes. The partial derivatives can be used in creating iterative search algorithms for fitting IIR filters to arbitrary targets. Examples of algorithms that use the derivatives are Steepest Descent or Conjugate Gradient.

It's not a trivial process though and there are lot of details to be worked out.

Poles, zeros and transfer function targets are all complex, so you either need to work with complex derivatives are use a suitable real-valued representation (real/imag, magnitude/phase, etc.).
It's important to formulate a suitable error criteria that properly
reflects the requirements and trade-offs for your specific
application,
The search algorithms often get stuck in local minima
Log spaced frequency grids may have to be pre-warped to straighten out the error surface a bit.
Many search algorithm have "heuristic" parameters that are hard to get right


Questions regarding OFDM, FFT, IFFT The questions below are for OFDM, IFFT, FFT operations.

Example 1

In OFDM modulation, a stream of parallel symbols are input to an IFFT operation of size $M$ for example and the output is time domain samples of size $M$.

My question is:

Let us take for example, the fact that the input symbols to the FFT operation are QAM symbols. Then, at the output of IFFT I receive $M$ samples. Each sample is weighted sum of modulated subcarriers.

Does it mean that each QAM symbol (input of IFFT) is sent out on ONE subcarrier?

Example 2

I have seen a paper that take the time domain samples of an impulse response of channel of length $T_{CH}$ puts them into  M-FFT operation and the result was $H(k)$ for $k\in[1:M]$ which the authors state are the frequency response of the $k$th subcarrier. Why is that the FFT generated the subcarrier frequency response? According to Example 1, each sample should be weighted sum of all input ?

Can anyone explain to me whats the connection between two examples?

Thanks and looking forward for your clarification.

To your first question, yes, each symbol is sent on a separate subcarrier; that's the basic idea of OFDM.

I'm not sure where your confusion lies with your second question. It's well-known that the Fourier transform of a system's impulse response gives its frequency response. In this case, the system is the communication channel. If you choose the sample rate of your impulse response measurement properly, and use the appropriately-sized FFT, you can get a frequency response measurement whose bins line up with the subcarrier frequency grid. This can be used for equalizing the channel.


Ideas on Matrix Factorization / Transformations for $ {L}_{1} $ Minimization I am starting with a typical $\ell_1$ basis pursuit problem:

$$
\min_{\mathbf{x}} \Vert \mathbf{x} \Vert_1 \quad \mathrm{s.t.} \quad \Vert \mathbf{ERx} - \mathbf{y} \Vert_2 \leq \epsilon,
$$

where $\mathbf{R}\in\mathbb{C}^{M \times P}$, and $\mathbf{E}$ is an $M' \times M$ row selection matrix that randomly chooses $M'<P$ rows from $\mathbf{R}$. It is also known that $\operatorname{rank}(\mathbf{ER})=M'$, so that the constraint system is full-rank but under-determined.

Since both $M$ and $P$ are on the order of $10^7$, the above problem is impractical because of memory and computational limitations. An algorithm exists that will let me efficiently compute $\mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{y}$, so I would like to see if I can combine this with a sparse approximation (described below) to solve the following equivalent problem:

$$
\min_{\mathbf{x}} \Vert \mathbf{x} \Vert_1 \quad \mathrm{s.t.} \quad \Vert \mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{ERx} - \mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{y} \Vert_2 \leq \epsilon,
$$

Since $\mathbf{E}$ is a row-selection matrix, the product $\mathbf{E^{\mathrm{H}}}\mathbf{E}$ is a diagonal matrix with a $1$ at element $(i,i)$ if row $i$ was selected and $0$ otherwise. Furthermore, given the nature of $\mathbf{R}$, even though the product $\mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{ER}$ is rather large at $P \times P$, it is also approximately sparse with only a relatively small number of terms with significant modulus per row/column. This makes me think that the above form is doable since I could compute just the sparse approximation and avoid the unrealistic memory and computation requirements.

However, $\mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{ER}$ has size $P \times P$ but is only of rank $M'$ and this is causing problems with the  basis pursuit algorithms as they expect the rows of the constraint matrix to be
linearly independent. Hence, my solicitation for input.

Are there any factorizations and/or transformations that could turn the constraints of the second form of the problem into under-determined yet full-rank system? Please do keep in mind the large sizes of these matrices. Any other suggestions/comments are also welcome. We can assume that $\mathbf{R}$ has rank $M$ or $P$, as there is freedom in its design.

Since $ \epsilon $ is a parameter you need to set why not trade it with another parameter you need to set to create an easily solvable problem (Relaxation of the Problem)?

You can transform the problem into the following form ($ {L}_{1} $ Regularized Least Squares):

$$ \arg \min_{x} \frac{1}{2} \left\| A x - z \right\|^{2} + \lambda \left\| x \right\|_{1} $$

There are many solvers for that even for large sparse matrices.
For any $ \epsilon $ from the original problem there is $ \lambda $ such that the solutions are the same.
Since $ \epsilon $ it to be arbitrarily set in the original problem you can trade it for $ \lambda $.

This can be solved with many solvers easily without the need for Matrix Inversion.


Sum of sines sound wave using microcontroller I need to produce a signal which contains only several frequencies, for example, 17KHz, 17.2 kHz,17.6 kHz. The components I'm using are microcontroller (PIC16F1783 by microchip) which works at 80 kHz frequency, and produces an analog signal, which is then fed into an audio amplifier and then into a speaker.

First, I'm trying to output a single frequency, say 17 kHz. under timing and low memory constraints my algorithm is as follows:

I have an array of 400 elements which contains samples of a sine wave at constant intervals
Every clock cycle(80,000^-1 sec) I output sine_array[index] to the audio amplifier (and then to the speaker) where index starts at 0 and increased each cycle by: (17,000/80,000)*400=85.

However, for some reason, besides seeing the pure 17 kHz frequency, I also hear (and see on the oscilloscope) parasitic 5 kHz frequency!
I've also noticed that for frequencies such as 16 kHz (index increased by 80), 8 kHz there are no parasitic frequencies.
So I though it may have something to do with the fact that 400/80 is an integer, whereas 400/85 isn't.

Can anyone please suggest a creative way to sort this issue? The final goal is, as stated, to produce a signal which is composed of several frequencies with a rather small difference (around 200 Hz).

To explain the problem (without solving it):

  However, for some reason, besides seeing the pure 17KHz frequency, I also hear (and see on the oscilloscope) parasitic 5KHz frequency!

Your sampling rate is 80 kHz, so the Nyquist frequency is 40 kHz.  Your generation method produces a distorted sine wave at 17 kHz, which produces harmonics at:

17×2 = 34 kHz
17×3 = 51 kHz
17×4 = 68 kHz
17×5 = 85 kHz
...

But everything above 40 kHz gets aliased, so

51 kHz becomes 40-(51-40) = 29 kHz
68 kHz becomes 40-(68-40) = 12 kHz
85 kHz becomes 40-(85-40) = -5 kHz, which then aliases off 0 Hz and becomes 5 kHz
...

  I've also noticed that for frequencies such as 16KHz (index increased by 80) ,8KHz there are no parasitic frequencies.

Now you've got harmonics at 32, 48, 64, 80, ...  These alias to 32 kHz, 32 kHz, 16 kHz, 0 kHz, etc.  So they aren't as obvious, but they're still there unless your wavetable is lining up perfectly with the samples.

If you look at your output using an FFT analyzer, you should be able to see these aliased harmonics more clearly:

http://www.music.mcgill.ca/~gary/307/week4/wavetables.html


Basic question about chirp signal Suppose I generate chirp signal with the following python code:

 t = np.linspace(0, 10, 512)
 f0 = 2
 f1 = 5
 fs = 2*f1
 w = chirp(t, f0=f0, f1=f1, t1=10, method='linear')
 result  = rfft(w)
 phase = np.zeros(257)
 mag = np.zeros(257)
 for binnum in range(len(result)):
   mag[binnum] = abs(result[binnum])
   phase[binnum] = np.arctan2(result[binnum].imag, result[binnum].real)
 plt.plot(fs * np.linspace(0, 0.5, 257),np.unwrap(phase))
 plt.show()
 plt.plot(fs * np.linspace(0, 0.5, 257),mag)
 plt.show()

And the phase plot is:

Magnitude plot is:

But it is very different from what is described here.

Well, looking at the phase graph I can see parabola from 0 to 0.8. The magnitude is definitly  not a constant.

What I'm doing wrong or misunderstood?

Thanks in advance.

I dont understand the definition of the chirp in the link you provided (perhaps chirp in frequency domain?), but maybe the following explanation could help.

In general, a linear chirp with a varying frequency given by $$f(t) = \alpha + \beta t$$ is defined by $$ x_1(t) = e^{2\pi j (\alpha t + 0.5\beta t^2)} $$

In that case, the phase of $x_1(t)$, the signal in time-domain, is indeed parabolic and the magnitude is constant and equals to 1.
Now, I am not familiar with python, but the chirp function in matlab produces a real signal rather than complex one. That is $$ x_2(t) = \sin(2\pi  (\alpha t + 0.5\beta t^2)) $$ In that case, you will get different results. However, you can use the Hilbert transform (I guess this is in python) which, given $x_1(t)$ will produce $x_2(t)$ using fft.


Matching masked images Let's say I have two images that are parts of the same bigger image, cut using two masks (for simplicity, let's assume that the masks are just rectangles of the same size). The images are presumed to overlap. The concrete examples are: image stabilization (two different views of the same background) and object tracking (matching the object on two images).

I've read this excellent article http://werner.yellowcouch.org/Papers/subimg/, but it deals with finding the displacement of a small image which is a part of a large one, while I want to find displacement between the images of the same size. My images intersect, but don't contain each other.

I want to find the coefficients of the translation transform which correctly overlaps the matching parts of the images.
I want to use something like FFT-based correlation matching. I have the following problems with simply using FFT-correlation:

The edges of the image frames are strong features. If the algorithm searches the displacement $[dx, dy]$ that maximizes something like $\sum_{x,y}{a[x, y]}b[x + dx, y + dy]$, it would just snap the images so that their frames match, disregarding the content. I heard about windowing, but I don't understand how to apply it here. If each image is tapered at the edges, the images would no longer match when positioned correctly.
If we use a correlation function that normalizes the result by dividing it by the area of intersection resulting from the displaced matching, the algorithm would just try to keep the intersection to a minimum (~1x1 pixel; the smaller the intersection area is, the easier it is to find a nearly perfect match.)

How to solve this matching problem using windowing or some other techniques?

Update: I know that the problem can be partially eliminated by discarding the low-frequency information (taking gradient or applying the Sobel operator). I'm trying to solve the problem without discarding the low frequency information.
Also, leaving only high-frequency features still doesn't solve the outlined problem: Imagine an image with a bright thin horizontal beam and a bleak square. These features are high-frequency and thus will remain after the filter is applied. Cross-correlation will than match the bright beam parts ignoring the square.

Here is the example which shows the matching problem:

Full image with frame borders:

Frames 1 and 2:

Correlation-based incorrect frame matching:

Ok, here's the simple method illustrated in Python:

# Correlation is convolution with one input reversed
corr = fftconvolve(img1, img2[::-1, ::-1])

# Get coordinates of peak in x, y pixels
peak = unravel_index(corr.argmax(), corr.shape)
# peak = (199, 189)
# Calculate offset from center point
print array(peak) - 199 # Offset is [  0 -10]

So it finds the vertical offset (10 pixels), but not the horizontal, as you've shown (bottom left output below).  This is because the cross-correlation output is stronger in the center when more of the input images are overlapping (gets dark toward the edges even if the inputs are identical).  (The found peak is shown with the red cross.)

And here's the method I described in the comments.  I have no idea if there's a mathematically rigorous foundation for this or if it has a name, but it produces the right results.  It may or may not be the same thing as normalized cross-correlation (I haven't bothered to decipher what Wikipedia is saying):

weighting = fftconvolve(ones(shape(img1)), ones(shape(img2)))
corr = corr / weighting # Normalize
peak = unravel_index(corr.argmax(), corr.shape)[::-1]
# peak = (224, 189)
print array(peak)-199 # Offset of [ 25 -10]

So this method works.  It normalizes every element of the output array by the number of overlapping pixels that were used to calculate that element, so that when only the thin edges of the images are overlapped, it still produces a strong output.  Now the offset from the dim square is not swamped by the white bar (lower right of image).

Note that while this works perfectly with your example images, this method doesn't work as well for real images, since the edges of the output tend to vary wildly compared to the middle, and can produce spurious peaks, but if you know the offset is constrained to less than the image width, you can look for the peak only in realistic locations near the center, ignore the edges, and it works.

Wikipedia suggests a method using windowed images which I thought might improve on this, but unless I'm doing it wrong, it doesn't:

h = hamming(shape(img1)[0]) # 1D Hamming window
ham2d = outer(h, h) # 2D Hamming window

h1 = img1 * ham2d
h2 = img2 * ham2d

corr = fftconvolve(h1, h2[::-1,::-1])
peak = unravel_index(corr.argmax(), corr.shape)[::-1]
# peak = (216, 189)
print array(peak) - 199 # Offset is [ 17 -10]

The windowing does too much, and causes it to find an erroneous peak that's in between the previous two results.


Phase angle operator, Implemention and drawbacks Let assume we have a signal as:
$\color{blue} {x(t)=M e^{(-\beta.t)} e^{j(\sigma .t+\delta )}}$
By applying phase angle operator we have:
$\color{blue} {Phase[x(t)]=Phase[M e^{-\beta.t} e^{j(\sigma .t+\delta )}]=\sigma.t+\delta}$
I have two questions:
 1. How can I implement the Phase angle operator, i.e. Phase, in practice ,like a chip.
 2. Does phase angle operator increase the effect of noise in a noisy signal?(for example some operators like derivative operator increases the noise effect in a noisy signal)
Kind Regards

As mentioned in my comment, you need the atan2 function. If it is not implemented in the library you're using, you can implement it yourself quite easily according to this definition. If your signal is $x(t)=x_r(t)+jx_i(t)$, where $x_r(t)$ and $x_i(t)$ are the real and imaginary parts, respectively, then for a given time $t$, the phase is computed as $\phi(t)=\text{atan2}(x_i(t),x_r(t))$.

The computation of the phase of a signal is indeed sensitive to noise. Note that - no matter how you implement it - there is always some direct or approximate computation of the ratio of the imaginary and the real part of the signal, and since both are corrupted by noise, the ratio can be very noisy. The computation of the ratio can also become numerically problematic if the denominator (i.e. the real part) becomes very small.


Semantic Distance Measure Between Images I would like to compare 2 generic images, with same sizes and normalized values. Which metric would be better than the baseline Eucledian distance?

One of the most known image similarity metric is the SSIM. You can take a look to the next links:

https://en.wikipedia.org/wiki/Structural_similarity
https://www.imatest.com/docs/ssim/
https://www.cns.nyu.edu/~lcv/ssim/

EDIT 1:
Okay, you refer to semantic similarity. I have not worked on that issue really, and probably there are new DNN solutions to that purpose. I have not checked the state-of-the-art related to this topic.
Even so, I have made a quick test using the skimage SSIM and the difference between images to give you some idea about what happens with these metrics on several images. Find in the next images some hints about their performance for your target application:

Please, notice the values of the SSIM, MSE and STD DIFF on the images. SSIM and MSE works very well to compare quality of the same image. But not for semantic comparison.
Cheers.


How to ensure anti-aliasing in image resampling at lower resolution? Suppose I have a 1D image with intensities $f_i$, $i=[0,n-1]$, where $n$ is the number of pixels. I want to downsample it to have $\frac{n}{2}$ pixels (let $n$ be an even number). For each new point I want to average the intensity over 2 neighboring pixels as

$$\tilde f_i = 0.5(f_{2i}+f_{2i+1}), \qquad i=[0,n/2-1]$$

Can somebody, please, answer the following questions

Will such transformation introduce aliasing (and if so, what is the correct way to resample with averaging over pixels?)
Is there some math procedure, using which I can answer questions like #1 myself?
How, in general, discrete signals are anti-aliased? From DSP books I know that for continuous signals you just apply low-pass filter to your input signal and then sample whatever you got after the low-pass transform. But here is the problem: for images from the very beginning I have only discrete signal and I cannot figure out how to use continuous coordinate low-pass filters for discrete signals. May be I should treat an image as a piece-wise continuous function, which is constant within each interval $[i,i+1)$ ?

A "1D image" is just a discrete signal. All the classic signal theory applies:

  Will such transformation introduce aliasing (and if so, what is the correct way to resample with averaging over pixels?)

So, what happens when you downsample is decimation. The aliases are the "folding" of existing spectrum over the conserved part of the spectrum; imagine you have a signal whose spectrum (for simplicity, assume that's its discrete Fourier transform) is

$[A\,|\,B\,|\,C\,|\,D\,|\,E\,|\,F\,|\,G\,|\,H]$

and you simply throw away half of the samples, the resulting spectrum is

$\frac 12 [A+E| B+F| C+G| D+H]$.

  Is there some math procedure, using which I can answer questions like #1 myself?

So what you'd need to do is suppress $[E\,|\,F\,|\,G\,|\,H]$ before throwing away the samples. Something that does that is usually called a low pass filter (LPF).

LPFs are often implemented as Finite Impulse Response Filter (FIR), where you take the current input sample, multiply it with a factor, add a previous input sample (also multiplied with a, possibly different, factor).

As it turns out, your averaging is indeed such an LPF FIR. It has the taps (coefficients) of

$[0.5, 0.5]$. To see this, consider your $2i+1$th sample to be your "current" one, and the $2i$th to be "one in the past".

Now, such a FIR can be analyzed; it's z-domain frequency response is this:

Notice how at the vertical bar, i.e. the point below which you want to preserve and above you want to suppress the spectrum, this only has a value of 0.5; in fact, about 22% of the total pixel "energy" that passes through your filter is in the "upper" spectral half. That's too much, and you'll get aliasing.

You simply might want to build a better filter. This implies you'll add more taps, i.e. more samples right and left of your current sample

There's plenty of tools out there to design filters. Scipy's filter design tools will absolutely do; just design a low pass filter.

For example, I've used GNU Radio's gr_filter_design tool, and specified the passband to end just below the half-band boundary (at frequency $0.45f_{sample}$), and the stopband to start just above that boundary ($0.55f_{sample}$)¹. Suppression should be 30dB, i.e. the magnitude in the $\le\frac1{1000}$ of the passband in the stopband. The resulting filter:

It has 13 taps, i.e. one coefficient for the "center" pixel and one for each of the six pixels left and right of it:

-0.00405147997662425, 0.009050002321600914, -0.02354922890663147, 0.04652633145451546, -0.07230162620544434, 0.09264583140611649, 0.9033603072166443, 0.09264583140611649, -0.07230162620544434, 0.04652633145451546, -0.02354922890663147, 0.009050002321600914, -0.00405147997662425

Notice the symmetry (typical for the class of linear-phase FIRs).

Original | Filtered with your averaging filter, then downsampled | filtered with the 13-tap filter, then downsampled (please view in full size; it doesn't make sense to compare things scaled by your browser)².

Visual comparison shows that the readability of the text suffers a bit more from the long filter response, while the reproduction of the fine structures of the dark cacao powder atop the tasty, juicy FFTiramisu is is better with the 13tap filter.
This is most probably an effect of the filter design (it's a design processed based on the Hamming window). Pre-scaling low pass filters employed by picture processing software often aren't even linear; that makes modeling significantly more complicated, but since sensory input in humans typically has logarithmic charactistics, I presume that hides "border effects" better.

  How, in general, discrete signals are anti-aliased? From DSP books I know that for continuous signals you just apply low-pass filter to your input signal and then sample whatever you got after the low-pass transform. But here is the problem: for images from the very beginning I have only discrete signal and I cannot figure out how to use continuous coordinate low-pass filters for discrete signals. May be I should treat an image as a piece-wise continuous function, which is constant within each interval [i,i+1)[i,i+1) ?

Well, in DSP downsampling is a very common thing to do, so you do with digital filters, as explained above. Don't consider your signal to be continuous; it's not.

¹ FIRs can't have infinitely steep responses. As a general rule, the closer the passband edge is to the stopband edge, the more taps you will need to realize the filter.

² Processing was done using Python, numpy, scipy.signal.lfilter, scipy.misc.im{read,save}.


How to smoothen signal with missing values before differentiation? I want to differentiate a noisy signal with many, randomly located, missing data values. Which smoothing techniques should be applied before differentiating the signal?

I have a velocity signal at 200 Hz, but due to wireless transmission about half the data values are missing (NaN), up to 8 in a row (40 ms between data points). I would like to get the derivative signal (acceleration), but without smoothening, it is too noisy. I'm mostly interested in a rough (!) estimation of maximum acceleration. What is the best techniques to handle the missing data? Linear interpolation between known data points?

The best tool for this job is normalized convolution. It can deal with missing samples as well as uncertainty.

The paper describing the method is "Normalized and Differential Convolution -- Methods for Interpolation and Filtering of Incomplete and Uncertain Data" by Hans Knutsson Carl-Fredrik Westin. There is a PDF on Semantic Scholar.

Normalized convolution basically takes a convolution with a Gaussian of the data, with 0 for the missing samples, and divides sample-wise with the result of the same convolution on a weight signal, which is 1 for known samples and 0 for the samples you set to zero. More formally:

$$
g = (( f m ) * h ) / (m * h)
$$

(With $h$ the convolution kernel and $m$ the weight signal). You can derive the analytic derivative of that equation, leading to a derivative operator that is robust to missing samples. See the paper for more details.


Synchronize device to audio stream - Approach I am working on a device which's actions should be synchronized to the audio played. The audio stream is known in beforehand and can be altered to contain signals for the synchronization. But the device and the PC playing the music do not have any other connection but the audio.
The synchronization should be in the regime 100ms.

The transfer function between the device and the PC is not known and can vary. Especially the gain can vary. Thererfore an approach not to sensitive to overmodulation would be preferrabel, if such an approach exists.

The current signal in the audio-stream is a click noise. A crosscorrelation is working quite well on the original data, but the signal gets less clear when a transfer function is introduced. And vanishes for an overmodulated signal, while for the human ear the signal is still present.

What would be a suitable approach to synchronize the device to the audio stream?

EDIT: It is an acoustical link.

One problem with an acoustic link is frequency dependent multi-path fading/reinforcement, which can vary widely as the source and mic positions change, or as the room acoustic vary (moving hard objects, etc.)

So I might try a "click" that has a wide and known pattern of several non-harmonically related spectra, and check the Hamming distance between the transmitted and received spectral pattern peaks.

Use a non-harmonically related set of frequencies so that the harmonics of any clipped spectra will be mostly ignored (depending on window duration), and because multi-path fading can act as a comb filter notch.


Beginner's guide to OFDM As the title suggests, I'm a bit of a dummy, and utterly perplexed by the concept of OFDM, overlapping the carriers at the clock nulls and somehow magically not creating an unusable amount of cross carrier interference. What's the science (or sorcery, as is so often the case with wireless communications) behind this concept?

In time domain, OFDM breaks one serial fast bit stream into many parallel slow bit streams.
Then, these parallel slow bit streams are multiplied with orthogonal sinusoids, where orthogonality between two sinusoids is defined with a summation over a certain time interval as

$$\sum \limits _{n=0} ^{N-1} \cos 2\pi k f_0 n \cdot \cos 2\pi k' f_0 n = 0$$ when $k \neq k'$.

This process is illustrated in detail in Figure below with the help of an example.

Now for each such individual signal, the multi-path components arrive in not
too distant future (just like a low rate stream). Hence, the
equalizer design is easy having less spread paths and consequently less interference with future symbols, provided that we find a way to separate the subcarriers at the Rx.

Separating these bits at the Rx is easy: we can correlate (multiply sample by sample and sum) the composite signal with just one subcarrier.  Utilizing their orthogonality property, contribution from all other subcarriers will cancel out to zero, while the contribution from the subcarrier with that frequency will pop out, scaled in amplitude by our modulation signal. The same procedure can be repeated for all other subcarriers.

Essentially, this is an operation of Discrete Fourier Transform (DFT) as

$$X[k] =\sum \limits _{n=0} ^{N-1} s[n] e^{-j2\pi \frac{k}{N}n}$$

which, for each $k$, will generate our modulate data.

With respect to frequency domain, the signal at the Rx is a product of the spectra of the Tx signal and the wireless channel. Thus, the channel will allow some frequencies to pass through unharmed while suppressing some others. The low data rate signal needs less manipulation by the Rx to get the original data back due to a narrower spectrum. Essentially for this kind of signal, the channel acts just as a single multiplier that can be equalized through estimating that number and dividing the Rx signal by that number. So the equalization reduces to a single division operation. On the other hand, a high data rate signal needs a lot of Rx processing to equalize the channel.

To solve this problem, what OFDM does in frequency domain is fairly simple. It just segments the available bandwidth into many parallel almost flat channels through utilizing those sinusoidal subcarriers. This slicing of the spectrum is drawn in Figure below, just like slicing a bread.

Hence, equalization for each narrow slice requires just a single division operation, rendering the computational load of the equalizer to a total of $N$ divisions.


Equivalent of Fourier transform for Gaussians? I have a long signal (million of samples) containing a lot of Gaussian peaks, whose standard deviation is random and about $5$ to $50$ samples wide. Sometimes, these peaks overlap, but not often. The signal also contains high frequency noise.

I would like to compute the distribution of standard deviations, in a similar way one can extract the "distribution of frequencies" using FFT. Currently, we do a first pass for detecting the peaks, then do a classical function fit in the peak region. This is not very robust and the distribution fluctuates a lot depending on the fitting algorithm fine tuning.

Is there a more robust approach which would do a kind of a Gaussian transform?

The (continuous) Fourier transform of a Gaussian is a Gaussian:

$$ \mathscr{F}\Big\{ x(t) \Big\} \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, dt $$

$$ \mathscr{F}\left\{ e^{-\pi t^2} \right\} = e^{-\pi f^2} $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha t^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} \qquad \alpha > 0 $$

$$ \mathscr{F}\left\{ e^{-\pi \alpha (t-\tau)^2} \right\} = \frac{1}{\sqrt{\alpha}} e^{-\frac{\pi}{\alpha} f^2} e^{-j 2 \pi f \tau} $$

$$ \mathscr{F}\left\{ \sum_m c_m e^{-\pi \alpha_m (t-\tau_m)^2} \right\} = \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\frac{\pi}{\alpha_m} f^2} e^{-j 2 \pi f \tau_m} $$

note these gaussians have "tails" that go on forever, but gaussians usually die off to very close to zero very rapidly.

assuming no time or frequency aliasing, let $x[n] \triangleq x\big(n/f_\text{s}\big)$ where $f_\text{s}$ is the sample rate.  Then the DTFT is

$$\begin{align}
X\big(e^{j\omega}\big) &\triangleq \sum_{n=-\infty}^{\infty} x[n] \, e^{-j \omega n} \\
&= f_\text{s} \sum_{n=-\infty}^{\infty} x\big(n/f_\text{s}\big) \, e^{-j f_\text{s} \omega (n/f_\text{s})} \frac{1}{f_\text{s}} \\
& \approx f_\text{s} \int\limits_{-\infty}^{+\infty} x(t) \, e^{-j 2 \pi \frac{f_\text{s} \omega}{2 \pi} t} \, dt \\
&= f_\text{s} \ \mathscr{F}\Big\{ x(t) \Big\} \Bigg|_{f=f_\text{s}\frac{\omega}{2 \pi}} \\
\end{align}$$

so if we let $$x(t) = \sum_m c_m e^{-\pi \alpha_m (t-\tau_m)^2}$$

$$\begin{align}
x[n] &= \sum_m c_m e^{-\pi \alpha_m ((n/f_\text{s})-\tau_m)^2} \\
&= \sum_m c_m e^{-\pi (\alpha_m/f_\text{s}^2) (n-f_\text{s}\tau_m)^2} \\
&= \sum_m c_m e^{-\pi \widehat{\alpha_m} (n-\widehat{\tau_m})^2} \\
\end{align}$$

then the DTFT is

$$\begin{align}
X\big(e^{j\omega}\big) & \approx f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\frac{\pi}{\alpha_m} f^2} e^{-j2\pi f\tau_m} \Bigg|_{f=f_\text{s}\frac{\omega}{2 \pi}} \\
&= f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\frac{\pi}{\alpha_m} \left(f_\text{s}\frac{\omega}{2 \pi}\right)^2} e^{-j2\pi \left( f_\text{s}\frac{\omega}{2 \pi} \right)\tau_m} \\
&= f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\pi\frac{f_\text{s}^2}{4 \pi^2 \alpha_m} \omega^2} e^{-j \omega (f_\text{s}\tau_m)} \\
\end{align}$$

Then, assuming you choose your DFT length $N$ to be sufficiently large (you're saying in the millions) to cover all of the gaussians with their specific delays $f_\text{s}\tau_m$ in samples and widths of ca. $\frac{f_\text{s}}{\sqrt{\alpha_m}}$.  Then the $N$-point DFT will sample this DTFT at $N$ equally spaced points around the unit circle.

$$\begin{align}
\mathcal{DFT:} \quad X[k] &= \sum_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
&= X\big(e^{j\omega}\big) \Bigg|_{\omega = 2\pi\frac{k}{N}} \quad \mathcal{:DTFT}\\
&= f_\text{s} \ \sum_m \frac{c_m}{\sqrt{\alpha_m}} e^{-\pi\frac{f_\text{s}^2}{4 \pi^2 \alpha_m} \omega^2} e^{-j \omega (f_\text{s}\tau_m)} \Bigg|_{\omega = 2\pi\frac{k}{N}}\\
&=  \ \sum_m \frac{c_m}{\sqrt{\alpha_m/f_\text{s}^2}} e^{-\pi\frac{f_\text{s}^2/\alpha_m}{4 \pi^2} \left(2\pi\frac{k}{N}\right)^2} e^{-j 2\pi\frac{k}{N} (f_\text{s}\tau_m)} \\
&= \sum_m \frac{c_m}{\sqrt{\widehat{\alpha_m}}} e^{-\frac{\pi}{\widehat{\alpha_m}} \left(\frac{k}{N}\right)^2} e^{-j 2\pi\frac{k}{N} \widehat{\tau_m}} \\
\end{align} $$

So summing up and I will remove the little hats in the discrete-time domain, if you have a collection of gaussian pulses with heights of $c_m$, peak width parameters $\frac{1}{\sqrt{\alpha_m}}$ (in units of samples), and peak positions at $\tau_m$ (also in units of samples), then your input should look like:

$$ x[n] = \sum_m c_m e^{-\pi \alpha_m (n-\tau_m)^2} $$

and the result of the $N$-point DFT should look like:

$$ X[k] \approx \sum_m \frac{c_m}{\sqrt{\alpha_m}} \ e^{-\pi/(N^2\alpha_m) k^2} \ e^{-j 2\pi(\tau_m/N) k} $$

for $-\tfrac{N}{2} < k < \tfrac{N}{2}$ .  Remember with the DFT, $X[k+N]=X[k]$ for all integer $k$.


Various size arrays interpolation - resampling sorry if it's stupid question but it’s third day that I try to figure it out how to do that. And just can’t find the solution.

Let’s say I have array $A$ size 4, and $B$ size 8. And I need to put/spread evenly all 8 values from $B$ to 4 places in $A$.

In that case it’s quite easy to imagine the solution, but let’s say I have $A$ size 44100, and $B$ size 48000. How to express 48000 values of $B$ by 44100 values of $A$.

I know there are a lot of free and ready to use resampling algorithms, by I don’t need to use them. I need to understand mathematical basis of those algorithms. And not only for audio. But for other purposes also. I need it for my study. Could anyone help me?
Best regards.

Re-sampling is annoyingly messy & complicated.

The math is straight forward enough. Let's start with "rational" sampling rate conversion. Rational means the ratio of the two sample rates can be expressed as a fraction of two (preferably small) integers. Example: going from 48 kHz to 40 kHz has a ratio of 5/6 (output/rate over input rate) and the frequency that both divide is 240 kHz

The steps are the following

Up-sample by the numerator to the common divisor frequency. This is done by inserting zeros between the samples of the input signal.
Up-sampling results in periodic repetition of the spectrum in the frequency domain. It creates so-called mirror spectra. In our example we have up sampled from 48 kHz to 240 kHz.
In order to get rid of the mirror spectra you need to lowpass filter the signal.  That lowpass filter needs to be chosen so there is no aliasing, i.e. below the smaller of the two Nyquist frequencies.
Downsample the low passed signal. This is simply done by throwing away the samples you don't need. In our example, you would only keep every sixth sample.

So the basic operations are: up-sampling, low-pass filtering, down-sampling.

This process tends to be expensive: you are running a very steep low-pass filter at a very high sampling rate. There are a few ways to make it a lot more efficient, especially if the lowpass filter is an FIR filter. There are two properties you can exploit:

You don't need to calculate the samples that you are going to throw away anyway
After up-sampling, most of your input samples into the low-pass filter are zero, so you can eliminate these from the convolution calculation

When you do both of these things, you end up with a so-called poly-phase filter. See for example http://www.ws.binghamton.edu/fowler/fowler%20personal%20page/EE521_files/IV-05%20Polyphase%20FIlters%20Revised.pdf

In essence you calculate each output sample by filtering the appropriate input samples with an FIR filter. The FIR filter is time-varying, i.e. it's different from one sample to the next. You can also interpret the whole process as filtering with a time-variant fractional delay: the filter you need is a function of how much time offset there is between a specific output sample and the input sequence. In our example you would need to sample at the times $0,1.2,2.4,3.6,4.8,6.0,7.2 $ etc

The number of different filters you need is the numerator in the ratio fraction. That's okay for 48 kHz to 40 kHz but less convenient for 48 kHz to 44.1 kHz where the integer ratio can only be reduced to 137/160. In many cases the ratio cannot be reasonably expressed as an integer ratio and often it's also time variant (different master clocks). In this case we need to do irrational sample rate conversion.

Irrational sample rate conversion typically also uses a poly phase filter with some reasonable number of phases (say 32 or thereabouts). For each output sample you then calculate the exact phase (or fractional delay) and then linearly interpolate between the two closest phases that you have.

Now comes the tricky part: The design and choice of the low-pass filter needs to be carefully matched to your specific application requirements. We can't use an ideal low-pass filter since it's infinitely long and non-causal. So when designing this filter you need to know things like

What signal to noise ratio are you shooting for?
What overall delay can you tolerate ?
What is the highest "good" frequency you want?
Are you sensitive to phase distortion, phase modulation ?
How much aliasing can you tolerate?
How much pre-ringing is acceptable?
How much CPU and memory can spent on this?

Then you need to feed all of this into the filter design process, which is a non-trivial problem.


Combining simultaneous measurements of the same signal for better SNR? I have multiple data sets of the same signal measured using the same hardware and I want to know what sort of techniques I can use to improve the SNR.

I was thinking if I adjust and ensure they have the same phase, i could coherently add them, but I didn't really see an improvement in SNR. I'm guessing the reason for this is because the noise gets added together as well.

It makes me wonder two things:

What are some valid ways to combine like signals, to improve SNR
Is there any way to exploit the similarity in signals to reduce noise?

Thank you

Coherent Gain happens when you have copies of the same signal but different noise. If the noise is independent in each copy, the gain should be $10log N$ dB where $N$ is the number of copies. For a well crafted measurement system, additive electronic noise should be independent. Things like ambient noise will have some correlation provided a reasonable distance between sensors and $5logN$ is a reasonable lower expected gain.

If you don’t see at least $5 log N$ dB, then there is probably something wrong. Coherent Gain works. If it didn’t, many here would not have jobs.

You might try cross correlating pairs of your copies as a diagnostic. Without knowing your data, I cant say with certainty that you should see a well defined peak but the absence of a well defined peak would be a cause for concern