how to calculate the maximum obtainable directivity for a microphone array? Introduction

The relationship between the maximum obtainable directivity and aperture of an antenna is described by:

$$ D = \frac{4\pi A_e}{\lambda^2} $$

Where $D$ is the directivity, $A_e$ the effective aperture and $\lambda$ the used wavelength. The derivation can be found in the following document: pdf.

However, if the antenna does not scan as a sphere but in $x, y$ but not $z$. Following the same path as the above-mentioned document the equation becomes:

$$ D = \frac{2\pi A_e}{\lambda^2}  $$

The aperture of a lossless isotropic antenna is then also $A_e =  \frac{\lambda²}{2 \pi}$ instead of $\frac{\lambda²}{4 \pi}$

Problem

I cannot manage to use that equation for a microphone array using delay-and-sum beamforming.

For example, if I have an array of 25 microphones spaced as a grid and that my setup has the following properties:

Smallest distance $d$ between them is $0.01$ meters.
Amount of orientations $o$ that it looks at is $64$.
Frequency $f$ to be detected is 8 kHz.

The wavelength $\lambda = \frac{v}{f} = \frac{343}{8000} = 0.0429 $

The aperture $A$ is $A = \frac{\lambda^2}{precision} = \frac{0.0429^2}{\frac{2\pi}{64}} = 0.0187 $

Thus:

$$ D = \frac{2\pi A_e}{\lambda^2} =  \frac{2\pi A}{\lambda^2} = \frac{2 \pi \cdot 0.0184}{0.0429^2} = 62 $$

I rounded the numbers for the simplicity of formulating my problem, if you do not round then you get $64$ instead of $62$ which is the number of orientations that it looks at furthermore the results cannot be true because:

It does not involve the number of microphones used.
It does not involve the smallest spacing between microphones.
I probably made a mistake somewhere.

How can I calculate the maximum obtainable directivity for an array of sensors?

When using a microphone array for baseband acoustic signals, Array Gain is a preferred performance metric over Directivity Index.

In geophysics, the term “stacking” is often preferred over the equivalent “delay-and-sum beam former. It is certainly more descriptive because you essentially time align signals like a stack and add them coherently.

If the noise is independent at each microphone, for the aligned stack, the upper bound on coherent gain is $10 Log N$.

Directivity is a function of the beam pattern which is a function of a single frequency.  For an RF system with an obvious carrier, the pattern at that nominal frequency is an often used way to characterize gain.

The stacked array gain is baseband. It is the upper bound on gain.


Separate Complicated Signal into Exponential Components This morning, I asked a similar question about separating simple signals - with constant frequency - into exponential components. User @AndreasH suggested that a hilbert transform can do this as follows:

$$\cos(\omega_1 t) + j ~\textrm{Hilbert}[\cos(\omega_1 t)]  = e^{j\omega_1 t}$$

This works great when $\omega_1$ is a real constant, and when time is linear. Is there a way to do it for a more complicated signal? For example,

$$x(t)=\cos(\omega_o t^2 ) $$

When I perform a hilbert transform on this one, I get something with a lot of error functions.

How can I get $y(t)=e^{j\omega_o t^2}$ from $x(t)$?

Hilbert transform method from your earlier question and Jason R's answer should still work. What you are interested in is extracting the instantaneous phase of your function. For a signal $x(t)$ this is given by $\phi(t) = \angle x_a(t)$ where $x_a(t)$ is the analytic signal of $x(t)$.

Digitally you can do this using a Hilbert transform too. Here's a code snippet in Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert
plt.style.use('presentation')

if __name__=='__main__':

  t = np.linspace(0,3,100)
  x = np.cos(10*t*t)

  y = np.imag(hilbert(x))

  x_a = x + 1j*y
  inst_phase = np.angle(x_a)

  true_phase = 10*t*t
  true_phase_wrapped = (true_phase-np.pi)%(2*np.pi)-np.pi

  plt.figure(1)
  plt.plot(t, x)
  plt.plot(t, inst_phase)
  plt.plot(t, true_phase_wrapped)
  plt.legend((r'$x(t)=\cos(10t^2)$', r'$\angle x_a(t)$',
      r'$10 t^2$ wrapped to $[-\pi,\pi]$'))

  plt.figure(2)
  plt.plot(t, np.unwrap(inst_phase))
  plt.plot(t, true_phase)
  plt.legend(('unwrapped inst. phase', 'true phase=$10t^2$'))
  plt.show(block=False)


What does the "smoothing optical flow" mean? I have this algorithm for Correlation Based Pyramids of images:

• Create two Gaussian pyramids from the two input
images

• Compute optical flow using “5×5” regions on
smallest pyramid level.

• Smooth optical flow, and use as initial guess for
higher resolution.

• Continue with next level. Search close to guess
from higher resolution.

What is the third step asking me to do?

Optical ﬂow is an approximation of the local image motion based upon local derivatives in a given sequence of images. That is, it speciﬁes how much each image pixel moves between adjacent images. (The same concept extends on 3D as well).

The concept of Optical flow derives from the fact that the moving patterns cause temporal varieties of the image brightness. It is assumed that all temporal intensity changes are due to motion only. This is usually true but there are many exceptions (see below). These variations are measured using derivatives. The optical flow methods try to calculate the motion between two image frames which are taken at times $t$ and $\delta t$ at every pixel position. The following constraint equation can be used to derive the optical flow - which is called as Image constraint equation or optical flow equation:

$$
{{\partial I} \over {\partial x} }\Delta x + {{\partial I} \over {\partial y}} \Delta y + {{ \partial I} \over {\partial t}} \Delta t = 0
$$

Refer to wiki page, or this lecture note for more details.

Thus the computation of diﬀerential optical ﬂow is, essentially, a two-step procedure:

measure the spatio-temporal intensity derivatives (which is equivalent to measuring the velocities normal to the local intensity structures) and
integrate normal velocities into full velocities, for example, either locally via a least squares calculation or globally via a regularization.

However, the above equation alone is not sufficient for computation since it involves two variables with only one equation/constraint. BKP Horn, introduced another constraint in their seminal paper.

BKP Horn, Brian Schunck "Determining  Optical  Flow" Artificial Intelligence, Vol. 16, No. 1–3, August 1981, pp. 185–203.

Accordingly, in a rigid body motion (with non deformation assumption), neighboring pixels have similar velocity. Based on this they derived another constraint based on velocity:

$$
\nabla ^2 u = {{\partial ^2 u} \over {\partial x^2}} + { {\partial ^2 u} \over {\partial y^2}} = 0
$$

and
$$
\nabla ^2 v = {{\partial ^2 v} \over {\partial x^2}} + { {\partial ^2 v} \over {\partial y^2}} = 0
$$
where
$$
u = {{\delta x} \over {\delta t}} \text{ and  } v = {{\delta y} \over {\delta t}}
$$

The above two equations are considered as Smoothing Constraints of optical flow.
Read this: cv-online page and another reference here.


Frequency response of FM modulation/demodulation chain with phase derivative demodulation Frequency modulating a carrier by white noise and then demodulating the complex signal using discrete derivative of phase it appears that the discriminator is acting as a low-pass filter.

How do I calculate the rolloff vs modulating frequency for a given phase-based discriminator's impulse response?

I know about Octave's freqz function, but applying it to the discriminator's impulse response gives the frequency response of the discriminator, which in my understanding would translate to amplitude distortion of the demodulated signal, so this is not directly the frequency response of the modulate->demodulate chain which I'm looking for.

pkg load signal;
orig_mono_240Khz = randn(1,150000);
orig_mono_240Khz = orig_mono_240Khz - mean(orig_mono_240Khz);
max_val = max([abs(max(orig_mono_240Khz)), abs(min(orig_mono_240Khz))])*1.2;

subplot(2,1,1);
plot(abs(fftshift(fft(orig_mono_240Khz))));
axis("tight");

phase_changes = j*75000*2*pi.*orig_mono_240Khz/max_val;
signal = resample(exp(cumsum(phase_changes/240000)), 10, 1);
signal = signal - min(signal);
signal = round(signal/max([max(imag(signal)),max(real(signal))])*255);

signal = signal - mean(signal);
rcv_240 = decimate(signal, 10);
unwrp_phase = unwrap(angle(rcv_240));
unwrp_phase = unwrp_phase - mean(unwrp_phase);
disc_resp = [1, 0, -1];
phase_drv = conv(unwrp_phase, disc_resp)(100:end-100);
phase_drv = phase_drv-mean(phase_drv);

subplot(2,1,2);
plot(abs(fftshift(fft(phase_drv))));
axis("tight");

You are essentially seeing a Cascade-Integrator-Comb (CIC) response which is identical to a moving average filter (Aliased Sinc function magnitude response) as seen with CIC filter structures.

Consider what is happening in units of phase:

You start with a white noise signal which is translated from magnitude directly to units of frequency in the FM modulation process.

To demodulate you convert to phase (which is an integration of the instantaneous frequency of your waveform, since phase is the integral of frequency) so if you evaluated the frequency response of the resulting phase versus time signal you would expect to see the low-pass frequency response of an integrator, rolling off at 1/f.   (You aren't actually doing any integration, and extracting the phase itself from your waveform is not filtering, the phase is already filtered as if you had integrated your original modulation waveform: For example, the phase for white-FM rolls off at 1/f in frequency). Also be aware that for constant envelope waveforms (no AM) the power spectral density is directly proportional to the phase fluctuations-- for small angles, the power of the sidebands is given directly by the phase versus time.

You then compute the differences in phase between samples over N samples to recover the modulating signal-- this is a digital approximation of differentiation and is the classic "Comb Filter", as with N large the magnitude response represents the teeth of a comb.

$$1-z^{-N}$$

So essentially the process matches the cascade of an accumulator (integrator) with a comb, which is identical to a moving average filter over N samples.

To see how these are equivalent:

Moving average as a geometric series:
$$ \begin{align}
& 1+z^{-1}+z^{-2}+z^{-3} + ... + z^{-N-1} \\
& = \sum_{n=0}^{N-1}z^{-n} \\
& = \frac{1-z^{-N}}{1-z^{-1}}
\end{align}$$

accumulator: $\frac{1}{1-z^{-1}}$

comb: $1-z^{-N}$

Also interesting is the poles and zeros of each which also explain how the moving average is identical to a CIC: The comb filter has N total zeros starting at z = 1 and equally spaced on the unit circle (the roots of $1-z^{-N}$). The integrator has a pole at z=1. A moving average filter has N-1 zeros, all at the same locations as the comb filter except for z=1.

From your spectrum it looks like you are using N=1. As you increase N you would expect the envelope of the result to approach a Sinc function consistent with a moving average filter.


transfer function and 'causal' signal - evaluate transfer function or use z-transform of input? From my studying difference equations and transfer functions, I understand that when a complex exponential input $x[n]=z_1^n$ is applied to an LTI system with transfer function $H(z)$, determining the output involves evaluating the transfer function at $z=z_1$, i.e. $H(z_1)$ which is a complex quantity that simply scales the input signal to produce the output.  The textbook I have uses a sinusoid input to show the effect (the transfer function causes a magnitude and phase change to the input).

The answers to this question I asked a while ago made the point of recognizing the difference between $a^n$ and $a^n \cdot u[n]$.  All of that was in the context of two sided signals and the bilateral z-transform.

In something I found online, in the context of causal signals and the unilateral z-transform, there was a question posed: "if $x[n]=(1/2)^n \cdot u[n] \to y[n]=\delta [n-2],$ what is the output for $x[n]=cos(\frac{\pi}{3}n)?$",  where the author determined the output by evaluating $H(z)$ determined from the first relationship at $z=e^{jn \pi /3}$.  Since this is a 'causal' context, isn't that cosine really $x[n]=cos(n\frac{\pi}{3}) \cdot u[n]$ and hence requires taking the z-transform of the input?

When can you just evaluate $H(z)$ at the value of interest, and when do you have to take the (unilateral) z-transform of the input?

For a discrete-time LTI system with transfer function $H(z)$, the response to $x[n]=z_1^n$ equals $y[n]=H(z_1)z_1^n$ if $z_1$ is inside the region of convergence of $H(z)$. This relationship holds whether or not the system is causal or stable.

From this relationship it follows that if the system is real-valued and stable, i.e., if its impulse response is real-valued and if $H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}$ exists, the response to a sinusoidal input $x[n]=\cos(n\omega_0)$ is given by $y[n]=|H(e^{j\omega_0})|\cos[n\omega_0+\phi(\omega_0)]$.

I think your confusion comes from assuming something you call a "causal context". It doesn't matter if a system is causal or not for above relationships to hold. You can apply a sinusoidal input signal to a causal system, and as long as the system is real-valued and stable, the output will be given by the relationship mentioned above. If the input is a sinusoid switched on at a finite $n=n_0$, the system's response can be computed by solving the convolution sum or, equivalently, by using the $\mathcal{Z}$-transform.


Calculating colour temperature from two sources of varying luminous intensity I'd like to calculate the resultant colour temperature if I were to mix two different temperatures, each at a different intensity in lumens.

I have series of LEDs in two different temperatures, and would like to adjust the combined output using the two sources, but I can't find any information on this. Is it as simple as finding the mean, scaling for a higher luminosity in one source? I'm assuming that the latter portion won't be a simple linear relationship due to human perception.

LEDs do not emit black-body (Planckian) radiation which has a specific spectrum, so instead of color temperature which concerns black-body radiators we have to resort to correlated color temperature (CCT). It is the temperature of a black body that emits light of a color that is perceptually closest to the color of light of the same brightness from, in this case, the LEDs. The CIE 1960 $(u,v)$ chromaticity space ignores brightness and is designed to be perceptually uniform. It is used to calculate correlated color temperature and is based on the CIE 1931 $(X, Y, Z)$ color space.

Your LED datasheet gives for the two LED types their nominal fluxes in lumens and also gives their CIE 1931 $(x, y)$ chromaticities. Flux in lumens is calculated from a light spectrum using the CIE 1924 photopic (normal color vision, not low-light) luminosity function which is identical to the CIE 1931 spectral efficiency for the $Y$ color component, by definition. So we can work with normalized $(X, Y, Z)$ color where $Y$ equals the flux in lumens. This normalization is compatible between LED types of different chromaticities. Knowing $Y$ and the chromaticity $(x, y)$, we can easily calculate the rest of the $(X, Y, Z)$ triplet describing a LED type's light color, using the CIE 1931 definitions of $x$ and $y$:

  $$x = \frac{X}{X+Y+Z},\quad y = \frac{Y}{X+Y+Z}\tag{1}$$

$$\Rightarrow\quad X = \frac{Y\,x}{y},\quad Z = \frac{(1 - x - y)\,Y}{y}.\tag{2}$$

The $(X, Y, Z)$ color space is a linear-intensity color space so the $(X, Y, Z)$ colors of the LEDs sum additively in additive mixing of their lights. The correlated color temperature is then calculated from the summed $(X, Y, Z)$ color by going back to $(x, y)$ using the above definitions and by converting to CIE 1960 $(u, v)$ by the definition of that chromaticity space in relation to CIE 1931 chromaticity:

  $$u = \frac{4x}{12y - 2x + 3},\quad v = \frac{6y}{12y - 2x + 3}.\tag{3}$$

These steps could also be combined to a single equation. What remains to do is to find the temperature of a black-body that gives $(u, v)$ coordinates (Planckian locus) nearest to this point by straight-line distance. Krystek, M. (1985), An algorithm to calculate correlated colour temperature. Color Res. Appl., 10: 38-40 recommends using a bisection search together with the following rational Chebyshev approximation valid for $1000 < T < 15000$ where $T$ is temperature in kelvins:

  $$\begin{eqnarray}u = \frac{0.860117757+1.54118254 \times 10^{-4}T + 1.28641212 \times 10^{-7} T^2}{1+8.42420235 \times 10^{-4}T + 7.08145163 \times 10^{-7}T^2},\\
v=\frac{0.317398726+4.22806245 \times 10^{-5}T + 4.20481691 \times 10^{-8} T^2}{1-2.89741816 \times 10^{-5}T+1.61456053 \times 10^{-7}T^2}.\end{eqnarray}\tag{4}$$

The Colour package (install using pip install colour-science) for Python includes the Krystek 1985 method but I just found a bug in their implementation in version 0.3.14 which will be fixed in version 0.3.15. They also have other similar methods in function uv_to_CCT and in other related functions. Example:

import numpy as np
import colour

# Subscript 0 denotes 1800K LED, 1 denotes 4000K LED.
# LED chromaticity center point:
xy_0 = np.array([0.5656, 0.4136])
xy_1 = np.array([0.3860, 0.3870])
# LED nominal flux (lumens) at nominal forward current:
flux_0 = 825
flux_1 = 1370
# Mixing weight (lumens)
weight_0 = flux_0
weight_1 = flux_1

# Proportional XYZ color at flux of 1 lumen (setting Y = 1):
XYZ_0 = np.array([xy_0[0]/xy_0[1], 1, (1 - xy_0[0] - xy_0[1])/xy_0[1]])
XYZ_1 = np.array([xy_1[0]/xy_1[1], 1, (1 - xy_1[0] - xy_1[1])/xy_1[1]])

# Correlated color temperatures of the two LED types
cct_0 = colour.xy_to_CCT(colour.XYZ_to_xy(XYZ_0), method='McCamy 1992')
cct_1 = colour.xy_to_CCT(colour.XYZ_to_xy(XYZ_1), method='McCamy 1992')

# Correlated color temperature of a mixture of the two LED types each at nominal current
colour.xy_to_CCT(colour.XYZ_to_xy(weight_0*XYZ_0 + weight_1*XYZ_1), method='McCamy 1992')

The obtained correlated color temperatures of the two LED types are 1746 and 3934 kelvins, close to their nominal CCTs of 1800 K and 4000 K. At nominal forward current for each LED type, the CCT of the mixture is 2715 kelvins. So the script seems to be working as intended.

We can also get a color temperature curve as function of mixing ratio, and see that it is not exactly linear, although this depends on the definition of a mixing ratio. Continuing in Python:

import matplotlib.pyplot as plt

ccts = np.zeros(101)
for percent_1 in range(101):
    weight_0 = (100-percent_1)*flux_0/100
    weight_1 = percent_1*flux_1/100
    ccts[percent_1] = colour.xy_to_CCT(colour.XYZ_to_xy(weight_0*XYZ_0 + weight_1*XYZ_1), method='McCamy 1992')

plt.plot(range(101), ccts, "-", range(101), cct_0*(100 - np.arange(101))/100 + cct_1*np.arange(101)/100, "--")
plt.ylabel("correlated color temperature (K)")
plt.xlabel("4000 K LED flux / nominal flux (%)")
plt.grid()
plt.show()

Figure 1. Correlated color temperature (solid blue) of a mixture of 1800 K and 4000 K LED types at different mixing ratios. 0 % and 100 % mean a single LED type operating at its nominal flux at its nominal forward current and the other LED type off. 50 % means each LED type is operating at half of its nominal flux. Calculations use manufacturer's LED data. A linear relation (dashed orange) is shown for comparison.


What Are the Differences in Four Different Types of Third Order Image Moments (Skewness)? I am attempting to analyze images using the third order image moment (aka skewness), but I am having trouble figuring out how to do so. As it turns out, there are four different types of third order image moments:
$$
\begin{array}{rcl}
\mu_{21} &=& M_{21} + 2\bar{x}M_{11} - \bar{y}M_{20} + 2\bar{x}^2M_{01}\\
\mu_{12} &=& M_{12} + 2\bar{y}M_{11} - \bar{x}M_{02} + 2\bar{y}^2M_{10}\\
\mu_{30} &=& M_{30} - 3\bar{x}M_{20} + 2\bar{x}^2M_{10}\\
\mu_{03} &=& M_{03} - 3\bar{y}M_{02} + 2\bar{y}^2M_{01}
\end{array}
$$
$$
\begin{array}{rcl}
M_{p,q} &=& \mbox{Raw moment of order}\ (p+q)\\
\mu_{p,q} &=& \mbox{Central moment of order}\ (p+q)\\
\end{array}
$$
Formulas pulled from: https://en.wikipedia.org/wiki/Image_moment
It's been very difficult to find out what the differences are between these 4 (if there are any) on the internet. If there's anyone here that could explain in layman's terms how they differ, I would greatly appreciate it.
This question is most likely asking something that's more fundamental to the calculation of moments themselves than specifically with regard to image moments. However, I have seen several questions about them here, so I figured it would be a good place to start.

It it based on the definition of Central Moments on 2D grid given in the Wikipedia Page.

$$ {\mu}_{pq} = \int\limits_{-\infty}^{\infty} \int\limits_{-\infty}^{\infty} (x - \bar{x})^p(y - \bar{y})^q f(x,y) \, dx \, dy $$

Since we're dealing with 2 random variables, a multiplication which the powers sum to 3 can be achieved in 4 different methods.

It's up to you to decide how to calculate the third moment.

Usually I'd use the Rotation Invariant Moments.
You can find the definition on the same page.
Enforcing the Rotation (And Scale) Invariance yields single definition.


Why are the basis functions for DFT so? When you get a DFT of a signal, you use the basis functions as:

$e^{-j2\pi kn/N}$

Why is it so? Why don't we use the conjugate, $e^{j2\pi kn/N}$, or any other function?

If you look at the idea of Continuous Time Fourier Series (CTFT), it says any periodic signal can be constructed by summation of infinite number of complex exponentials. But in descrete case , only 'N' different  complex exponentials are enough , because there are only 'N' distinct complex  exponentials exists (N is the period of discrete signal).

DFT is nothing but DFS (Discrete Fourier Series).So You can extend the
same idea to get the answer.


how to compute a discrete fourier transform on fragmented data I have a list of collected data points that I need to take a DFT of, however with the problem that about one quarter of the data points in the middle are missing, and so even though the existing data points are evenly spaced their is still a large time gap that needs to be accounted for. I looked up some things on a variable sampling rate DFT, but it was first of all over my head and secondly I am not sure if it is applicable. My question then is how would I account for that time gap in my calculation. It would be best to have everything explained in as close to laymen's terms as possible as I am out of my depth in this.

Thank you for any help.

Basically, it is impossible to calculate the full DFT of your data. Take a look at the DFT matrix.

All of its elements are non-zeroes. Thus, each and every signal element is important for any frequency.

Take a look at this another way, had it been possible to do what you want, you could find the DFT of the signal, do an inverse DFT transform, and get the full signal back! Obviously, this is wrong. You don't have this information.

Nevertheless, it is possible to compute the frequencies of your data approximately. Recall that one of the possible explanations of DFT transform is that each of its coefficients is the best fit to a sum of shifted sine waves, $A* sin(w (x-x_0))$, or $a sin(wx)+b cos(wx)$, where $w$ is dependent on the coefficient index. Thus, you should do a least-squares fit to $a_i,b_i$, which are the coefficients of $a+bi$ in DFT.


State-Space Representation of Forward and Backward Filters In [1], the author shows an efficient way of implementing the forward and backward filter using matrices. One can also implement this using filtfilt command in MATLAB. However, I am interested in implementing this using matrices.

Background

Any digital filter can be represented in its state-space format as follows:
\begin{align}
x(t+1) &= F x(t) + G u(t) \notag \\
y(t) &= H x(t) + D u(t) \notag
\end{align}

This filter can be expressed as
\begin{align}
\mathcal{Y} = \mathcal{H} U + \mathcal{O}x_0
\end{align}
where $\mathcal{Y} = [y_0, y_1,\ldots,y_{N-1}]^T$ is the vector of outputs, $U = [u_0,u_1,\ldots,u_{N-1}]^T$ is the vector of inputs, $x_0$ is the initial state,
\begin{align}
\mathcal{H} = \begin{bmatrix} D  & 0 & \ldots & 0 \\
                              HG & D & 0 & \vdots \\
                              \vdots & & \ddots & \\
                              HF^{N-2}G & \ldots & HG & D
\end{bmatrix},
\end{align}
and
\begin{align}
\mathcal{O} = \begin{bmatrix} H \\ HF \\ \vdots \\ HF^{N-1} \end{bmatrix}.
\end{align}

Problem Statement

In forward-backward filtering proposed in [1], the author mentions that the forward and backward filters are different (generally speaking). However, if we look at the implementation of the forward-backward filter derivation, the same filter is used, i.e., only the time-reversal of the input and outputs are performed but the filter transfer function does not change. What should be an efficient representation of the state-space filter while applying the forward filter and backward filter, i.e., determine $(F_f, G_f, H_f, D_f)$ and $(F_b, G_b, H_b, D_b)$.

Reference

[1] F. Gustafsson, "Determining the initial states in forward-backward filtering," in IEEE Transactions on Signal Processing, vol. 44, no. 4, pp. 988-992, Apr 1996. http://www.diva-portal.org/smash/get/diva2:315708/FULLTEXT02

The backwards filter, as the name suggests filters the input signal with the same filter $(F,G,H,D)$ backwards in time

$$
\begin{align}
x(t-1) &= F\,x(t) + G\,u(t) \\
y(t) &= H\,x(t) + D\,u(t)
\end{align} \tag{1}
$$

But in order to define a filter using a state space model, the update law does have to go forwards in time, which can be done by solving the first equation for $x(t)$

$$
x(t) = F^{-1}\,\left(x(t-1) - G\,u(t)\right) \tag{2}
$$

however $x(t)$ can only a function of $x(t-1)$ and $u(t-1)$. This can be solved by using a new definition for the state $z(t) = x(t) + F^{-1}\,G\,u(t)$, so

$$
z(t) = x(t) + F^{-1}\,G\,u(t) = F^{-1}\,x(t-1) = F^{-1}\,\left(z(t-1) - F^{-1}\,G\,u(t-1)\right) \tag{3}
$$

which can be rewritten in the following more standard form

$$
z(t+1) = F^{-1}\,z(t) - F^{-2}\,G\,u(t). \tag{4}
$$

Substituting the definition of $z(t)$ in to the expression for $y(t)$ gives

$$
y(t) = H\,x(t) + D\,u(t) = H\,\left(z(t) - F^{-1}\,G\,u(t) \right) + D\,u(t) \tag{5}
$$

which can be simplified down to

$$
y(t) = H\,z(t) + \left(D - H\,F^{-1}\,G\right)\,u(t). \tag{6}
$$

Combining equation $(4)$ and $(6)$ gives the following (forwards) state space model, which is equivalent to the (backwards) state space model in equation $(1)$

$$
\begin{align}
z(t+1) &= F^{-1}\,z(t) - F^{-2}\,G\,u(t) \\
y(t) &= H\,z(t) + \left(D - H\,F^{-1}\,G\right)\,u(t)
\end{align} \tag{7}
$$

So the backwards filter based on the forwards filter $(F,G,H,D)$ can be written as $(F^{-1},-F^{-2}\,G,H,D - H\,F^{-1}\,G)$.


About Discrete Fourier Transform vs. Discrete Fourier Series I am new to the field of signal processing. I am wondering what is the difference between DFS(Fourier Series) vs. DFT(Fourier Transform).

For common applications, usually we get a segment(length N) of digital waveform(like a audio segment), and then we apply FFT(DFT) and then do post-analysis with it.

I am wondering if we can use DFS(thus not using DFT at all) all the time and just assume the waveform segment is repeated with period N.  Would this naive thinking/approach cause any problems?

okay, i'm gonna expound a little.
quoting (except for any typos that may result) from the 1989 text of O&S (Introduction to Chapter 8, The Discrete Fourier Transform, p 514):

Although several points of view can be taken toward the derivation and interpretation of the DFT representation of a finite-duration sequence, we have chosen to base our presentation on the relationship between periodic sequences and finite-length sequences. We will begin by considering the Fourier series representation of periodic sequences.  While this representation is important in its own right, we are most often interested in the application of Fourier series results to the representation of finite-length sequences. We accomplish this by constructing a periodic sequence for which each period is identical to the finite-length sequence. As we will see, the Fourier series representation of the periodic sequence corresponds to the DFT of the finite-length sequence. Thus our approach is to define the Fourier series representation for periodic sequences and to study the properties of such representations. Then we repeat essentially the same derivations assuming that the sequence to be represented is a finite-length sequence. This approach to the DFT emphasizes the fundamental inherent periodicity of the DFT representation and ensures that this periodicity is not overlooked in applications of the DFT.

section 8.1, p 516 on the DFS:

Eq. (8.11) $\quad \tilde{X}[k] = \sum\limits^{N-1}_{n=0} \tilde{x}[n] \ e^{-j2\pi n k/N} $

Eq. (8.12) $\quad \tilde{x}[n] = \frac{1}{N} \sum\limits^{N-1}_{k=0} \tilde{X}[k] \ e^{+j2\pi n k/N} $

regarding the DFS, $\tilde{x}[n]$ (with the tilde) is defined to be periodic with period $N$ such that $$ \tilde{x}[n+N] = \tilde{x}[n] \quad \forall n $$ and $\tilde{X}[k]$ turns out to also be periodic with period $N$ (so $ \tilde{X}[k+N] = \tilde{X}[k] \quad \forall k $)
later, in section 8.6, p 532 on the DFT:

Eq. (8.59) $\quad X[k] = \begin{cases} \sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\ 0, & \text{otherwise} \end{cases} $

Eq. (8.60) $\quad x[n] = \begin{cases} \frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ 0, & \text{otherwise} \end{cases} $

Generally the DFT analysis and synthesis equations are written as

Eq. (8.61) $\quad X[k] = \sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N} $

Eq. (8.62) $\quad x[n] = \frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N} $

In recasting Eqs. (8.11) and (8.12) in the form of Eqs. (8.61) and (8.62) for the finite-duration sequences, we have not eliminated the inherent periodicity. As with the DFS, the DFT $X[k]$ is equal to samples of the periodic Fourier transform $X(e^{j\omega})$, and if Eq. (8.62) is evaluated for values of $n$ outside the interval $0 \le n \le N-1$, the result will not be zero but rather a periodic extension of $x[n]$. The inherent periodicity is always present. Sometimes it causes us difficulty and sometimes we can exploit it, but to totally ignore it is to invite trouble.

so the first obvious thing i would say is that the tildes used for the DFS (to explicitly depict a periodic sequence) are symbols and still do not change any mathematical fact.  the direct relationship between the periodic $\tilde{x}[n]$ and the "finite-length" $x[n]$ is
$$ \tilde{x}[n] = x[n \bmod N] \qquad \forall n \in \mathbb{Z}, \ N \in \mathbb{Z}>0$$
where $ \qquad\qquad\qquad n \bmod N = n - N \left\lfloor \frac{n}{N} \right\rfloor  $
now i know some folks will point to the Eqs. (8.59) and (8.60) definition of the DFT that has truncated (to $0$) values outside of the interval $0 \le n,k \le N-1$.
however, that definition is contrived.  it could just as well be expressed as
$\quad X[k] = \begin{cases}
\sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\
5, & \text{otherwise}
\end{cases} $
$\quad x[n] = \begin{cases}
\frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ 5, & \text{otherwise}
\end{cases} $
or
$\quad X[k] = \begin{cases}
\sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\
5000, & \text{otherwise}
\end{cases} $
$\quad x[n] = \begin{cases}
\frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ 5000, & \text{otherwise}
\end{cases} $
or
$\quad X[k] = \begin{cases}
\sum\limits^{N-1}_{n=0} x[n] \ e^{-j2\pi n k/N}, & 0 \le k \le N-1 \\
\text{the man on the moon}, & \text{otherwise}
\end{cases} $
$\quad x[n] = \begin{cases}
\frac{1}{N} \sum\limits^{N-1}_{k=0} X[k] \ e^{+j2\pi n k/N}, & 0 \le n \le N-1 \\ \text{and his hot girlfriend}, & \text{otherwise}
\end{cases} $
because that $0$ in that contrived DFT definition will never ever be used in any theorems regarding the DFT.  when that contrived definition is used for the DFT, then when using any DFT theorems to do any real work (other than the linearity and scaling by constant theorems), then one must use modulo arithmetic in the arguments of $x[n]$ or $X[k]$.  and using that modulo arithmetic is explicitly periodically extending the sequence.
so (sorta responding to hotpaw) there are two or three processes that you should think about when using the DFT on a real signal.

the sampling process. what happens to the spectrum of $x(t)$ when you sample it with a "dirac comb" or whatever you want to call the sampling function?

windowing to finite length. what happens when you window either $x(t)$ or the sampled version, $x[n]$, with a rectangular window of length $N$?

periodic extension.  what happens when you periodically extend it by repeatedly shifting the windowed $x[n]$ by $N$ samples and overlap and add it?

deal with each step by itself.


Is the DCT prone to spectral leakage like the DFT? My understanding is that the DCT assumes the input to be an even periodic signal.
This avoids the discontinuity that can occur in the periodic extension assumed by the DFT (see image).
Does this lack of discontinuity prevent spectral leakage from occurring in the DCT?
If so, why are input signals still often windowed prior to DCT?
Thanks

A window function other than rectangular can be applied to suppress sidelobes also with the discrete cosine transform (DCT). Window functions are also sometimes used together with some flavors of DCT as a lapped transform, with pre and post windowing providing some protection against reflection artifacts arising from processing (such as quantization) of the DCT data.
Here is tested to multiply a real sinusoidal input by 1) a rectangular window or 2) a Hamming window, before DCT analysis. The tested inputs are cosines with different phase shifts. Phase shifts that are multiples of 90 degrees give identical results in the magnitude of DCT. Each plot shows the magnitude of DFT in decibel scale for multiple different fractional offsets of the sinusoid frequency from a DCT bin frequency, in units of bin width. Octave's dct was used for the calculations.
DCT of cosine with 0 degree phase shift:
Rectangular window:

Hamming window:

DCT of cosine with +45 degree phase shift:
Rectangular window:

Hamming window:

DCT of cosine with -45 degree phase shift:
Rectangular window:

Hamming window:

Note that depending on the phase of the real sinusoid, its DCT, and DCT of the sinusoid multiplied by a window function, will be more or less spiky, which is another phenomenon than the window's much lower sidelobes.
Octave source:
pkg load signal;
phase = 2*pi/8; # Phase shift of the cosine
N = 256; # Sequence length
M = 8; # Number of fractional shifts that are in range 0..1 of DCT bin width
n = [0:N-1]'; # Zero-based index
omega = 2*pi*(0.25 + 0.5*[0:M]/M + 0.5*round(N/M)*[0:M]); # Cosine frequency
x = 2*real(exp(i*(n.*omega/N + phase))); # Cosine with phase shift
plot(n/(M*round(N/M)), 20*log10(abs(dct(x)/sqrt(N))));
xlim([0, 1]);
ylim([-60, 0]);
title(["rectangular window, cosine phase shift = " num2str(360*phase/(2*pi)) " deg"]);
xlabel("bin (not labeled), for each peak: fractional offset of cosine from bin center");
ylabel("|DCT(x)[k]| (dB)");

y = x .* hamming(N);
plot(n/(M*round(N/M)), 20*log10(abs(dct(y)/sqrt(N))));
xlim([0, 1]);
ylim([-60, 0]);
title(["Hamming window, cosine phase shift = " num2str(360*phase/(2*pi)) " deg"]);
xlabel("k (not labeled); for each peak: fractional offset of cosine from bin center");
ylabel("|DCT(x)[k]| (dB)");


Signal remains constant when implementing FSK I am trying to implement FSK Modulation. For bit 0 and 1, I defined a frequency wave (for instance, for bit 1, c = cos(2pift), with t = 0 : 1/fs : (1/fs)*N - 1/fs, where fs is the data rate). When I'm simulating, c is remaining constant (only value 1 for cosinus).

Why ? Because If a use a sinus, it is not remaining constant.

Thanks for avance.

Regards.

Code :

N = 1000; % number of bits
R = 1000; % bit rate
fs = 1/R;  % sampling rate

t = 0 : fs : N*fs - fs;
bits = round(rand(1,N));

f1 = R * 6;
f2 = R * 4;

c1 = A*cos(2*pi*f1*t); % bit 1
c2 = A*cos(2*pi*f2*t); % bit 0

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.


FIR design using firpm (Parks-Mclellan) Can you please help me with some ideas on how to implement an FIR filter in Matlab using Parks-Mclellan algorithm?

My specifications are the passband cut-off frequency, the transition bandwidth. The I also have the passband ripple, which is an alternating value $ \delta_p=\pm0.2$ and the stopband ripple $\delta_s$. This because the weigth in the stopband frequency range comes from: $ \rho(\omega)=\delta_p/\delta_s$. The design using firpm should meet these specifications.

I have two troubles: The first one is that I have to find the optimal minimum order of the filter N to get the desired  design without using any Matlab builtin function (I saw there is one called firpmord, but again I can't use it), so I imagine I need some kind of iteration but don't know how to put this in code. Also can't find a way to account for the alternating ripple to use it in my calculations. I'm just looking for some advice on how to achieve this.

I'm just starting, trying to catch it:

wp=0.1;
w_tr=0.3;
w_stop=1;
Gp=1;
dp=0.2; % this should be plus minus
ds=0.05;
N=100
h=firpm(N,[0 wp w_tr w_stop],[Gp Gp 0 0],[1,dp/ds])
[h_m,w]=freqz(h,1);
semilogy(w,abs(h_m));

First of all make sure that your $\delta_p$ value is linear (and not in dB). The function firpm needs a linear value. Don't worry about the passband error ripple having both signs; the value of $\delta_p$ used to compute the weighting function is just the maximum error magnitude, which is $0.2$ in your case. So you correctly defined the stopband weight by $0.2/0.05=4$.

The question is if you need to write an algorithm that automatically finds the necessary filter order, or if you can just try to find the order manually. If you need an automatic procedure you could use the second output argument of [h,err]=firpm(...), which is the maximum approximation error. You could design filters with increasing orders in a loop until the error returned as a second argument from firpm is smaller than the desired value for $\delta_p$ (because the approximation error equals weight times deviation, and since you chose the passband weight to be $1$, the error equals $\delta_p$).

Finally, note that the stopband edge equals the passband edge plus the transition bandwidth. So if the transition bandwidth is really $0.3$, as in your example, then the stopband edge would be the cut-off frequency (passband edge) plus $0.3$ (which equals $0.4$ in your example). These specifications are extremely mild and can be satisfied by a filter of very low order.


How to find out if a transform matrix is separable? In image processing, when we have a series of basis images, how could we know if the transformation is separable or not?
For example, I know that following bases are separable and transformation can be done with matrix $A$, but I don't know how to find $A$ from $A_1$, $A_2$, $A_3$ or find separability.
\begin{equation}
A_1 = \frac{1}{2}\begin{bmatrix}
    1 & 1 \\
    1 & 1
\end{bmatrix}
A_2 = \frac{1}{2}\begin{bmatrix}
    1 & -1 \\
    1 & -1
\end{bmatrix}
A_3 = \frac{1}{2}\begin{bmatrix}
    1 & -1 \\
    -1 & 1
\end{bmatrix} \\
A = \frac{1}{\sqrt 2}\begin{bmatrix}
    1 & 1 \\
    1 & -1
\end{bmatrix}
\end{equation}

I admit I did not really thought about it before. I hope my notations won't be too sloppy.

I assume that given an operator matrix $A(u,v)$, you can apply this operator as a transform on an image $I$, to obtain an image in a novel domain $J(u,v)$. For instance a Fourier kernel would give $$a_{m,n}(u,v) = \exp^{-2\pi \imath \left(um/M+vn/N\right)}\,.$$

For each choice of a $(u,v)$ pair, you get a specific kernel, for instance a basis element.

The transform is separable, in the common sense, when you can write it as a element-wise product, for instance:
$$a_{m,n}(u,v) = b_{m}(u) . c_{n}(v)  $$
where $b(u)$ and $c(v)$ are 1D  vectors:
$$A(u,v) = b(u) . c^T(v)\,,  $$
one acting on  rows, the other on columns.

Thus, the vectors being of rank at most one (the rank of a matrix is lower than its minimum size), $A$ is at most of rank $1$ as well. I take it as sufficient condition. If I rule out uninteresting matrices and vectors of rank $0$, that should be a necessary condition. More details for instance at A rank-one matrix is the product of two vectors.

So:

if $A$ is of rank $>1$, it does not appear to be separable (at least with this definition)
if $A$ is of rank $1$, take a non-zero vector for the span, and then multiply it by appropriate coordinates to get the matrix back.

For $A_3$, the determinant is zero, so the matrix is rank $1$. You will get vector $[1 -1]^T$ in the image, and quickly:

$$2 A_3 = [1 -1]^T [1 -1]\,.$$

Note, the separation is not unique.


Getting a part of the invert of Si function I'm working on an algorithmus that sould  need the inverse of the SI function.

since this last once isn't oviouse to calculate, I'm trying to get only the parobla part from the SI function, to be clair , say I have this figure :

I'm only interessted in the part from the max of each curve tell the point of intersection with the following one. So my question what is the fomula for the "parabola part " in si function.

the ploted function are all : si(3pi X /4)

thanks for any help.

There are two problems, I think: fitting a parabola to the higher part, and moving the parabola to the next part of the sequence. This will address fitting the parabola.

I'd just use a least squares technique to find the $a$, $b$, and $c$ parameters.

This shows the part of the paper of interest:


Confused with using Fourier Transform Properties I'm seeking the Fourier transform of $t \cdot (u(t+1)-u(t-1))$

Given the fact that (The first minus here should be a plus)

Using properties.
My initial approach was:

but I wasn't sure how I should go around nesting properties, so I decided that I can perhaps:

However, checking Wolfram, we have:

I don't understand what's going on, my questions are:
1- What went wrong in my second attempt?
2- How would I carry on if I chose to go with my first attempt?
Sorry if these questions seem naive or anything, I've just been so recently introduced to the topic.

Your first approach doesn't work because the original signal is simply not the integral of the signal you sketched below it. Just graphically integrate that signal to see this. Your second approach doesn't work for the very same reason. If you integrate the rectangular signal  you don't get a signal that is zero for $t>1$. The value of the integral would remain at its final value and wouldn't change anymore because the original signal doesn't change for $t>1$.
That exercise can be solved in several ways, and you should try them all to get some practice. First, just compute the Fourier integral without trying any smart tricks. The signal is simple enough for this to be a good option, even if just as a reference.
Second, convince yourself that if $X(j\omega)$ is the Fourier transform of $x(t)$, then $jX'(j\omega)$ is the Fourier transform of $tx(t)$, where  $X'(j\omega)$ is the derivative of $X(j\omega)$ with respect to $\omega$. Since the signal is defined as $t$ times some other signal (whose Fourier transform is very basic), using the above-mentioned property is quite obvious.
A third approach could be to define a signal $x_1(t)$ which satisfies
$$x(t)=\int_{-\infty}^tx_1(\tau)d\tau$$
Then compute the Fourier transform of $x_1(t)$ and use the integration property of the Fourier transform to derive $X(j\omega)$. This is what you actually tried to do, but your definition of $x_1(t)$ was wrong. As a hint, note that for the correct definition of $x_1(t)$ you need two Dirac delta impulses.
Obviously, if you do things right, all three results should agree.


Please help me understand this FM carrier frequency inequality The link of the statement is evidenced in this book:

I don't understand how equation 85 come? I try to use Carson rule to see this inequality but turns out I have $2 \Delta f$ instead of $3 \Delta f$.

So someone could help me to understand this?

Checking the book you linked to, this equation refers to FM modulation over a (weakly) non-linear channel, modeled as a memoryless third-order non-linearity (Eq. 82). The third-order non-linearity produces signals at twice and at three times the carrier frequency $f_c$ (Eq. 84). Note that also the frequency deviations of these additional signal components at $2f_c$ and at $3f_c$ are doubled and tripled, respectively.

If $\Delta f$ is the frequency deviation of the original FM signal, the signal component at $2f_c$ has a frequency deviation of $2\Delta f$, and according to Carson's bandwidth rule, that signal component extends down to approximately

$$f_1=2f_c-2\Delta f-W\tag{1}$$

where $W$ is the bandwidth (highest frequency) of the message signal. The original FM signal centered at $f_c$ extends up to approximately

$$f_2=f_c+\Delta f+W\tag{2}$$

In order to be able to separate the original FM signal from the signal component at $2f_c$ created by the non-linearity, we require $f_2<f_1$, i.e.,

$$f_c+\Delta f+W<2f_c-2\Delta f-W\tag{3}$$

or, equivalently,

$$f_c>3\Delta f+2W\tag{4}$$

which is the same as Eq. $(85)$ in the book.


Why do we need to conjugate complex signals in autocorrelation and cross correlation Why is it necessary to conjugate $f(t)$ while performing auto correlation or cross correlation with respect to $g(t)$, if $f(t)$ and $g(t)$ are complex signals?

$\underline{Prologue :}$

Let me ask you another question, How will you compare two complex numbers $U$ (a+jb) and $V$ (c+jd)?
By comparing magnitude? Subtract them and take real part? Multiply them and compare?

Since any complex number involves two entities ( one for magnitude |$z$| and other for argument $\theta$ ) any comparison involves comparison of those two entities.
Suppose we represent $U$ and $V$ in a complex argument plane as $M e^{j \theta}$ and $N e^{j \phi  }$  if we multiply both of them together, we get

$U$ x $V =  M e^{j \theta}   *  N e^{j \phi } = MN e^{j ( \theta +  \phi )}  $

If we seek maximum value out of the above multiplication then

$|U$x$V|_{max} =  |MN|_{max} * |e^{j ( \theta +  \phi )}|_{max} $

Suppose our comparison is only about the relative location of $V$ with respect to $U$ in the complex argument plane, for the time being, forget about the magnitudes $M$ and $N$.

maximum value of the second term $\;e ^{j ( \theta +  \phi )}$  occurs at 1.

for that to happen
$\theta+\phi = 2\pi k \;\;\;\;$ where $\:k\:\epsilon\;( 0,1,2...$

if $k = 0$,

$\theta + \phi = 0$

$\phi =  - \theta $

ie, when two complex numbers with bounded magnitudes get multiplied, their maximum value occurs, when the second number ($V$) is in the negative argument direction of the first number ($U$), ie in the direction of its conjugate ($U^*$).

ie
maximum of $U \times V$ occurs when   $V$ points in the  direction of $U^*$ { the orange dashed line in the picture }

$\textit{-or conversely-}$

If $V$ is getting multiplied by $U^*$ then their maximum occurs when  $V \textit{points in the direction of } U$

In another words $U^* \times V$   gives nearness of $V$ w.r.t $U$ in terms of angle ( complex argument ) and it will decrease as a function of $cos \:\alpha$, where $ \alpha$ is the angle between $U$ and $V$.

$\underline{Epilogue:}$

Considering two signals $f$ and $g$, correlation is given as

$\hspace{1cm} (f \otimes g)(\tau) = \int_{-\infty}^\infty f^*(t) g(t + \tau) dt$

At any instant '$t$', a signal is just a point in the complex plane. So at any moment $f$ and $g$ are just two points in the complex plane. Then the job of comparison of two signal becomes a mere comparison of two complex points.
With this logic, to compare two signals with a given lag '$\tau$' between them, just multiply point by point in '$t$' by taking conjugate of the other signal. If you get a big number that means both reference and compared signals are just looking in the same direction in the complex plane.
When integrated on the whole '$t$', the value of the integration shows the how much similarity is between signal $f$ and signal $g$.
that is which is known as the cross-correlation between two signals.

{ I didn't talk about magnitudes $M$ and $N$, right? The method in which everything is put in a single benchmark by avoiding magnitudes $M$ and $N$ is known as Normalized cross-correlation }


Histogram Intersection with two different bin sizes? I have two images: A and B. Image A is a screenshot of image B, the original image. While there is noise in image A, the original features are still distinguishable.

Histogram intersection tells us that:

$$K(a,b) = \sum_{i=1}^{n}\min(a_i,b_i)$$

where if $K$ is large, $a$ and $b$ are similar.

If we process image A to create a histogram with an optimal amount of ($n \pm x$) bins, how do we apply histogram intersection in this case where image A and B don't have the same amount of histogram bins?

I don't see any other way than to adapt one of the histograms such that both histograms have the same number of bins before performing histogram intersection. This will always involve some guesswork since the histograms themselves don't give you enough information to do this accurately. I would adapt the histogram with the larger number of bins to the number of bins of the other histogram. This is a simple example:

h1 = [12 17]
h2 = [9 10 11]

So in this case you want to change h2 such that it has two bins like h1. Since the interval boundary between the two intervals of h1 is right in the middle of bin 2 of h2, it is reasonable to divide bin 2 of h2 evenly between the two other bins. So you get

h2_new = [14 16]

If the new boundary falls somewhere else inside a bin interval just divide the interval by the respective percentage. This step of course involves guesswork and it would be much better if you could avoid it by making sure that both histograms are computed using the same number of bins.

As a final note, usually histogram intersection is normalized such that its maximum value is $1$:

$$K(a,b)=\frac{\sum_{i=1}^n\min (a_i,b_i)}{\min\left(\sum_i a_i,\sum_i b_i\right)}$$


Smoothing power spectrum by convolution with boxcar function I am trying to smoothing a signal's power spectrum by convolving the spectrum with a boxcar function in frequency domain. However, the result is obviously not what I expected: original two frequency spikes become three and the frequencies are different.

Could anyone possibly point out my problem? Is the problem in the idea of convolving with a boxcar or in the MATLAB code?
Thanks...
t = 0:.001:.25;
x = sin(2*pi*50*t) + sin(2*pi*120*t);

f = 1000/251*(-floor(length(t)/2):length(t)/2);

% power spec
Fy = fftshift(fft(x,length(t)));
Pyy = Fy.*conj(Fy)/length(t);
plot(f(length(t)/2:length(t)*3/4),Pyy(length(t)/2:length(t)*3/4))
title('Power spectral density')
xlabel('Frequency (Hz)')
figure

% rectangle function
rectLen = 5;
rw = ones(rectLen, 1);
plot(rw)

% convolution in frequency domain to smooth
Fy_conv_rw = conv(Fy, rw, 'same');
Pyy_conv_rw = Fy_conv_rw.*conj(Fy_conv_rw)/length(t) / sum(abs(rw));
plot(f(length(t)/2:length(t)*3/4),Pyy_conv_rw(length(t)/2:length(t)*3/4))
title('Power spectral density smoothed')
xlabel('Frequency (Hz)')

The spectrum is periodic so you need circular convoluioint, not linear convolution
You are smoothing the  complex linear spectrum. Due to the phase of the complex ammplitude you will get significant interference in your smoothing window. Smoothing is typically in the energy domain and the phase is discarded. Complex smoothing is fairly difficult to do.


Why does the control effort increase if the bandwidth of the system is higher then the bandwidth of the plant? I am studying control systems and I am focusing on the control effort.
I have seen that if I place a pole at higher and higher frequencies, if the bandwidth of the system goes above the bandwidth of the plant, the control effort increases at high frequencies. In particular, at higher frequencies than the bandwidth of the plant, the control effort behaves as the controller, while at lower frequencies it behaves as the plant, and the frequency response of the control effort increases until the frequency of the open loop bandwidth.
The definition of control sensitivity function is:
$Q(j\omega)=\frac{C(j\omega)}{1+C(j\omega)P(j\omega)}$
by looking at this definition, we can make some approximations, so if $|L(j\omega)|>>1$ then $Q(j\omega)=\frac{1}{|P(j\omega)|}$ and if $|L(j\omega)|<<1$ then $Q(j\omega)=|C(j\omega)|$ where $L(j\omega)=C(j\omega)P(\omega)$. From the plot below, this means that for high frequencies, the control effort behaves like the controller, while for low frequencies it behaves like the plant.
This is the plot I have created using Matlab.

where by control effort I meant the frequency response of the control sensitivity function and the open loop is $L(j\omega) = C(j\omega)P(j\omega)$. Moreover, I have plotted the frequency response of the plant and of the controller. Just to be more clear as possible, the code for plotting this plot is:
s = tf('s');
P = 1/[(1+s)*(1+0.05*s)^2];
C = (s+1)/s;

tau_3 = 0.001;
CF_3 = [(1+s)*(1+0.05*s)^2]/((1+tau_3*s)^3);

Q3 = (C+CF_3)/(1+P*C);

C3 = (C+CF_3);

L3 = C3*P;

figure;
bodemag(Q3,'g',C3,'r',P,'b',L3,'y'),grid
legend('control effort Q(jw)','controller C(jw)','plant P(jw)','open loop L(jw)')
title('behaviour of the control sensitivity function')

In this case I am considering a two degrees of freedom controller, which is why I have used C+CF_3 in my code.
Why does this happen?

Intuitively (and I leave it as an exercise for you to work it out with math), it's because you have to push harder on the plant to get it to move fast.

Consider a bowl of honey and a spoon -- moving the spoon around in the bowl slowly takes little effort, but moving it fast requires a lot of effort, because honey is viscous.

It may be most informative to take a simple plant model, like a 1st-order low-pass filter, and without ever designing a closed-loop controller, design a signal for it to follow that has significant energy content above it's rolloff frequency.  Then look at the relative amplitudes of the components of that drive signal compared to the resulting components of the motion of the plant.

Even in open loop your control effort is higher at higher frequencies, for reasons very much analogous to that spoon in honey.


Where did the length of time’s period disappear in Periodic Fourier Series Discrete Time In continuous time Periodic Fourier Series has smallest n as possible, since it is an integral and a length of the repeating time (period time) which is T0. In discrete time however we don’t have units of time at all. Explain in your own words the intuition of avoiding time units in Discrete Time Fourier Series of a signal periodic in time.

Time T0 should be devided by N and the time intervals between N0 N1 and N2 should all be T/N. But, that T0 could be more than 2pi! Could you please explain the 2pi/N? What are the units of 2pi/N? The units of 2pi/T in contniuous time were 1/second. N is an integer, it shall have no units, it shall be the numbers of parts we devide some time T0 by. No? Yes? I'm confused. Could you help? Thanks ahead.

A DFT (Discrete Fourier Transform) is applied to an interval of a signal with uniformly sampled values on the underlying dimension, which is often time.  $N$ is a variable that is conventionally used to define the number of sample points.  The set of points in the interval is called the DFT frame.  Thus $N$ has units of $\frac{samples}{frame}$.

The DFT does not care about the sampling rate.  The sampling rate is used to convert frequencies in the DFT of cycles per frame to cycles per the underlying dimension's unit.

$2\pi$ has units of $\frac{radians}{cycle}$.  The $kth$ bin in a DFT represents a frequency of $k$ $\frac{cycles}{frame}$.

Therefore, the expression

$$ \frac{2\pi}{N}k $$ has units of

$$ \frac{\frac{radians}{cycle}}{\frac{samples}{frame}} \cdot \frac{cycles}{frame} = \frac{radians}{sample} $$

Now, look at the unnormalized definition of the DFT:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i \left(\frac{2\pi}{N}k\right)n } $$

The units of $n$ is samples, so the units of the exponent is:

$$ \frac{radians}{sample} \cdot sample = radians $$

Which is what it should be.

$F_s$ is the usual symbol for the sampling rate, also called the sampling frequency and sometimes the sampling density.  With seconds as the unit of time, it has units of $\frac{samples}{second}$.  I do not like the use of $Hz$ as a unit here, it causes confusion with $Hz$ being a unit of the frequency domain.

$$ T = \frac{N}{F_s} $$

$$ \frac{seconds}{frame} = \frac{ \frac{samples}{frame} }{ \frac{samples}{second} } $$

The frequency of the $k$th bin.

$$ f_k = \frac{k}{T} = k \frac{F_s}{N} = \frac{k}{N} F_s $$

$$
\begin{aligned}
\frac{cycles}{second} &=  \frac{ \frac{cycles}{frame} }{  \frac{seconds}{frame} }  \\
&= \frac{cycles}{frame}   \frac{ \frac{samples}{second} }{ \frac{samples}{frame} }  \\
&= \frac{ \frac{cycles}{frame} }{ \frac{samples}{frame} } \frac{samples}{second} \\
\end{aligned}
$$

Frequency spacing on the DFT is also known as the bin width.  The difference in frequency between any two adjacent bins (except DC and Nyquist) is going to be $\frac{F_s}{N}$ in units of $\frac{cycles}{second}$ also called $Hz$.


Relation between constellation SNR and baseband SNR Let's say I have an I/Q modulation, ie. QPSK.

I can simulate an the effects of an AWGN channel by adding to my QPSK symbols the realizations of a complex normal.

$\mathcal{CN}(0,\frac{1}{\mathrm{SNR}}) = \mathcal{N}(0,\frac{1}{2 \mathrm{SNR}}) + j \quad \mathcal{N}(0,\frac{1}{2 \mathrm{SNR}})$

with
$\mathrm{SNR} = 10^\frac{\mathrm{SNR_{db}}}{10}$

I should also be able to generate a baseband waveform with my QPSK symbols, ie using a 15kHz carrier, add gaussian noise to my baseband signal and then make a soft demodulation.

Using the same variance on both channel models, I observe a lower SNR on the soft decoded symbols than if I directly apply the noise realizations over the I/Q constellation.

Why are both different? Shouldn't I expect even better SNR when applying AWGN to the baseband waveform, as the noise isn't band limited?

What's the relation between the noise variance applied to the constellation and the noise variance applied to the baseband signal?

Edit:

Defining the quadrature carriers as:

$\phi_1=\sqrt{\frac{2}{T}}\cos(2 \pi f_c t),\quad 0 \le t \le T$

$\phi_2=\sqrt{\frac{2}{T}}\sin(2 \pi f_c t),\quad 0 \le t \le T$

With $T$ the symbol interval and $f_c$ the  carrier frequency.

The baseband waveform for an I/Q symbol $s_i=[s_{i1} s_{i2}]$ is:

$s(t)=s_{i1} \phi_1 + s_{i2} \phi_2$

The received waveform at the output of the AWGN channel is

$x(t)=s(t)+ w(t)$

With $w(t)$ a continuous-time wide-sense-stationary white noise process.

Assuming perfect coherence, the demodulation performed by a correlator is:

$\hat{s}_{i1}=\int_0^T x(t) \phi_1 dt = s_{i1}+\mathrm{w_{i1}}$

$\hat{s}_{i2}=\int_0^T x(t) \phi_2 dt = s_{i2}+\mathrm{w_{i2}}$

I'm looking for the relation between $w(t)$ and the I/Q noise $[\mathrm{w_{i1}},  \mathrm{w_{i2}}]$.

My waveform simulation will sample $x(t)$.

I think that the question by OP xvan has several misconceptions in it and these propagate into the answer that he has provided for his own question.

Define
\begin{align}\phi_1(t) &= \sqrt{\frac 2T}\cos(2\pi f_ct),~ 0 \leq t \leq T,\\
\phi_2(t) &= -\sqrt{\frac 2T}\sin(2\pi f_ct),~ 0 \leq t \leq T,
\end{align}
where $f_c$ is the carrier frequency are two orthogonal unit energy RF signals and thus $$s(t) = s_{I}\phi_1(t)+s_{Q}\phi_1(t)\tag{1}$$ is a RF QPSK signal with rectangular baseband pulses (unless, of course, $s_I$ and $s_Q$ are themselves some suitably chosen pulse shapes (e.g. raised-cosine) whose dependence on $t$ has been suppressed by the OP.  (Note that I have inverted $\phi_2(t)$ for reasons that will become apparent soon).  $s(t)$ as defined in $(1)$ is a real-valued signal and one does not add complex Gaussian noise to it; one must add (real-valued continuous-time) white Gaussian noise to $s(t)$ to get $x(t)$.  It is also helpful to assume that $f_cT$ is an integer, and if not, at least that $f_T \gg 1$.
The complex-valued baseband symbol $[s_1, s_2]$ (or $\mathbf s = s_1+js_2$ if one likes) takes on values $\left[\pm \sqrt{\mathcal E_b}, \pm \sqrt{\mathcal E_b}\right]$ where $\mathcal E_b$ is the signal energy per bit. These are the constellation points in the complex baseband equivalent of the RF signal $(1)$. Thus,
\begin{align}
s(t) &= \pm\sqrt{\frac{2\mathcal E_b}{T}}\cos(2\pi f_ct)
\mp \sqrt{\frac{2\mathcal E_b}{T}}\cos(2\pi f_ct)\\
&= \Re\left(\mathbf s \sqrt{\frac 2T}\exp(j2\pi f_c t)\right) \tag{2}
\end{align}
showing that the QPSK signal is the real part of the result of (amplitude-)modulating the complex baseband signal $\mathbf s$ onto the complex carrier $\sqrt{\frac 2T}\exp(j2\pi f_c t)$.
Note also that if I had stuck with the OP's definition of $\phi_2(t)$, $(2)$ would have been
$$s(t) = \Re\left(\mathbf s \sqrt{\frac 2T}\exp(-j2\pi f_c t)\right)$$
which would have been a disaster since I don't believe in negative frequencies on Mondays, Wednesdays and Fridays.
To the signal in $(1)$, we add continuous-time white Gaussian noise $W(t)$ which is a hypothetical zero-mean process with autocorrelation function $R_W(t) = \frac{N_0}{2}\delta(t)$ where $\delta(t)$ is the Dirac delta, that is, the power spectral density (PSD) $S_W(f)$ has value $\frac{N_0}{2}$ for all $f, -\infty < f < \infty$. Hypothetical because the noise has infinite power. More practically, the result of passing $W(t)$ through a BIBO filter with transfer function $H(f)$ results in a zero-mean WSS Gaussian process with PSD $\frac{N_0}{2}|H(f)|^2$ while
$\int_{0}^{T} \phi_i(t)W(t) \,\mathrm dt$ is a zero-mean Gaussian random variable with variance
\begin{align}
\sigma^2 &= E\left[\int_{0}^{T} \phi_i(t)W(t) \,\mathrm dt\int_{0}^{T} \phi_i(u)W(u) \,\mathrm du\right]\\
&= E\left[\int_{0}^{T} \int_{0}^{T}\phi_i(t)\phi_i(u)W(t) W(u) \,\mathrm dt \,\mathrm du\right]\\
&= \int_{0}^{T} \int_{0}^{T}\phi_i(t)\phi_i(u)E[W(t)W(u)] \,\mathrm dt \,\mathrm du\\
&= \int_{0}^{T} \int_{0}^{T}\phi_i(t)\phi_i(u)\frac{N_0}{2}\delta(t-u) \,\mathrm dt \,\mathrm du\\
&= \int_{0}^{T} \phi_i(u)\left[\int_{0}^{T}\phi_i(t)\frac{N_0}{2}\delta(t-u) \,\mathrm dt\right] \,\mathrm du\\
&= \frac{N_0}{2}\int_{0}^{T} \phi_i^2(u)\,\mathrm du\\
&= \frac{N_0}{2}
\end{align}
since $\phi_1(t)$ and $\phi_2(t)$ are unit-energy signals. A calculation similar to the one above shows that $\int_{0}^{T} \phi_1(t)W(t) \,\mathrm dt$ and $\int_{0}^{T} \phi_2(t)W(t) \,\mathrm dt$ are uncorrelated random variables and hence independent random variables since they are jointly Gaussian random variables.  In short, for $i=1,2$,
\begin{align}
\hat{s}_i &= \int_{0}^{T} x(t)\phi_i(t) \,\mathrm dt\\
&= \int_{0}^{T} [s(t) + W(t)] \phi_i(t) \,\mathrm dt\\
&= s_i + W_i
\end{align}
where $s_i$ has value $\pm \sqrt{\mathcal E_b}$ while $W_1$ and $W_2$ are independent Gaussian random variables with identical variance $\frac{N_0}{2}$.

Now, the OP wants to sample $x(t)$ a large number of times (perhaps multiple samples per period of the carrier) and use the samples presumably to simulate the operations of the correlators in the receivers to make sure that the correlators are doing what the math says that they are doing. This is a whole different kettle of fish since one cannot sample white noise in any meaningful sense and what the variance of the samples turns out to be depends on the details of the actual sampling device and not any model of the device that one might devise in MATLAB. If the "sampler" effectively grabs a very very short segment of $x(t)$ and returns the average value of $x(t)$ over that segment, then model that, and not that the sampler puts out the instantaneous value of $x(t)$.

But, let's take a more realistic/charitable view and assume that the noise process is band-limited white noise whose PSD has value $\frac{N_0}{2}$ only for $|f| \in \left[f_c-\frac 1T, f_c+\frac 1T\right]$ which is the nominal null-to-null passband for the QPSK signal. The noise autocorrelation function $R_W(t)$ is thus $\frac{N_0}{2}\operatorname{sinc}\left(\frac tT\right)\cos(2\pi f_ct)$ from which the most important take-away is that samples from $W(t)$ that are spaced $T$ seconds apart are independent, but samples spaced less than $\frac{1}{2f_c}$ seconds apart are very definitely correlated.  The samples are indeed Gaussian random variables and thus the Riemann sums that be used are indeed Gaussian random variables, but we can't just add up the variances of the variables in the sum; they are correlated random variables and so we must take the covariances into account as well. This messes up the calculations and getting it all straightened out is beyond what one should expect in an answer in a forum such as dsp.SE. What the OP has written in his own answer is not quite right, but is not correctable in a few words.


A State space model for discrete Sine wave Using kalma filter I'm looking to apply an optimal LQR filter to a discrete signal of the form
$$
x[n]=A\sin(\omega_0n+\phi)+v[n]
$$
The amplitude $A$ and the phase $\phi$ are unknown variables I want to estimate using the filter, and $v[n]$ is an uncorrelated noise signal of variance $\sigma_v^2$.

I don't know how to build a state model to generate this sine wave, and proceed from there.

You could use a nonlinear Kalman filter, such as the extended Kalman filter (EKF), and track the phase and frequency as your state variables.

In this case, your Kalman filter is essentially acting like a phase-locked loop (PLL).

Example reference


How to estimate the autocorrelation from nonuniformly spaced data Assume a continues-time random process $X(t)$ sampled nonuniformely in time to acquire discrete signal $x[n]$. The sampling times are known but the autocorrelation is not. Is there an accurate approach to estimate the autocorrelation function? This problem is challenging since at some lags there is no sample to average.

Any suggestion even in a special case, e.g. a simple autocorrelation function with an exponential decay (continuous-time AR process of order one) would be appreciated.

Edit

I have already read some papers, so I would like either an analytical answer or a peace of algorithm/code to address this problem for instance for an AR(1) process assuming non-uniformly spaced samples(e.g. in uniformly distributed random times).

As suggested by  Marcus Müller, interpolation in the time domain could be a solution. I never had to perform such a task, and the outcomes may depend in the nonuniformity of your sampling. I propose a use of the Wiener-Khinchin theorem:

  The Fourier transform of the autocorrelation function is the power
  spectrum, or equivalently, the autocorrelation is the inverse Fourier
  transform of the power spectrum

You may use for instance:

a non-uniform DFT (for instance NUFFT),
variants of the Lomb-Scargle periodogram (aka Least-squares spectral analysis),

and get a regularly sampled power spectrum, which I hope you could inverse Fourier back to an estimate of the autocorrelation. As I have no practical experience, if you intend to try this, I'd really appreciate you could share the outcomes, good or bad, along with the specificities of your signal:

what kind of nonuniformity (additive random, jitter, lacunary)?
type of signal and noise?

Additional references:

Power Spectral Density of Unevenly Sampled Data by Least-Square Analysis: Performance and Application to Heart Rate Signals, Laguna et al., 1998
Spectral estimates from nonuniform samples, L. Shaw, 1971


Guitar pitch detection I try to develop a guitar game in android platform.

And I need to do the real-time pitch detection of guitar chord.

I will get the input from the microphone, and then analyze the input (the input playing which kind of guitar string/chord).

I find two kinds of method that I can use, one is YIN, another one is FFT.

Which method can get better performance and exact result?

Your question is ambiguous.  Chord, pitch and frequency detection/estimation are very different things.  Bare FFT magnitudes are really bad at guitar pitch modeling.  YIN is a better monophonic pitch estimator.  Polyphonic chord estimation is vastly more difficult, and may or may not even be possible, depending on your input data or time/pitch accuracy requirements (exactly what is your "exact result"?) plus any computational limitations ("real time" latency requirements, CPU performance/cores/GPUs available, etc.)


Binomial Approximation of Gaussian Distribution It is said that we may use the binomial coefficients ( a layer from Pascal's triangle) to approximate the 1-D Gaussian kernel with certain $\sigma$, where $\frac{n}{4} = \sigma$ and $n$ is the index of the layer.

This works really nice when we want to generate a gaussian quickly. But is there a way to prove this? Or just coincidence?

And how do we decide the variance of certain layer of Pascal's triangle? Say the third layer is [1 3 3 1] (consider [1 1] as the first layer), and it is used to approximate gaussian kernel of  $\sigma = \sqrt{\frac{3}{4}} = \frac{\sqrt{3}}{2}$. But how do we prove this?

Let $X \sim B(n, p)$. At the third layer $\begin{bmatrix}1 & 3 & 3 & 1\end{bmatrix}$ of the Pascal's triangle, we have $X \sim B(3, \frac{1}{2})$. Which means $P(X = k) = {3 \choose k} (\frac{1}{2})^k (\frac{1}{2})^{3-k} = \frac{{3 \choose k}}{2^3}$. The variance is given by $\sigma^2 = n p q = 3 \cdot \frac{1}{2} \frac{1}{2} = \frac{3}{4}$ and the std by $\sigma = \sqrt{\frac{3}{4}}$. You can derive the variance by looking at $n$ independent Bernoulli trials.

In general, a 1d binomial filter/kernel of size $n$ with $p = \frac{1}{2}$ is given by $\begin{bmatrix}P(X = 0) & P(X = 1) & \cdots & P(X = n-1)\end{bmatrix}$. As $n \to \infty$, this kernel will approximate the normal distribution with mean $\frac{n}{2}$ and variance $\sigma^2 = \frac{n}{4}$. If you want a different variance/std, you should change $p \neq 0.5$.

As somebody already wrote in the comments, you should look at the De Moivre-Laplace theorem to prove this. However, this is only a special case of the more general central limit theorem. By showing that the binomial distribution is asymptotically equivalent to the normal distribution as $n\to\infty$, I would argue that we have shown already both directions.


How can it be that there is a series of integrals in Fourier Series if it’s a projection on a continuous basis? If the process of finding Fourier coefficient is finding the projection of a signal on a member from an orthonormal basis, basis which is continuous in frequency. How can it be that Fourier coefficient is an integral and signal representation, which is the sum of these coefficients multiplied by the projections, is a series? The coefficient basis members are continuous, how can they be discrete in the signal representation while being continuous in the coefficient representation? The logic says: we took projection in continuous space, then we shall continue in continuous space while we multiply by that projection. Why do we have a mix of integral and a series then?

  How can it be that there is a “series of integrals” in Fourier Series if Fourier is just an operation of projecting a continuous basis and multiplying the projection on the (continuous) basis members once again?

The short answer is that the sinusoidal components are only orthogonal if they are integer multiples of each other. Any other choice of frequency would not result in $<s1, s2> = 0$, which can be easily confirmed by comparing such cases (integer and non-integer related sinusoids). Thus the choice of basis signals is discrete in units of their frequency (the frequency domain), but the time domain waveforms that each of these discrete frequency values represent is continuous.

The frequency domain signal can still be a continuous signal if we allow zero-value signals, as it is just the non-zero values that would necessarily be discrete. That said, we can use discrete signals to describe the waveform in frequency without any loss of information.

The intuitive view starts from understanding that the Fourier Series Expansion states that we can represent a continuous time single-valued function (with no discontinuities) of duration $T$ as a series of sinusoidal components (sines and cosines) each with a frequency of $1/T$ and integer multiples.

Consider such an arbitrary function below:

Mathematically we need not limit time to $T$ in the the reconstruction; time can extend to $±\infty$ as long as we allow the signal in time to be periodic over period T, we will achieve the same result as derived in the Fourier Series Expansion. Since the waveform repeats exactly in time over every interval $T$, we need only express the result over period $T$ but extending it as such provides further insight into the Fourier Series Expansion, such as the inability to converge if there was a discontinuity between the end and the beginning of the interval. Here in the graphic I purposely made that boundary continuous so that the repeated waveform has no discontinuity at each transition, thus this is an example that the expansion will converge everywhere from $0$ to $T$ where in general such a continuity would not exist and you will typically see severe ringing at the start and end of the reconstructed waveform.

Thus the waveform components in the decomposition MUST also be cyclical over the interval $T$, otherwise they cannot possibly converge to the same waveform over the next period of $T$, meaning the signal cannot be periodic over $T$.  (As demonstrated with the possible first two components of the illustrated waveform as shown below.)

It is periodicity specifically in the time domain that is associated with discrete values of frequency in the frequency domain since the other frequencies if they did exist would destroy the periodicity over that interval: In order to be discrete in frequency, periodicity must exist in the time domain. Consider a sine-wave itself which is periodic it time and exists at one frequency. As demonstrated with the expansion an arbitrary periodic waveform can have multiple frequency components but they all must be integer multiples of the fundamental component at frequency $1/T$. I like to explain it like this: If I could call out your name exactly once per second, and exactly the same way each time--the resulting frequencies of my voice would only exist at 1 Hz and integer multiples of 1 Hz, all the way out to the 4KHz or so where spectral content exists.

This also works in reverse: periodicity in frequency is associated with discrete values in time (sampling). When we sample a time domain signal, we end up with multiple copies of the signal spectrum in frequency, spaced by the sampling rate. Same thing- in this case the frequency spectrum is continuous and periodic while the time domain signal is discrete. The DFT is an example that is discrete both in time and frequency, and therefore there is an implied periodicity both in time and in frequency!


Segmenting and labeling chessboards through computer vision With almost no experience/knowledge in computer vision, I am currently trying to segment and label the separate squares on a chessboard (i.e. a1 ... h8). My goal is to use the information to physically play chess with an AI using a webcam and a robotic arm.

I am writing the vision program in C# (using Emgu CV) for simpler programming and GUI implementation. As my first step, I attempted to detect squares by following this tutorial.

As we can see, not all of the squares were detected and I believe some of the squares were detected more than once (as 91 boxes were created).

Questions

What are the factors to consider in determining the suitable values for Canny threshold and threshold linking? i.e gray.Canny(cannyThreshold, cannyThresholdLinking)
The squares in the original image appear as straight lines to us but to the computer, it is slightly tilted. What is the cause of this? Is it because of aliasing at the edges?
Why are most of the squares not detected? Could it be the lighting?
What are the proper techniques/methods I should investigate to solve this problem? I am only familiar with the keywords: Canny, Hough (for circles I believe) and Sobel at the moment.

Code

Image<Bgr, Byte> image = new Image<Bgr, Byte>(Emgu_CV.Properties.Resources.Chessboard).Resize(0.15, INTER.CV_INTER_LINEAR);
Image<Gray, Byte> gray = image.Convert<Gray, Byte>(); //.PyrDown().PyrUp() worsens detection

imageBox.Image = image;

double cannyThreshold = 180.0;
double cannyThresholdLinking = 150.0;

Image<Gray, Byte> cannyEdges = gray.Canny(cannyThreshold, cannyThresholdLinking);

List<MCvBox2D> boxList = new List<MCvBox2D>();

using (MemStorage storage = new MemStorage())
    for (Contour<Point> contours = cannyEdges.FindContours(); contours != null; contours = contours.HNext)
    {
        Contour<Point> currentContour = contours.ApproxPoly(contours.Perimeter * 0.05, storage); // don't quite understand the purpose of *0.05
        if (contours.Area > 250 && contours.Area < 10000)
        {
            if (currentContour.Total == 4)
            {
                bool isSquare = true;
                Point[] pts = currentContour.ToArray();
                LineSegment2D[] edges = PointCollection.PolyLine(pts, true);

                for (int i = 0; i < edges.Length; i++)
                {
                    double angle = Math.Abs(edges[(i + 1) % edges.Length].GetExteriorAngleDegree(edges[i]));

                    if (angle < 80 || angle > 95)
                    {
                         isSquare = false;
                         break;
                    }

                    if (isSquare) boxList.Add(currentContour.GetMinAreaRect());
                }
            }
        }
    }

Image<Bgr, Byte> squareImage = image.CopyBlank();
foreach (MCvBox2D box in boxList)
    squareImage.Draw(box, new Bgr(Color.Yellow), 1);
outputBox.Image = squareImage;

I have once done something similar. You should give a try to OpenCV's checkerboard detection. There are also works on Sudoku square detection, which might give you a hint.

If these don't work for you, then I would suggest the following: Forget about detection the checkerboard. Initially, assume that every piece is in place and take the snapshot of the board (you can even manually place an initial board, doesn't matter). Then, whenever a piece is moved, grab a new image and apply a differencing of two images (this is somewhat like a background subtraction). Note that this way, you easily get the location of the moved piece. And since you know the piece and what was in the previous location beforehand, you get the complete information of the move.

Use a downward looking camera, just like you showed.


GNU Radio | Packet Encoder/Decoder | Unwanted bytes added to the transmitted stream (A follow up question can be found here: link)

In an attempt to understand how the packet encoding blocks work I tried to mimic the experiment performed by Aaron, published in this post, but reduced the flowgraph to the components that interested me the most, namely the Packet Encoder, the Packed to Unpacked and the Packet Decoder. I have also used a text file in_strings.txt instead of the cat_in.png.

Follows the flowgraph I have  implemented:

My expectation: I was hoping that out.txt (File Sink) would be equal to in_strings.txt (File Source).

What really happened: out.txt is very similar to in_strings.txt but it contains additional unwanted x characters that show up after 65535 expected characters.

Follows a preview of the diff side by side between the input and the output:

Follows the in_strings.txt file: https://pastebin.com/NWa7NEww

I wonder why those x's show up in the output?
Thanks in advance!

Packet Encoder and Decoder are broken; they drop data.

That's why they are in the deprecated category (for years now!). We've removed them, because as a project, GNU Radio has not been able to fix them (and also, they were terrible from an architecture point of view).

So there's exactly one solution: don't use packet encoder / decoder.


Correlation in DSBSC modulation Suppose we have a DSBSC signal given by :

  $$s(t) = 10 \cos (2\pi f t + \phi ) \cdot m(t) $$

  where carrier is random signal with phase ($\phi$) uniformly
  distributed in $[0,2\pi]$.

  If the power spectral density of the message signal $m(t)$ is :

  Q.How can I express the power spectral density of $s(t)$ in terms of given PSD of $m(t)$?

My attempt:

By Wiener–Khinchin relations,

Autocorrelation of $m(t) = 3 \mathrm{sinc}^2(3000t)$

Now to get the autocorrelation of $s(t)$,

$$R_{ss} (T) = E[ x(t)x^*(t+T)\cdot m(t) m^*(t+T) ]$$

where $x(t)$ is the carrier.

If I assume they are independent, then :

$$R_{ss}(T) = R_{mm}(T) R_{xx}(T)$$

which gives me

$$R_{ss}(T) = 150\cos(2\pi T)\cdot\mathrm{sinc}^2(3000T)$$

The random phase $\phi$ makes the carrier signal $c(t)=A\cos(2\pi f_0+\phi)$ wide-sense stationary with the autocorrelation function

$$R_c(\tau)=\frac{A^2}{2}\cos(2\pi f_0\tau)\tag{1}$$

The Fourier transform of $(1)$ is the power spectrum of the carrier:

$$S_c(f)=\frac{A^2}{4}\left[\delta(f-f_0)+\delta(f+f_0)\right]\tag{2}$$

where $\delta(f)$ is the Dirac impulse.

Since the carrier signal $c(t)$ and the message signal $m(t)$ are assumed to be uncorrelated, the autocorrelation function of the modulated signal $s(t)=c(t)m(t)$ equals the product of the autocorrelation functions of $c(t)$ and $m(t)$:

$$R_s(\tau)=R_c(\tau)R_m(\tau)\tag{3}$$

From $(3)$ it follows that the power spectrum $S_s(f)$ of $s(t)$ is the convolution of the power spectra of $c(t)$ and $m(t)$:

$$S_s(f)=S_c(f)\star S_m(f)=\frac{A^2}{4}\left[S_m(f-f_0)+S_m(f+f_0)\right]\tag{4}$$


In a periodic square wave, how to change in a particular condition? The periodic square wave,
$$x(t) = 1, |t| < T_1 \\
x(t) = 0, T_1 < |t| < T/2$$

What I wonder is ..

(1) If $T_1$ is constant and $T \to \infty$, is x(t) a constant function?

(2) If $T_1$ is constant and $T \to 0$, is this non-sense?

(3) If $T$ is constant and $T_1 \to 2/T$, is x(t) a constant function?

(4) If $T$ is constant and $T_1 \to 0$, is x(t) a impulse function?

Actually, I am a undergraduate and it is my first class about signal. I don't know impulse function correctly. Could you explain impulse function in terms of (4) ?

The above questions are my intutitions. I don't know answers.

If $T_1$ remains constant and $T\rightarrow\infty$ then you get a single square pulse of length $2T_1$. There won't be another square pulse because the period is infinite.

$T_1$ constant and $T\rightarrow 0$ is indeed nonsensical because the definition of the signal assumes $T_1 < T/2$.

If $T$ is constant, you can get a constant function if $T_1=T/2$ (and not $2/T$).

If $T$ is constant and $T_1\rightarrow 0$ your function becomes 0. You will only get an impulse if you keep the area under the rectangle constant, i.e. if its height increases as its width decreases.


Is my answer incomplete? Checking the stability of a system Yesterday, during my exam, I had the following exercise:

  Given $$H(s) = \frac{1}{s^2+2s+4}$$ check if it's stable.

which was supposed to be the hardest (since it was the last one). From my knowledge, I quickly found out the poles, checked the real part of the pole (which was $-1$) and said that the system is stable. Is my argument incomplete? Everyone in the classroom seemed to write a lot more at that final exercise.

Did you mention the region of convergence?  The system can be stable, but anti-causal depending on the ROC.


MATLAB Adaptive Filters In the audio domain, I am currently attempting to use MATLAB to distil:

$\textrm{signal}$ from $\textrm{noise + signal}$
$\textrm{noise}$ from $\textrm{noise + signal}$
using two adaptive filters $\rightarrow$ two results.

I get the first answer quite effectively using:

hFDAF1 = adaptfilt.fdaf(AdaptFiltLength,StepSize,Leakage,Delta,Lambda)

[Errors1, Adapt_Out_Audio] = filter(hFDAF1, NoisePlusAudio, Audio)

where NoisePlusAudio is a microphone in the room, and Audio is music playing through speakers in the room.

How do I get the adaptive filter to remove the Audio from NoisePlusAudio, just leaving the background noise of the room, much like an echo canceller?
The following does not work, it just lets all the sound through.

[Errors2, Adapt_Out_Noise] = filter(hFDAF2, Audio, NoisePlusAudio)

I'd love to know the answer ....

Reason for it might be the continuous adaptation, even in non speech regions. So, you can try Voice Activity Detector (VAD) for knowing the exact regions for adaptation of filter, and processing will be frame wise (20 ms).


Why MATLAB fft cos makes imaginary parts? The cos fourier transform has no imaginary parts, but in this code it has imaginary parts that little big.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(7);

fs=128;

t=-4*pi:1/fs:4*pi-1/fs;

x=1*cos(0.5*t);

X=fft(x);

N=length(x);

n=0:N-1;

f=n*fs/N;

phs=angle(fftshift(X));

subplot(2,1,1); plot(f,abs(fftshift(X))/N,'LineWidth',1.5);

title('FFT'); ylabel('Magnitude'); xlabel('Frequency(Hz)');

subplot(2,1,2); plot(f,phs*180/pi,'-o'); xlabel('Frequency(Hz)');

ylabel('Phase(Degree)');

figure(8);

Re=real(X); Im=imag(X);

subplot(2,1,1); stem(f,Re); ylabel('Re(G(f))'); xlabel('Frequency(Hz)');

subplot(2,1,2); stem(f,Im); ylabel('Im(G(f))'); xlabel('Frequency(Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

If I make time range a multiple of 0.5, the results of imaginary parts has very little value that is e-14~15.

Let me know why this problem happened...

Those $10^{-15}$ terms after the FFT are due to numerical roundoff errors... Apart from that, you should properly frame the cosine into a fully periodic length, as it seems you did, then you can expect the imaginary part to be as close to zero as the numerical format allows...
The following line of OCTAVE / MATLAB code shows the effect of frame length (aka aperture size) on the FFT computation:
L = 16   ; period of cosine
N = 128  ; frame length of sample block (integer multiple of L)
figure,stem(imag(fft(cos(2*pi*[0:N-1]./L),N))); title('imaginary part of the FFT');

Now let us change the period of the cosine so that the frame does not include an exact number of periods:
L = 15  ; period of cosine
N = 128  ; frame length of sample block (non-integer multiple of L)
figure,stem(imag(fft(cos(2*pi*[0:N-1]./L),N))); title('imaginary part of the FFT');

Now if your result display more deviations than those found on the first plot, then your frame length is probably the problem. If you, on the other hand, observe very close to zero (probably randomly distributed) values as in the first plot, then that's probably due to the numerical roundoff errors involved in the FFT computation.


Impact of modulation index in FM If I have the spectrum of a FM modulation, does the frequency of the carrier disappear when the modulation index is too big?
Does the amplitude of any side band depend from the modulation index?
And does the total number of side bands depend on the modulation
index?

Assuming an FM modulation of a sine wave of frequency $f_m$. Theoretically, the modulated signal is an infinite sum of cosines at the frequencies $f_c + nf_m$ ($n\in \mathbb{Z}$) which the amplitudes depend on the coefficients $J_n(\beta)$ (of course in addition to the amplitude of the carrier wave), where $\beta$ is the modulation index and $f_c$ is the carrier frequency. At this stage, the answer to the question Q.2. is yes. For the other items:

Q.1. The PSD of the modulated signal in $f_c$ is determined for $n=0$, then the PSD at that frequency depends on the coefficient $J_0(\beta)$. Knowing that $J_n(\beta)$ trends to zero as $\beta$ trends to infinity, so the PSD of the modulated signal in $f_c$ trends to zero as the modulation index trends to infinity.
Q. 3. Theoretically the bandwidth of an FM modulated signal is infinite. But when considering Carson bandwidth rule ($B_C$), the answer to your question is yes since this rule considers the modulation index. $B_C=2(1+\beta)W$ where $W$ is the baseband bandwidth of the modulating signal (assuming a sine wave, $W=f_m$).


From Mathematics to DSP I'm currently working on multiple mini pet-projects in my spare time, grounding myself in the basics of audio signal processing, and trying to develop algorithms for polyphonic pitch detection. As I'm having a lot of joy (even in the failure) in doing so, I am considering moving into a career in signal processing (I'm aware this is a vast area!) My background is in mathematics (I have a  bachelors degree in applied mathematics, so the theory and practice excites me) and I have no grounding whatsoever in electrical or electronic engineering.

I'm therefore wondering the following:

Assuming I can develop a project portfolio demonstrating basic competence, what specific careers would be open to someone of my background in something highly geared to signal processing? I would be particularly interested in sound engineering, but would also be interested in compression algorithms and coding theory.
Is a supplementary MSc. necessary to get my foot in the door to job positions? I am open to both industry and academic positions.
Would an MSc. in something involving machine learning and neural networks be a.) relevant, and b.) worth the dough? I am considering getting out of industry for a year to work more on my projects in the context of such a degree programme. I am getting the impression that combined with a knowledge of signal processing, some knowledge of A.I. and self-learning systems would go a long way.

Thanks in advance for your answers.

In general, Maths to audio signal processing is a very natural progression. I seem to recall my Master's and PhD supervisors both did their undergraduate degrees in Maths.

To your questions:

This is the hardest question to answer. Depends what's in your "portfolio of basic competence". In addition to the understanding of signal processing theory that you're developing you need to show you're able to implement these algorithms probably either in hardware or software. Given little EE grounding, I'd suggest software. Look at developing your software engineering skills. They will serve you well in either academia or industry.
MSc will definitely help but is not absolutely necessary. The types of algorithms you are describing could conceivably be the subject of an MSc level project. I work at a company that makes hardware for audio signal processing. The statistic they gave me is that around 60% of staff have a Master's or higher degree.
Yes, combining signal processing and machine learning is a very wise idea if you want to work towards the research end of the field. Several problems in audio signal processing - e.g. audio source separation and speech recognition- are increasingly finding their state-of-the-art solutions come from machine learning.


Grayscale conversion before applying sobel gradient operator I am currently replacing a hand-made image processing library by OpenCV. The programmer who worked on this library had implemented the sobel gradient function. This implementation is right except that he first converts his image to grayscale with the formula thereafter (the arithmetic mean):

Y = 0.333 * R + 0.333 * G + 0.333 * B

In OpenCV, I see that they use the formula:

Y = 0.299 * R + 0.587 * G + 0.114 * B

To convert an image to grayscale.

Is the hand-made implementation wrong or can it be motivated by some purpose. Should I always use the OpenCV conversion version?

The values used by OpenCV are values used by Luma coding in video systems. Your original code just used an equal weight for each channel.

This may or may not be important for your application (probably not really).


Compressive Sensing - Location of Non Zero Elements I am new in compressive sensing and I would like to know if there already exists a deterministic or probabilistic approach to estimate the locations of the $k$ non-zero entries in the $k$-sparse vectore $x$.

Actually some greedy algorithms like MP, exactly do that. They iteratively find signal supports (location of non-zero elements) using measurements and measurement matrix, and reconstruct original signal.
However, in case you are asking for a method to find signal supports from measurements all at once, I think there is no such thing (may be someone discover it someday!).
If there was a such a thing, the computational cost of the CS recovery algorithms wouldn't be so high.

Fore further information on greedy recovery algorithms, take a look at these:

http://www.mathworks.com/help/wavelet/ug/matching-pursuit-algorithms.html

https://en.wikipedia.org/wiki/Matching_pursuit


Spectral Content of two waveforms I have found this interesting question here:

Which of these two waveforms has a richer spectral content, and why?

Intuitively speaking, I think that the first one is richer, but I cant really justify it! Cheers for any input.

Your intuition is correct. The reason is that the rectangular pulse train has  discontinuities, whereas the rectified sine is continuous (but its derivative has discontinuities). The consequence of this is that the Fourier coefficients of the rectangular pulse train decay as $1/n$, whereas the coefficients of the rectified sine decay as $1/n^2$ (where $n$ is the index of the Fourier coefficients).


Problem with a code of FIR Lowpass filter I'm making an interpolator (filter) to remove aliasing but this doesn't work.

I see also some strange behaviour: x[i - j], but what if j > i? I simply followed FIR definition but still don't understand why this doesn't work.

const size_t M = 200;
double h[M + 1]; // h is filled with coeff

inline void process(double* x, double* y, size_t len)
{
    for (size_t i = 0; i < len; ++i)
    {
        y[i] = 0.;
        for (size_t j = 0; j < M + 1; ++j)
        {
            y[i] = y[i] + x[i - j] * h[j];
        }
    }
}

The problem is, as you've already pointed out, the choice of the summation limits. Obviously, the array index shouldn't become negative. Look at the definition of convolution:

$$y(n)=\sum_{k=-\infty}^{\infty}h(k)x(n-k)\tag{1}$$

Now assume, as in your case, that the filter length (i.e. the number of coefficients) equals $M+1$, and that the input signal is a finite length sequence in the range $0\le n<N$. Then you have two conditions on the summation index $k$:

$$0\le k \le M\quad\text{and}\quad 0\le n-k \le N-1$$

With these conditions the sum in (1) becomes

$$y(n)=\sum_{k=\max\{0,n-N+1\}}^{\min\{M,n\}}h(k)x(n-k)\tag{2}$$

It is the sum in (2) that you need to implement.

Also note that if you want to implement convolution, then the output signal will be longer than the input signal, so you would need to define two different variables len. But if you only want to compute the output signal as long as there is a non-zero input signal present, then your implementation is fine.

Another improvement would be to have a routine that just takes a scalar input value (the current input sample) and computes the corresponding output sample. In this way you don't need to know the length of the input signal and you don't need to store the complete input signal in an array because you just process each incoming sample directly. Of course you need an internal buffer storing the last $M$ input values. Your version is appropriate for offline-processing of not too long signals.


Two different definitions of cross-correlation I have come across two different definitions of the cross-correlation function:

$$
R_{XY}(\tau) = \int_{-\infty}^{+\infty} x(t) y(t+ \tau) dt
$$

and

$$
R_{XY}(\tau) = \int_{-\infty}^{+\infty} x(t) y(t- \tau) dt
$$

Are these two definitions equal?

Also is $ R_{XY}(\tau) = R_{YX}(\tau)$?

Start with:

$$ R_{XY}(\tau) = \int_{-\infty}^{+\infty} x(t) y(t+ \tau) dt $$

Then it follows:

$$ R_{YX}(\tau) = \int_{-\infty}^{+\infty} y(t) x(t+ \tau) dt $$

Let $s = t+ \tau $:

$$ R_{YX}(\tau) = \int_{-\infty}^{+\infty} y(s -\tau ) x(s ) ds $$

Reverse the product and substitute $t$ for $s$:

$$ R_{YX}(\tau) = \int_{-\infty}^{+\infty}  x(t ) y(t -\tau ) dt $$

Now it is pretty clear that

$$ R_{XY}(\tau) = R_{YX}(-\tau) $$

The question of odd or even is irrelevant because they are not the same function.

Note:  Your second function is assumed to be $ R_{YX}(\tau) $, otherwise your are simply giving contradictory definitions to the same function.