Applying Nyquist's sampling theorem to a real signal

I do not see specifically what you think is incorrect. I follow your thinking and do not see a flaw but let's review the key characteristics below in case that fills in a hole in the thought process:

For periodic systems we can consider a time span that extends over one period (0 to T) or from $-\infty$ to $+\infty$. Both are equivalent as long as we agree that the $\infty$ case must be periodic.

Similarly for sampled systems we can consider a frequency span that extends over one frequency range, 0 to $F_s$, where $F_s$ represents the sampling frequency, or from $-\infty$ to $+\infty$. Both are equivalent as long as we agree that the $\infty$ case must be similarly periodic in frequency.  (Note the base frequency that repeats can also be given as $-F_s/2$ to $+F_s/2$ but the range 0 to $F_s$ is most common when referring to the DFT.)

So the DFT itself is both periodic in time and in frequency, representing a sampled system that is periodic in time.

For any other systems; aperiodic or not sampled, we cannot use the DFT and call that equivalent, we can only approximate those systems with the DFT.

Below is a summary of the various Fourier Transforms that we can use in those other cases:

CASE 1: Continuous and Aperiodic in time: The Fourier Transform

First case we have is the Fourier Transform, where the time function is by definition continuous and extending from $-\infty$ to $+\infty$. Since it is continuous in time (not sampled) the frequency result will be aperiodic (not repeating). Since it is aperiodic in time, the frequency result will be continuous.

CASE 2: Continuous and periodic in time: The Fourier Series Expansion

Note that x(t) can still be periodic in the case above; there is nothing that excludes taking the Fourier Transform as written in the graphic above on a periodic function. This will result in the continuous Fourier Transform result as being 0 except for the frequencies at multiples of 1/T, including zero, which may have non-zero values at only these locations.  Knowing this, we can solve for just those frequencies in this case, resulting in the Fourier Series Expansion:

CASE 3: Discrete and aperiodic in time: The Discrete Time Fourier Transform

The third case is to now consider the sampled version of our first case, where we have an aperiodic time domain function again, but this time it is sampled in time. This is the Discrete Time Fourier Transform. We have a similar result to the first case, the Fourier Transform, with the only change that it is now repeating in frequency. As previously described the repetition boundaries occur at multiples of $F_s =  1/T$, so if our frequency spectrum in our first case extended beyond $F_s$, we would get frequency folding (aliasing) as a result and therefore would not be able to represent our time domain waveform exactly (the sampling rate was not high enough!).

(Side note: Zero padding prior to taking the DFT approximates the DTFT, the more zeros you add in, the more the frequency result approximates a continuous function, and that function is samples of the DTFT!).

CASE 4: Discrete and periodic in time: The Discrete Fourier Transform

Finally, we come to the DFT (the FFT is an algorithm to compute the DFT exactly). In this case, when considering the same time axis as above that extends from $-\infty$ to $+\infty$, the time domain waveform must be repeating (periodic).  The periodicity in time results in discrete frequencies that can only exist at integer multiples of 1/T Hz where T represents one period in time in seconds. Note that like the Fourier Series Expansion, the periodicity is implied; the actual transform considers N samples in time and provides N samples in frequency, but the periodicity described is evident when both axis are extended to infinity.

I hope this answered your question. If not please comment on the remaining points of confusion and I will try and update accordingly.


Derive the Forward Euler substitution for transfer function

To the extent you can factor the transfer function into individual integrator sections of the general form $\frac{1}{s}$ you can make this substitution, which is an approximation of the Matched-$z$ Transform where you substitute every $s$ for $s=\frac{\ln(z)}{T}$. (map from $s$ to $z$ using $z =e^{sT}$).

This results in first order forms given by $H_\mathrm{a}(s) = k_\mathrm{a}\frac{s-b}{s-a}$ being replaced with $H_\mathrm{d}(z) = k_\mathrm{d}\frac{z-e^{bT}}{z-e^{aT}}$ where $k_\mathrm{a}$ and $k_\mathrm{d}$ are fixed gain constants. $k_\mathrm{d}$ can be determined by using $H_\mathrm{a}(s=0) = H_\mathrm{d}(z=1)$ for low-pass forms or $H_\mathrm{a}(s=j \pi/T ) = H_\mathrm{d}(z=\pi)$ for high-pass forms.

For all pole forms this is also equivalent to what you would get using the mapping from $s$ to $z$ using the method of Impulse Invariance. The Matched-$z$ Transform is not recommended for higher order polynomials with zeros and poles since neither the time or frequency domain response is preserved (The method of Impulse Invarance or Bilinear Transform would more often be used, each with their own advantages).

This post summarizes various mapping techniques further: Are there alternatives to the bilinear transform?


What is the algorithm to generate sine waves of arbitrary frequency in the STFT domain?

The FFT of a non-integer multiple frequency is very complex, you don't want to go there!
There are various oscillator schemes that are a few MAC instructions per sample such as https://vicanek.de/articles/QuadOsc.pdf
Checkout a review here: https://www.njohnson.co.uk/pdf/drdes/Chap7.pdf


It's part of a feedback loop -- presumably setting the sample rate in some real-time audio thingie, since I see keywords in the code like 'resample' and 'jack' and 'PLAY'.
// Run loop filter and set resample ratio.
So, this tells you what it does!
Elsewhere in the code, err is computed as the timing error between the input stream and the output stream.  So, this is a feedback loop.
_z1 += _w0 * (_w1 * err - _z1);
_z2 += _w0 * (_z1 - _z2);

Low-pass filter the error with a pair of first-order low-pass filters.  (_z1 is err filtered, _z2 is _z1 filtered).  Note that for this to work, _z1 and _z2 need to be persistent.  In C/C++, this means that they should be at least at file scope -- and in my humble opinion, this means they should have longer, descriptive names, even to the point of declareing a struct called "filter_state" -- it never hurts to hand your reader a clue-by-four.
I'm not going to dredge through all the loop analysis, but the cutoff frequencies of these filters should be significantly higher than the target loop bandwidth -- and the loop bandwidth is determined, in part, by the gain of the timing error detection process.  So you can't just slap this into your code willy-nilly.
_z3 += _w2 * _z2;
This is an integrator.
_rcorr = 1 - _z2 - _z3;

if (_rcorr > 1.05) _rcorr = 1.05;
if (_rcorr < 0.95) _rcorr = 0.95;

This does the final calculation of the command variable to the resampler.  Note, however, that we have a controller with an integrator, and we're limiting the output of the integrator, but we're not limiting the state of the integrator.  If you implement this code, search on "integrator anti-windup" for why this is a mistake, and why _z3 should be limited in magnitude, somewhere between $\pm 0$ and $\pm 0.1$, depending on your desired recovery characteristics from a large excursion.
_resamp.set_rratio (_rcorr);
This is actually setting the resample ratio, which in turn affects the value of err in the next iteration of the loop.


Demodulating upper sideband (USB) signals

Apart from scale factors, your USB signal is

$$f(t)=m(t)\cos(2\pi\nu_ct)-\hat{m}(t)\sin(2\pi\nu_ct)\tag{1}$$

(as you've correctly stated in your question). Now if you multiply with a (coherent) carrier $2\cos(2\pi\nu_ct)$, you get

$$\begin{align}g(t)&=2\big[m(t)\cos(2\pi\nu_ct)-\hat{m}(t)\sin(2\pi\nu_ct)\big]\cos(2\pi\nu_ct)\\
&=2m(t)\cos^2(2\pi\nu_ct)-2\hat{m}(t)\sin(2\pi\nu_ct)\cos(2\pi\nu_ct)\\&=m(t) (1+\cos(4\pi\nu_ct))-\hat{m}(t)\sin(4\pi\nu_ct)\tag{2}\end{align}$$

Applying a low pass filter to $(2)$ will filter out all terms at $2\nu_c$, leaving you with the message signal $m(t)$.


Definitely yes. But after IFFT operation data block is CP extended before transmission.
Each sample of OFDM symbol from the output of IFFT is weighted sum of all QAM/PSK modulated subcarriers from the input of IFFT.
In the MIMO case all spatial-time algorithms are applied in the frequency domain. So in transmitter you should do encoding with QAM subcarriers before IFFT and do the distinct IFFT for each spatial channel created. While in receiver you should do decoding after FFT applied separately to each spatial channel again. By doing this you exploit the fact OFDM is multiplexing of the large number (actually M) of narrow band channels and by appropriate processing (actually FFT/IFFT pair and CP insertion technique) we can treat each subcarrier as narrow band. It is very useful in the case of MIMO because MIMO decoding algorithms are much too simpler when dealing with narrow band (so flat fading channel is the case) and for each subcarrier we can define simple MIMO matrix of weights $[W_{ij}]$. That's why OFDM + MIMO is so popular combination today.


Let's assume that $A(\omega)$ is complex valued:

$$A(\omega)=B(\omega)e^{j\phi_A(\omega)}\tag{1}$$

where $B(\omega)$ is a (possibly bipolar) real-valued function, and $\phi_A(\omega)$ is a (real-valued) phase function.

The phase of

$$H(e^{j\omega})=A(\omega)e^{-j(\alpha\omega-\beta)}=B(\omega)e^{j\phi_A(\omega)}e^{-j(\alpha\omega-\beta)}=B(\omega)e^{j\phi(\omega)}\tag{2}$$

is then

$$\phi(\omega)=\phi_A(\omega)-(\alpha\omega-\beta)\tag{3}$$

(where I define the phase such that jumps by $\pi$ are taken care of by the sign of $B(\omega)$, not by $\phi(\omega)$). For the phase $\phi(\omega)$ to be linear, $\phi_A(\omega)$ must be constant or linear:

$$\phi_A(\omega)=a\omega+b\tag{4}$$

With ($4$) the total phase becomes

$$\phi(\omega)=(a-\alpha)\omega+(b-\beta)=-\alpha'\omega+\beta'\tag{5}$$

Combining $(2)$ and $(5)$ we get

$$H(e^{j\omega})=B(\omega)e^{j\phi(\omega)}=B(\omega)e^{-j(\alpha'\omega-\beta')}\tag{6}$$

which has the same form as the equation in your question. Note that in $(1)$ we have defined $B(\omega)$ to be real-valued.

This shows that $(6)$ (or, equivalently, the equation in your question) is the most general form of the frequency response of a linear phase filter. If $A(\omega)$ were chosen to be complex-valued, then $A(\omega)$ must have a linear phase (otherwise $H(e^{j\omega})$ wouldn't have a linear phase), and that linear phase can be absorbed by the constants $\alpha$ and $\beta$, such that the expression for the frequency response can always be reduced to the form $(6)$ with a real-valued function $B(\omega)$.


Interpretation of the eigenvalues of the inverse Hessian in a KLT tracker

Think of them as 2D smoothness terms.
The smoother the patch, the lower the matrix rank and the closer the matrix is to being singular.

On an straight edge (not a corner), just one eigenvalue will be large.
On a corner both will be large.

Using eigenvalues means that the angle of the edge is not a factor, and at any angle, an edge will give just one large e.v.


Why does the general formula for DFT of real sinusoids not apply here?

It's not a contradiction. It's just that for an alternating signal the two contributions at $\pm k_0$ coincide and add up. Note that in your example with $N=4$, you have $k_0=\pm 2$, but these two values of $k_0$ are in fact identical due to the periodicity of $N=4$: $X[-2]=X[2]$.

The same would be true for $k_0=0$, i.e. a constant signal (with frequency zero and amplitude $A_0$). In that case you'd get $X[0]=A_0N$, just like in the case of a signal at Nyquist (maximum frequency, as in your example) where you get $X[N/2]=A_0N$. For all other values of $k_0$ you get two contributions with amplitudes $A_0N/2$.

Regarding your second question, this is an identity that is very important in DSP. First of all note that due to periodicity of the complex exponential function you have the equality

$$\sum_{n=-N/2}^{N/2-1}e^{-j2\pi (k-k_0)n/N}=\sum_{n=0}^{N-1}e^{-j2\pi (k-k_0)n/N}\tag{1}$$

For $k=k_0$ the complex exponential in $(1)$ becomes $1$, so the sum equals $N\cdot 1=N$. For $k\neq k_0$ you can use the formula for the geometric series to see that

$$\sum_{n=0}^{N-1}e^{-j2\pi (k-k_0)n/N}=\frac{1-e^{-j2\pi (k-k_0)}}{1-e^{-j2\pi (k-k_0)/N}}=0\tag{2}$$

because for $k\neq k_0$ we have $e^{-j2\pi (k-k_0)}=1$.


FIR difference equation vs convolution summation

I should have read further before asking this question! Anyway, I've found out that there are different ways to implement the FIR system as mathematically described in Equation 7.2.1 in the post above. Among these realizations/structures is the direct-form, which is the simplest to begin with. Equation 7.2.4 shows the transformation of the general equation in 7.2.1 in accordance with the direct-form.


Yes image scaling can also be performed in frequency domain using DFT/FFT.

For example, given an $N \times M$ image, you can implement the interpolation by integer $K$ in the frequency domain by enlarging the $N \times M$ point DFT/FFT of the original image into $K \cdot N \times K \cdot M$ size where new DFT/FFT samples are all zero; except at the four corners; effectiveley you are zero padding the original DFT/FFT in frequency domain. The effect in time domain wil be interpolation.

The following is an example of 2x image zooming using FFT.

K = 2;
I = im2double(imread('Cameraman.tif');
S = size(I);

F = fft2(I);
F2 = repmat( K*K*F , K , K);

Mask = zeros(K*S);
Mask(1:S(1)/2+1 , 1:S(2)/2+1) = 1;
Mask(end-S(1)/2 : end , 1:S(2)/2+1) = 1;
Mask(1:S(1)/2+1 , end-S(2)/2 :end ) = 1;
Mask(end-S(1)/2 : end , end-S(2)/2 :end ) = 1;

F2 = Mask.*F2;

I2 = real( ifft2( F2 ));
figure,imshow(I);title('Original');
figure,ishow(I2);title('2x Zoomed');

You should adjust the mask carefully for even and odd length dimensions...


In framing of audio samples', what is need of frame shift while giving frame size??

By performing the windowing with overlap we are artificially increasing our time resolution (larger granularity of features in time). This is especially useful when frame duration is long (bad time resolution, very good frequency resolution), thus yielding kind of extra 'time resolution'.
Usually no one is using the rectangular window, but other types such as Hamming, Hann, etc. As you probably know, windows tend to zero samples at both of their ends. By doing data you are kind of loosing some data. Let's imagine that you have some very short event you want to detect, at point where two windows are joined. If there is no overlap, then these samples will be almost removed. But if you perform 50% of overlap, there will be another window, right in the centre, and you won't lose anything.


Flat fading criterion of OFDM subcarrier spacing

Your mathematical derivation is correct, your $H[k]$ is the single-tap equalizer (i.e. one tap for each subcarrier, and the subcarriers do not mix with each other. That's the orthogonal in OFDM).

Let me try to explain this a bit more general, without going into coherence bandwidth and flat fading. To my understanding, explaining it with $B_c$ and flat fading is a bit superficial and undergraduate explanation, such that you understand the idea, but can't really prove it works. Even, the coherence bandwidth is not strictly defined, how can one design a system according to that?

Assuming that the channel is time-invariant over all times, the channel is an LTI system. The eigenfunctions of LTI systems are complex exponentials (of all frequencies). Note that these complex exponential range over all times (i.e. they are infinitely long). So, in principle, one could transmit data on all frequencies (infinitely small apart), but one would have to wait infinitely long until the signals arrived. The principle bases on the convolution theorem, i.e. your channel in time-domain does a convolution, so in frequency domain, it's doing elementwise multiplication. That's the continuous-time principle of OFDM. Clearly, not very attractive for practical implementations.

What can you do? We've know that two complex exponentials of different frequency are orthogonal to each other when considering the whole time axis, i.e. $\int_\mathbb{R}\exp(j2\pi f_1 t)\exp(j2\pi f_2 t)dt=\delta(f_1-f_2)$. We want to have orthogonality between the signals, (because we want to have single-tap equalization), but we dont want to wait infinitely long. Unfortunately, the orthogonality between to arbitrary frequency complex exponentials is not given, when the time interval is shorter: $\int_T\exp(j2\pi f_1 t)\exp(j2\pi f_2 t)dt\neq\delta(f_1-f_2)$. However, there are frequencies $f_1, f_2$ such that the orthogonality holds, namely $\int_T\exp(j2\pi f_1 t)\exp(j2\pi (f_1+\frac{n}{T})t)dt=\delta_n$ (here $\delta_n$ is the Kronecker symbol, i.e. 1 for $n=0$, else zero.). I.e. two exponentials that are $n/T$ apart in frequency are orthogonal over a time interval of $T$.

So, in principle, we could use time-intervals of $T$ and use only frequencies that are $1/T$ apart from each other to transmit orthogonal signals. But, there's more. If the orthogonality should hold at the receiver, the received signal needs to be a linear combination of complex exponentials. We know, that complex exponentials are eigenfunctions to the LTI system, so everything should be fine, right? No! Only infinitely long complex exponentials are eigenfunctions. If we transmit $\exp(j2\pi f_1t)\text{rect}(t/T)$, we will not get a complex exponential at the output, so it wont work. So, what to do? Essentially, we'd have to transmit infinitely long complex exponentials with frequency distance $1/T$ to get signals at the receiver, which are orthogonal in the interval $T$. Hm, still doesn't sound very appealing.

Here's the solution: The impulse response of the LTI system (i.e. the channel) is usually not infinitely long, but has a length $T_C$ (at least, we approximate/assume that it's zero for $t>T_C$). If the channel has length $T_C$, then for a given point in time $t$, only the signal times from $t-T_C$ to $t$ have an influence onto the output signal at the current point in time. So, we dont need to transmit infinitely long complex exponentials, but just of length $T+T_C$, and we get orthogonality in the time interval $[T_C,T+T_C]$.

You see where this is leading to, right? The extra time in the beginning is what we denote as the Cyclic prefix. Why do we call it cyclic? Because, the exponentials with frequency distance $1/T$ are periodic with period $T$. So, the signal in time $[0,...,T_C]$ is exactly equal to the signal at $[T,...,T+T_C]$.

To summarize: Transmitting complex exponentials of duration $T+T_C$ over an LTI system with impulse duration $T_C$ yields orthogonal signals over the interval $[T_C,...,T_C+T]$. Orthogonal means single-tap equalization is possible. What is the coefficient for the equalizer? Due to the convolution theorem, it is just the value of the frequency response of the channel at the carrier frequency.

This is the time-continuous explanation, which can become intuitive. However, I personally like more the discrete-time explanation, using linear algebra.

Let $F$ be the N-point Fourier transform matrix, i.e. $FF^H=F^HF=I$. So, the transmitted signal without CP is $y=F^Hx$ where $x$ is the data. Adding a CP is done by $y_{cp}=Cy$ with $C$ being a $N+N_{CP}\times N$ matrix. Then, we the received signal becomes

$$z=HCF^Hx$$ where $H$ is the Toeplitz matrix that does the convolution (I ignore noise here). $z$ is of length $N+N_{CP}+L-1$, where $L$ is the channel length. Then, at the receiver, the CP is removed, by a matrix $D$ of dimension $N\times N+N_{CP}+L-1$, i.e.

$$w=DHCF^Hx$$

Finally, we go to the frequency domain:

$$W=FDHCF^Hx$$

Let's have a look at the $N\times N$ matrix $DHC$: It's a matrix that performs circular convolution with the impulse response of your channel. In finite-discrete time the Convolution theorem states, that circular convolution equals elementwise multiplication of the DFTs. So, $FDHCF^H$ is a diagonal matrix, containing the frequency response of the channel on the diagonal. There you have your single-tap equalizer again (inversion by diagonal matrix is just division by the diagonal).


Using the traditional instruments for measuring high frequency signal

Not sure to understand what you are asking, Microwaves wavelengths range between 300 MHz (1 m) and 300 GHz (1 mm) which is indeed of the same order of magnitude as the probe.

1GHz corresponds to a wavelengths of about 30 cm  and is also in the Microwave range.

Are you saying that it is possible to measure a frequency of 1GHz with a probe while it is usually said that Microwaves can't be measured with it?


How the increasing of sampling frequency in OFDM didn't cause increase on the required channel bandwidth

You can say the null subcarriers have their own bandwidth if you define an alphabet including "zeros" and use null carriers to transmit these "zeros". As the comment of MBaz,

  If we start with the discrete-time OFDM symbol s[n], and increase the
  sampling rate (defined by the time interval between successive samples
  of s[n]), then the OFDM symbol becomes "shorter" and therefore its
  bandwidth will increase.

But it is not the case. The OFDM symbol $T_u = N_{dft} T_s = 1/\Delta f$.

In this bandwidth interpretation, I think the term "virtual subcarrier" is more appropriate.

What you send over by electromagetic wave is data in $N$ active subcarriers $< N_{dft}$, or the frequency positions of null/virtual subcarriers can be used for other systems; because when you revert back the time domain signal to the frequency domain, what you need is the frequency positions of these active subcarriers.

It is easier to think OFDM as FDM. In this case you have $N$ narrowbands of bandwidth $B_{nb}$ and a sample rate $F_s > N B_{nb}$ is required for practical implementations, such as easily-realizable anti-aliasing filters. OFDM is simply one implementation of this system by using $\Delta f = B_{nb}$ and $F_s = N_{dft} \Delta f = N_{dft} B_{nb} > N B_{nb}$.


It's really hard to say which one is the best unless we define more precisely what we mean by "best." Some algorithms measure performance using MSE, others use more fancy metrics like entropy.

Since you know the angle and direction of blur, it seems you have a good idea about the blurring function (point spread function). I would start off with a minimum mean squared error type filter (Wiener filter).


The reason the Radix-4 FFT is of interest is in the simplicity of multiplying by $\pm j$ in actual implementation. Below shows the Radix-4 4 point DFT core processing element as part of the Radix-4 FFT Butterfly in comparision to the Radix-2 FFT butterfly (with 2 point DFT core processing element) and the resulting decrease in number of operations, applicable when the input signal is of a length that is a power of 4 (or for portions of the signal that are).

To further see how the radix-4 algorithm reduces multiplications consider the number of real multipliers and additions in a complex multiplier:
$$V_1V_2 =  I_{out}+jQ_{out}= (I_1+jQ_1)(I_2+jQ_2) = (I_1 I_2 - Q_1 Q_2) + j(I_1Q_2+I_2Q_1)$$
Where we see from above there are 4 real multipliers and 2 additions (Where the real and imaginary components are maintained on separate paths so not actually added).
Now consider if $V_2$ above was $j$, then $I_2 = 0$ and $Q_2 = 1$:
$$V_1V_2 =  V_1j = (I_1 0 - Q_1 1) + j(I_1 1 +0 Q_1) = -Q_1 + jI_1$$
We see that multiplying a complex number $I_1+jQ_1$ by $j$ is done by simply changing the sign of $Q$ and then swapping $I$ and $Q$. Similarly multiply by $-j$ would be done by changing the sign of $I$ and then swapping $I$ and $Q$.
The radix-4 allows for an increased number of multiplications by $\pm j$ where the radix-2 solution minimizes all complex multiplications including multiplying by $\pm j$.


Assess quality of a bandpass filter without human intervention

The Remez exchange algorithm not only returns the optimal filter coefficients but also the approximation error $\delta$. This error determines the maximum amplitude error in the pass bands and in the stop bands. The actual amplitude error in band $i$ is given by $\delta/W_i$, where $W_i$ is the error weight in band $i$.

If the algorithm converges, then you know for sure that the maximum errors $\delta/W_i$ occur at the band edges of band $i$. And you know the band edges because they are part of the specs.


Implementing MVDR beamformer in the STFT domain

In this context $\Phi_y$ often describes the (estimated) power spectral density matrix  of $y(f)$, which is
$$\Phi_y(f) = E\{y(f) y^H(f)\} =
\begin{bmatrix}
 y_1 \cdot y_1^*(f),& y_1 \cdot y_2^*(f),& \dots    \\
 \vdots        & \ddots    &          \\
 y_M \cdot y_1^*(f),&               & y_M \cdot y_M^*(f)  \\
   \end{bmatrix},$$
where $M$ is the number of channels, and each element of the matrix corresponds to the cross-power spectral density of two channels for this frequency. Consistent estimates of PSDs in blockwise STFT processing can be obtained, for example, using Welch's method (pwelch in MATLAB).

I am also unsure about your steering vectors and complex conjugate transposes (which I assume is the $\cdot^H$ in your notation). They need to make the steering vector compensate the delay for a specific direction.

From my memory, I think the formulas should be

$$ h_\text{MVDR}(f) = \frac{\Phi_y^{-1}(f)d(f)}{d^H(f)\Phi_y^{-1}d(f)}$$

and

$$ Z(f) = h_\text{MVDR}(f)^T y(f).$$

I think the orginal MVDR paper is this one by Capon. You might want to double-check notations since a complex conjugate does change the beamformer's output quite a bit.

Edit: Thanks to Stanley Pawlukiewicz for pointing out my error!


@omar: I don't know how you obtained your expression for the 'h' coefficients but it's wildly incorrect.  Study the MATLAB documentation regarding the 'firpm()' and 'fir1()' filter design commands.  Experiment, experiment, and experiment again until you understand those commands. Then when you try to use one of those commands to design your desired filter you'll find that your filter needs a truly HUGE number of taps.


What you're seeing is acutance, as described in the wikipedia article on bicubic interpolation.

However, for downsampling an image, cubic interpolation is perhaps not the ideal choice, but makes more sense for upsampling.

Note that if you're downsampling by more than a facor of 4, there will be many pixels in the original image that are simply never used in the downsampled version.

The easiest way to get around this is to downsample by iteratively taking the mean value of 4 neighbouring pixels, until you're within a factor of 2 of the desired image size. Then you can do the last part by cubic interpolation.

If you're want high quality resampling, try taking a look at the ImageMagik documentation for inspiration. In particular they recommend using windowed filters for downsampling, such as Lanczos or Sinc.


What are the main advantages and disadvantages of the lifting scheme?

Unlike the convolution, the lifting scheme can

map integers to integers (CDF 5/3 transform in case of JPEG 2000),
asymptotically reduces the computational complexity by a factor two,
can be computed in-place,
allows to simply treat signal boundaries (periodic symmetric extension in case of JPEG 2000).


Not a definite answer, but here's some analysis on the data:

strct=load('/tmp/mozilla_mmatthe0/fhss_1.mat');
FHSS_1=strct.fhss_1;

start = 14653;
Ns = 3650;
Ncp = Ns/8;
N_total = round(Ns + Ncp);

symind = 10;

X1 = FHSS_1(start+symind*N_total+(1:Ns));
subplot(211);
plot(mag2db(abs(fft(X1))))

fstart = 2593;
fend = 2701;
xlim([fstart-50, fend+50]);

QAM = X1(fstart:fend);

subplot(223);
plot(QAM, 'x-')
axis equal;

subplot(224);
delta = QAM(1:end-1) ./ QAM(2:end);
hold off;
plot(abs(delta)); hold on;
plot(angle(delta), 'k-');
plot(abs(QAM), 'r-');
legend({'abs(delta)', 'angle(delta)', 'abs(QAM)'});

My observations:

the SNR seems to be fairly high (the noise floor in the spectrum is much below the signal part)
There seems to be a strong correlation between adjacent symbols along the subcarriers. I mean, when drawing the constellation symbols along the subcarriers, and connecting the adjacent subcarriers with a line (i.e. a normal "plot", I get a somehow spiral (bottom left figure).

My conclusion: The QAM symbols are highly structured. This does not seem like a data transmission. Instead, it appears to be something like a preamble/synchronization sequence, and these dont need to follow a given constellation. So, the answer might be, that there is no constellation behind.

This answer is to start some discussion on the data. It's not definite.


Can you turn a square wave into a sine wave using a low pass filter?

To be able to analyze what a low pass filter does first you would need to understand what a Fourier transform is, hence some theory first.
The Fourier transform essentially represents the time-domain information stored in some function (a square wave in your case) in terms of frequencies.
A simple example would be a sine wave $sin(2\pi f_{c} t)$ which contains just a single frequency, $f_{c}$. In the frequency domain, after taking the Fourier transform, it produces 2 frequencies at $f_{c}$ and $-f_{c}$ (the negative value occurs because of the mathematical representation of the Fourier transform using complex exponentials, but nevertheless the positive frequency $f_c$ is what one would expect).
Sine wave frequency response http://www.gaussianwaves.com/gaussianwaves/wp-content/uploads/2013/12/Sine_wave_power_calculation_matlab_2.png?02772a
If you are clear with that, then let's move on to what a filter is. A filter modifies the input and generates some output from it and is often represented in the frequency domain => A filter modifies the Fourier transform (aka frequency spectrum) of an input. Basically, it may attenuate some parts of the frequency spectrum and magnify other parts, and may add a "phase change" (frequencies manifest themselves as complex numbers, and phase changes become relevant then)
Now, lets see what exactly a low pass filter does to a square wave.
A low pass filter can be understood as a filter that cuts-off frequencies above a certain threshold when applied to some input. High frequencies which correspond to rapid transformations in the time domain are removed from the input and only the lower frequencies remain.

(source: www.sfu.ca)
In the image attached above, the function on the left is a square wave and the one on the right is the magnitude plot of its corresponding Fourier transform (as mentioned before, the frequencies are actually complex numbers and the entire Fourier transform would be complete only with the addition of a "phase plot" as well).
Observing the magnitude plot, if one were to pick a low pass filter with cut-off frequency between that corresponding to bin #1 (1st harmonic: 100 Hz) and bin #2 (3th harmonic: 500 Hz), the other higher frequencies would get cut out and there would only be 2 frequencies that remain: -100 Hz and 100 Hz which would correspondingly represent a sine wave (as mentioned above).
Answer:
It is therefore viable to use a low pass filter to create a sine wave out of a square wave. However one is restricted to a sine wave with the frequency the same as the period ($1/f_o$) of the input square wave. Higher harmonics cannot be generated (a band-pass filter would be required for the same).
Image sources: google
EDIT: Thanks to @Matt for pointing out that the question says square wave and not rectangular function. Edited answer to reflect the same.


Finding delay where squared error between signal and data is minimal

Just doing something like:
$$
e_\delta = \sum_{n=0}^{N-1} \left| x[n-\delta] - x_\delta[n]\right|^2
$$
for several values of $\delta$ and choosing
$$
\hat{\delta} = \arg \min_{\delta} e_\delta
$$
(i.e. the value of delta that minimizes $e_\delta$) will probably be good enough.
Example signals below:

gives the following error plot

The error plot seems to give the right value of about 5.5 samples for this contrived example data.

Code Only Below
from numpy import log10, asarray, polyfit, ceil, arange, exp, sin, pi, log, random, sum, diff
import matplotlib.pyplot as plt
from scipy.signal import kaiserord, lfilter, firwin, freqz

T = 100
time_period = list(arange(1,T))

K1 = 10
K2 = -10
tau = 200
measurement = [K1 + K2*exp(-x/tau) for x in time_period]

fir_filter = [0,0,0,0,0,0.5,0.5,0,0,0,0]
channel_1 = lfilter(fir_filter,1, measurement)

max_delay = 20

def error_calculation(signal, delayed_signal, delay):
    error = [ (signal[time-delay]-delayed_signal[time])**2 for time in list(arange(delay,T-delay)) ]
    return sum(error)

plt.figure(0)
plt.figure(figsize=(20,20))
plt.plot(time_period, measurement)
plt.plot(time_period, channel_1)
error = []
for delay in arange(1,max_delay):
    error.append(error_calculation(measurement,channel_1,delay))

plt.figure(1)
plt.figure(figsize=(20,20))
plt.plot(arange(1,max_delay), error)


Minimum sampling frequency, quantization, and bitrate calculation

It is not clear in which part of the problem statement you experience issues.
However, a few hints may help you get the gist of such problems

The sampling frequency, it is intuitive that in order for a system to monitor the state of another system it must poll it at higher rates than the highest rate in which the monitored signal is able to change its state. But what is the minimum sampling frequency that gives our sampling system the ability to correctly approximate the continuous signal in the discrete time domain?

This frequency threshold is at two times the highest frequency of the continuous signal:
$$
f_{sampling} \gt 2\cdot f_{bandwidth}
$$
The theorem that the problem statement suggests is the Shannon-Nyquist sampling theorem, the theorem that the above inequality comes from.
In real sampling systems, the sampling frequency is never taken exactly on this threshold, usually the threshold is taken a little bit higher than $2\cdot f_{bandwidth}$, for example at $2.2\cdot f_{bandwidth}$.
However, in theory it is safe to assume a sampling frequency at the Nyquist rate (i.e. the aforementioned frequency threshold).
Your analogue signal's highest frequency is $14.5\, Hz$, ranging from very low frequencies so it is safe to assume, this actually the bandwidth frequency.
All the above summarized in the following equation:
$$ f_{sampling} = 2\cdot f_{bandwidth}$$

Bitrate, in order to calculate the bitrate of your sampling system it is necessary to know the sampling frequency. Fortunately, this was the first the step of the exercise. The second component for the bitrate calculation is the actual length of each sample in bits. The minimum number of bits that are necessary in order to represent a quantized value in x distinct levels is given by:

$$
2^{bitlength}\, =\, x
$$
or
$$
bitlength = \lceil\log_2{x}\rceil
$$
You could actually deduce the bitrate equation by performing a dimensional analysis, because it is stated that it is measured at kbps(kilobits per second). For instance, we know that $Hz$ is measured in $\frac{samples}{sec}$ and kb stands for kilobits meaning $1000\, \frac{bits}{sample}$, so:
$$
bitrate = bitlength\cdot f_{sampling} \left[\dfrac{bits}{sample}\dfrac{samples}{sec} \right]
$$
or in kbps
$$
bitrate = bitlength\cdot f_{sampling}\cdot 10^{-3} \left[kbps\right]
$$


How do I optimize the window lengths in STFT?

This is the classic "uncertainty principle" of the Fourier Transform. You can either have high resolution in time or high resolution in frequency but not both at the same time. The window lengths allow you to trade off between the two.

If you want to detect "events" in your EEG signal with a resolution of say 10ms, then this should be your window length. This will give you a frequency resolution of about 100 Hz.


Given MFCC values for some wav’s, how can I easily tell if all speakers are the same?

In general, you would deal with this kind of problem using some form of Clustering.

The fundamental premise behind this idea is that the features of similar data items would have similar values.

How you assess similarity, how you assess "a cluster" and what this means for the particular application you are dealing with are very much context specific issues.

In this particular example, the features are the MFCC coefficients. If these are similar between two frames we would tend to think that the frames contain more or less signals produced from the same source.

The simplest clustering to apply would be k-means clustering but this is not going to help too much here, especially if you don't know how many speakers you might get in a recording. So, long story short, you will most likely gravitate to something like hierarchical clustering.

You don't have to extract the means of each feature to cluster them. Doing so would be equivalent to training and using a minimum distance classifier. The only thing that is needed here to complete that is the "ground truth" knowledge that will tell you which centroids belong to which speakers.

Note here that if you have two speakers talking in a mostly silent room, the "feature space" of the clustering is not necessarily expected to be divided into two distinct "point clouds". It could, for example, contain separate regions for background silence, one speaker talking, the other speaker talking, both speakers talking at the same time. This last class might contain further grading depending on how loudly was each speaker talking while interrupting the other.

Hierarchical clustering allows you to focus on details with varying degrees of "attention", provided that the MFCC clusters the feature space with sufficient accuracy.

The other thing you can do, is a Principal Component Analysis on the MFCC matrix which (hopefully) will return one component for each one of these "major events" happening in the file.

Hope this helps.

EDIT:
Forgot to add, if you are looking for the terminology that could return to you some useful results in varying levels of complexity, this is "Speaker Diarisation". See for example here.


Spatial receive diversity: a different view

multiple options

So, we're talking about Selection Combining (SC); there's other combining methods, too!

independent receptions

Good that you're stating that, this will be handy in a minute.

Is there a mathematical derivation that tackles this subject from this point of view?

The very definition of a CDF tackles that, but only if the channel realizations have a distribution that have a finite magnitude maximum, which is usually not the case.
Let $|h_i|$ be the magnitude of the $i$th of the $N$ individual channels. Then, the cumulative density function of the "best of all $N$ channels" is simply defined as $\eqref{def}$
\begin{align}
F_{\text{SC},N}(x) &= P\left(\max_{i=1,\ldots,N}|h_i| < x\right) \label{def}\tag{1}\\
&=P\left(|h_1|<x \,\wedge\, |h_2|<x\,\wedge\,\ldots\,\wedge\,|h_N|<x\right)\tag2{}\\
&=\prod_{i=1}^N P\left(|h_i|<x\right)\tag3&\text{due to independence}\\
&=\prod_{i=1}^N F_{|h_i|}\left(x\right)\tag3 & \text{definition of "CDF"}\\
&=\prod_{i=1}^N F_{|h_1|}(x)\tag4 &\text{all identically distr.}\\
&= \left(F_{|h_1|}(x)\right)^N\tag5\label{ord}
\end{align}
I'm almost certain you've seen $\eqref{ord}$ already – you're looking for the concept of diversity order.
Now, let's see about your infinite number of antennas. Two ways to say the same thing:

You know that this doesn't make sense. If I'm sampling any continuous distribution, such as an absolute normal distribution, an infinite amount of times, I will surely hit any non-zero-mass interval. So, even from basic definition of probability, the answer is already clear: you can assume that you get an arbitrarily large channel coefficient! Congratulations!
$\lim_{N\to\infty} k^N = 0\quad\forall |k|<1$, no matter what $k$ actually is. But if $k=F(x)$ is a cumulative probability function (CDF), it's always $<1$ (unless you hit the maximum value), by definition of CDF. In other words: you surely never get anything but 0; that means, the probability of getting a value out of your diversity combining that is smaller than any value that's not the maximum of your possible channel realizations is 0, which means you always get the maximum possible channel.

This makes no sense physically, because it especially means that with finite transmit energy, one of your infinitely many receive antennas gets infinitely much receive energy, if you assume any fat-tailed distribution of channel coefficient magnitudes (e.g. normal). That's (one of the many reasons) why your considerations don't matter for communication engineering.


Sine Wave spectrum showing multiple frequencies - Audacity

Following could be the reason for extra frequencies -

Non Ideal characteristics of speaker used to play the generated sine wave.
Non Ideal characteristics of microphones used to record the played sine wave.

Also one needs to make sure that given sampling rate(48kHz) is supported by the speaker/sound card and the recording microphone otherwise extra frequency components will be introduced.


Upsample first by 12 and then downsample by 13, you'll get very close to 37Mhz, 36.92MHz to be precise, should be good enough.

If you want exactly 37Mhz, then upsample by 37 and then downsample by 40, but this is more calculation intensive with not much gain, so better opt for the first strategy. Also antialsiaing demands are less for the first case.


To build on Laurent's answer, here is an example.  Top frame shows an example transient signal: a damped sine wave.  As the signal decays very quickly, the first 0.1s of the signal is the most interesting.  Second frame shows hann window.  Hann window is almost zero over that first 0.1s.  Third frame shows what happens when you apply hann window to the signal, the most interesting part of it is completely thrown away.


MATLAB has a lot of ways to design filters. You can try designfilt as a starting point: https://uk.mathworks.com/help/signal/ref/designfilt.html. You'd choose the options to set up a FIR LPF and give it your desired cutoff and other parameters (stopband attenuation, ripple allowance, method of filter design).
Choosing the cutoff frequency to be a function of the bit rate is not what you want though. Doing this will distort your signal which you are interested in. Instead, for GFSK, I'd recommend simply computing the 99, 95, or 90 percent bandwidth and setting the cutoff frequency to be half of that. MATLAB has an function called obw which gives the 99% bandwidth but you could write a function myObw(signal, percBandwidth) to calculate any percentage you'd like.
For something like this, it is a good idea to see if what you did actually makes sense. For example, plotting your filter's frequency response over the top of the signal's is a good place to start. The picture below illustrates what I mean and directly shows the problem with choosing the cutoff frequency to be half of your bit rate.


Inverse system of a system with integral

It seems that we have the same homework. You probably are Greek.

After having the same question, I came to the conclusion that in order for this system to be invertible, you have to prove that for any given $$x_1(t),  \, x_2(t)$$ the following sentence is correct: $$ x_1(t) \ne x_2(t) \,\, {\Rightarrow} \,\, y_1(t) \ne y_2(t) $$

As like you proved when a math function is $"1-1"$.

For this very example, if you have the following the following two signals:

$$
\
  x_1(t) = δ(t-4)
\
$$

and

$$ x_2(t)=0 $$

it turns out that while initially $ x_1(t) \ne x_2(t) $, we have $ cos(πt/8)x_1(t) = cos(πt/8)x_2(t) $, for any $t$. This means that

$$ \int_{t-1}^{t+1}\cos\left(\frac{\pi\tau}{8}\right)x_1(\tau)d\tau = \int_{t-1}^{t+1}\cos\left(\frac{\pi\tau}{8}\right)x_2(\tau)d\tau $$

and the system is not invertible.


Resolution of Discrete Fourier Transform is 1/T - Mathematical proof?

$\frac{Fs}{N}$ is not equivalent to $\frac{1}{T}$, since $Fs = \frac{1}{T}$,  $\frac{1}{T N} \neq \frac{1}{T}$ for $N>1$ (as it is implied here of course).

The point about the resolution of the Discrete Fourier Transform follows from its definition where an $N$ point transform, "divides the circle" in $N$ discrete points. As you can see $k n$ (the discrete frequency and sample "instance" variables) are multiplied by the constant $\frac{i 2 \pi}{N}$.

Here, it is still possible to "sample" a non-integer harmonic (say $k = 1.267$) but that would not improve the situation if you only have a limited number of true $n$ (samples).

Hope this helps.


Why can't this system be adequately represented using a z-domain transfer function?

A transfer function describes an LTI system. As such, the given system can be described by a transfer function. However, if there are non-zero initial conditions, the system is no longer linear because there's a contribution in the output that does not depend on the input signal but only on the initial conditions. Consequently, the transfer function cannot directly be used to compute the system's response if there are non-zero initial conditions.
Nevertheless, the (unilateral) $\mathcal{Z}$-transform can still be used to compute the system's response, even with non-zero initial conditions by transforming the difference equation and using
$$\mathcal{Z}\big\{ y[n-k]\big\}=z^{-k}Y(z)+\sum_{m=0}^{k-1}z^{-m}y[m-k],\qquad k>0\tag{1}$$
EXAMPLE:
Let's use a simple example with a similar pole-zero cancellation as in the original problem to illustrate the point. Consider a system described by
$$y[n]-y[n-1]=\alpha \big(x[n]-x[n-1]\big)\tag{2}$$
The corresponding transfer function is
$$H(z)=\frac{Y(z)}{X(z)}=\frac{\alpha(1-z^{-1})}{1-z^{-1}}=\alpha\tag{3}$$
Clearly, $y[n]=\alpha x[n]$ is a solution of $(2)$. It is also the only solution if we require the system to be linear. However, it is not the only solution if we allow non-linear systems because there are infinitely many solutions of the form
$$y[n]=\alpha x[n]+c\tag{4}$$
with an arbitrary constant $c$. Note that these solutions cannot be inferred from the transfer function $(3)$.
Let's now use the $\mathcal{Z}$-transform to solve $(2)$ with initial conditions $y[-1]\neq 0$ and $x[-1]=0$. Transforming $(2)$ using $(1)$ gives
$$Y(z)(1-z^{-1})-y[-1]=\alpha X(z)(1-z^{-1})$$
which results in the following $\mathcal{Z}$-transform of the output:
$$Y(z)=\alpha X(z)+\frac{y[-1]}{1-z^{-1}}\tag{5}$$
In the time domain this becomes
$$y[n]=\alpha x[n]+y[-1]u[n]\tag{6}$$
where $u[n]$ is the unit step. Eq. $(6)$ is just a causal version of $(4)$.
This shows that the $\mathcal{Z}$-transform can be used to compute the system's response with non-zero initial conditions, even though the transfer function alone is inadequate for solving the problem.


Someone's written a php interface to OpenCV - https://github.com/mgdm/OpenCV-for-PHP, though I haven't tested it.

As you've tagged face detection I'll also mention that someone has written a php interface specifically to the face detection functionality of openCV too; http://www.xarg.org/project/php-facedetect/

Alternatively you could try (and this would probably offer more flexibility);

[Note: this is just a short list; searching e.g. "python image processing" will give more options]

Python

OpenCV bindings
scipy.ndimage
Scikit-image

Java

OpenCV bindings
ImageJ api

C++

OpenCV
CImg

The top answer to this question also has a good discussion, though I'd be careful about using Matlab in this instance as you might fall foul of their licence restrictions.


Confusion understanding digital signal processor?

The staple processor instruction of specialized DSP processors is a multiply-accumulate instruction. The vector multiplication and summation are performed in a single command fetch, it's a sort of SIMD (but not exactly). Only advanced MCU's have this instruction in their command set. Hence the advantage of DSP processors over general-purpose MCU's w.r.t. ALU performance on typical signal-processing operations, even when comparable technological process are used in their manufacturing. The same for FPGA devices with internal DSP blocks. In short, less command fetch cycles, higher performance.


Making a 3D point cloud from multiple RGB-D images

Simple answer: Use LiveScan3D. Should be easy enough.
Longer answer: You are looking for a way to 'register' multiple point clouds. This is an extremely well studied problem in computer vision. In your case I would consider two scenarios:

Target-based (photogrammetry-like): Place some targets in the world such as markers (see OpenCV Aruco markers). Estimate the marker pose in each camera image (the camera center is set to be the origin). Once the marker pose is computed in all views, you could transform each 3D cloud into the coordinate frame of one of the cameras (choice is arbitrary).

Markerless: The 3D geometry itself is usually rich enough to relate the views together. This is called registration. You are indeed looking into a multiview registration problem. The typical pipelines follow two main steps: (i) pairwise registration; you can use for example FGR for this. Then these individual pairwise transformation estimates have to be linked together. If your view-graph in which there exist an edge for each pair of cameras seeing common structure, is simple, you could compute the minimum spanning tree and chain the transformations along the tree to compute all the coordinate frames with respect to the first one (free to choose).

Note that, with the advances in deep learning we can now solve this problem in an end-to-end manner, efficiently:
Gojcic, Z., Zhou, C., Wegner, J. D., Guibas, L. J., & Birdal, T. (2020). Learning multiview 3D point cloud registration. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (pp. 1759-1769). https://arxiv.org/pdf/2001.05119.pdf
The first option requires you to physically place fiducial tags in the scene but is easier to implement. The second option is more automatic, but might be harder to get to work.


Let's assume we have $ m $ points in the Fourier Domain (Discrete Fourier Domain) for our reference frequency response. We'll designate their values by $ y \in \mathbb{C}^{m} $.

Also we need a filter of $ n $ coefficients to approximate the reference response as good as it can. We'll use $ h $ for that.

We need now the Fourier Matrix, $ F \in \mathbb{C}^{m \times n} $ in order to write the following linear complex least squares problem:

$$ \arg \min_{h} \frac{1}{2} {\left\| F h - y \right\|}_{2}^{2} $$

You could easily write this also in a Sum Form.
The nice thing about it is the solution is give by the classic Least Squares solution - $ \hat{h} = {\left( {F}^{H} F \right)}^{-1} {F}^{H} y $.

Of course, in practice, you need many point in frequency domain to add more and more restrictions. You could also write a continuous version of this problem using the DTFT.


There is a huge gap between the formalism and construction for the continuous wavelet analysis (Mexican hat or Morlet wavelet, the latter being only approximate or complex Gaussian wavelet) and the discrete  wavelet transforms (DWT). Let us make a long story short.

Continuous wavelet transforms (CWT) hinge on a continuous wavelet $\psi$, for which for any (nice) function $f(x)$ can be represented as (or analyzed) a weighted ($c_{a,b}$) sum of translates and dilations  $\psi\left(\frac{x-b}{a}\right)$. Whether the synthesis formula is equal to $f(x)$ requires $\psi$ to satisfy a moderately mild admissibility condition.

Then a question arose: how can we discretize all of the above to get tractable algorithms? In other words, are there $a_0$ and $b_0$ such that projections onto enumerable series

$$\psi\left(\frac{x-k b_0 a_0^j}{a_0^j}\right)$$

with $j$ and $k$ integers can still be perfect? The answer is a big yes, provided that $a_0b_0 < C_\psi$, a constant depending on $\psi$ (this is a king of Pauli-Heisenberg-Gabor inequality, result by Ingrid Daubechies, 1984). This becomes a wavelet frame: a redundant set of vectors that exactly recovers signals.

Now, with DWT, we want a base, i.e. a non-redundant frame. Something that keeps the size of the initial data. From the above harmonic analysis techniques, if you choose an admissible wavelet, there is very little chance that a basis could be derived from it. Haar or Frankling bases were known examples, and Yves Meyer in 1985 derived a non-trivial example, the Meyer basis.

To build discrete wavelet bases, it was discovered that a more efficient method was needed.
On a discrete space, as a given resolution, you except that there exists an orthogonal basis from an orthogonal template of a translated scaling function or father wavelet. Then, you embedded different resolutions,  and you get a 2-scale linear  equation linking  $\psi(t)$ and $\psi(t/2)$. Using orthogonal projections, you can build another  2-scale linear equation linking a mother wavelet at different dilations. Those equaltions tell you than you can implement the above with convolution and subsampling. Consequently, there are polynomial relations between scale. Solving such nonlinear equations yield a series of weights/linear coefficients providing a hierarchy of discrete wavelets that, plugged into the formalism, are genuine CWT. However, in most cases, they don't have nice close formulae like Gaussians or Morlet. The latter in turn have no discrete counterpart.

In Stéphane Mallat "A wavelet tour", chapter one on sparse representations says:

  Continuous versus Discrete and Finite Originally, many signal processing engineers were wondering what is the point of considering
  wavelets and signals as functions,since all computations are performed
  over discrete signals with conjugate mirror filters.Why bother with
  the convergence of infinite convolution cascades if in practice we
  only compute a finite number of convolutions? Answering these
  important questions is necessary in order to understand why this book
  alternates between theorems on continuous time functions and discrete
  algorithms applied to finite sequences. A short answer would be
  "simplicity." In L2(R), a wavelet basis is constructed by dilating and
  translating a single function $\psi$. Several important theorems
  relate the amplitude of wavelet coefficients to the local regularity
  of the signal $f$. Dilations are not defined over discrete sequences,
  and discrete wavelet bases are therefore more complex to describe.The
  regularity of a discrete sequence is not well defined either, which
  makes it more difficult to interpret the amplitude of wavelet
  coefficients. A theory of continuous-time functions gives asymptotic
  results for discrete sequences with sampling intervals decreasing to
  zero.This theory is useful because these asymptotic results are
  precise enough to understand the behavior of discrete algorithms. But
  continuous time or space models are not sufficient for elaborating
  discrete signal-processing algorithms.The transition between
  continuous and discrete signals must be done with great care to
  maintain important properties such as orthogonality. Restricting the
  constructions to finite discrete signals adds another layer of
  complexity because of border problems. How these border issues affect
  numerical implementations is carefully addressed once the properties
  of the bases are thoroughly understood.


The sampling theorem states that $f_\mathrm{S} \geq 2f_\mathrm{max}$, where $f_\mathrm{S}$ and $f_\mathrm{max}$ are the sampling and maximum signal freuqency, respectively. But there's an additional condition: The equal sign only holds if the signal spectrum does not contain a dirac impulse at $f_\mathrm{S}/2$ which is clearly the case in your example. Therefore $f_\mathrm{S} > 2f_\mathrm{max}$ has to be fulfilled.
A reference for this statement can be found here, for example:

If a signal contains a component at exactly B hertz, then samples spaced at exactly 1/(2B) seconds do not completely determine the signal, Shannon's statement notwithstanding. This sufficient condition can be weakened, as discussed at Sampling of non-baseband signals below.
More recent statements of the theorem are sometimes careful to exclude the equality condition; that is, the condition is if x(t) contains no frequencies higher than or equal to B; this condition is equivalent to Shannon's except when the function includes a steady sinusoidal component at exactly frequency B.

The above is an excerpt of an earlier version of the according Wikipedia article.


I would like to direct you to 3 references:

  C. Steger: “Extracting Curvilinear Structures: A Differential
  Geometric Approach”. In B. Buxton, R. Cipolla, eds., “Fourth European
  Conference on Computer Vision”, Lecture Notes in Computer Science,
  Volume 1064, Springer Verlag, pp. 630-641, 1996.

  C. Steger: “Extraction of Curved Lines from Images”. In “13th
  International Conference on Pattern Recognition”, Volume II, pp.
  251-255, 1996.

  C. Steger: “An Unbiased Detector of Curvilinear Structures”. IEEE
  Transactions on Pattern Analysis and Machine Intelligence, vol. 20,
  no. 2, pp. 113-125, 1998.

In these works, Steger develops a subpixel curvilinear structure extraction algorithm, which is found to work very well for such images. I have also run the algorithm on your images and here you go:

Lines visualized on original image

Lines only

I set the line thickness to 2 for better visualization, normal contours are at sub-pixel level.


i think it's pretty obvious that only the component with angular frequency of $4 \pi $ gets through the LTI system to the output.  so the output should be $\frac{1}{4} \cos(4 \pi t) $.


Names of system functions in frequency domain

$\Omega$ and $\omega$ are frequencies (in radians). Usually, $\omega$ is used, but if you need to deal with continuous-time as well as discrete-time systems, it's common to use one for discrete time and one for continuous time in order to distinguish the two. There is no real standard as to which one is which. I've seen $\Omega$ used for continuous time as well as for discrete time. Having said that, the standard text book Discrete-Time Signal Processing by Oppenheim and Schafer uses $\Omega$ for continuous-time systems and $\omega$ for the normalized frequency of discrete-time systems.
$H(\Omega)$ and $H(\omega)$ are (continuous-time and discrete-time) frequency responses. Note that in discrete time it's also common to use $e^{j\omega}$ (instead of $\omega$) as the argument of the frequency response. This is discussed in the answers to this question. $H(s)$ is the transfer function in continuous time, and $H(z)$ is the transfer function in discrete time. $H(s)$ is the Laplace transform of the continuous-time impulse response $h(t)$, and $H(z)$ is the $\mathcal{Z}$-transform of the discrete-time impulse response $h[n]$. So with transfer functions it's always clear that the variable $s$ (or $p$ in some older texts) is used in continuous time and $z$ is used in discrete time.
The magnitudes $|H(\Omega)|$ and $|H(\omega)|$ are simply called what they are: magnitudes of frequency responses, or sometimes magnitude responses.
The magnitudes $|H(s)|$ and $|H(z)|$ don't have any special names as far as I know, and they are also rarely used with general complex arguments.


FIR filter compensator when using a CIC decimation filter

There was a similar question, https://dsp.stackexchange.com/a/1551/306, and the following is a subset of the answer from the other post.

Generally, to compensate a CIC filter the inverse of the CIC filters response
can be used to generate the compensation filter.  The CIC has
a response of 2

$$
  H(\omega) = \left| \frac{sin(\omega D/2)}{sin(\omega M/2)} \right|^N
$$

Where D is the differentiate delay, M is the decimation rate, and N is the
filter order (number of cascaded filters).  The inverse can be specified as

$$
  H(\omega) = \left| \frac{sin(\omega M/2)}{sin(\omega D/2)} \right|^N
$$

Once we have the
frequency response of the compensation filter, we can simply choose the length
of FIR filter that we desire.  The length of the FIR is application
specific.  Obviously the longer the FIR filter the better compensation.

The following are plots of this straight forward compensation.

The following is the Python code to create the frequency responses and plots.

    import numpy as np
    from numpy import sin, abs, pi
    import pylab

D = 1; M = 7; N = 3

Hfunc = lambda w : abs( (sin((w*M)/2)) / (sin((w*D)/2.)) )**N
HfuncC = lambda w : abs( (sin((w*D)/2.)) / (sin((w*M)/2.)) )**N

w = np.arange(1024) * pi/1024

G = (M*D)**N
H = np.array(map(Hfunc, w))
Hc = np.array(map(HfuncC, w))
# only use the inverse (compensation) roughly to the first null.
Hc[int(1024*pi/M/2):] = 1e-8
plot(w, 20*log10(H/G))
plot(w, 20*log10(Hc*G))
grid('on')

See, Altera, "Understanding CIC compensation filters", for other approaches and $sinc^{-1}$ approximation 1.

1 Altera, "Understanding CIC compensation filters"

2 R. Lyons, "Understanding Digital Signal Processing",
    2nd ed., Prentice Hall, Upper Saddle River, New Jersey, 2004


Question on the definition of one-sided PSD in AWGN channel

The relation between one-sided and two-sided noise PSD is correct. But note that the power of a complex random variable $n$ is given by $E\{|n|^2\}$, and not by $E\{n^2\}$. So in your case you obtain $E\{|n|^2\}=\sigma^2_x+\sigma^2_y$.


If the low pass filter has a linear phase response, then the group delay it introduces is indeed constant, and it equals $\tau_g=(N-1)/2$ samples, with $N$ being the number of taps. If $N$ is odd, then the group delay is an integer number of samples, and in this case you just have to shift the output signal $\tau_g$ samples to the left as compared to the original signal:

$$y[n]=x[n]-x_f[n+\tau],\quad \tau\in\mathbb{N}$$

This is exactly what the Matlab code in your question does. Don't forget to make sure that the number of filter taps $N$ is odd!

Alternatively, you could delay the original signal $x[n]$ by the filter's group delay before subtracting from it the filtered signal.


Is it possible to extract peak locations in the time domain using help from fourier/wavelet analysis?

You can figure where the peaks would be if the other tone wasn't there by an accurate estimate of the phase.  This location will be a rougher estimate in the presence of the other tone.  When you add a slanted signal (say a line) to a mode (a peak), it will shift the location of the peak somewhat.  Since you are interested in the higher frequency tone, the lower frequency one is going to look like a broad hill in comparison which is going to be less interfering than if you were doing it the other way around.
(I'm using the word "tone" broadly, meaning a pure sinusoidal signal.)


If you refer to the documentation for wavedec, a signal $x$ is decomposed on one level into two sets of coefficients: $cA_1$ and $cD_1$. They correspond to a low-pass and a high-pass filter applied to $x$, following by a downsampling. As a result, if you reconstruct $x_1$ with waverec from $cA_1$ only (setting $cD_1$ to zero), and $x_2$ from $cD_1$ only (setting $cA_1$ to zero), $x_1$ will mostly correspond to the lower half of $x$ spectrum, and $x_2$ the upper half.

The same reasoning works on several levels: if your signal has a range of frequency in $[0\,,f]$, $cD_1$ gathers coefficients mostly from $[f/2\,,f]$, $cD_2$ gathers coefficients mostly from $[f/4\,,f/2]$, etc.

So for a sampling frequency of $44100$ Hertz, the bands would be:

$cD_1$: $11025 \to 22050$
$cD_2$: $5512.5\to 11025 $
$cD_3$: $2756.25 \to 5512.5$
$cD_4$: $1378.125  \to 2756.25$
$cD_5$: $689.0625 \to 1378.125 $

If you want to filter out a frequency band, you can zero wavelet coefficient whose spectrum intersect that frequency band, and reconstruct the data. This is a form of thresholding in the wavelet domain. As you guessed, it can work in a way similar  to Fourier.

Indeed, thresholding and shrinkage are very effective with wavelets, possibly more than with a Fourier transform, for denoising. In the wavelet domain, you can design the shrinkage to preserve specific time intervals, allow smooth transitions, etc.

But the wavelet filters are imperfect filters. And downsampling cause aliasing, causing a not-so-clean filtering. To perfect a pure band-pass filter, I would not recommend the DWT (discrete wavelet transform), unless the wavelet is of quite high order.


Look at the equation for the frequency magnitude response of a sinc filter. It will be of the form:
$$\frac{\sin(\alpha)}{ \alpha}$$
Look carefully at the $\alpha$ term and determine for what frequencies the $\sin(\alpha)$ numerator will be equal to zero. The magnitude response notches occur at those frequencies.


What is the estimator used in this paper for system identification?

Interesting paper!

Q1: yes, seems like OLS, non-blind. Usually OLS is non-blind. Some sort of regulated least squares can be used in a bayesian ML-setting, blind. This paper looks interesting: http://www.cs.berkeley.edu/~jordan/papers/lindsten-etal-sysid12.pdf

Q2: PBRS: https://en.wikipedia.org/wiki/Pseudorandom_binary_sequence. There is a lot of literature on this, e.g. Coleman Brosilow, Babu Joseph, "Techniques of Model-Based Control"

Q3:they say it is semi-blind since the dynamics of the sequence generator is known.


Why is there a difference between MATLAB's FFT of a zero-padded signal compared to the non-padded case

If you write

FFT(data)

you basically sample a continous "spectral curve" at length(data) many points. If you do zero padding

FFT(data,len)

where length(data)<len you effectively sample the curve at more points (len) that are more closely spaced. But it is still is the same curve. It was just sampled at different frequencies.

If this is really the only difference you made, then you're probably witnessing a "sampling issue". By that I mean the curves appear to be different because you've sampled it at different points and there is some oscillations going on. Chances are that you havn't applied any windowing which would smooth out the curve a bit and make this effect (otherwise known as spectral leakage) less likely.

You could try one of the following things:

do even more zero padding for an increased sample rate (samples/Hz) and check out the shape
apply a window function to smooth out the curve before the FFT
use PWELCH

PWELCH does more than simply an FFT. It includes windowing and averaging of overlapping but smaller windows. This is probably what you should be doing given your "non-periodic" hint.


What is the right algorithm to detect segmentations of a line chart?

Basically, it looks like you need to know some threshold to detect peaks:

detect bottom edge (Xbe)
calculate the moving average (until the last edge)
detect peak itself (current value - average value > threshold)
start of peak found: X1=Xbe
detect peak end (current value - average value < threshold)
detect bottom age (Xbe)
end of peak found: X2=Xbe

Bottom edge detection in the simplest way is when you make moving subtract for two consecutive values and sign changes from negative to positive.


Controllability's dependence on B matrix

There are different ways of expressing a control system. For example,

Observable Canonical Form
Controllable Canonical Form

If you express your system in controllable canonical form, then $B$ take the form of a column matrix with all zeros except one in the top position.

Link: https://en.wikibooks.org/wiki/Control_Systems/Standard_Forms


Precise 5th and 7th harmonics of a sampled sine wave

This answer discusses the harmonic spectra of the quantized sequence in five cases:

limit $f/f_s \to 0$,
synchronous sampling of a cosine with rational $f/f_s$,
synchronous sampling of a sinusoid of arbitrary phase with rational $f/f_s$,
asynchronous sampling of a sinusoid with rational $f/f_s$, and
asynchronous sampling of a sinusoid with rational $f/f_s$ in the limit $L/A \to
0$,

where $f$ is the frequency of the sinusoid, $f_s$ is the sampling frequency, $L$ is the fundamental period of the quantized sequence, and $A$ is the amplitude of the sinusoid. Decibels full-scale (dBFS) shall be referenced to the sinusoid amplitude $A$ at 0 dBFS.

Limit $f/f_s \to 0$

Let the zero-to-peak amplitude of the continuous-time prototypical sinusoid, which is being sampled and quantized, be a positive real number $A$. In the limit $f/f_s \to 0,$ where $f$ is the sinusoid frequency and $f_s$ is the sampling frequency, the amplitudes of the harmonics are found in the Fourier series representation of the continuous-time piece-wise constant quantized waveform. The terms of the Fourier series are directly the harmonics of the fundamental frequency. We can construct the waveform as a sum of components like the red curve here, illustrated for sinusoid amplitude $A = 7$:

Figure 1. prototypical sinusoid, of amplitude $A=7$, and its quantization in blue, with a component highlighted in red.

We are working with cosine instead of a sine because the math is nicer this way. Because of the symmetries, each component has only odd harmonics and only cosine terms in its Fourier series. The non-zero coefficients of the Fourier series of a component of integer amplitude $k \in 1\ldots \operatorname{round}(A)$ are given by:

$$a_n = \frac{4}{\pi}\int_{x_0}^{x_1}\cos(nx)\,k\,dx = \frac{4k}{\pi n}\big(\sin(x_1n) - \sin(x_0n)\big),\quad n\text{ odd},\tag{1}$$

with the component having the said amplitude over the interval $x_0 < x < x_1$ in the first quarter-period of the cosine. I did not bother to explicitly write the contributions from the symmetrically arranged non-zero pieces of the component in the other quarters, because they will contribute identically to the odd harmonic cosine terms in the Fourier series. Instead I simply included the implicit factor 4 in the equation.

The Fourier series of the full piece-wise constant waveform is a sum of the Fourier series of the components. The boundaries of the pieces that we need to include in the sum are:

$$\begin{gather}0 < x < \operatorname{acos}\left(\frac{\operatorname{round}(A) - 0.5}{A}\right),\quad \text{if }k = \operatorname{round}(A),\\
\operatorname{acos}\left(\frac{k + 0.5}{A}\right) < x < \operatorname{acos}\left(\frac{k - 0.5}{A}\right),\quad \text{if }k \in 1\ldots \operatorname{round}(A)-1,\end{gather}\tag{2}$$

where $\operatorname{acos}\left(\frac{s}{A}\right)$ comes from solving $A\cos(x) = s$ in the first quarter-period. The non-zero coefficients of the Fourier series of the full piece-wise constant waveform are then given by:

$$\begin{gather}a_n = \frac{4\operatorname{round}(A)}{\pi n}\,\sin\bigg(\operatorname{acos}\left(\tfrac{\operatorname{round}(A) - 0.5}{A}\right)n\bigg) +\\\frac{4}{\pi n}\sum_{k=1}^{\operatorname{round}(A)-1} \Bigg(k\sin\bigg(\operatorname{acos}\left(\tfrac{k - 0.5}{A}\right)n\bigg) - k\sin\bigg(\operatorname{acos}\left(\tfrac{k + 0.5}{A}\right)n\bigg)\Bigg),\quad n\text{ odd}.\end{gather}\tag{3}$$

The rounding convention does not matter as long as it is to a nearest integer.

8-bit quantization

It is not out of reach to compute Eq. 3 for $A=2^{23} - 1$ and $n \in {5, 7},$ as parameterized in the question, using something like Python's mpmath. Let's try that but with a 8-bit sinusoid first, because we can then evaluate $a_n$ at multiple $n$ in a reasonable time and perhaps learn something from the results. In Python:

import numpy as np
import mpmath as mp
def a_n(A, n):
    if int(n) & 1:
        return 4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
    else:
        return 0

A = 2**7-1  # amplitude (positive real)

for n in np.arange(1, 10**4, 2):
    print(str(n)+","+str(20*mp.log10(mp.fabs(a_n(A, n))/A)))  # n, a_n (dBFS)

Figure 2. Amplitudes of harmonics of the quantized prototypical sinusoid of amplitude $A=2^{7}-1$. Only odd harmonics are present.

It took some minutes on my PC to get the results (Fig. 2). The graph looks jagged. I think this is because the waveform is similar to a sawtooth riding on a sinusoid, and near the zero crossing the sawtooth frequency stays approximately the same over a large portion of the waveform, being $2\pi A \approx$ 797.964534 times the fundamental frequency at the zero crossing as an asymptotic approximation in the limit $A\to\infty$ of the frequency $\pi/\operatorname{asin}\big((k + 0.5)/A\big)$ calculated based on the distance between the discontinuities of the saw wave just around the zero crossing. Indeed, the largest amplitude is found nearby at $|a_{787}| =$ -68.30 dBFS. This is much higher than the amplitude of a "random" early harmonic. Such peaks are repeated later at approximate harmonics of the sawtooth frequency but less clearly and with a roll off arising from the factor $1/n$ in Eq. 3, making the asymptotic decay -20 dB/decade (~-6 dB/octave) as $n\to\infty$. The amplitude of the fundamental frequency is 0.00066 dBFS, I think because the rounding at the top of the sinusoid is "pulling it up" from the 0 dBFS amplitude of the prototypical cosine.

24-bit quantization

The following continuation of the earlier Python script computes in high precision your desired numbers for $A=2^{23} - 1$, in the limiting case $f/f_s \to 0$:

A = 2**23-1 # amplitude (positive real)

n = 5  # number of the harmonic (odd positive integer)
mp.mp.prec = 70  # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS
mp.mp.prec = 140  # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dBFS

n = 7  # number of the harmonic (odd positive integer)
mp.mp.prec = 70  # Precision 1
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB
mp.mp.prec = 140  # Precision 2
20*mp.log10(mp.fabs(a_n(A, n))/A) # in dB

mp.mp.prec = 53  # Restore default precision

I'm only reporting as many digits as are agreed about by repeated computation using two different mpmath precision settings, including also the digits that got changed by $\pm1$ due to larger changes in later digits. After a total of a few hours of computation on my PC, the result is that with $A=2^{23} - 1$ and $f/f_s \to 0,$ the amplitude of the 5th harmonic is -226.91150085 dBFS and the amplitude of the 7th harmonic is -226.9115030 dBFS. I also computed the amplitude of the fundamental frequency, 3.9195785E-11 dBFS. I also computed in a full day the amplitudes of a number of higher harmonics, in Python:

A = 2**23-1 # amplitude (positive real)

for m in range(1, 20):
  n = mp.mpf(10)**m+1 # number of the harmonic
  mp.mp.prec = 70  # precision
  a_n_low_prec = a_n(A, n)
  mp.mp.prec = 90  # precision
  a_n_high_prec = a_n(A, n)
  print(str(n)+","+str(20*mp.log10(mp.fabs(a_n_low_prec)/A))+","+str(20*mp.log10(mp.fabs(a_n_high_prec)/A)))  # n, a_n

Figure 3. Amplitudes of a selection of odd harmonics of the quantized prototypical sinusoid, with $A=2^{23}-1$.

The corner frequency before the roll-off is again approximately the sawtooth frequency at the zero crossing of the quantized prototypical sinusoid, $2\pi A = 2\pi(2^{23}-1) \approx$ 5.270717225E7 times the fundamental frequency $f$.

Synchronous sampling of a cosine with rational $f/f_s$

Let the prototypical sinusoid, which is being sampled and quantized, be a cosine without a phase shift. That is to say, one of the samples is at the peak of the sinusoid. Let $f/f_s$ be equal to a rational number $c/d$, with $c$ and $d$ integer. The fundamental period $L$ of the periodic sequence of quantized samples is:

$$L = \frac{d}{\operatorname{gcd}(c, d)},\tag{4}$$

where gcd denotes the greatest common divisor. A real discrete Fourier transform (real DFT) of length $L$ has bin frequencies that are the harmonic frequencies of the period-$L$ frequency. The period-$L$ frequency is equivalent to the frequency of the prototypical sinusoid only if the numerator of the irreducible fraction representation of $f/f_s$ equals 1. Otherwise the real DFT bin frequencies contain also other frequencies than those equivalent to in-band harmonics of the sinusoid.

What we have learned from the results so far is that the amplitudes of the harmonics of the quantized prototypical sinusoid do not start to roll off until after the approximately $2\pi A$th harmonic. This means that with rational $f/f_s$, the amplitude of an in-band harmonic frequency of the quantized sequence is typically almost fully determined by aliased frequencies. As an example of such aliasing, with $f_s =$ 48 kHz and $f =$ 1000 Hz, the 5th harmonic gets summed with the aliases of the 43th, 53th, 91th, 101th, 139th, 149th, etc. harmonic:

$$\begin{eqnarray}
&& \ldots\\
&=& 48000\text{ Hz}\times -3 + 1000\text{ Hz}\times 149\\
&=& 48000\text{ Hz}\times -2 + 1000\text{ Hz}\times 101\\
&=& 48000\text{ Hz}\times -1 + 1000\text{ Hz}\times 53\\
&=& 1000\text{ Hz}\times 5\\
&=& 48000\text{ Hz}\times 1 - 1000\text{ Hz}\times 43\\
&=& 48000\text{ Hz}\times 2 - 1000\text{ Hz}\times 91\\
&=& 48000\text{ Hz}\times 3 - 1000\text{ Hz}\times 139\\
&=& \ldots\end{eqnarray}\tag{5}$$

The effect of summation of the harmonics aliased to a frequency depends also on the phase of the sinusoid being quantized, because the phases of the aliased harmonics control whether there is constructive or destructive interference. Aliasing does not change the phase of a harmonic. A phase shift of the quantized prototypical sinusoid is equivalent to shifting it in time, which accounts to a phase shift of each unaliased harmonic proportional to its number. In this answer we will only be sampling and quantizing a prototypical cosine without a phase shift, so to compute a single bin of the real DFT of the discrete-time periodic sequence of samples, it would suffices to sum those $a$ coefficients of the Fourier series of the quantized prototypical sinusoid that alias to the equivalent frequency of the bin. We'll use the normalization convention that $|\operatorname{bin}_1| = 1$ for the real DFT of a sinusoid of amplitude $1$ and of angular frequency $2\pi/L$, where $L$ is the length of the real DFT in time domain. If the numerator of $f/f_s$ expressed as an irreducible fraction equals 1, we have:

$$\operatorname{bin}_n = a_n + \sum_{k=1}^{\infty}\left(a_{kL-n} + a_{kL+n}\right),\tag{6}$$

where $a_n$ is given by Eq. 3 for odd $n$ and $a_n$ is zero for even $n$. With $f = $ 1000 Hz, $f_s = $ 48 kHz, on my PC it would take perhaps thousands of years to compute Eq. 6 to reasonable accuracy for $A=2^{23}-1$. I also have difficulties evaluating the series even with $A=2^3-1$, with the sequence of partial sums oscillating. I think this difficulty arises from that with $A = 2^3 - 1$ and $f/f_s = 1/48$, the sequence before quantization contains values $-3.5$ and $3.5$ exactly at a discontinuity in the quantized prototypical sinusoid. In any case, with rational $f/f_s$, the real-DFT-based approach to determine the amplitudes of the harmonics is much more practical.

24-bit quantization

The following Python script computes using the real-DFT-based method the amplitudes of the harmonics of a synchronously sampled cosine for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz:

import numpy as np
from sympy import Rational
def bins(A, c, d):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
    return np.fft.rfft(waveform)/(A*L/2)

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)]  # real DFT magnitudes (dBFS)

While not required by the script, in this case $c/d = 1/48,$ so each real DFT bin corresponds to a harmonic of the prototypical sinusoid. The fundamental frequency's amplitude measures as -1.75E-7 dBFS, the 5th harmonic as -160.90 dBFS and the 7th as -160.75 dBFS. Unfortunately I could not find a suitable multiple precision Fast Fourier Transform (FFT) library, but I think those three numbers came out in sufficient precision using NumPy's rfft. The precision is not sufficient to calculate the amplitude of say the 3rd harmonic seemingly at about -333 dBFS (Fig. 4).

Figure 4. Magnitudes of the real DFT bins of the quantized sequence of a synchronously sampled cosine with $A=2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz. These results are only applicable when one of the samples is at the peak of the sinusoid, so not applicable to arbitrary phase shifts or asynchronous sampling.

With this $f/f_s$ the fundamental period is $L = 48$. The computed amplitudes of the harmonics are also correct for phase shifts of the prototypical cosine that are multiples of $2\pi/L$, because the cosine is sampled also at those phases, or such a multiple plus $\pi$, because of the symmetry properties of cosine. In the presented case the compatible phase shifts of the cosine are multiples of $7.5°$.

We have encountered here a troublesome combination of $L = 48$ and $A = 2^{23}-1$. We have samples $A\cos(8\times2\pi/L) =$ $A\cos(40\times2\pi/L) =$ $A/2$ and $A\cos(16\times2\pi/L) =$ $A\cos(32\times2\pi/L) =$ $-A/2$, and as $A$ happens to be an odd integer, those samples will always have a fractional part of $1/2$. The above reported amplitudes of the harmonics may thus be dependent on the rounding mode, or, the computation of the phase and/or evaluation of the cosine may have given a small numerical error which has determined the direction of rounding. This may even break the symmetry properties of the sequence, affecting which amplitudes of the harmonics of the period-$L$ frequency are zero.

Method validation

We can validate the two approaches against each other by choosing a pair of $A$ and $f/f_s$ that does not result in sampling of the quantized prototypical cosine at any of its discontinuities. Firstly, this removes having to choose arbitrarily which way a number exactly half-way between integers should be rounded, and secondly, should help in convergence of the infinite series in Eq. 6. Let's try $A = 8$ and $f/f_s = 1/48$. First, the real-DFT-based method, in Python:

import numpy as np
from sympy import Rational
def bins(A, c, d):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d)))
    return np.fft.rfft(waveform)/(A*L/2)

A = 8  # amplitude (positive reals)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
[20*np.log10(float(abs(x))) for x in bins(A, c, d)]  # real DFT magnitudes (dBFS)

This results in the following amplitudes of the harmonics, printed in precision beyond their accuracy, in dBFS:

-inf, 0.005622747208892056, -inf, -331.19944825653926, -inf, -50.41376796795221, -inf, -35.41599672115829, -inf, -327.2878768807577, -inf, -38.14783193548751, -inf, -45.42857299685606, -inf, -324.8188873214554, -inf, -40.41625528549065, -inf, -36.12873038111221, -inf, -337.55682735433885, -inf, -33.04816436790989, -inf

Then, aliasing the harmonics of the quantized prototypical cosine by evaluating Eq. 6, with the series truncated to the first approximately $10^4\pi A$ terms, continuing with the imports and constants of the previous Python script:

import mpmath as mp
def a_n(A, n):
    if int(n) & 1:
        return 4*mp.floor(A + 0.5)/(mp.pi*n)*mp.sin(mp.acos((mp.floor(A + 0.5)-0.5)/A)*n) + 4/(mp.pi*n)*mp.nsum(lambda k: k*( mp.sin(mp.acos((k-0.5)/A)*n) - mp.sin(mp.acos((k+0.5)/A)*n) ), [1, mp.floor(A + 0.5)-1])
    else:
        return 0

mp.mp.prec = 100
L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
if c == np.gcd(c, d):
    for n in range(L/2):
        bin = a_n(A, n) + mp.nsum(lambda k: a_n(A, k*L - n) + a_n(A, k*L + n), [1, int(10000*mp.pi*A)])
        float(20*mp.log10(mp.fabs(bin)/A)) # harmonic amplitude (dBFS)
else:
    print("Error: fundamental frequency must be in bin_1")

mp.mp.prec = 53

which results in the following amplitudes, printed in precision beyond their accuracy, in dBFS:

-inf, 0.005622452922818016, -inf, -165.35811090601203, -inf, -50.41370472648242, -inf, -35.41599394709451, -inf, -155.34167605811237, -inf, -38.14782298030069, -inf, -45.42859376902997, -inf, -164.6586753088176, -inf, -40.41623886311205, -inf, -36.12872396082294, -inf, -149.78722963456818, -inf, -33.04815743325275

The results seem to agree within their numerical accuracy (Fig 6.).

Figure 6. Comparison of amplitudes, in dBFS, of harmonics of a quantized cosine sampled synchronously, with $A = 8$ and $f/f_s = 1/48$, calculated by the real-DFT-based method (blue squares) and by summing coefficients $a_n$ of harmonics that alias to the same frequency (orange crosses).

The 3rd harmonic amplitude calculated by truncating the sum in Eq. 6 converges very slowly as function of the number of terms (Fig. 6). This is probably the case also with the other low-amplitude harmonics. It does not seem feasible to reach the accuracy of the real-DFT-based method using the harmonic aliasing method.

Figure 7. $\operatorname{bin}_3$ as function of the number of terms in truncation of the sum in Eq. 6.

Synchronous sampling of a sinusoid of arbitrary phase with rational $f/f_s$

24-bit quantization

The real-DFT-based method is simple to modify for arbitrary phase shifts of the prototypical cosine. For $A = 2^{23}-1$, $f =$ 1000 Hz, $f_s =$ 48 kHz, and a phase shift $\phi = 0.123$ in radians, the amplitudes of the harmonics of the fundamental period can be computed by the following Python script:

import numpy as np
from sympy import Rational
def bins(A, c, d, phi):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
    return np.fft.rfft(waveform)/(A*L/2)

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
phi = 0.123  # phase shift of the cosine (radians)
[20*np.log10(float(abs(x))) for x in bins(A, c, d, phi)]  # real DFT magnitudes (dBFS)

Figure 8. Magnitudes of the real DFT bins of the quantized sequence of a synchronously sampled cosine with phase shift $\phi = 0.123$, $A=2^{23}−1$, $f=$ 1000 Hz, and $f_s$= 48 kHz.

Asynchronous sampling of a sinusoid with rational $f/f_s$

Asynchronous sampling with rational $f/f_s$ can be modeled by assuming that $f/f_s$ is very close to a rational number. Then, locally the sequence of samples appears periodic with a computable harmonic spectrum. After sufficient time has elapsed, the phase of the prototypical sinusoid will have drifted to a different value and with a different local harmonic spectrum. For a random such sufficient time, the drift is a random variable between $0$ and $2\pi$. A reasonable probabilistic analog of the amplitude of a real DFT bin of the periodic sequence is the expected value of the square of its absolute value, $\operatorname{E}[|\operatorname{bin}_n|^2]$, which can be expressed in the same dBFS scale. We can estimate its value using uniform statistical sampling over phase shifts $\phi$ in a non-redundant range $0 < \phi < 2\pi/L$, where $L$ is the fundamental period of the locally periodic sequence. The end points of the range shall be deliberately left out so that we don't have the previously encountered problem with rounding modes and such.

24-bit quantization

$\operatorname{E}[|\operatorname{bin}_n|^2]$ is estimated for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz, by the following Python script, which also estimates $\max|\operatorname{bin}_n|$ over random phase shifts, representing the worst-case non-fundamental harmonic amplitudes in the local spectra:

import numpy as np
from sympy import Rational
def E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N):
    L = Rational(d, np.gcd(c, d))  # Fundamental period (samples)
    num_bins = int((L + 2)/2)
    sum_sq_bins = np.zeros(num_bins)
    max_sq_bins = np.zeros(num_bins)
    for n in range(N):
        phi = float((0.5 + n/N)*2*np.pi/L)
        waveform = np.around(A*np.cos(np.dot(range(L), 2*np.pi*c/d) + phi))
        sq_bins = np.square(np.absolute(np.fft.rfft(waveform)/(A*L/2)))
        sum_sq_bins = sum_sq_bins + sq_bins
        max_sq_bins = np.maximum(max_sq_bins, sq_bins)
    return [sum_sq_bins/N, max_sq_bins]

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
N = 2**18 # size of statistical sample (positive integer)
results = E_sq_abs_bins_and_E_max_sq_abs_bins(A, c, d, N)
[10*np.log10(float(x)) for x in results[0]]  # est. expected value of real DFT squared magnitudes (dBFS)
[10*np.log10(float(x)) for x in results[1]]  # est. max real DFT squared magnitudes (dBFS)

The results convergence quite slowly as the size of the statistical sample N is increased, and the accuracy requirements stated in the question are almost but not quite fulfilled by a day's calculation on my PC for $\operatorname{E}[|\operatorname{bin}_n|^2]$ and not at all for $\max|\operatorname{bin}_n|$. In the results obtained (Fig. 9) $\operatorname{E}[|\operatorname{bin}_n|^2]$ is about -157.1 dB and $\max|\operatorname{bin}_n|$ about -147.2 to -146.0 dBFS for non-fundamental harmonics. For the fundamental frequency, $\operatorname{E}[|\operatorname{bin}_1|^2]$ is in the order of -0.95E-11 dBFS, and $\max|\operatorname{bin}_1|$ is in the order of 3.4E-7 dBFS or higher.

Figure 9. Estimated $\operatorname{E}[|\operatorname{bin}_n|^2]$ (blue squares) and $\max|\operatorname{bin}_n|$ (red lines) for local $\operatorname{bin}_n$ in slow-drifting asynchronous sampling, for $A = 2^{23}-1$, $f =$ 1000 Hz, and $f_s =$ 48 kHz, estimated by a uniform statistical sample of phase shifts with a sample size of $2^{18}$.

The estimates show that asynchronous sampling gives a smoothing effect on the statistics of the local amplitudes of the harmonics, compared to synchronous sampling with a fixed alignment of the prototypical cosine with respect to the sampling grid. However, with asynchronous sampling, sometimes the worst case local harmonic peaks occur, whereas with synchronous sampling the harmonic amplitudes are well controlled and the worst cases are avoided by a suitable choice of phase shift.

Asynchronous sampling of a sinusoid with rational $f/f_s$ in the limit $L/A\to 0$

In this section of the answer, we will consider only the question's case that $L$ is even and $f/f_s = 1/L$, that is to say that $\operatorname{bin}_1$ corresponds to the frequency of the prototypical sinusoid.

If the length of the fundamental period of the locally periodic sequence of samples $L$ calculated by Eq. 4 is very small compared to the amplitude $A$ of the prototypical sinusoid, then for different phase shifts of the sinusoid, the quantization error can perhaps be considered a sequence of independent pseudorandom numbers distributed uniformly from -0.5 to 0.5, with an antiperiodicity constraint if $L$ is even (and some other constraint if $L$ is odd). Independence of the quantization errors would require that at a number of certain small time lags the autocorrelation of the quantized prototypical sinusoid approaches zero, but this is not something we know to be a fact. If it was true, then in the limit $L/A \to 0$ any sequence consisting of ±0.5's could be found as a quantization error sequence at some phase shift $\phi$ of the prototypical cosine, and within those sequences there would be sequences that maximize $|\operatorname{bin}_n|$ of the real DFT of the quantization error sequence. Due to orthogonality of the real DFT basis vectors, if $n$ does not correspond to the fundamental frequency of the sinusoid, then $\operatorname{bin}_n$ of the quantized sequence equals $\operatorname{bin}_n$ of the quantization error. If $n$ corresponds to the fundamental frequency, then $\operatorname{bin}_n$ of the quantized sequence equals $\operatorname{bin}_n$ of the quantization error plus $Ae^{i\phi}$.

Even if we do not assume independence of quantization errors, we can use the line of thinking of the previous paragraph to find an upper bound of $\max|\operatorname{bin}_n|$, the maximum absolute bin value for a harmonic of the prototypical sinusoid. A worst-case sequence of ±0.5's that maximizes $|\operatorname{bin}_n|$ could be found by an exhaustive search in the binary sequences, but can also obtained directly. A worst-case sequence is any sequence $\operatorname{seq}_k = 0.5\operatorname{sgn}\big(\cos(2\pi nk/L + \chi)\big)$, where $\operatorname{sgn}$ is the sign function, and with a phase shift $\chi$ that does not result in a zero crossing of the expression's sinusoid being at any sample point (the integers). One worst-case sequence is found in the limit of approaching from one side the transition between two worst-case sequences:

$$\begin{eqnarray}\operatorname{seq}_k &=& 0.5\times\begin{cases}\operatorname{sgn}\big(\sin(2\pi nk/L)\big),&\text{if }\sin(2\pi nk/L) \ne 0,\\
\cos(2\pi nk/L),&\text{if }\sin(2\pi nk/L) = 0,\end{cases}\\
&=& 0.5\times\begin{cases}\operatorname{sgn}\big(\sin(2\pi nk/L)\big),&\text{if }2kn \text{ mod } L \ne 0,\\
\cos(2\pi nk/L),&\text{if }2kn \text{ mod } L = 0,\end{cases}\end{eqnarray}\tag{7}$$

where $\text{ mod }$ denotes remainder in integer division in a numerically safe expression. An upper bound of $\max|\operatorname{bin}_n|$ is then found by real DFT of the worst-case sequence.

24-bit quantization

For odd $n$ the sequence given by Eq. 7 satisfies the antiperiodicity constraint for even $L$, and we can calculate an upper bound of $\max|\operatorname{bin}_n|$ of the maxima of local harmonic amplitudes $|\operatorname{bin}_n|$ for $A = 2^{23}-1$, $f =$ 1000 Hz and $f_s =$ 48 kHz, in Python:

import numpy as np
from sympy import Rational
def bins_max_err_seq(A, c, d, n):
    L = int(Rational(d, np.gcd(c, d)))  # Fundamental period (samples)
    waveform = 0.5*np.sign(np.sin(np.dot(range(L), 2*np.pi*n/L)))
    for k in range(L):
        if 2*k*n % L == 0:
            waveform[k] = 0.5*np.cos(k*2*np.pi*n/L)
    return np.fft.rfft(waveform)/(A*L/2)

A = 2**23-1  # amplitude (positive real)
c = 1000  # numerator of f/f_s (positive integer)
d = 48000  # denominator of f/f_s (positive integer)
if c == np.gcd(c, d) and int(Rational(d, np.gcd(c, d))) % 1 == 0:
    for n in np.arange(1, int(d/c/2 + 1), 2):
        print(str(n) + "," + str(20*np.log10(float(abs(bins_max_err_seq(A, c, d, n)[n])))))  # real DFT magnitude of bin n for the maximum quantization error sequence (dBFS)
else:
    print("Error: Fundamental frequency must be in bin 1 and L must be even\n")

For the odd harmonics, for $A = 2^{23}-1$, $f =$ 1000 Hz and $f_s =$ 48 kHz, we get locally $\max|\operatorname{bin}_n| \le $ -142.34 dBFS when the number of the harmonic $n$ is divisible by 3, and $\max|\operatorname{bin}_n| \le $ -142.39 dBFS otherwise which applies to the 5th and to the 7th harmonic. This seems compatible with Fig. 9, with an about 4 dB higher upper bound of the maximum local amplitudes of harmonics compared to the maximum local harmonic amplitude found by numerical search in actual quantization error sequences. A willing eye can also see an elevation in $\max|\operatorname{bin}_n|$ of Fig. 9 when $n$ is divisible by 3.

This answer has reached the maximum allowed length of 30000 characters.


How to represent the nonlinear model as a state space in Unscented Kalman Filter

Rather than write

$$
x_t = Ax_{t-1} + f(u_{t-1},\mathbf{w})
$$

as the state update equation, I'd write:

$$
\xi = \left[ \begin{array}{c}
x_t\\
u_t
\end{array}
\right]
$$

and then
$$
\xi_t = g(\xi_{t-1},A,\mathbf{w})
$$

so that
$$
y_t = \left[ 1\ \ \ 0 \right] \xi_t + v_t
$$

Then you could apply the non-additive noise formulation of the EKF to get your filter equations.


What type of signal do you need to demodulate (BPSK, FSK, OFDM, etc.)? While a lot of practical systems have them, you don't always need an equalizer or phase-locked loop, especially if you're just making a toy simulation where you control all the conditions. If you've never done this before, I wouldn't recommend jumping in with all of that included. Each of those topics that you brought up is a large area of study in itself, so I can't really cover all of the theory in a single answer. I would operate in the following order:

Generate a modulated signal with no noise. Build your receiver and show that you can extract the information data out successfully (since there are no real impairments, you should be able to do this perfectly in simulation).
Add white Gaussian noise to the input signal that your receiver operates on. Test the receiver to calculate its bit error rate at its output as a function of the signal-to-noise ratio in the receiver input signal. Compare these results to theoretical results found in a digital communications textbook (like this or this) to verify that you're operating correctly.
Add a frequency and/or phase offset to the receiver's input signal. Experiment with different techniques for correcting the offset in your receiver (as you noted, a phase-locked loop is one such way). Compare the receiver's bit error rate with and without the phase/frequency synchronizer in place.
Add channel effects to your simulation (i.e. multipath). This can be simulated by convolving the transmitted signal with a simulated channel impulse response. There are a lot of different channel models out there, depending upon the type of scenario that you want to test. At the receiver, add an equalizer that attempts to correct for the channel effects. The easiest-to-implement adaptive equalizer is probably one based on the LMS algorithm. Compare the bit error rate at the receiver with and without the equalizer in place.

This is very high-level and each step is complicated in and of itself, but you're certain to run into trouble if you try to do everything at once with no background.


The phrase "zero mean" is primarily a phrase from the mathematical field of statistics.  For us, a "zero mean" sequence means a sequence of numbers whose mathematical average is zero.

By the way, MATLAB's randn(1,N) command does not generate exactly a zero mean sequence. It generates an N-length sequence of random numbers that fluctuate randomly above and below an amplitude of zero, but the sequence's mean is not guaranteed to be zero.  The larger is N the closer the average (mean) of randn(1,N) will be to zero.  Try these:

mean(randn(1,5))  and  mean(randn(1,500))

I have no idea what "the amplitude interval (0,0 1,0)" means.


Correct magnitude spectra of a cosine DFT?

you're missing the frequency component at $k=N-k_0$.