I don't know if it qualifies as quite "visually stunning", but you might want to check out my blog article:  DFT Graphical Interpretation: Centroids of Weighted Roots of Unity

The concept of the $1/N$ normalized DFT as a center of mass calculation was a major "aha moment" for me.  It is a good answer for "What does the DFT really mean?"

By request, here is one of the figures from my article:

A little explanation is in order.  The top graph is a time domain representation and the polar graphs on the bottom are on the complex plane.  The left most circle is for bin zero, aka the DC bin, the second is bin one, and so on.  The little blue circle is the center of mass and is also the bin value as as complex number.

$$ A = \frac{2 \pi n}{N} $$

This sample has 3 cycles per frame with a phase of 3.  Bin three (the fourth polar graph) clearly shows the bin value has magnitude of $1/2$ and the phase value of 3 is almost $\pi$ and therefore almost halfway around the circle.

There are many more examples, and more thorough explanations with the math in the article.


At which point in the signal processing chain to use Hann Window?

You should apply window in time domain. The purpose is to make your signal go to zero at boundaries of your frame/sample block. It is necessary because DFT assumes signal to be periodic, and zeroing it at the ends makes it so. Obviously if length of your frame is constant, then precompute window values.

Enough words:

It just occured to me that you want to get the frequency resolution of 0.2 Hz - is that correct? You must keep in mind that frequency resolution is given by: $\Delta f = \frac{f_s}{N} $. So unless your sampling frequency $f_s$ is not lower than approx. 410 Hz (unlikely) you won't achieve the resolution of 200 mHz. You need to extend your observation window. Also padding your signal with zeros will not help too much - this will be only an interpolation of your spectrum points.

EDIT:

If you perform any FIR filtering, then apply the window afterwards. Windowing should be done just before the FFT calculation to make your signal 'disappear' at the ends (especially when filtering produces artifacts for first output samples). Regarding your frequency resolution, with 2kHz of $f_s$ and $N$ equal to 2048 you get approximately 0.98 Hz of a resolution. If  you want more than that, then you might think about extending your sample block or decreasing the sampling frequency. Just keep in mind equation I mentioned above.

EDIT 2:

I presume that the function you pasted (and one I moved to your original question) is only for calculation of window. Do I understand you correctly? Or is it function that is performing windowing of your signal (variable ADC_Values[i] is present). Check the following things:

I corrected your equation for window function - notice that it is
dependent only on n - there are no signal samples.
If you want to perform windowing within your function then code
should look more like:

ADC_Values[i] = (0.5 - (0.5 * cos ( (2.0 * PI * i) / (INPUT_SAMPLES - 1)))*ADC_Values[i];

Window calculation itself is given by:

`win[i] = (0.5 - (0.5 * cos ( (2.0 * PI * i) / (INPUT_SAMPLES - 1)));`

Therefore I suggest you to precompute the win once and for all, passing it later to your function and performing ADC_Values[i] *= win[i];. This will save some time.

I am assuming that you making typecasting correctly. If you have any doubts, then extract your signal block before and after filtering, together with window samples and windowed signal. This should help you to debug the code.

EDIT 3:

I am afraid that you still do not understand the principle. Changing of window can only remove leakage effects - kind of noise covering your harmonics. It will not increase your frequency resolution. In the matter of fact, it will make it worse - as rectangular window has narrowest main-lobe. Let me repeat once again - if you want to increase your frequency resolution then you must use the appropriate buffer size and sampling frequency, according to equation $\Delta f = \frac{f_s}{N} $.


Is there an intensity estimation algorithm working solely in time domain

Particle velocity and sound pressure are related by the field impedance. The field impedance is a function of the media and the boundary conditions of the field.

For a reasonably shaped wave in a reasonably well behaved medium, that impedance is simply $\rho \cdot c$, i.e. density of the medium times speed of sound. For air that is about 400 Pa/(m/s). There is no phase shift and not a lot of frequency dependency.

In this case the intensity signal can be calculated as

$$I(t) = \frac{p^2(t)}{\rho \cdot c}$$

This will only get you the magnitude. Keep in mind that intensity is a vector with the same direction as the particle velocity. If you want to integrate Intensity to get actual power you need to take into account the angle between the Intensity vector and the normal vector of your integration surface.

If you only have the pressure signal, you have to make some assumptions around the direction of the vector. Pressure is a scalar and there is no spatial information contained.


Kalman Filter Process Noise - Model Where It Vanished

The process noise $w(t)$, which is typically assumed to be a zero mean, white Gaussian noise with (power) variance $\sigma^2$, is used to account for any mismatches between the assumed dynamic model of the states and the actual truth.

If your assumed constant state (DC) model perfectly matches with the underlying nature of the observed signal, then you should set the process noise to zero; i.e., your model is exact.

Then you only have measurement noise due to the sensor's intrinsic noise.

If you do not set your $\sigma^2$ to zero (adding non-zero process noise to a constant dynamic model) then your state estimation will never reach steady-state, but wander around it; as you are telling the Kalman filter that the true signal is not a constant...


Implementing an L-System conversion from R2 to R1 for symmetric, square matrices

I've defined the following function in R

State 0 indicates that one is located within a right-angle region

State 1 indicates that one is located within a Hilbert curve

State 2 indicates that one is located within a wedge region

rot <- function(n, x, y, rx, ry) {
  if (ry == 1) {
    if (rx == 0) {
      x <<- n-1 - x
      y <<- n-1 - y
    }

    temp <- x
    x <<- y
    y <<- temp
  }
}

xy2d <- function(n, x, y) {
  # rx = x region (right = 1, left equals 0)
  # ry = y region (top = 1, bottom = 0)
  # s = half the dimension
  # d = linear coordinate value
  d=0L
  state = 0L
  for (s in 2L^(log2(n/2L):0)) {
    rx = bitwAnd(x, s) > 0
    ry = bitwAnd(y, s) > 0
    if (state == 0) {
      # *--2  | 0  ->  *--2   | 1 -> 0--3   | 2 -> *--2   |
      # 0--1  |        0--1   |      1--2   |      1--0   |

      # if x on right region, add lower left triangle (region 0)
      #      (s * (s + 1) / 2) * rx
      # if y on upper region, add lower right square as well (both regions 0 and 1)
      #     (s * s) * ry
      d = d + (s * (s + 1) / 2) * rx + (s * s) * ry
      d
      # change the state
      #  rx  ry | (rx + ry)
      #  0   0  |  0
      #  1   0  |  1
      #  1   1  |  2
      state = rx + ry
      state
    }else if (state == 1) {
      # 0--3  | 0  ->  0--1   | 1,2 -> 0--3   | 3 -> 2--3   |
      # 1--2  |        3--2   |        1--2   |      1--0   |

      # Determine how many squares to add
      # equals number in s-by-s square times it's order in the curve
      #  rx  ry | (3 * rx) ^ (1 - ry)
      #  0   1  |  0 ^ 0 = [00] ^ [00] = 00 = 0
      #  0   0  |  0 ^ 1 = [00] ^ [01] = 01 = 1
      #  1   0  |  3 ^ 1 = [11] ^ [01] = 10 = 2
      #  1   1  |  3 ^ 0 = [11] ^ [00] = 11 = 3
      d = d + s * s * bitwXor(3 * rx, 1 - ry)
      d
      # ROTATE
      rot(n, x, y, rx, ry)
      x
      y

    }else if (state == 2) {
      # *--2  | 0  ->  2--1   | 1 -> *--2   | 3 -> *--2   |
      # 1--0  |        3--0   |      1--0   |      0--1   |

      d = d + (s * s) * (rx == ry) + (s * (s + 1) / 2) * (rx == 1 && ry == 1)

      #  rx  ry |  Target | state - (rx + ry)
      #_________|_________|__________________
      #  1   0  |  1      | 2 - (1 + 0) = 1
      #  0   0  |  2      | 2 - (0 + 0) = 2
      #  1   1  |  0      | 2 - (1 + 1) = 0

      # new_state = current_state - (rx + ry)
      state = 2 - (rx + ry)
      if (state == 1) {
        # If the state goes from state-2 to state-1, then rotate the hilbert
        # curve 180 degrees (i.e. reflect horizonally and vertically)
        x <- n-1 - x
        y <- n-1 - y
      }
    }
  }
  return(d)
}

xy2d(8,5,0)

[1] 16
```


Which one of these signals will be least distorted by this communication channel

Normally, a signal is composed of several frequencies, which strengths are given by the Fourier Transform of the signal.

$x_1(t)=\delta(t) \Leftrightarrow \mathcal{F}(x_1(t))=C$, i.e. the Fourier transform is a constant (value of constant depends on the actual normalization of the Fourier Transform). Hence, it has all frequencies contained, with equal power.
$x_2(t)=5 \Leftrightarrow \mathcal{F}(x_2(t))=\delta(f)$, i.e. the Fourier Transform of a constant is a Dirac impulse at frequency zero. Hence, this signal has only one frequency component, the zero-frequency.
$x_3(t)=10e^{j1000t} \Leftrightarrow \mathcal{F}(x_3(t))=C\delta(f-1000/2\pi)$ i.e. the Fourier Transform of a complex exponential is a Dirac impulse at the position of the frequency of the exponential. Hence, it has only one component, at frequency $1000/2\pi$
$x_4(t)=1/t \Leftrightarrow \mathcal{F}(x_4(t))=C\text{sign}(f)$, i.e. the Fourier transform of the inverse is a sign function. Hence, it contains all frequencies.

Accordingly, the answer to your question would be signal 2, as it only contains the DC frequency ($f=0$).


“The Fourier transform cannot measure two phases at the same frequency.” Why not?

It's because the simultaneous presence of two sinusoidal signals with the same frequency and different phases is actualy equivalent to  a single sinusoidal at the same frequency, but, with a new phase and amplitude as follows:
Let the two sinusodial components be summed like this :
$$ x(t) = a \cos(\omega_0 t + \phi) + b \cos(\omega_0 t + \theta) $$
Then, by trigonometric manipulations it can be shown that :
$$ x(t) = A \cos(\omega_0 t + \Phi) $$
where
$$A = \sqrt{ a^2 + b^2 + 2 a b \cos(\theta-\phi) } $$ and
$$ \Phi = \tan^{-1}\left(\frac{ a \sin(\phi) + b\sin(\theta) } { a \cos(\phi) + b\cos(\theta) } \right) $$
hence you actually have a single sinusoidal (with a new phase and amplitude), and therefore nothing to distinguish indeed...


Reduce the number of frequencies after FFT

Apply a sliding window function of 500 samples length to the signal, for example a Hamming window. Choose a suitable overlap between windows e.g. 75%. Calculate the spectrum of each window and then average over all windows. Then the frequency bins lie exactly at the natural numbers as you wanted.

But beware: each bin does not only represent one exact frequency, but the energy in a frequency band, in this case  $[f-0.5,f+0.5]$Hz


Group delay of $H(\omega)= 1- re^{j \theta}e^{ - j \omega} $

The book's formula is right.

Let $$H(w) = 1 - r e^{j(\theta - w)} = [1-r \cos(\theta - w)] + j [-r \sin(\theta - w)]$$ Since the group delay $\tau$ is the negative of the derivative of the phase of $H(w)$, we first define the phase as:
$$\phi(w) = \tan^{-1}\left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)$$

Using the derivative rule for the inverse tangent as:
$$\frac {d}{dw} \left( \tan^{-1}(u)\right) = \frac{u'}{1+u^2}$$

We shall compute the derivative of the phase as;
$$ \phi^{'}(w) = \frac{ \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^{'} }{ 1 + \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^2 }$$

which continues as:

$$\begin{align}
\phi^{'}(w) &= \frac{ r \cos(\theta-w) \cdot[1 -r \cos(\theta-w)] + r \sin(\theta-w) \cdot [-r \sin(\theta-w)]}{(1-r \cos(\theta-w))^2 \cdot \left( 1 + \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^2 \right) }\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 \cos^2(\theta-w) - r^2 \sin^2(\theta-w) }{(1-r \cos(\theta-w))^2 \cdot \left( 1 + \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^2 \right) }\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 } { (1-r \cos(\theta-w))^2 \cdot \left( \frac{ (1-r \cos(\theta-w))^2 + r^2 \sin^2(\theta - w)} { (1-r \cos(\theta-w))^2} \right) }\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 } { (1-r \cos(\theta-w))^2 + r^2 \sin^2(\theta - w)}\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 } { 1 - 2r\cos(\theta-w) + r^2 \cos^2(\theta - w) + r^2 \sin^2(\theta - w)}\\
\end{align}
$$

Which finally simplifes to
$$ \phi^{'}(w) = \frac{ r \cos(\theta-w) -r^2 } { 1-2r \cos(\theta-w) + r^2}$$

Then since the group delay is the negative of this we finally have:
$$ \tau(w) = -\phi^{'}(w) = \frac{ r^2 - r \cos(\theta-w) } { 1 + r^2 -2r \cos(\theta-w) }$$

$$ \tau(w) = -\phi^{'}(w) = \frac{ r^2 - r \cos(w-\theta) } { 1 + r^2 -2r \cos(w -\theta) }$$

Which is the same as the book's formula.

Coming to Matlab's output, you may try simplifications. But they are known to produce more verbose outputs than possible with simplifications.


Why Should the Sum of Masks / Kernel / Filter in Image Processing be One?

When you're trying to create a blurring kernel (namely a weighted mean of the data) you keep the mask sum as one.
This preserves the mean of the filtered image (think about the DFT, if the sum mask is 1 its DC gain is 1 which means it keeps the sum of the image which means it keeps the mean value of the image.)

When you're trying to apply a high pass filter you want to have gain of zero at the DC coefficient.
This is done by having the sum of the mask as zero (same reason as above.)


Angle of arrival (AOA) estimation using FMCW radar using MUSIC algorithm

Yes, it is proper/practical way to obtain correlation matrix. just multiply your sample with its conjugate to get auto-correlation and with conjugate of other antenna elements sample to get cross-correlation. arrange them in matrix to get correlation matrix($N$$X$$N$).
you can also look out for MVDR implementations in matlab which uses similar method for obtaining noise correlation matrix.


how to understand the MIMO zero forcing as a bank of decorrelators?

It seems like you understand your question already but are unsure it is correct? Let the decorrelator be $\mathbf{D}=(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T$, then the symbols post-decorrelator are $\mathbf{z}=\mathbf{D}\mathbf{y}=\mathbf{D}(\mathbf{Hx}+\mathbf{n})$. The first stream, $z_1$, is computed by taking the first row of $\mathbf{D}$, call the row vector $\mathbf{d}_1$ (first row of $\mathbf{D}$), and performing the operation: $z_1=\mathbf{d}_1(\mathbf{Hx}+\mathbf{n})$. You get $z_1=x_1 + \text{filtered noise}$ because $\mathbf{d}_1\mathbf{H}=[1, 0, ..., 0]$.


Opening is the morphological operation used for example in removal of small particles or some noise in binary image. Additionally it is widely used in hand-writing recognition where you want only the 'skeleton' of your letter.

Mathematically it is a dilation of the erosion. Erosion sets current pixel to 0 if any pixel in neighbourhood is 0. Dilation on the other hand sets pixel to 1 if any neighbour is 1.

You can think of that in following way: first you do the erosion ($N$ steps), which is shrinking boundaries of your objects. After that you perform dilation (also $N$ steps), which is expanding boundaries of your objects. Because small ones were removed in erosion step you will bring back biggest elements.

Here is some image showing this operation:


How to use Robust PCA for image inpainting?

As far as I can remember, one of the main interests of using R-PCA for image inpainting is that a pixel accurate defective pixel masks is not necessary.
Intuitively, the main idea is that the observed image has low rank, but that the defects (given by the mask);

disturb this low-rank property
form a sparse distribution.

Thus, image inpainting using R-PCA is given in this context by:

$X = L + S,$

where $L$ is the restored (inpainted) image and $S$ is the defective pixel mask.

This approach was implemented for example for background subtaction (motion detection task): the authors perform a low rank reconstruction of some shopping mall hallway. The moving pedestrians end up in the sparse component (=> mask) and the background is given by the low rank part. The backgournd does not have any "holes" in place of the pedestrians since a hole breaks the low rank constraint.

For more insight about this approach, see Yi Ma's low rank page.


First, maybe you could give us more context (is that homework?).

I think that what is meant by "housekeeping operations" are the data transfers between registers / state variables ("delay propagation") that naturally arise when naively implementing filters. For example, a naive implementation of a 4-tap FIR would look like this:

output = h[0] * input + h[1] * x[0] + h[2] * x[1] + h[3] * x[2]
// "Housekeeping" operations
x[2] = x[1]
x[1] = x[0]
x[0] = input

The data transfers operations can be avoided on some hardware architectures (especially DSP) by using modulo addressing. I think this is what is meant by "suitable addressing strategies".


What is a Kravchuk transform and how is it related to Fourier transforms?

Transliterations of  Ukrainian names have different avatars in English (and in others languages as well). You can find Kravchuk polynomials, and other papers like  On Krawtchouk Transforms or Krawtchouk polynomials and Krawtchouk matrices. You can find as well Kravchuk orthogonal polynomials.

As they form an orthogonal basis of polynomials (as well as many others, listed on polpak: Bernoulli, Bernstein, Tchebychev, Hermite, Laguerre, Legendre, Zernike), they are candidates for a transform. Derived moments are used in image processing, and the following paper seems to have a wide audience:

Image analysis by Krawtchouk moments, IEEE Transactions on Image Processing 2003

  A new set of orthogonal moments based on the discrete classical
  Krawtchouk polynomials is introduced. The Krawtchouk polynomials are
  scaled to ensure numerical stability, thus creating a set of weighted
  Krawtchouk polynomials. The set of proposed Krawtchouk moments is then
  derived from the weighted Krawtchouk polynomials. The orthogonality of
  the proposed moments ensures minimal information redundancy. No
  numerical approximation is involved in deriving the moments, since the
  weighted Krawtchouk polynomials are discrete. These properties make
  the Krawtchouk moments well suited as pattern features in the analysis
  of two-dimensional images. It is shown that the Krawtchouk moments can
  be employed to extract local features of an image, unlike other
  orthogonal moments, which generally capture the global features. The
  computational aspects of the moments using the recursive and symmetry
  properties are discussed. The theoretical framework is validated by an
  experiment on image reconstruction using Krawtchouk moments and the
  results are compared to that of Zernike, pseudo-Zernike, Legendre, and
  Tchebyscheff moments. Krawtchouk moment invariants are constructed
  using a linear combination of geometric moment invariants; an object
  recognition experiment shows Krawtchouk moment invariants perform
  significantly better than Hu's moment invariants in both noise-free
  and noisy conditions.

Later, you can read:

Image Analysis Using Hahn Moments, IEEE Transactions on Pattern Analysis and Machine, 2007, where Hahn moments generalize Chebyshev and Krawtchouk moments:

  This paper shows how Hahn moments provide a unified understanding of
  the recently introduced Chebyshev and Krawtchouk moments. The two
  latter moments can be obtained as particular cases of Hahn moments
  with the appropriate parameter settings and this fact implies that
  Hahn moments encompass all their properties. The aim of this paper is
  twofold: (1) To show how Hahn moments, as a generalization of
  Chebyshev and Krawtchouk moments, can be used for global and local
  feature extraction and (2) to show how Hahn moments can be
  incorporated into the framework of normalized convolution to analyze
  local structures of irregularly sampled signals.

In Wikipedia's Discrete Fourier transform we find:

  The choice of eigenvectors of the DFT matrix has become important in
  recent years in order to define a discrete analogue of the fractional
  Fourier transform—the DFT matrix can be taken to fractional powers by
  exponentiating the eigenvalues (e.g., Rubio and Santhanam, 2005). For
  the continuous Fourier transform, the natural orthogonal
  eigenfunctions are the Hermite functions, so various discrete
  analogues of these have been employed as the eigenvectors of the DFT,
  such as the Kravchuk polynomials (Atakishiyev and Wolf, 1997). The
  "best" choice of eigenvectors to define a fractional discrete Fourier
  transform remains an open question, however.


Minimizing Time Sidelobes with Pulse Compression

You could define a weight sequence $w(n)$ that weighs the error differently for different time indices $n$:

$$\epsilon^2=\sum_{n=0}^{N+K-2} w(n)\big|d(n)-z(n)\big|^2\tag{1}$$

Samples that are too large get a higher weight, whereas samples that are smaller than necessary can get a lower weight.
Of course, you cannot know in advance which samples should get which weight, so you need to iteratively solve the RMS problem minimizing $(1)$, and adapt the weight sequence $w(n)$ in each iteration depending on the magnitude of the error $e(n)=d(n)-z(n)$.

A possible update rule for the weights $w(n)$ is

$$w^{(k)}(n)=\frac{\left|e^{(k)}(n)\right|}{\displaystyle\frac{1}{N+K-1}\sum_n\left|e^{(k)}(n)\right|}\cdot w^{(k-1)}(n)\tag{2}$$

where $k$ is the iteration index. You probably need to experiment a bit with the update equation to adapt it to your desired result, but the general idea is to increase the weight where the resulting error is too large and decrease it where it is smaller than needed.

This type of algorithm is called iteratively reweighted least squares.


Soft question: Why do we need to reconstruct a signal?

The concept of reconstruction has nothing to do with the application, rather it has to do with the question: did I get the same signal that is really there. If you cannot recreate the signal back, that means the conversion process is loosing/modifying underlying information, which in most cases you do not want to happen.

So the confidence on the Fourier Transform, or to use the digitized data in the controller, comes from the fact that the loss of information is within the tolerance limits of your application. You can measure the quality of your digitized signal, by actually re-synthesizing the signal and comparing it to the original signal.


Why is the impulse response function of this system 0?

This system
$$ y(t) = t^2 x(t) $$

is not LTI and therefore does not have an impulse response of the form $h(t) = \mathcal{T}\{\delta(t)\}$.

So your statement $h(t) = t^2 \delta(t)$ is not correct... Hope this solves your confusion.


Every system will be different, so it’s impossible to make a blanket statement that applies to everything, but IFF systems can be used to accomplish the task. The signal processing here can be as basic as looking for a predetermined pulse pattern, or require demodulation/decryption of some “special” signal. The Wikipedia article I linked covers the topic in pretty great detail.


Shifted output of input vs output of shifted input

This is all about consequently substituting variables.
Shifting by $t_0$ means replacing the variable $t$ by $t-t_0$. If $y(t)=x(2t)$, then $y(t-t_0)=x(2(t-t_0))=x(2t-2t_0)$.

If you have a shifted input $x_2(t)=x(t-t_0)$, then the output is $y(t)=x_2(2t)$, i.e., the variable $t$ is multiplied by $2$. Replacing $t$ by $2t$ in $x(t-t_0)$ results in $x(2t-t_0)$.

So the slide is correct.

EDIT: Concerning shifting the output, imagine a triangular input $x(t)$ with its apex at $t=0$. The output is a compressed version of the input, i.e., it is also triangular with its apex at $t=0$. If the output is shifted by $t_0$ the apex is shifted to $t=t_0$. The function $x(2t-t_0)$ has its apex at $t$ satisfying $2t-t_0=0$, i.e., at $t=t_0/2$, so the correct formula for the shifted output is $x(2t-2t_0)$.


Your conclusion is wrong: the system is time variant.

Regarding your question: for a system to be time invariant, it needs to be time invariant for any delay.

Even if there are specific delays when the system behaves as time-invariant, in general the system still is classified as time-variant.


Using all-pole filter to model the Room Impulse Response

It's not because of reverberation.

When you want to model the Frequency Response of the room, it's common to simplify your approximation by using either all-pole or all-zero models. You don't want to use the full zero-pole model.

To get some intuition:

zeros correspond to time delays and antiresonances
poles correspond to resonances of your Room Response

In practice all-zero models are not being used due to various reasons, such as:

required filter length is comparable to the IR length, and almost 40x the length of corresponding all-pole filter
filter will be valid only for specific distances and positions between the source and receiver (remember: time-delays).

That is why the all-pole model is used instead. As mentioned above, poles correspond to the resonances, i.e. standing waves, which are:

independent of the source location (quite intuitive)
independent of the receiver location (except of the nodes)

Additionally the required filter length is way less than in case of all-zero models. According to Mourjopoulos, for $RT_{60}\approx0.5 \mathrm{s}$, the required order is within the range of $50 < N < 500$.  The same author, concludes that all-pole filters are easier to manipulate than all-zero filters, due to their filter length. Morevover, author mentioned that all-pole filters are sufficient approximation than using raw impulse response data.

Here is some literature:

  Mourjopoulos J. - On the Variation and Invertibility of Room Impulse
  Response
  Functions

  Mourjopoulos J., Paraskevas M. A. - Pole and Zero Modelling of Room
  Transfer
  Functions


I have a pressure signal and want to do SPL analysis on it

Steps are quite simple:

Apply the window to your signal
Calculate the DFT of your signal
Scale the DFT properly (half of spectra and coherent gain)
Convert the result to decibel scale with a reference acoustic pressure of $2\cdot 10^{-5} \mathrm{Pa}$

I don't have MATLAB with me right now so here is some Python code.
Considering 1 kHz sinusoid with amplitude of 1 Pascal (94 dB SPL). Resulting plot:

Source code:

#!/usr/bin/env python

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':

    plt.close('all')

    # Generate 94 dB (Peak) SPL sinusoid
    fs = 48000
    duration = 1
    npts = int(fs*duration)
    t = np.linspace(0, 2*np.pi, npts)
    f = 500
    ref = 2.0e-5 # 20 uPa
    amp = 1      # 1 Pa -> 94 dB SPL
    s = amp * np.sin((f*duration) * t)

    # Window signal
    win = np.hamming(npts)
    signal = s * win

    sp = np.fft.fft(signal)
    freq = np.fft.fftfreq(npts, 1.0/fs)

    # Scale the magnitude of FFT by window energy and factor of 2,
    # because we are using half of FFT.
    # To obtain RMS values, divide by sqrt(2)
    sp_mag = np.abs(sp) * 2 / np.sum(win)

    # Shift both vectors to have DC at center
    freq   = np.fft.fftshift(freq)
    sp_mag = np.fft.fftshift(sp_mag)

    # Convert to decibel scale
    sp_db = 20 * np.log10( sp_mag/ref )

    plt.subplot(2,1,1)
    plt.plot(t, s)
    plt.xlabel('Time [s]')
    plt.ylabel('Acoustic pressure [Pa]')
    plt.grid(True)

    plt.subplot(2,1,2)
    plt.semilogx(freq, sp_db)
    plt.xlim( (0, fs/2)  )
    plt.ylim( (0, 100))
    plt.grid(True)
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('SPL [dB]')
    plt.show()


You are removing the scaling effect yourself when you perform the following operations:

Pd = abs(Yd/Ld)
P = abs(Y/L)

If you don't divide by the length, then you'll see that one has half the amplitude of the other.


Eigenfunction of LTI causal system Z-transform

Hint: $y[n] = H(e^{j (\pi / 2)}) x[n]$.

Since $x[n]$ is a complex exponential, it is an eigenfunction of a LTI system with eigenvalue given by the frequency response at the frequency of $x[n]$ (which in this case is $H(e^{j (\pi/2)})$. In less fancy words, if you put a complex exponential into a LTI system, you get the same complex exponential out, except scaled by the transfer function evaluated at the frequency of the complex exponential.


When you use machine learning algorithms on data sets, you use one part of the data (the training set) to train your algorithm (i.e., feature extraction). Once the training is completed, you'll need to evaluate the performance of the trained algorithm and you do this by applying it to new data, that is, the second part of your original data (the test data).

With this strategy you will be able to evaluate whether your trained algorithm indeed extracts valid features of your data (i.e., it works well on the training data and on the test data), or if it was "overlearned" or "overfitted" (i.e., it performs only well on the training data and bad on the test data).

For further information see wikipedia on training set, test set, and cross validation .


Why not use the same "standard" exponentials for both continuous and discrete time

This has to do with the way the Laplace transform and the $\mathcal{Z}$-transform are defined:
$$\mathcal{L}\big\{x(t)\big\}=\int_{-\infty}^{\infty}x(t)e^{-st}dt\tag{1}$$
$$\mathcal{Z}\big\{x[n]\big\}=\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{2}$$
Of course, the question remains why they are defined that way. The answer lies in the fact that in continuous time we usually have to deal with differentiation and integration, which in the Laplace transform domain correspond to multiplication with $s$ and $1/s$, respectively. Whereas in discrete time, we usually have to deal with difference equations formulated in terms of time delays and time advances, which in the $\mathcal{Z}$-transform domain correspond to multiplications by $z^{-1}$ and $z$, respectively.
Consequently, the definitions $(1)$ and $(2)$ lead to transfer functions of linear time-invariant (LTI) systems that are ratios of polynomials in $s$ for continuous time, and polynomials in $z$ for discrete time. In this way, the choices of $s$ and $z$ in $(1)$ and $(2)$ and the corresponding definitions of exponentials as $e^{st}$ and $z^n$ appear natural.


Realisation of a passive RC low pass filter

Note that your desired passive RC network is really just the cascade of multiple first-order RC filters, each of which have the transfer function $\frac{1}{1+RCs}$.  By choosing the R and C values to use in each stage, you can place pole locations where you want them in the $s$-plane.

Here's the problem: recall that the pole locations occur where the denominator of the transfer function is equal to zero. Therefore, for a particular RC stage, its corresponding pole location is:

$$
s = -\frac{1}{RC}
$$

Note that this value is purely real. This makes such a topology infeasible for implementing a prototype such as the Butterworth filter. Its poles are located in a semicircle in the left half of the (complex) $s$-plane. Since you can only place poles on the real axis using your RC stage topology, you cannot implement a Butterworth network. You end up having the same problems with other types like Chebyshev or elliptic structures.

The solution is to use additional elements like inductors or amplifiers; these provide more degrees of freedom in the resulting transfer function to place poles at other places in the $s$-plane. That's why all the references you can find use them!


How can I extract information about timbre from a .wav file. In which form can I study those features?

I think the time history of the wave file is enough to obtain all sorts of features. There are mathematical formulae to extract the features from a given time series.

Some important features are:

Time Domain Audio features

(a) Short-term Energy, (b) Zero Crossing Rate, (c) Entropy of Energy, (d) Basic Statistics (short term mean and standard deviations).

Frequency Domain Audio Features

(a) Spectral Centroid and Spread, (b) Spectral Entropy, (c) Spectral Flux, (d) Spectral Rolloff, (e) Mel-Frequency Cepstrum Coefficients, (f) Chroma Vector

More details regarding audio analysis can be found in this book
https://www.elsevier.com/books/introduction-to-audio-analysis/giannakopoulos/978-0-08-099388-1

If you have enough data, using these features you can build a machine learning model to correctly classify an unknown sound.


Signal Processing in Amateur Radio Astronomy

The signal processing operation that measures how "alike" two signals are is called cross-correlation.  For cross-correlation to work both signals must be sampled at the same rate.  If they aren't, you must resample one of them to the rate of the other.


Allowing for differences between Matlab and Octave ( I only have Gnu Octave to play with,) and correcting the errors noted in the comments, I get the following:
Needed fixes:

Values for Wp and Ws divided by half of sampling rate - the cheb1ord function requires frequencies be given relative to the Nyquist frequency.
Removed 's' (not allowed for Octave anyway.)

To make the code work for Octave, I had to make the following changes:

Remove duty cycle value from square
Remove 's' from cheb1ord and cheby1

Code for Octave (should run under Matlab as well:)

%Making the square pulses
Fs = 1000000;
nyq=Fs/2;
t = 0:1/Fs:50e-4;
f = 2000;

%square pulse with 50% duty cycle
xt = square(2*pi*f*t);

%Making the filter
Wp = [1000/nyq 3000/nyq];
Ws = [100/nyq 3900/nyq];
[n,Wn] = cheb1ord(Wp,Ws,3,20);
[b,a] = cheby1(n,3,Wn,'s');

%Filtering with the coefficients
xt_filtered = filter(b,a,xt);

%Plotting the filtered signal
subplot (211);
plot(t,xt);
subplot(212);
plot(t,xt_filtered);


Finding optimal decision threshold for binary comm system

Let P(0) ( P(1) ) is a probability of transmitting bit zero (one); P(e|0) ( P(e|1) ) is a probability of error when detecting bit zero (one).
The probability of erroneous detection is
$$
P(e) = P(e|0)P(0) + P(e|1)P(1) \tag {1}
$$

Let $V_0$ ( $V_1$ ) be a nominal signal voltage of bit zero ( one ) signal at the transmitter.
$$
P(e|0) = \int_T^{\infty}{{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(\nu-V_0)^2/2 {\sigma}^2\right) d\nu } \\
P(e|1) = \int_{-\infty}^T{{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(\nu-V_1)^2/2 {\sigma}^2\right) d\nu }
$$
where $T$ is a detection threshold (OP's $\mu$).

Differentiating $P(e)$ of eq.1 w.r.t. $T$, we arrive at
$$
-P(0){{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(T-V_0)^2/2 {\sigma}^2\right)} + P(1){{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(T-V_1)^2/2 {\sigma}^2\right)} \tag {2}
$$
To find an optimal threshold, we equate the expression in eq.2 to zero:
$$
P(0)\exp\left(-{\frac {(T-V_0)^2} {2 {\sigma}^2}}\right) = P(1)\exp\left(-{\frac {(T-V_1)^2} {2 {\sigma}^2}}\right) \\
T = {\frac {V_0+V_1} {2}} + {\sigma}^2{\ln{(\frac {P(1)} {P(0)})}}
$$
that is, exactly the OP's  result.
The OP's reference source may compute the optimal threshold for the system in which the probabilities of zero/one bit transmissions are equal. The variance ${\sigma} ^2$ is multipled by $\ln{(\frac {P(1)} {P(0)})}$, and it disappears from the expression for optimal threshold, when the zero/one probabilities in transmission are equal ($P(0)=P(1)=1/2$).


Algorithm/Code for finding surface shape of an image?

The 3x3 linear system appears to be that for photometric stereo, where the $i$-th row is $$l_i\cdot n = b_i,$$
with

$l_i$ is the 3D light source vector, e.g., $l_1$ the first row of the matrix, (0.2425, 0, -0.9701)
$n = (x_1, x_2, x_3)$ is the unknown surface normal vector, or "shape"
$b_i$ is the brightness observed when the surface was illuminated with $l_i$

Solving this 3x3 linear system gives the normal $n$ for a single surface point ("middle" according to "Example 1").  (You may wish to divide $n$ by its length $\sqrt{x_1^2+x_2^2+x_3^2}$ so that it is unit length.)

Usually, by surface shape in this context of photometric stereo one means the height function $z(x,y)$, where we associate a height to each pixel $(x,y)$.  To obtain this, we first need to solve a related 3x3 system at each pixel location $(x,y)$ (not just the "middle"), so that we have a surface normal $n(x,y)$ everywhere.  Then we solve for that function $z$ whose normal function is that $n$ we obtained via many 3x3 solutions above.  Finding $z$ requires we impose integrability on vector field $n$ (because modeling error and noise may mean there exists no $z$ having normals $n$).  However, I am digressing from your original "Example 1", which seeks a single surface normal.

See also: http://pages.cs.wisc.edu/~lizhang/courses/cs766-2008f/syllabus/10-09-shading/shading.pdf and https://en.wikipedia.org/wiki/Photometric_stereo

Note that the above assumes diffuse reflection (no surface glossiness), surface smoothness, and no shadows nor occlusion.  Light sources must be distant if we keep them constant over image.


Does the DFT calculate spectral components up to half the sampling frequency, $f_s/2$?

Depends a bit on the indexing convention, but in typically you would interpret the frequency interval as $[-f_s/2, f_s/2]$ and not as $[0,f_s]$. The DFT is periodic in N so you you have

$$X(N-1) = X(-1) $$

Keep in mind that that the sampling theorem requires the signal to be band limited to $f_s/2$ so assuming that you have actual independent information about frequencies higher than $f_s/2$ is misleading.

For real signals, it's conjugate symmetric anyway, i.e. $f_{-k} = f^*_k$ so there is only independent information on half of the spectrum anyway.


Why $x(t) = \cos(\frac{\pi}{3}t)$ and $y(t) = \sin(\frac{\pi}{3}t)$ are considered identical signals?

There is obviously an error in the problem formulation. From the given solution, the problem should read

  With $x(t) = \cos(\omega_x(t + \tau_x) + \theta_x)$ and
  $y(t) = \cos(\omega_y(t + \tau_y)+ \theta_y)$, determine for which of the
  following combinations $x(t)$ and $y(t)$ are identically equal for all
  $t$.

In this case the given solution $\omega_x=\omega_y$ and $\omega_x\tau_x+\theta_x=\omega_y\tau_y+\theta_y+2\pi k$ makes sense. The solution to the problem as it is given (with $y(t)=\sin(\omega_y(t + \tau_y)+ \theta_y)$), would be

$$\omega_x\tau_x+\theta_x=\omega_y\tau_y+\theta_y-\frac{\pi}{2}+2\pi k$$


Correlation between real and imaginary parts of a Fourier transform of zero mean Gaussian

Recall the following properties of the Fourier transform:

If $x(t)$ is an even function, then its Fourier transform $X(\omega)$ is purely real.
If $x(t)$ is an odd function, then its Fourier transform $X(\omega)$ is purely imaginary.

Thus, we can think of the real and imaginary parts of the Fourier transform of a zero-mean Gaussian random process $x(t)$ as the Fourier transforms of two separate inputs $x_e(t)$ and $x_o(t)$: the components of the process that have even and odd symmetry, respectively. We split the process into these components as follows:

$$
x_e(t) = \frac{x(t)+x(-t)}{2}
$$

$$
x_o(t) = \frac{x(t)-x(-t)}{2}
$$

Note that $x(t) = x_e(t) + x_o(t)$. Moving through this step by step,

Your observation was that $\text{Re}\{X(\omega)\}$ and $\text{Im}\{X(\omega)\}$ (the real and imaginary parts of the process's DFT) are uncorrelated.
Using the Fourier transform properties mentioned above, we can deduce that $\text{Re}\{X(\omega)\} = X_e(\omega)$ and $\text{Im}\{X(\omega)\} = X_o(\omega)$; the real and imaginary parts of $X(\omega)$ are none other than the Fourier transforms of $x_e(t)$ and $x_o(t)$, respectively.
Therefore, your observation is equivalent to saying that $X_e(\omega)$ and $X_o(\omega)$ are uncorrelated.
Since the Fourier transform is a one-to-one mapping between the time and frequency domains, I posit that the lack of correlation between $X_e(\omega)$ and $X_o(\omega)$ would imply a lack of correlation between $x_e(t)$ and $x_o(t)$ as well.

What is the correlation between $x_e(t)$ and $x_o(t)$? Simple:

$$
\begin{align}
\mathbb{E}(x_e(t)x_o(t)) &= \mathbb{E}\left(\left(\frac{x(t)+x(-t)}{2}\right)\left(\frac{x(t)-x(-t)}{2}\right)\right) \\
&= \mathbb{E}\left(\frac{1}{4}\left(x^2(t) - x^2(-t)\right)\right) \\
&= \frac{1}{4} \left(\mathbb{E}(x^2(t)) - \mathbb{E}(x^2(-t))\right) \\
&= \frac{1}{4} \left(\sigma^2 - \sigma^2\right) \\
&= 0
\end{align}
$$

As expected, the even and odd components are uncorrelated. So, to summarize, I would say the following:

The real and imaginary components of a Fourier transform correspond to the individual Fourier transforms even and odd components of the input function.
For a zero-mean Gaussian random process, these even and odd components are uncorrelated.
Therefore, their Fourier transforms (the real and imaginary components that you asked about) are also uncorrelated.

Edit: To address your followup:

If $x(t)$ is Gaussian, then its even and odd components $x_e(t)$ and $x_o(t)$ are as well, due to the property that any weighted sum of Gaussian random variables is also Gaussian.
If $x_e(t)$ and $x_o(t)$ are Gaussian random processes, then their Fourier transforms $X_e(\omega)$ and $X_o(\omega$) are as well. This follows from the same property as the previous statement; if you look at the transform, you're computing a weighted sum of a bunch of Gaussian random variables.
If $X_e(\omega)$ and $X_o(\omega)$ are Gaussian, and they are uncorrelated with one another (as described above), then they are also independent. This is a property of the Gaussian distribution.


How to model Tape Saturation (Audio DSP)?

Saturation can be modeled to first order as a symmetrical distortion of the signal, but it also has history dependence--hence the term hysteresis in describing magnetic response to a field that is applied.  Bias (an AC field applied to the tape at the same time as the signal to be recorded) linearizes the curve and reduces hysteresis effects by eradicating the memory of the medium.  If your objective is to simulate the sound of a tape recorder, you should consider hysteresis effects, and the fact that high frequency signals are self-biasing, so that they effectively saturate at a different (lower) level than low-frequency ones.  (This was the insight behind Dolby HX processing, which really made cassette tapes sound a lot better.)


Unless I'm missing something here it would just be the difference in frequency between the highest frequency sine wave and the lowest:

BW = Fmax - Fmin


How do I sum the discrete cosine transform cosine waves?

You used the DCT-II with a certain scaling convention:

$$X[k]=s[k]\sum_{n=0}^{N-1}x[n]\cos\left(\frac{\pi}{2N} (2n+1)k\right),\qquad k=0,1,\ldots,N-1\tag{1}$$

where $s[0]=1/\sqrt{N}$, and $s[k]=\sqrt{\frac{2}{N}}$ for $k=1,\ldots,N-1$. Its inverse is given by an appropriately scaled DTC-III:

$$x[n]=\sum_{k=0}^{N-1}s[k]X[k]\cos\left(\frac{\pi}{2N} (2n+1)k\right),\qquad n=0,1,\ldots,N-1\tag{2}$$

where $s[k]$ are the same scaling coefficients as used in $(1)$.

So, given the scaled DCT-II coefficients $X[k]$, the original signal can be written as

$$x[n]=\sqrt{\frac{2}{N}}\left[\frac{1}{\sqrt{2}}X[0]+X[1]\cos\left(\frac{\pi(2n+1)}{2N}\right)+X[2]\cos\left(\frac{\pi(2n+1)}{N}\right)+\ldots\right]$$


Does it make sense to have complex numbers representing real-world audio signals?

So it seems like real-world (discrete) audio signal might have complex values when being represented digitally,

No, you misunderstood that. The discrete audio time signal doesn't have non-real values. The Fourier transform can have such.

  but does this make sense?

It doesn't need to make sense. It's just math. I sometimes remind myself of that – it helps to occasionally get some distance from the very human desire to put sense into things.

However, assume any signal has Fourier transform $S(f)$, and that transform is purely real. Now you shift that signal in time by $\Delta t$. From knowing the Fourier transform properties, or by just applying the integral (continuous FT) or the sum (discrete FT), you instantly see that this means multiplication with a $e^{j2\pi \Delta t f}$. In the general case, this produces non-zero imaginary parts where there were none before.

In other words: even if you have a signal whose FT is purely real, simply by shifting it in time a bit, you can always get a non-zero imaginary part in the spectrum. That applies to all signals, be it audio, or I/Q baseband of a radio receiver.

  If yes, then how can we interpret this?

I think it makes sense for you to just continue to follow these lectures, I don't want to spoil too much, but:

only time-symmetrical (to be exact, time-hermitian, but since audio signals are real, that's identical to symmetry) signals have a real spectrum. If you see a non-real spectrum, you can simply tell the signal wasn't hermitian (symmetric).

Also, it tells you that if you want to deal with the spectra of signals, you always have to consider both, amplitude and phase, to represent a signal. And neither of these two aspects can be defined on the real part of the spectrum alone! The first is the root of the sum of squares of the real and the imaginary part (so the imaginary part is important), and the second is the arcustangens of the ratio of imaginary and real part. Phase "contains" the time-shift, measured in full signal periods, for harmonic signals, by the way, which is a class of signals very important for Audio.

  Isn’t the value of the input supposed to represent change in air pressure?

yes. But that the time-domain signal, not the frequency-domain signal.

  So if for an input signal x, say x[0] was 2 + 3j, what would this tell us about the pressure on the diaphragm at that point of time?

You're confusing a real time-domain signal with its Fourier transform. A single scalar value, like the momentary air pressure, is real. there's no $+3j$; there can't be, in time domain. However, in frequency domain, there might be an imaginary part.

Make sure that you understand this:

A signal is just a changing value over time, or some other axis.
The most intuitive representation hence is one along that axis, for example, an audio signal might be a change in pressure over time.

However, there's other ways of expressing the same signal. The Fourier Transform of a signal describes the same signal, but is a different function than the time signal. We often speak of "frequency domain" (after FT) and "time domain".

I don't know your background, but early on in my studies, we were taught of vector spaces, bases and base transforms.

"Time steps" is just one basis with which you can describe a signal; "Fourier Coefficients" is another one. The same signal has different coefficients (ie. representations as vectors of numbers) regarding these two bases. And the discrete FT is really, really, really nothing but a plain, old, boring base transform matrix with a few nice properties!


How to calculate the required bandwidth from a spectrum diagram? [FM]

There are many definitions of bandwidth; see for example section 2.2 in this book (PDF linked at the end of the page). In your case, it is easy to use the "absolute bandwidth", which is defined as $B=f_{max}-f_{min}$, where $f_{max}$ is the largest frequency in the spectrum, and $f_{min}$ is the smallest.

According to this definition, you'd have $B=(100.028-99.972)\mathrm{\,MHz}=56\mathrm{\,kHz}$.


Compressing plotted data with variable sample rate?

Here is what I did. The basic idea is to mimic what you would see if you were plotting every single sample, just do it much faster

Calculate roughly how many data points will fall in one horizontal pixel. There is no point in plotting more points than the are pixels on the screen.
Divvy up the data in blocks so that you roughly get one block per pixel
Calculate the min and the max for each block
For each block just add the max and the min in the correct order to your plot data

This way you get a plot array that has about twice the number of data points than there are pixels which is about the highest number that makes sense. You have to recalculate whenever you zoom, change the data range or re-size the window, but that tends to be fairly quick. If you don't want to handle window resizing, you can simply just calculate for the worst case (full screen).

Here is a code snippet that shows that in Matlab (it omits the ordering for brevity). The hi-res and the low-res plot look almost the same.

%% reduced sample plot example
n = 200000; % lots of samples
fs = 44100;
% noise with exponential decay
x = randn(n,1).*exp(-(0:n-1)'/fs/3);
% block approach
npix = 2000; % emulate an HD display
x1 = reshape(x,n/npix,npix); % chop into blocks
x2 = [min(x1); max(x1)];  % cal min and max
x2 = x2(:); % interleave
figure(1); plot(x); % plot with 200000 data points
figure(2); plot(x2); % plot with 2000 data points


I thought my network was experiencing outages. I collected bandwidth (Mbps) measurements every 30 min.
. . .
What would be some ways to analyze the signal to see the time scale of outages?

An "outage" is a random event. Therefore, looking for regularity in outages might be hinting at something more widely already known.
The answers to the questions in the comments are very much dependent to the objectives. But I will try to make this a bit more constructive.
For example:
The CDF presented is over 1400 instances at 30 minute intervals and it therefore gives you the probability of sustained download at a particular speed over 1400*30 minutes or about 30 days (If the horizontal axis depicts the index of the measurement). If you were asking if your line would sustain a download speed over the next 15 minutes, the answer would be "I don't know", because you have not measured anything more frequently than 15 minutes.
On their own, wavelets and discrete fourier analysis cannot tell you much. They can tell you that a particular segment of the signal suddenly requires many more components to describe it and this usually happens at discontinuities but this helps in recognising a "spike" which you can already do here by thresholding (for example).
Rather, if you suspect that there is some regularity in the outages and you want to estimate it, then what you can do is some form of time series decomposition.
You can run autocorrelation for example and this will already give you a handle on periodicity. Or a full blown discrete fourier transform (or even wavelet transform) to determine periodic components more accurately.
By determining the Trend and Seasonality components of your signal, you can then create "filters" to discover areas where it veers off its usual course. (And then of course, having discovered the timestamps of those events, you can analyse (or try to predict) the time that passess between them and see if there is anything useful in there).
But notice here, we always need a time reference. For example, when you assess Trend, a trend over what time scale? This will then determine how often you need to observe the signal.
Hope this helps.


Initial Rest Condition for LCCDE causal LTI systems

If you are given an input $x(t)$ with $x(t)=0$ for $t<t_0$, and you specify an initial condition $y(t_1)=0$ for $t_1>t_0$, then the resulting system is generally non-causal, because we already know the system's response at $t_1>t_0$, regardless of the input signal in the interval $[t_0,t_1]$.

For a system to be linear (in the sense used by Oppenheim) and causal, we need to specify the initial condition at the point $t_0^-$, i.e., we require that the left-sided limit $y(t_0^-)$ satisfies $y(t_0^-)=0$.

For the given example, if an impulse is applied at $t=0$, we require $y(0^-)=0$ for the system to be linear and causal. Note that there is no requirement on the right-sided limit $y(0^+)$.

Clearly, the response $y(t)=e^{-at}u(t)$ satisfies the initial condition $y(0^-)=0$.


Unable to understand how the paper simplifies the covariance matrix - Kalman filter

In order to derive equation (20) you can use the following steps:

From substitution of  equation (16) into equation (18)

$$  P_n = P_{n|n-1} - K_nP^t_{x_ny_n}$$

Now plugging right side equation of (19) into the last equation

$$  P_n = P_{n|n-1} - K_n(P_{y_ny_n} - R_n)=$$
$$ =P_{n|n-1} - K_nP_{y_ny_n} + K_nR_n $$

Now Pluging (16) again we know that
$$ K_nP_{y_ny_n} = P_{x_ny_n}$$
putting this back we get

$$ P_n=P_{n|n-1} - P_{x_ny_n} + K_nR_n $$

Using left side equation of (19) we get the desired result
$$P_n = K_nR_n$$


The central notion is that  the considered algebras are generated by a shift operator $\mathcal{S}$. Every element (filter) of the algebra should be (by definition) generated by a (finite) number of algebraic operations (addition, product, scalar multiplication) involving the generator and its inverse. This is explained in more details in Section II-A, around Equation (30), and footnote 19. So if the shift commutes with any filter, since any filter is made of a combination of shifts and inverse shifts, then any filter commutes with any filter.
To me, this is really related to the subset of shift generated algebras they consider (and not on the module of signals).
For instance, with standard discrete filters this unit element would be the Kronecker $\delta_0$ (the "unit" algebra element relative to the multiplicative $\cdot$ operation of the algebra, here the convolution). A (finite support) filter would be written as sums and multiplies of shifts, with scalars corresponding to filter taps: $$h_0S^0(\delta_0)+h_1S(\delta_0)+h_2S^2(\delta_0)+\cdots$$
and commute with basic shifts.


What is the use of the Welch lower bound when k > 1?

Given any collection of $m$ unit-norm (or unit-Euclidean-length)
vectors in a real or complex $n$-dimensional vector space $\mathbb R^n$ or
$\mathbb C^n$, let $c_{\max}$ be the maximum magnitude of the
inner product between two vectors in the set. How small can $c_{\max}$ be?
One answer to this question is given by the
Welch bounds which are a collection of lower bounds on the minimum possible
value of $c_{\max}$. These bounds apply to all collections of $m$ unit-norm
vectors in $n$-dimensional space,
and they even hold for the case when the collection is a multiset (and thus
two vectors are identical, and so $c_{\max}$ necessarily equals $1$).

The Welch bounds are of the form
$$c_{\max}^2 \geq \sqrt[k]{C_k(m,n)}, ~k = 1, 2, 3, \ldots \tag{1}$$
or more correctly
$$c_{\max}^2 \geq \sqrt[k]{\max\{C_k(m,n), 0\}},~ k = 1, 2, 3, \ldots \tag{2}$$
because the function $C_k(m,n)$ can have negative value. Note for
example, that for $k = 1$, $(1)$ becomes
$$c_{\max}^2 \geq \frac{m-n}{n(m-1)}$$
where the right side is negative when $m < n$ and $0$ when $m = n$.
This is as it should be because we know that in a $n$-dimensional
space, we can always find $m \leq n$ orthonormal vectors and so
$c_\max^2 = 0$ for this set.
Thus, $(2)$ which becomes
$$c_{\max}^2 \geq \max\left\{\frac{m-n}{n(m-1)}, 0\right\}\tag{3}$$
is a much better way of expressing the bound.
When $m \gg n$ and so $(m-n)/(m-1)$ is approximately $1$, then
(as Welch pointed out as also have you in your question)
the right side of $(3)$ approaches $n^{-1}$ as $m \to \infty$ and so
we have that
$$c_\max  \geq \frac{1}{\sqrt{n}}\tag{4}$$
as a limiting form of the lower bound.

Turning to the question asked about Welch bounds for $k \geq 2$,
these are useful when $m \gg n$ and give larger lower bounds than
the right side of $(4)$. If I remember correctly, the Welch bound
for $k = 2$ gives something like
$$c_\max \geq \frac{2^{1/4}}{\sqrt{n}}$$
when $m \geq \frac{n^3+n^2}{2}$ which is an improvement on $(4)$.

For a discussion of some of these points (as applied to
crosscorrelations and autocorrelations instead of inner products),
I refer you
to my paper "Bounds on crosscorrelation and autocorrelation of sequences"
in IEEE Transactions on Information Theory, November 1979 which is
behind IEEE's paywall unless you have institutional access through your
employer or have paid personally for IEEExplore access.


Auto-correlation function, an inverse problem

Let's look at the case $x[n] \in \mathbb{R}$, where $x[n]$ is real.

Autocorrelation is basically convolution of the signal with it's time inverse. This can be easily expressed in the frequency domain.

$$ \mathscr{F}\Big\{ r_{xx}[n] \Big\} = \mathscr{F}\Big\{ x[n]  \Big\} \cdot  \mathscr{F}\Big\{ x[-n] \Big\} $$

$$R_{xx}(\omega) = X(\omega)\cdot X^*(\omega) = \Big| X(\omega) \Big|^2
$$

So it's easy to see that the Fourier Transform of the auto correlation is simply the magnitude squared of the Fourier Transform of the input signal. That's sometimes referred to as the Power Spectrum.

It's also easy to see that information gets lost in the process. There are $N$ unique values going in but because of the symmetry properties of the auto correlation there are only $\frac{N}{2}$ unique (independent) values coming out. Looking in the frequency domain, we can see that the phase is lost.

  If yes, is it unique? and what is the method to find $x[n]$?

No, it's not unique

  Can we find the corresponding $x[n]$ ?

There is an infinite number of $x[n]$.

Take the Fourier Transform of the autocorrelation
Take the square root
Add an aribtiraty (but odd-symmetric) phase function
Do an inverse Fourier Transform

Any signal derived this way will have the same original auto correlation function.

  if no, what other constraint properties should we add to $C[n]$, in order to make it yes?

You can't make it a yes, since it's not unique. No matter what auto correlation you choose, there will be infinite $x[n]$ that will have it as an autocorrelation.


Does Zero Padding Distort the Spectrum of a Signal?

DFT vs DTFT, Fourier Transform
The problem appears rooted in viewing DFT as a 'special case' of the continuous Fourier Transform, and of its input as some signal with legitimate frequency contents. That's a fallacy.
The DFT is a standalone mathematical transformation which does not condition itself upon anything but the input being finitely long, and having finite values. It's an information encoding algorithm. It operates on any and all finite information - be it signals, Shakesperean sonnets, or a cup of water - as long as represented numerically.
DFT can be shown to be equivalent to solving $N$ simultaneous equations defining an $N$ long input; it's a complete, lossless transform.

What the spectrum "should be" is not what we get

Should be: amounts of sinusoidal frequencies comprising the input data.$^1$
Is: amounts of sinusoidal frequencies that sum to the input data.$^1$

See here. The input can be a pulse, or a single nonzero value and the rest zeros; there's nothing of periodic nature about either of these, nor do they need to stem from a physical source which is sinusoidal. The only truth is that the coefficients we find give us the sinusoids that sum to the input; with another transform, it can be triangle waves.
1: for a real-valued input. For complex, something similar may hold, but it's still a sum of bases (discussion here).

DFT coefficients are not same as DTFT or FT coefficients
Both the forward and inverse transforms to acquire respective coefficients imply very different things for what the coefficients "mean". Under certain assumptions they can be treated equivalently, but not fundamentally.
This is clear from the inverse-transform perspective, or reconstruction; assume for now a rectangular-windowed DTFT:

DFT: sum of $N$ bases, each $N$-long, finite in value.
DFT coefficient: describes one such basis. $N$ coefficients in total.
DTFT: infinitesimal sum (integral) of infinitely many bases, each $N$-long, finite in value.
DTFT coefficient: describes one such basis. infinitely many coefficients in total.

From forward-transform:

DFT coefficient: obtained by adding (discretized complex phasor at integer frequency) * (input), pointwise products.
DTFT coefficient: obtained the same way. Except not really; see below.

Result? If you use one transform's coefficients in another transform, you get nonsense.

Generator vs Data
A generator is a complete descriptor - example: function, continuous. For describing a 1D signal, the signal is infinitely resolved; we can determine every possible intermediate value. From a generator we can get data, but not vice versa; even a truly continuous segment has infinite interpolations.
The DTFT (no windows) works with generators, not data. Its coefficients are described by a function (continuous), which are also determined by an operation between functions (one describing an infinitely long but discrete input, another the complex phasor at a continuous $\omega$).
The rectangular-windowed DTFT, on the other hand, works with both, but only "with permission" from generators - that is, the resulting spectrum happens to be the same whether we compute for $n=0$ to $N-1$ or $-\infty$ to $\infty$.
In contrast, the DFT works with data. No generators. But we can use generators if we wish - it's now the reverse case, generators "with permission" from data.

Zero-padding: DFT coeffs treated as DTFT coeffs -- all code.
Before the meat, last piece of background. Take 1Hz over two seconds; what's the DFT (which freq is nonzero)? 2. Now take 0.1Hz over 10 seconds; what's the DFT? 1. The DFT correlates against basis functions spanning the length of the input signal, with no regard to "actual time". $k=1$ corresponds to 1 cycle over the length of the signal.
We can use this to interpret "treat DFT coeffs as DTFT coeffs". Suppose
$$s[n] = \cos(2 \pi n / N),\ N = 128,\ n=[0, ..., N - 1].$$

Pre-pad, the spectrum is just what we expect: nonzero everywhere but at $k=1$ and $k=128$. Post-pad - a madhouse. What's the deal?
Recall how "0.1Hz over 10 seconds" $\rightarrow k=1$? By the same token, the original $f=1$ is now "seen" as $k=2$ - i.e., it correlates with the $k=2$ basis. Then what does $k=1$ "see"? Let's see:

(Only real component of basis shown for clarity; can add the imaginary component for completeness)
Put differently, the "normalized" frequency of the DFT is divided by two as a result of doubling the original signal's length with the padding, so $k_p=2 \leftrightarrow k = 1$, and $k_p=1 \leftrightarrow k = 0.5$. Or, integer frequencies with a padded signal's spectrum correspond to fractional frequencies in the unpadded's spectrum.
Is there more evidence? Yes; first, see this. If above is true, then we
can correlate the original signal with a $k=0.5$ basis directly, and yield the same coefficient as with $k_p=1$. And indeed:

Why one window for both padded and unpadded? ... because they look identical - check the code. Mean absolute difference: 6e-14 (within float).

In what sense do we "sample the DTFT"? What's the role of zero-padding?
Former should be clear now; the DTFT evaluates itself for every $\omega$. For any given $\omega$, this  exactly the same as taking DFT with an equivalent $k$ basis. In other words, DTFT correlates with neighboring fractional frequencies, just like DFT correlates with integer frequencies. Simpler, "DTFT is DFT with $k$ spanning continuously" (but not quite; see caveats in above sections).
So why does zero-padding "result in" the DFT sampling DTFT's spectrum? The answer is closer if we reformulate the question: Why does zero-padding result in the DFT correlating with fractional frequency bases?
Yes, they're indeed one and the same question, which was answered in the previous section. A DFT with $2N$ coefficients is equivalent to a DFT with $N$ coefficients at integer frequency bases, and another DFT with $N$ coefficients at fractional frequency bases$^1$.
But why zero-padding? Why not one-padding, or pi-padding? Because it's the only way for an $M>N$-DFT with integer frequencies to be equivalent to an $M$-DFT with fractional frequencies$^1$. Refer to the blue/orange figure; note how the sum of products of the orange with the blue for $k_p=2$ is the same as for $k=1$ with the unpadded blue (not shown; just imagine). If the zeros were anything else, the values might have been proportional, or somehow similar, but not equal. In technical terms, zero is orthogonal to all basis functions over any interval.
1: if $M=3N$, then the integer-DFT will correlate $k=[0, 1, ..., 3N-1]$ bases with the padded signal, and the "mixed" DFT will correlate $k=[0, 0.33, 0.66, 1, ..., N-1, N - .66, N - .33]$ with the original signal. Note that if $M$ is not an integer multiple of $N$, then the padded's integer-DFT may not have any integer equivalents with the unpadded's DFT.

Does zero-padding distort the spectrum?
The grand-prize. The answer to which is... yes and no:

No: zero padding is equivalent to correlating with neighboring, non-integer frequencies; there is no "distortion" here, the coefficients describe exactly what they intend to.
Yes: if "the spectrum" is supposed to describe the frequencies in the input signal (e.g. $f=1$), then correlating with bases neighboring $1$ yields non-zero coefficients, implying non-zero frequencies around $1$, whereas the signal has no such frequencies. Further, only integer-ratio frequencies are completely decorrelated (e.g. 1 and 2, 1.5 and 3), so we'll get many high and low non-zero frequencies, whereas the input has only one frequency.

So what's the deal? The deal is - DFT does NOT describe input's spectrum. This is the entire point of this answer, which was so warmly-received because of a technicality I was quick to retract with a NOTE: on top. Ignore the downvotes, they mean nothing without explanation.
Since it's not the spectrum we expect, it's meaningless to speak of a "distortion" as if it ever gave the actual spectrum. The input could be $f=1.1$, and then the unpadded DFT would be "distortion", and the only true spectrum would stem from incrementing the bases by $1.1$.
This also explains why we can't do much better than the DFT; we don't know what the "actual" frequency of a signal is without its generator (function). The DFT is simply our "best guess", one which must be interpreted carefully.

Possible to pad without distorting?
Yes, but in a specific sense: if the number of zeros added is an integer-multiple of the signal's length, then the DFT will contain the unpadded's spectrum. For example, pad by $2N$; then, the DFT, from original's point of view, is doing $k=[0,.33,.66,1,...,N-1,N-.66,N-.33]$. Notice how $[0, 1, ..., N-1]$ are all captured.
In other words, we can pad without losing the original spectrum.

Re: other answers

The fact that the input is discrete or finite has virtually nothing to do with this - which this upvoted answer contradicts. Though, a little more work is due to show this; short version is, imagine the "winding machine" winding forever, and why that yields an impulse, and how that relates to the finite spike in DFT. Whatever the answer, it's way besides the point.


Let's say you want to filter a signal $x[n]$ through a Gaussian filter with impulse response $h_g[n]$ and a Moving Average (the "sqwave") filter with impulse response $h_s[n]$.  Then the resulting operation is

$$y[n] = h_s[n] * (h_g[n] * x[n]) = (h_s[n] * h_g[n]) * x[n]$$

So the a priori convolution of the two filter tap sets into one filter taps set, just reduces the two filtering operations into one filtering operation.

The Moving Average filter is just a low pass filter with a poor transition band and poor sidelobe supression.  (The one used here isn't even normalized to unity gain.)

I can't definitely answer "why?", except to say that it looks like someone wanted to add a cheap low pass filter as a sort of channel filter.  That moving average filter has been there since 2007, when GNURadio's git history starts.  The old CVS repositories aren't on-line anymore, so the original rationale is probably lost to history.


Two different answers while doing Inverse Z-transform

Q1: You should ask yourself if

$$\frac{z}{3z^2 - 4z + 1}\stackrel{?}{=}\frac{\frac 3 2}{z-1} - \frac{\frac 1 2 }{z-\frac 1 3}$$

really holds. You'll find out that you forgot to scale correctly.

Q2: There is no reason to consider $X(z)/z$ instead of $X(z)$ in this case.

And concerning your question about the anti-causal sequence, a multiplication by $z^{-1}$ always corresponds to replacing $n$ by $n-1$, no matter if a sequence is right-sided or left-sided.


Why does FFT not have an effect on my smoothed signal?


Keypoint orientation with the oFAST from Rublee's ORB

The point descriptor usually requires samples at some specific locations from around the detected point.

If you have point location refined and orientation assigned, you just shift and rotate the locations before computing point descriptor.

Some simple descriptors require samples at discrete locations (e.g. 15x15 pixel image patch) and applying shift+rotation will make the sample locations non-discrete (the patch is arbitrarily rotated and shifted) and you have to use bilinear or bicubic interpolation to take the pixel values.

Here is a sample of detected points - notice the corresponding patches are rotated:

To get e.g. 64x64 matrix of values for each patch, one has to sample points rotated accordingly.


Wiener filter optimal coefficients and the input signal

For conveniece I'm going to drop the time indexes. Let $e=d-\mathbf{a}^T\mathbf{x}$, where $d$ is the desired signal and $\mathbf{a}$ weights the elements of the observations/input vector $\mathbf{x}$. The goal is to minimize expected squared error with respect to the vector $\mathbf{a}$.

$$\min J=\min E[e^2],$$

where $E$ is the expectation operator.
To find the minimum we take the derivative of the objective function w.r.t. $\mathbf{a}$ and set it equal to zero. Therefore
$$
\mathbf{g}=\frac{\partial E[e^2]}{\partial \mathbf{a}}=E[2e\frac{\partial e}{\partial \mathbf{a}}]=0,$$

but from the defintion of $e$ earlier, we have
$$ \frac{\partial e}{\partial \mathbf{a}}= -\mathbf{x}$$
Therefore
$$ E[2e\frac{\partial e}{\partial \mathbf{a}}]= -E[2e\mathbf{x}]=0. $$ The factor 2 can be removed and thus $E[e\mathbf{x}]=0$, which means the error $e$ is orthogonal to the input signal $\mathbf{x}$.


How to create a sine wave that transitions through frequencies and returns to phase

For the first continuous $f_0=300\mbox{Hz}$ signal, after $2N$ samples (for a duration of $2N/f_s$ at a sampling rate $f_s$) the phase will have changed by:
$$
\begin{align}
 \Delta \phi_1
   &= 2\pi \frac{f_0}{f_s} (2N-1)
\end{align}
$$

On the other hand, as the second signal ramps up linearly (in small quantized steps as per your phasor implementation) from $f_0=300\mbox{Hz}$ to $f_1=305\mbox{Hz}$ then down from 305Hz to 300Hz over the same $2N$ samples, the phase will have changed by:
$$
\begin{align}
 \Delta \phi_2
   &= \sum_{n=1}^{N} 2\pi \left(\frac{f_0}{f_s} + \left(\frac{f_1}{f_s}-\frac{f_0}{f_s}\right)\frac{n}{N}\right)
      + \sum_{n=1}^{N-1} 2\pi \left(\frac{f_1}{f_s} + \left(\frac{f_0}{f_s}-\frac{f_1}{f_s}\right)\frac{n}{N}\right) \\
%   &= 2\pi \frac{f_1}{f_s}
%      + 2\pi \sum_{n=1}^{N-1} \left(\frac{f_0}{f_s} + \left(\frac{f_1}{f_s}-\frac{f_0}{f_s}\right)\frac{n}{N} + \frac{f_1}{f_s} + \left(\frac{f_0}{f_s}-\frac{f_1}{f_s}\right)\frac{n}{N}\right) \\
%   &= 2\pi \frac{f_1}{f_s}
%      + 2\pi \sum_{n=1}^{N-1} \left(\frac{f_0}{f_s} + \frac{f_1}{f_s}\right) \\
   &= 2\pi \frac{f_1}{f_s} + 2\pi \left(\frac{f_0}{f_s}+\frac{f_1}{f_s}\right) \left(N - 1\right)
\end{align}
$$

For those two signals to be in phase after $2N$ samples, their phase difference must be a multiple of $2\pi$. Thus the constraint is that
$$
\begin{align}
  \Delta \phi_2
    &= \Delta \phi_1 + 2\pi k \\
  2\pi \frac{f_1}{f_s} + 2\pi \left(\frac{f_0}{f_s}+\frac{f_1}{f_s}\right) \left(N - 1\right)
    &= 2\pi \frac{f_0}{f_s} \left(2N-1\right) + 2\pi k
\end{align}
$$
for some arbitrary integer $k$, or equivalently:
$$
\begin{align}
  N &= k \frac{f_s}{f_1-f_0}
\end{align}
$$
Which, provided that $kf_s/(f_1-f_0)$ is an integer greater than 0, yields a signal duration of
\begin{align}
  \frac{2k}{f_1-f_0}
\end{align}

In your specific case for a frequency excursion of $5\mbox{Hz}$, the signals will return in phase for modulation durations of 0.4, 0.8, 1.2, ... seconds.

Note that if you need to sustain the signal for some $M$ samples at $f_1$ (instead of immediately going back down after the rise), then the constraint becomes:
$$
\begin{align}
  N+M &= k \frac{f_s}{f_1-f_0}
\end{align}
$$
In other words the duration of the rise and plateau together would, in your case, need to be a multiple of 0.2 seconds. So for example if you want to have a plateau of 1 second at 305Hz, you might choose a rise + plateau of 1.2 seconds, giving 0.2 seconds for the ramp up and another 0.2 seconds for the ramp down (for a total signal duration of 1.4 seconds)


Calculate phase shift between two signals

Do this procedure on both signals, and subtract the results to get what you want (it would be a nearly constant sequence whose value is phase shift in radians)

Perform FFT on the signal to obtain its spectrum.
Find the dominant harmonic of the signal : $f_c$ (where the peak of FFT occurs).
Filter the signal with a narrow-band linear phase filter (center = $f_c$, bandwidth as less as possible).
Again, perform FFT on the signal and extract its phase.

Note : It is assumed that the signal has a dominant frequency otherwise a phase shift is pointless


How to invert FFT these magnitude equations to get FIR filters?

As @RBJ commented the first one is a Gaussian filter. The Fourier Transform of a Gaussian is a Gaussian so this one can be easily created by sampling the impulse response given by the inverse Fourier Transform.

For case (1), an FIR filter with a magnitude response given by:

$$|H(f)| = e^{-cf^2}\tag{1}\label{1}$$

Given the Fourier Transform of a Gaussian:

$$ e^{-\pi t^2} \leftrightarrow e^{-\pi f^2}\tag{2}\label{2}$$

And using the Fourier Transform scaling property:

$$\frac{1}{|a|}x(t/a) \leftrightarrow X(af) \tag{3}\label{3}$$

We can rewrite $\ref{1}$ as:

$$X(af) = |H(af)| = e^{-\pi (af)^2}\tag{4}\label{4}$$

with $a = \sqrt{c/\pi}$

If we restrict $c$ to be a positive real number, then from $\ref{3}$ we get:

$$\mathscr{F}^{-1}\{ e^{\pi(af)^2}\}=\frac{1}{|a|}e^{-\pi(t/a)^2} = \sqrt{\frac{\pi}{c}}\left( e^{\frac{(\pi t)^2}{c}}\right)\tag{5}\label{5}$$

This is demonstrated with an increasing number of coefficients showing the excellent match that is achieved.

The MATLAB/Octave code used to create the filter coefficients was:

fs = 1;   # sampling rate
c = 100;
N = 27;   # number of taps, must be odd
t = [-(N-1)/2: (N-1)/2]/fs;   # time axis
coeff = sqrt(pi/c)*exp(-pi^2*t.^2/c); # sampled imp resp

Case (2) appears to be a high-pass filter which could be solved using the approach above by sampling the inverse Fourier Transform for $e^{-a/f}$ and $e^{-a/\sqrt{f}}$.

I couldn't easily proceed from this point and posted on the mathematics site for a possible solution:

https://math.stackexchange.com/questions/3666703/inverse-fourier-transform-for-e-a-f-and-e-a-sqrtf

A very simple but sub-optimum approach is to derive the filter coefficients from the inverse DFT. This is the frequency sampling approach to FIR design and provides an exact match at the frequency sampling locations, but results in significantly more ripple than optimized algorithms such as least-square and Parks-McLellan. The ripple can be significantly reduced by windowing the resulting coefficients, but for this application any windowing will significantly reduce the null at DC. However, in the spirit of "needing something that will work one way or another" the Frequency Sampling approach can be considered if a larger number of taps is feasible. Regardless of approach, the number of taps for an FIR filter is proportional to the inverse of the width of the transition band of the filter, so if the parameters are such that the filter is a tight DC-notch filter, a large number of taps with any traditional FIR approach will be necessary (an IIR filter may make a lot more sense, or a modification of Richard Lyons clever implementation of a linear=phase DC notch filter with CIC structure: https://www.dsprelated.com/showarticle/58.php ).

The example plots below show the result for a 275 and 501 tap FIR filter using the Frequency Sampling approach with the following example parameters as given in the code below to implement an FIR for Case (2):

fs = 1;
a = 500; b = 50; d = 50000; g = 100;
N = 275;
f = linspace(-fs/2, fs/2, N);
magC2 = e.^(-(a+b.*(sqrt(g*abs(f))))./(d*abs(f)));
coeffC2 = fftshift(ifft(ifftshift(magC2)));

The above plot appears to show an excellent match in both cases, however zooming in on the stop band shows that the frequency sampling approach only maintains a match (with the tighter parameters used) down to around -25 dB with 275 taps, while it maintains a match down to around -65 dB with 501 taps. Pursuing either an approach of sampling the computed impulse response from the inverse Fourier Transform, or a least squares algorithm to match the magnitude response will allow for a solution with significantly fewer taps. Also if the notch needed is actually not as tight as I have demonstrated here, the number of taps required will be significantly reduced as well.

Also observe how we can estimate the number of taps needed for any approach from the ideal in the plot above: if we want/need to maintain high accuracy out to -80 dB rejection, we observe for this case with the particular parameters chosen that the steepest transition the slope is approximately 65 dB over a frequency of 0.005 samples/cycle. From fred harris' rule of thumb this would be:

$$N \approx \frac{A}{22}\frac{1}{\Delta f} = (65/22)(1/0.005) = 590 $$

We can also see how increasingly challenging it will be to maintain a tight match at deeper rejections given the slope is continuously increasing.


How can I realize a complex communication system in a real time application?

The output of qammod in your program is a sequence of amplitudes. If $a_k$ is the real part of the sequence, and $b_k$ the imaginary part, then the QAM baseband signal is:
\begin{align*}
s_{BB}(t)&=\sum_k a_k p(t-kT_p) + j\sum b_k p(t-kT_p) \\
&= A(t) + jB(t).
\end{align*} In this equation, $T_p$ is the pulse (baud) rate, and $p(t)$ is the pulse shape.

In order to convert this to a real, physical passband signal with carrier frequency $f_c$, you upconvert the baseband signal to frequency $f_c$: $$s_{PB}(t)=\text{Re}(s_{BB}(t)e^{j2\pi f_ct}).$$ While the previous equation is mathematically nice, in practice this is implemented as: $$s_{PB}(t)=A(t)\cos(2\pi f_ct)+B(t)\sin(2\pi f_ct)$$ with possible scale factors to keep the signal energy constant.

Once you have $s_{PB}(t)$, all that remains is amplify it, possibly bandpass filter it to eliminate out-of-band emissions, and connect an antenna.

Note that Matlab works in the discrete-time domain and the signals I've described are continuous time. At some point in your system, you will need to do a digital to analog conversion. Exactly where that will happen depends on many factors. I recommend that you take a look at some common RF architectures, such as heterodyne, homodyne, zero-IF and so on. They can give you ideas on how to implement the analog part of your system. One good reference is "RF System Design of Transceivers for Wireless Communications" by Q. Gu. Another very good book is Steven Tretter's "Communication System Design Using DSP Algorithms: With Laboratory Experiments".


How to find average and power of a signal

No solutions to homework problems, but here are two hints:

what is $e^{j2\pi nk}$ for integer $k$ and $n$?
do the values of $x[n]$ depend on $n$?


DAFX is a pretty good one.

Harmony Central used to have a surprisingly good knowledgeable base for algorithms but that seem to have disappeared since Guitar Center took them over.

A lot of source code can be found at http://www.musicdsp.org/, but it's fairly varied in quality.


Possible to Derive Receive Signal Power from RSSI?

I also read this from a response to a USRP user's question about RSSI measurements:

    [The] Received Signal Strength [Indicator is] always relative to some signal model, incorporating considered bandwidth, assumptions on the modulation scheme, duration of transmission, generally: It's a estimation of received signal strength based on some property of the received signal.

Yes, I wrote that, and I still think it applies, so:

  What is the difference between RSSI measurements and the dBm values that I record from the spectrum analyzer?

Well, the RSSI is calculated based on some signal-model based estimator, whereas the dBm value on your spectrum analyzer is calculated with another metric – namely, the power measured going through a filter with the filter bandwidth you configured at your spectrum analyzer.

These two might, or might not be related, taking into account that even a signal with very low power spectral density values might be sufficiently good for something that has e.g. matched filters to extract the signal from noise, or any other form of processing gain.

So, your only option is to understand both estimators, with the Spectrum Analyzer estimator being relatively simple, and the flex estimator being free and open source (so you can just read the code ), and relate the two; so basically, you'll have to find, given a (class of) signal(s) $s(t)$, the functions $f_{SA}$, $f_{flex}$ and finally $g$ from

$$\begin{align}
P_\text{spectrum analyzer}\left(s(t)\right) &= f_{SA}\left(s(t)\right)\\
P_\text{flex frame sync RSSI}\left(s(t)\right) &= f_{flex}\left(s(t)\right)\\
\leadsto
P_\text{flex}'\left(P_\text{spectrum analyzer}\right) &= g(f_{SA}\left(s(t)\right))\text{ .}
\end{align}$$

I've just skimmed the liquiddsp code for the RSSI estimation, and the interesting part seems to happen in qdetector_cccf_execute_align.

However, this can never work without making a lot of assumptions. For example, if your signal happens to have a very narrow, but very strong single tone interferer, it will probably not break flex reception, and if that RSSI estimator is any good, not change RSSI very much, but it will drastically change the dBm value displayed by your spectrum analyzer. So if the question was

  Can I find an estimator that always gives me the same reading as the liquidDSP flex RSSI based on the dBm value displayed by a spectrum analyzer?

then the answer is a

  No. See above proof by counterexample (single tone interferer).

However:

  If I cannot relate these two measurements by a correction factor, is there a way to obtain the receive signal power measurements (in dBm) programmatically (GNUradio perhaps)?

Sure. Why not? You got the raw samples, don't you? So do a magnitude squared of them, and find a linear correction factor. That will work (have done so multiple times), because now you're actually measuring to things that relate to each other (signal amplitude as seen by your SDR frontend's ADC and signal power as seen by the spectrum analyzer simply have a quadratic relation).

Note that you will have to re-calibrate for different frequencies, sampling rates, if settable, analog frontend widths (only applies to B2xx, E3xx, and a few very old daughterboards for your N2xx currently; by the way, the TVRX* boards have always-on AGC, so they can't be used for this).

The relationship between input voltage and digital number amplitude is linear for quite a large range for the Ettus frontends – however, as soon as your input amplifier and mixers start reaching saturation, this will no longer apply, so you did good to use a few dB of attenuation.

If you happen to have a WBX, SBX, CBX or UBX daughterboard, there's IP2/IP3 measurements (along with effective vs set gain) available that will help you estimate the region of linearity on http://files.ettus.com/performance_data/.


Raised-Cosine Filter - calculating magnitude response

The frequency $f$ is no parameter but it's the independent variable. The frequency response $Z(f)$ is a function of frequency. The only real parameter is $\beta$, defining the excess bandwidth, but you already know that $0\le \beta\le 1$. You also know that $T_s$ is the symbol period, so $1/T_s$ is indeed the baud rate, which is part of the system's specifications.

From the formula you can see that the minimum bandwidth is half the baud rate (for $\beta=0$, i.e. zero excess bandwidth). For $\beta=1$ (i.e., $100\%$ excess bandwidth), the bandwidth is twice as large, i.e. it equals the baud rate $1/T_s$.


How to detect, using matlab, the center of the biggest circle?

Would the built-in function imfindcircles work for you?  The radius and center are both returned. The examples are included in the same page. The function is based on Hough transform so it might be an overkill.

Another much faster method is regionprops. You may use this function only after binarizing the image first, and that means there must be a contrast difference between the region inside and outside the circle. You may try something like this (copied from this post):

st = regionprops( ~BW, 'Area', 'Centroid', 'PixelIdxList' );

sel = [st.Area] > numel(BW)*0.5; % exclude the small regions
st = st(sel);

cntr = .5 * [size(BW,2) size(BW,1)]; % Compute region distance to center of image
d = sqrt( sum( bsxfun(@minus,vertcat( st.Centroid ), cntr ).^2, 2 ) );

[mn idx] = min(d); % pick the region closest to center

res = false(size(BW));
res( st(idx).PixelIdxList ) = true; % the region outside the circle is
                                      set to false, thus excluded

Note that in the code above, cntr = .5 * [size(BW,2) size(BW,1)]; assumes the circle is in the center of the whole image. You need to estimate the rough position of the circle in your case as a prior knowledge and change this line.

Several assumptions mentioned above make the regionprops option more tricky than imfindcircles.


This can't be the way to go: If you have an ideal LPF, one-bit DAC is always sufficient. (Hint: Draw the spectrum of a square wave. Draw the pass band of your LPF on top.)

You don't have an ideal LPF; or you might want to produce a range of sine frequencies (not only a single one) spanning more than an octave, so that you can't just use a fixed LPF + 1 bit DAC.

The easy way to calculate THD is realizing that it's a relation between powers. So, you need to know how much "error power" there is in a quantized signal. There's this term quantization noise that captures that – and it's well-known for sinusoidal signals!


Image processing - Is mean filter shift-invariant?

The 1D mean filter defined by the following I/O relation is time/shift invariant;

$$y[n] = T\{x[n]\} = \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d_1}^{n+d_2}{x[k]}$$

For $d_1$ and $d_2$ being positive integers, the filter computes the mean (or the average) of its input $x[n]$ in the range of its summation index. To show that it is shift invariant we shall simply show the truth of the following statement:  $y[n-d] = T\{x[n-d]\}$ for and $d$ and any input $x[n]$

$$ y[n-d] = \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d-d_1}^{n-d+d_2}{x[k]} ~~~\overset{?}{=}~~~ T\{x[n-d]\} = \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d_1}^{n+d_2}{x[k-d]}$$

Which is easily shown to be true by making the simple change of variable in the second sum as $k-d=k'$ and then optionally replacing $k'$ with $k$ again:
\begin{align}
T\{x[n-d]\} &= \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d_1}^{n+d_2}{x[k-d]} ~~~~ &\scriptstyle{\text{replace k-d with k'}}\\
T\{x[n-d]\} &= \frac{1}{d_1 + d_2 + 1} \sum_{k'=n-d-d_1}^{n-d+d_2}{x[k']} &\scriptstyle{\text{replace k' with k again}}\\
T\{x[n-d]\} &= \frac{1}{d_1 + d_2 + 1} \sum_{k=n-d-d_1}^{n-d+d_2}{x[k]} &\scriptstyle{\text{recognise that this is =  y[n-d]}}\\
T\{x[n-d]\} &= y[n-d] &\scriptstyle{\text{Hence the given System is LTI}}\\
\end{align}

This is extended to 2D in a straight-forward manner


How can i do that?

Generate a 370 ms long Hann window, take a 370 ms chunk of your audio, then multiply them together, sample-by-sample.  You're just taking a chunk of audio and fading it in and out at the ends to get rid of any glitches at the start and stop of the chunk.

  What is meant by 28/32?

I'm not sure.  I'd guess it means that the waveforms overlap by 28/32 = 87.5%.  This just means that the chunks of audio you read should overlap with chunks you've already read.

Presumably your sample rate is 44100 Hz, and the window length is 214 = 16384 samples = 371.5 ms?  So

Generate a 16384-point Hann window and store it in memory
Read audio from sample 0 to sample 16383
Multiply them together to get a windowed chunk
Do FFT or whatever other processing on the windowed chunk
Read audio from sample ~2048 to 2048 + 16384 (87.% overlap with previous chunk)
Multiply with window to get a windowed chunk
Do FFT or whatever
Read audio from sample ~4096 to 4096 + 16384
...

The numbers are probably off-by-one, but you get the idea.


How to tune PID parameters to track a family of reference signals?

As already mentioned, a linear system that tracks a step reference, will track it for arbitrary values and time of application. Due to the internal model principle, in a feedback control system, if the denominator of the product
$$
C(s) P(s)
$$
contains a pole at $s = 0$, then it will robustly track a step function.
The control law
$$
C(s) = \frac{k_{d} s^{2} + k_{p} s + k_{i}}{s}
$$
has a pole at $s = 0$. If the plant  $P(s)$  does not have a zero at $s = 0$, then your system will track an arbitrary step function $R(s) = \frac{A}{s}$, if it is closed-loop stable.

Since the complementary transfer-function, from reference to output, is $Y(s)/R(s) = T(s)$, and
$$
T(s) = \frac{C(s) P(s)}{1 + C(s) P(s)} = \frac{\frac{1}{s} \bar{C}(s) P(s)}{1 + \frac{1}{s} \bar{C}(s) P(s)} = \frac{\bar{C}(s) P(s)}{s + \bar{C}(s) P(s)} ,
$$
this can be seen by using the final value theorem:
$$
\lim_{t \rightarrow \infty} y(t) = \lim_{s \rightarrow 0} s \frac{A}{s}\frac{\bar{C}(s) P(s)}{s + \bar{C}(s) P(s)} = A
$$


Correct form of discrete-time Fourier series representation

Because the discrete-time basis functions $e^{j2\pi kn/N}$ are not only $N$-periodic in $n$ (as is the signal $s[n]$), but also in $k$. So adding terms for $k\ge N$ does not add anything new, simply because

$$e^{j2\pi (k+N)n/N}=e^{j2\pi kn/N}e^{j2\pi n}=e^{j2\pi kn/N}$$


As people have said in the comments, to have it sound right you will need to low pass filter before subsampling.

The matlab function decimate does this automatically and gives some control over the type of low pass filter used.


Signal remains constant when implementing FSK

The sampling rate you're using is too low. The cosine signals have frequencies 4,000 and 6,000 Hz, but your sampling frequency is 1000.

An FSK signal has a very large bandwidth -- in theory, it is infinite. This is because of the sharp transistions between pulses of different frequencies. My suggestion is to use a sampling frequency that is around 5 to 10 times larger than your carrier, or around 25,000 to 50,000 samples per second in your case.


Pulse Doppler Radar: time delay in down conversion

The delay is dealt with by knowing it. It is calculated and measured and any residual error is minuscule compared with any other uncertainty.

It’s dealt with design and calibration.


Removing white pixels from a b/w image in MATLAB?

So...you want to remove the fat white borders from around the black region??

You could try search for the first and last row,column that are 0 in the second image, using (e.g.) find:

% assuming BW2 is 2D, one channel.
% the reason the first has a BW2' and the second is BW
% is because of Matlab being column-major when it does `find`:
% have to make sure we actually get the first/last indices.
[col1,row1] = find(BW2'==0,1,'first'); % top-left corner
[row2,col2] = find(BW2==0,1,'last'); % bottom-right corner
BW3 = BW(row1:row2, col1:col2);
imshow(BW3)

This does assume that the top-left and bottom-right corner of the image proper won't be part of the outline - this is a very low chance.

If you want to be sure however, you could instead look for the first row that is not all 255 and first col that is not all 255, and similarly for the last row/col (ie locate margin by looking for entire rows/cols being 255) and crop like that:

cols = sum(abs(255-BW2),1); % work out difference from 255
rows = sum(abs(255-BW2),2);
col1 = find(cols,1,'first');
col2 = find(cols,1,'last');
row1 = find(rows,1,'first');
row2 = find(rows,1,'last');
BW3 = BW(row1:row2,col1:col2);

(Note this will fail on your jpeg because of the compression it uses -- the border is not uniquely 255, there are a few 251's in there. But it will work on your code if you don't save BW out as a jpeg first and then read it back in. If you're going to save images, use a lossless format like TIF or PNG).


Difference between 'spatially' and 'instantaneously' independent

The components are instantaneously uncorrelated if
$$
E[y_i(t_1) y_j(t_2)] = 0
$$
for all $i \not= j$ and for all $t_1,t_2$ and when $i = j$ :
$$
E[y_i(t_1) y_i(t_2)] = 0
$$
for all $t_1 \not= t_2$.

The components are spatially uncorrelated if
$$
E[y_i(t) y_j(t)] = 0
$$
for all $i \not= j$ (assuming a change in index corresponds to a spatial change).

Note that correlation is not quite the same thing as statistical independence.


Can we use Tiva C series TM4C123GXL for signal processing experiment purposes?

This question is apparently out of scope. But I would like to provide you a simplest guideline so that you can by yourself decide on whether a given microprocessor would be adequate for various DSP tasks.

1- The Floating Point Capabilities:
Nowadays any processor claiming high quality DSP computation must have strong floating point hardware support. Look for its FLOPs (floating point operations) throughput. Also look for single precision (32 bits) or double precision (64 bits) FP capabilities. For most multimedia purposes (speech, music, video, image) 32 bit is more than enough. For more accurate scientific computations 64 bits is necessary. This being said, the experts can use integer capabilities of a microprocessor to implement most DSP algorithms. But that would require a lower level handling of the mpu. Depending on what kind of an algorithm you want to implement, you should first determine number of floating point operations per second and then check for the MPU if it can support that.

2- Parallel Execution Capabilities (SIMD)
Many DSP algorithms (especially image and video prcessing types, or communication channel processing) will benefit from single instruction multiple data (SIMD) processing pipeline. It can be expensive to include in a simple MPU but if you have this capability, it's worth the money.

3- ADC / DAC capabilities are vital as you will eventually be inputting an analog signal, convert to digital, process it and cpnvert to analog before outputting. So look for for High Speed, High Accuracy ADC and DAC. Or if your application can tolerate low accuracy go with that one.  And there must be enough number of channels so that if you have multiple signals being input or output you can use that.

4- Memory: in addition to limited program code memory for several applications, data has to be buffered inside RAM. this ram can be on-chip or out of chip. First determine how much ram you will need and then look for wheter it can be on chip (very fast) or out of chip.


Why analog anti aliasing filter is used before analog to digital converter when there is already a digital filter after ADC?

Question 1:
The anti-aliasing filter before the ADC is exactly for the purpose of rejecting high frequencies, that will become lower frequencies (i.e. aliasing) after the ADC. The digital lowpass after the ADC cannot help here, as the aliasing has already happened. Consider this example:

Your ADC has a sampling frequency of Fs=100kHz.
Your input signal is a sum of two sine waves, with frequencies 10kHz and 220kHz.
After ADC you would find two sine waves: one at 10kHz, one at 20kHz (220kHz-2*Fs).
Hence, you have aliasing occured, and no digital lowpass can remove this aliasing.

Question 2:
Without more information on the system this cannot be answered. However, here are some thoughts:

filter of 5kHz requires a sampling frequency of at least 10kHz (ideally). You state you only need 2.5kHz. I think you mix something here.
in reality, no anti-aliasing filter is a perfect low-pass, hence its cutoff-frequency does not mean, that higher frequencies are perfectly blocked. Instead, they are more and more attenuated. To cope with non-ideal anti-aliasing filters, the sampling frequency should be higher than 2 times the cutoff. However, 200kHz from your example occurs still quite high for me.