Auto Detection of Rotation Angle on Arbitrary Image with Orthogonal Features I have a task at hand where I need to detect angle of an image like the following sample (part of microchip photograph). The image does contain orthogonal features, but they could have different size, with different resolution/sharpness. Image will be slightly imperfect due to some optical distortion and aberrations. Sub-pixel angle detection accuracy is required (i.e. it should be well under <0.1° error, something like 0.01° would be tolerable). For reference, for this image optimal angle is around 32.19°.

Currently I've tried 2 approaches:
Both do a brute-force search for a local minimum with 2° step, then gradient descend down to 0.0001° step size.

Merit function is sum(pow(img(x+1)-img(x-1), 2) + pow(img(y+1)-img(y-1)) calculated across the image. When horizontal/vertical lines are aligned - there is less change in horizontal/vertical directions. Precision was about 0.2°.
Merit function is (max-min) over some stripe width/height of the image. This stripe is also looped across the image, and merit function is accumulated. This approach also focuses on smaller change of brightness when horizontal/vertical lines are aligned, but it can detect smaller changes over larger base (stripe width - which could be around 100 pixels wide). This gives better precision, up to 0.01° - but has a lot of parameters to tweak (stripe width/height for example is quite sensitive) which could be unreliable in real world.

Edge detection filter did not helped much.

My concern is very small change in merit function in both cases between worst and best angles (<2x difference).

Do you have any better suggestions on writing merit function for angle detection?

Update: Full size sample image is uploaded here (51 MiB)

After all processing it will end up looking like this.

If I understand your method 1 correctly, with it, if you used a circularly symmetrical region and did the rotation about the center of the region, you would eliminate the region's dependency on the rotation angle and get a more fair comparison by the merit function between different rotation angles. I will suggest a method that is essentially equivalent to that, but uses the full image and does not require repeated image rotation, and will include low-pass filtering to remove pixel grid anisotropy and for denoising.

Gradient of isotropically low-pass filtered image

First, let's calculate a local gradient vector at each pixel for the green color channel in the full size sample image.

I derived horizontal and vertical differentiation kernels by differentiating the continuous-space impulse response of an ideal low-pass filter  with a flat circular frequency response that removes the effect of the choice of image axes by ensuring that there is no different level of detail diagonally compared to horizontally or vertically, by sampling the resulting function, and by applying a rotated cosine window:

$$\begin{gather}h_x[x, y] = \begin{cases}0&\text{if }x = y = 0,\\-\displaystyle\frac{\omega_c^2\,x\,J_2\left(\omega_c\sqrt{x^2 + y^2}\right)}{2 \pi\,(x^2 + y^2)}&\text{otherwise,}\end{cases}\\
h_y[x, y] = \begin{cases}0&\text{if }x = y = 0,\\-\displaystyle\frac{\omega_c^2\,y\,J_2\left(\omega_c\sqrt{x^2 + y^2}\right)}{2 \pi\,(x^2 + y^2)}&\text{otherwise,}\end{cases}\end{gather}\tag{1}$$

where $J_2$ is a 2nd order Bessel function of the first kind, and $\omega_c$ is the cut-off frequency in radians. Python source (does not have the minus signs of Eq. 1):

import matplotlib.pyplot as plt
import scipy
import scipy.special
import numpy as np

def rotatedCosineWindow(N):  # N = horizontal size of the targeted kernel, also its vertical size, must be odd.
  return np.fromfunction(lambda y, x: np.maximum(np.cos(np.pi/2*np.sqrt(((x - (N - 1)/2)/((N - 1)/2 + 1))**2 + ((y - (N - 1)/2)/((N - 1)/2 + 1))**2)), 0), [N, N])

def circularLowpassKernelX(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda y, x: omega_c**2*(x - (N - 1)/2)*scipy.special.jv(2, omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = 0
  return kernel

def circularLowpassKernelY(omega_c, N):  # omega = cutoff frequency in radians (pi is max), N = horizontal size of the kernel, also its vertical size, must be odd.
  kernel = np.fromfunction(lambda y, x: omega_c**2*(y - (N - 1)/2)*scipy.special.jv(2, omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = 0
  return kernel

N = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
window = rotatedCosineWindow(N)

# Optional window function plot
#plt.imshow(window, vmin=-np.max(window), vmax=np.max(window), cmap='bwr')
#plt.colorbar()
#plt.show()

omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernelX = circularLowpassKernelX(omega_c, N)*window
kernelY = circularLowpassKernelY(omega_c, N)*window

# Optional kernel plot
#plt.imshow(kernelX, vmin=-np.max(kernelX), vmax=np.max(kernelX), cmap='bwr')
#plt.colorbar()
#plt.show()

Figure 1. 2-d rotated cosine window.

Figure 2. Windowed horizontal isotropic-low-pass differentiation kernels, for different cut-off frequency $\omega_c$ settings. Top: omega_c = np.pi, middle: omega_c = np.pi/4, bottom: omega_c = np.pi/16. The minus sign of Eq. 1 was left out. Vertical kernels look the same but have been rotated 90 degrees. A weighted sum of the horizontal and the vertical kernels, with weights $\cos(\phi)$ and $\sin(\phi)$, respectively, gives an analysis kernel of the same type for gradient angle $\phi$.

Differentiation of the impulse response does not affect the bandwidth, as can be seen by its 2-d fast Fourier transform (FFT), in Python:

# Optional FFT plot
absF = np.abs(np.fft.fftshift(np.fft.fft2(circularLowpassKernelX(np.pi, N)*window)))
plt.imshow(absF, vmin=0, vmax=np.max(absF), cmap='Greys', extent=[-np.pi, np.pi, -np.pi, np.pi])
plt.colorbar()
plt.show()

Figure 3. Magnitude of the 2-d FFT of $h_x$. In frequency domain, differentiation appears as multiplication of the flat circular pass band by $\omega_x$, and by a 90 degree phase shift which is not visible in the magnitude.

To do the convolution for the green channel and to collect a 2-d gradient vector histogram, for visual inspection, in Python:

import scipy.ndimage

img = plt.imread('sample.tif').astype(float)
X = scipy.ndimage.convolve(img[:,:,1], kernelX)[(N - 1)//2:-(N - 1)//2, (N - 1)//2:-(N - 1)//2]  # Green channel only
Y = scipy.ndimage.convolve(img[:,:,1], kernelY)[(N - 1)//2:-(N - 1)//2, (N - 1)//2:-(N - 1)//2]  # ...

# Optional 2-d histogram
#hist2d, xEdges, yEdges = np.histogram2d(X.flatten(), Y.flatten(), bins=199)
#plt.imshow(hist2d**(1/2.2), vmin=0, cmap='Greys')
#plt.show()
#plt.imsave('hist2d.png', plt.cm.Greys(plt.Normalize(vmin=0, vmax=hist2d.max()**(1/2.2))(hist2d**(1/2.2))))  # To save the histogram image
#plt.imsave('histkey.png', plt.cm.Greys(np.repeat([(np.arange(200)/199)**(1/2.2)], 16, 0)))

This also crops the data, discarding (N - 1)//2 pixels from each edge that were contaminated by the rectangular image boundary, before the histogram analysis.

$\pi$
$\frac{\pi}{2}$
$\frac{\pi}{4}$
$\frac{\pi}{8}$
$\frac{\pi}{16}$
$\frac{\pi}{32}$
$\frac{\pi}{64}$
—$0$Figure 4. 2-d histograms of gradient vectors, for different low-pass filter cutoff frequency $\omega_c$ settings. In order: first with N=41: omega_c = np.pi, omega_c = np.pi/2, omega_c = np.pi/4 (same as in the Python listing), omega_c = np.pi/8, omega_c = np.pi/16, then: N=81: omega_c = np.pi/32, N=161: omega_c = np.pi/64 . Denoising by low-pass filtering sharpens the circuit trace edge gradient orientations in the histogram.

Vector length weighted circular mean direction

There is the Yamartino method of finding the "average" wind direction from multiple wind vector samples in one pass through the samples. It is based on the mean of circular quantities, which is calculated as the shift of a cosine that is a sum of cosines each shifted by a circular quantity of period $2\pi$. We can use a vector length weighted version of the same method, but first we need to bunch together all the directions that are equal modulo $\pi/2$. We can do this by multiplying the angle of each gradient vector $[X_k,Y_k]$ by 4, using a complex number representation:

$$Z_k = \frac{(X_k + Y_k i)^4}{\sqrt{X_k^2 + Y_k^2}^3} = \frac{X_k^4 - 6X_k^2Y_k^2 + Y_k^4 + (4X_k^3Y_k - 4X_kY_k^3)i}{\sqrt{X_k^2 + Y_k^2}^3},\tag{2}$$

satisfying $|Z_k| = \sqrt{X_k^2 + Y_k^2}$ and by later interpreting that the phases of $Z_k$ from $-\pi$ to $\pi$ represent angles from $-\pi/4$ to $\pi/4$, by dividing the calculated circular mean phase by 4:

$$\phi = \frac{1}{4}\operatorname{atan2}\left(\sum_k\operatorname{Im}(Z_k), \sum_k\operatorname{Re}(Z_k)\right)\tag{3}$$

where $\phi$ is the estimated image orientation.

The quality of the estimate can be assessed by doing another pass through the data and by calculating the mean weighted square circular distance, $\text{MSCD}$, between phases of the complex numbers $Z_k$ and the estimated circular mean phase $4\phi$, with $|Z_k|$ as the weight:

$$\begin{gather}\text{MSCD} = \frac{\sum_k|Z_k|\bigg(1 - \cos\Big(4\phi - \operatorname{atan2}\big(\operatorname{Im}(Z_k), \operatorname{Re}(Z_k)\big)\Big)\bigg)}{\sum_k|Z_k|}\\
= \frac{\sum_k\frac{|Z_k|}{2}\left(\left(\cos(4\phi) - \frac{\operatorname{Re}(Z_k)}{|Z_k|}\right)^2 + \left(\sin(4\phi) - \frac{\operatorname{Im}(Z_k)}{|Z_k|}\right)^2\right)}{\sum_k|Z_k|}\\
= \frac{\sum_k\big(|Z_k| - \operatorname{Re}(Z_k)\cos(4\phi) - \operatorname{Im}(Z_k)\sin(4\phi)\big)}{\sum_k|Z_k|},\end{gather}\tag{4}$$

which was minimized by $\phi$ calculated per Eq. 3. In Python:

absZ = np.sqrt(X**2 + Y**2)
reZ = (X**4 - 6*X**2*Y**2 + Y**4)/absZ**3
imZ = (4*X**3*Y - 4*X*Y**3)/absZ**3
phi = np.arctan2(np.sum(imZ), np.sum(reZ))/4

sumWeighted = np.sum(absZ - reZ*np.cos(4*phi) - imZ*np.sin(4*phi))
sumAbsZ = np.sum(absZ)
mscd = sumWeighted/sumAbsZ

print("rotate", -phi*180/np.pi, "deg, RMSCD =", np.arccos(1 - mscd)/4*180/np.pi, "deg equivalent (weight = length)")

Based on my mpmath experiments (not shown), I think we won't run out of numerical precison even for very large images. For different filter settings (annotated) the outputs are, as reported between -45 and 45 degrees:

rotate 32.29809399495655 deg, RMSCD = 17.057059965741338 deg equivalent (omega_c = np.pi)
rotate 32.07672617150525 deg, RMSCD = 16.699056648843566 deg equivalent (omega_c = np.pi/2)
rotate 32.13115293914797 deg, RMSCD = 15.217534399922902 deg equivalent (omega_c = np.pi/4, same as in the Python listing)
rotate 32.18444156018288 deg, RMSCD = 14.239347706786056 deg equivalent (omega_c = np.pi/8)
rotate 32.23705383489169 deg, RMSCD = 13.63694582160468 deg equivalent (omega_c = np.pi/16)

Strong low-pass filtering appear useful, reducing the root mean square circular distance (RMSCD) equivalent angle calculated as $\operatorname{acos}(1 - \text{MSCD})$. Without the 2-d rotated cosine window, some of the results would be off by a degree or so (not shown), which means that it is important to do proper windowing of the analysis filters. The RMSCD equivalent angle is not directly an estimate of the error in the angle estimate, which should be much less.

Alternative square-length weight function

Let's try square of the vector length as an alternative weight function, by:

$$Z_k = \frac{(X_k + Y_k i)^4}{\sqrt{X_k^2 + Y_k^2}^2} = \frac{X_k^4 - 6X_k^2Y_k^2 + Y_k^4 + (4X_k^3Y_k - 4X_kY_k^3)i}{X_k^2 + Y_k^2},\tag{5}$$

In Python:

absZ_alt = X**2 + Y**2
reZ_alt = (X**4 - 6*X**2*Y**2 + Y**4)/absZ_alt
imZ_alt = (4*X**3*Y - 4*X*Y**3)/absZ_alt
phi_alt = np.arctan2(np.sum(imZ_alt), np.sum(reZ_alt))/4

sumWeighted_alt = np.sum(absZ_alt - reZ_alt*np.cos(4*phi_alt) - imZ_alt*np.sin(4*phi_alt))
sumAbsZ_alt = np.sum(absZ_alt)
mscd_alt = sumWeighted_alt/sumAbsZ_alt

print("rotate", -phi_alt*180/np.pi, "deg, RMSCD =", np.arccos(1 - mscd_alt)/4*180/np.pi, "deg equivalent (weight = length^2)")

The square length weight reduces the RMSCD equivalent angle by about a degree:

rotate 32.264713568426764 deg, RMSCD = 16.06582418749094 deg equivalent (weight = length^2, omega_c = np.pi, N = 41)
rotate 32.03693157762725 deg, RMSCD = 15.839593856962486 deg equivalent (weight = length^2, omega_c = np.pi/2, N = 41)
rotate 32.11471435914187 deg, RMSCD = 14.315371970649874 deg equivalent (weight = length^2, omega_c = np.pi/4, N = 41)
rotate 32.16968341455537 deg, RMSCD = 13.624896827482049 deg equivalent (weight = length^2, omega_c = np.pi/8, N = 41)
rotate 32.22062839958777 deg, RMSCD = 12.495324176281466 deg equivalent (weight = length^2, omega_c = np.pi/16, N = 41)
rotate 32.22385477783647 deg, RMSCD = 13.629915935941973 deg equivalent (weight = length^2, omega_c = np.pi/32, N = 81)
rotate 32.284350817263906 deg, RMSCD = 12.308297934977746 deg equivalent (weight = length^2, omega_c = np.pi/64, N = 161)

This seems a slightly better weight function. I added also cutoffs $\omega_c = \pi/32$ and $\omega_c = \pi/64$. They use larger N resulting in a different cropping of the image and not strictly comparable MSCD values.

1-d histogram

The benefit of the square-length weight function is more apparent with a 1-d weighted histogram of $Z_k$ phases. Python script:

# Optional histogram
hist_plain, bin_edges = np.histogram(np.arctan2(imZ, reZ), weights=np.ones(absZ.shape)/absZ.size, bins=900)
hist, bin_edges = np.histogram(np.arctan2(imZ, reZ), weights=absZ/np.sum(absZ), bins=900)
hist_alt, bin_edges = np.histogram(np.arctan2(imZ_alt, reZ_alt), weights=absZ_alt/np.sum(absZ_alt), bins=900)
plt.plot((bin_edges[:-1]+(bin_edges[1]-bin_edges[0]))*45/np.pi, hist_plain, "black")
plt.plot((bin_edges[:-1]+(bin_edges[1]-bin_edges[0]))*45/np.pi, hist, "red")
plt.plot((bin_edges[:-1]+(bin_edges[1]-bin_edges[0]))*45/np.pi, hist_alt, "blue")
plt.xlabel("angle (degrees)")
plt.show()

Figure 5. Linearly interpolated weighted histogram of gradient vector angles, wrapped to $-\pi/4\ldots\pi/4$ and weighted by (in order from bottom to top at the peak): no weighting (black), gradient vector length (red), square of gradient vector length (blue). The bin width is 0.1 degrees. Filter cutoff was omega_c = np.pi/4, same as in the Python listing. The bottom figure is zoomed at the peaks.

Steerable filter math

We have seen that the approach works, but it would be good to have a better mathematical understanding. The $x$ and $y$ differentiation filter impulse responses given by Eq. 1 can be understood as the basis functions for forming the impulse response of a steerable differentiation filter that is sampled from a rotation of the right side of the equation for $h_x[x, y]$ (Eq. 1). This is more easily seen by converting Eq. 1 to polar coordinates:

$$\begin{align}h_x(r, \theta) = h_x[r\cos(\theta), r\sin(\theta)] &= \begin{cases}0&\text{if }r = 0,\\-\displaystyle\frac{\omega_c^2\,r\cos(\theta)\,J_2\left(\omega_c r\right)}{2 \pi\,r^2}&\text{otherwise}\end{cases}\\
&= \cos(\theta)f(r),\\
h_y(r, \theta) = h_y[r\cos(\theta), r\sin(\theta)] &= \begin{cases}0&\text{if }r = 0,\\-\displaystyle\frac{\omega_c^2\,r\sin(\theta)\,J_2\left(\omega_c r\right)}{2 \pi\,r^2}&\text{otherwise}\end{cases}\\
&= \sin(\theta)f(r),\\
f(r) &= \begin{cases}0&\text{if }r = 0,\\-\displaystyle\frac{\omega_c^2\,r\,J_2\left(\omega_c r\right)}{2 \pi\,r^2}&\text{otherwise,}\end{cases}\end{align}\tag{6}$$

where both the horizontal and the vertical differentiation filter impulse responses have the same radial factor function $f(r)$. Any rotated version $h(r, \theta, \phi)$ of $h_x(r, \theta)$ by steering angle $\phi$ is obtained by:

$$h(r, \theta, \phi) = h_x(r, \theta - \phi) = \cos(\theta - \phi)f(r)\tag{7}$$

The idea was that the steered kernel $h(r, \theta, \phi)$ can be constructed as a weighted sum of $h_x(r, \theta)$ and $h_x(r, \theta)$, with $\cos(\phi)$ and $\sin(\phi)$ as the weights, and that is indeed the case:

$$\cos(\phi) h_x(r, \theta) + \sin(\phi) h_y(r, \theta) = \cos(\phi) \cos(\theta) f(r) + \sin(\phi) \sin(\theta) f(r) = \cos(\theta - \phi) f(r) = h(r, \theta, \phi).\tag{8}$$

We will arrive at an equivalent conclusion if we think of the isotropically low-pass filtered signal as the input signal and construct a partial derivative operator with respect to the first of rotated coordinates $x_\phi$, $y_\phi$ rotated by angle $\phi$ from coordinates $x$, $y$. (Derivation can be considered a linear-time-invariant system.) We have:

$$\begin{gather}x = \cos(\phi)x_\phi - \sin(\phi)y_\phi,\\
y = \sin(\phi)x_\phi + \cos(\phi)y_\phi\end{gather}\tag{9}$$

Using the chain rule for partial derivatives, the partial derivative operator with respect to $x_\phi$ can be expressed as a cosine and sine weighted sum of partial derivatives with respect to $x$ and $y$:

$$\begin{gather}\frac{\partial}{\partial x_\phi} = \frac{\partial x}{\partial x_\phi}\frac{\partial}{\partial x} + \frac{\partial y}{\partial x_\phi}\frac{\partial}{\partial y} = \frac{\partial \big(\cos(\phi)x_\phi - \sin(\phi)y_\phi\big)}{\partial x_\phi}\frac{\partial}{\partial x} + \frac{\partial \big(\sin(\phi)x_\phi + \cos(\phi)y_\phi\big)}{\partial x_\phi}\frac{\partial}{\partial y} = \cos(\phi)\frac{\partial}{\partial x} + \sin(\phi)\frac{\partial}{\partial y}\end{gather}\tag{10}$$

A question that remains to be explored is how a suitably weighted circular mean of gradient vector angles is related to the angle $\phi$ of in some way the "most activated" steered differentiation filter.

Possible improvements

To possibly improve results further, the gradient can be calculated also for the red and blue color channels, to be included as additional data in the "average" calculation.

I have in mind possible extensions of this method:

1) Use a larger set of analysis filter kernels and detect edges rather than detecting gradients. This needs to be carefully crafted so that edges in all directions are treated equally, that is, an edge detector for any angle should be obtainable by a weighted sum of orthogonal kernels. A set of suitable kernels can (I think) be obtained by applying the differential operators of Eq. 11, Fig. 6 (see also my Mathematics Stack Exchange post) on the continuous-space impulse response of a circularly symmetric low-pass filter.

$$\begin{gather}\lim_{h\to 0}\frac{\sum_{N=0}^{4N + 1} (-1)^n f\bigg(x + h\cos\left(\frac{2\pi n}{4N + 2}\right), y + h\sin\left(\frac{2\pi n}{4N + 2}\right)\bigg)}{h^{2N + 1}},\\
\lim_{h\to 0}\frac{\sum_{N=0}^{4N + 1} (-1)^n f\bigg(x + h\sin\left(\frac{2\pi n}{4N + 2}\right), y + h\cos\left(\frac{2\pi n}{4N + 2}\right)\bigg)}{h^{2N + 1}}\end{gather}\tag{11}$$

Figure 6. Dirac delta relative locations in differential operators for construction of higher-order edge detectors.

2) The calculation of a (weighted) mean of circular quantities can be understood as summing of cosines of the same frequency shifted by samples of the quantity (and scaled by the weight), and finding the peak of the resulting function. If similarly shifted and scaled harmonics of the shifted cosine, with carefully chosen relative amplitudes, are added to the mix, forming a sharper smoothing kernel, then multiple peaks may appear in the total sum and the peak with the largest value can be reported. With a suitable mixture of harmonics, that would give a kind of local average that largely ignores outliers away from the main peak of the distribution.

Alternative approaches

It would also be possible to convolve the image by angle $\phi$ and angle $\phi + \pi/2$ rotated "long edge" kernels, and to calculate the mean square of the pixels of the two convolved images. The angle $\phi$ that maximizes the mean square would be reported. This approach might give a good final refinement for the image orientation finding, because it is risky to search the complete angle $\phi$ space at large steps.

Another approach is non-local methods, like cross-correlating distant similar regions, applicable if you know that there are long horizontal or vertical traces, or features that repeat many times horizontally or vertically.


Why does my zoom audio, when there are nearby microphones, seem to usually take an open loop format? The signal from the speaker on device A gets into the microphone from device B, then the amplifier and then the speaker from device B, back into device A. The result is a sort of runaway signal that is always reamplifying itself?
Why can't I create a sort of closed loop system when I have nearby microphones? What would I need to do? What gear and/or software would I need to create a closed feedback loop?

This effect is usually undesirable and is prevented by (acoustic) echo cancellation algorithms. These algorithms estimate the distortion of the incoming signal on its way to the microphone (via the loudspeaker and the room acoustics). That distortion is modeled by a filter and the output of that filter is subtracted from the microphone signal, such that the contribution of the incoming signal to the outgoing signal is minimized.
Zoom has an option to turn off echo cancellation, so you could try playing around with that.


PSOLA how to calculate synthesized Marks I have successfully extracts all Pitch Markers from waveforms, I used autocorrelation to find my pitch periods, The markers can be observed in the figure below:

I want change my signal (Pitch or Speed) using PSOLA, I need find the synthesized Marks, for example if I choose change my Pitch one octave UP(0.5), for it I need find the new markers to synthesize my signal, you can see in this figure below one example.

How can I calculate the new marks to be synthesized (as Synthesized markers)?

Let $t_a(n)$ be the time of the n-th analysis mark ("Pitch marker"), $P_a(n)$ the period at the n-th analysis mark (it can be defined as $t_a(n+1) - t_a(n)$, $t_s(n)$ be the time of the n-th synthesis mark, $n_s(k)$ is a number whose integer part is the index of analysis segment to be overlapped-added when processing the n-th synthesis mark (it can also represent a weighted interpolation between two consecutive analysis marks).

The synthesis marks will be defined by the recurrence:
$$n_s(0) = 0$$
$$t_s(0) = t_a(0)$$
$$t_s(k) = t_s(k-1) + \alpha(n) P_a([n_s(k)])$$
$$n_s(k) = n_s(k-1) + \alpha(n)$$

$\alpha(n)$ is a stretch factor. It is 0.5 if you want to double the frequency. If you want unvoiced frames to be processed naturally (no stretching for them), you can actually make it dependent on the voiced/unvoiced analysis. In this case: it is 1 if the analysis mark $[n_s(k)]$ is unvoiced, and 0.5 or whatever is your pitch-shift factor when it is voiced.


Names of system functions in frequency domain I was just trying to refresh my systems theory known from long ago, and I realized that I had forgetten the name of the basic functions.  Specifically, what are the names of these functions respectively:
$$H(\Omega),~ H(\omega),~ H(s),~ H(z),~ |H(\Omega)|,~ |H(\omega)|,~ |H(s)|,~ |H(z)|$$
I'm assuming that $\Omega$ is the continuous frequency in radians, and $\omega$ is defined as the discrete frequency in radians.
I'm guessing that $H(\Omega)$ is the frequency response, and that $H(s)$ is the transfer function... but i'm not sure about the other names.

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


Calculating 1/3 Octave Spectrum from FFT / DFT I am not often on this forum and I am not an expert on the subject.

I struggle with the theory of FFT / DFT and the 1/3 octave spectrum.

Assume I have a DFT analysis of a given signal.
It (the DFT analysis) consists of many equidistant frequency bins that contain the corresponding amplitude (RMS or Peak).
I now want to calculate a 1/3 octave spectrum which has different frequency bins that are not equidistant.
Assume that one of the 1/3 octave bins reaches from 100 Hz to 200 Hz and the bin size of my DFT is 10 Hz.
Assume the amplitude of the DFT is already in RMS.

Here is what I do at the moment (in Matlab) but it seems that it's  not correct:

Adding all DFT amplitudes within 100 Hz and 200 Hz in linear scale (Pascal in my case).
After adding the RMS amplitudes (within one 1/3 octave bin) I just calculate corresponding the Decibel (dB) value - that's it.

Is there more to it? Can you point me in the correct direction?

  I do not "want" (I have to mimic another software) to use so-called
  1/3-Octave-Band Filter Banks like in this Matlab example.

Related

1/3 octave spectra from fft

Dr. Manuel Kuehner,
  You are close. You need to take the square root of the linear values squared.

$$P_{\mbox{total_linear}}=\sqrt{p_1^2+p_2^2+...}$$

$$P_{\mbox{total_dB}}=20 log_{10}\left( P_{\mbox{total_linear/20E-6}} \right)$$

FYI: I wrote a MATLAB function to do exactly as you request. It is here

Looking to read? See page 16 of this book:
https://www.amazon.com/Noise-Control-Engineers-Harold-Lord/dp/0070387389


What is the precise definition of a nonlinear signal? I tried to find the the precise definition of a nonlinear signal. However it is not easy. How can we define a nonlinear signal?

In my opinion, the term nonlinear signal doesn't make much sense. If you refer to nonlinear signal processing then this is not about the processing of nonlinear signals, but about the nonlinear processing of signals, i.e., about the use of nonlinear systems.

A system $\mathcal{T}$ is nonlinear if it doesn't obey the superposition principle. I.e., if $y_1(t)=\mathcal{T}\{x_1(t)\}$ is the response to input $x_1(t)$, and $y_2(t)=\mathcal{T}\{x_2(t)\}$ is the response to input $x_2(t)$ then

$$\mathcal{T}\big\{\alpha x_1(t)+\beta x_2(t)\big\}=\alpha\mathcal{T}\big\{x_1(t)\big\}+\beta\mathcal{T}\big\{x_2(t)\big\}$$

is generally not satisfied.

Well-known examples of such nonlinear systems are median filters and Volterra filters, but there are many trivial operations, such as taking the magnitude of a signal, which must be classified as nonlinear processing.


Map time/scene differences between two similar videos Goal:
Given two videos that are episodes from series but in different format (one Bluray/one TV for example), map the time differences between the two (which scenes are not in one of the two) and be able to compensate for some subtitles lags between version.
Current state:
I've began something thanks to this question (Map time difference between two similar videos) and started building a python script so now I'm able to transform each frame of each video to a one floating point value.

The problem is now, I dont know how to detect those lags as you can see in the picture. I'm trying some cross correlation maybe with some window but I'm not sure of how I should handle the added scenes in both signals (which means here signal that cannot be correlated).
Questions:
Do you know any algorithm/method that can help me with this process ?
Do you know what the user in this thread (Map time difference between two similar videos) means by "added a Gaussian of that length with amplitude as the phase difference to the mapping (at the appropriate interval)" ?
Thanks for your help !

This was an interesting read and it looks like you're not too far off! I think that what you have left to do is divide one of the signals into intervals (let's call it signal 1) and cross-correlate each segment of signal 1 with the entirety of signal 2. The peak in the cross-correlation will indicate where there is the strongest match for that segment in signal 2. You then calculate the difference in frames between where the segment was in signal 1 and where its strongest match is in signal 2.
The third step in the post you mentioned appears to be Kernel Density Estimation. The author placed a Gaussian at the center of each segment from signal 1. The standard deviation of the Gaussian, I assume is one-sixth the segment length, and the amplitude of the Gaussian is the frame difference. The Gaussian should extend from time 0 to the length of signal 2. He then summed all the Gaussians to obtained a continuous function.
I think you can achieve the same thing if you do the following: once you have the frame difference as a function of time, you convolve that function with a Gaussian whose standard deviation is one-sixth the segment length you used to divide signal 1.
I hope this helped!


Definition of the DFT / FFT Bin Size When reconstructing the frequency domain for an FFT, what is the most self-consistent way to do this -- i.e. how is it best to define the bin width, $\Delta f$?

For example, previously I thought it a simple matter of saying,
$$\Delta f = span / N_{points,FFT}$$ where the span is the size of the window you want to look at and $N_{points,FFT}$ is the number of points that should appear in the final FFT result. This makes very clear sense especially when you deal with some FFT spectrum analyzer where you set the span and the number of points in the final FFT (some times called the number of FFT lines) can also be set but tends to have some upper value for that device.

However I have also seen definitions attached to the sampling rate for the incoming transient signal to be FFT'd, as $$\Delta f = f_{sr} / N_{points,FFT}$$ or even as $$\Delta f = 2 f_{sr} / N_{points,FFT}$$ depending on how you define $N_{points,FFT}$. I.e. the number of points that goes into making the FFT or the number of points that will appear in the final FFT result because half the spectrum is thrown away due to mirroring.

Maybe it doesn't matter and it just depends on the way you choose to define things. But if anyone has any insight that would be helpful!

I am adding a little more to clarify my confusion.

Suppose you have some FFT spectrum analyzer that can process at maximum 2048 points in a time transient. That then means The final FFT would have 1024 points -- as we discard the mirror image. I then (again this is all in the context of making an FFT with some FFT spectrum analyzer device) say I want to set my span of this FFT to $50\ \rm{Hz}$. This would give a resultant FFT with a bin width of $\Delta f = 50 / 1024 = 0.049\ \rm{Hz}$. To me this must imply that the sample rate, $f_{sr}$, must vary according to the desired span.

So is this is all question of what your free parameters are, whether the maximum number of points in a transient signal is fixed, or the maximum sample rate you can take, and so on?

Let's assume we have a sample rate of $f_s=10 kHz$ and FFT size of $N=1000$

Your bin spacing is $\delta f = 10 Hz$. It's simply the sample rate divided by the FFT size. That's all there is to it.
Keep in mind that the FFT requires both the time domain and the frequency domain signal to be perodic. Most signals are not, so there needs to be some sort "accommodation" that's specific to the needs for your application.
A real-valued sine wave with a frequency that is an integer multiple of the bin spacing will only show up in two bins symmetric around zero. I.e. a 100 Hz sine wave will show up in bins +10 and -10 (assuming we index the bins from $[-N/2,N/2-1]$
A complex exponential  with a frequency that is  integer multiple of the bin spacing will only show up in one bin. I.e. a +100 Hz complex exponential will show up bin +10.
A cosine wave with a frequency that's NOT an integer multiple of the bin spacing (e.g. 105Hz) will spread out over all bins. The bulk of the energy will be in the neighborhood of 105 Hz but there will be energy in ALL bins. This is a direct consequence of bullet 2): the FFT assumes periodicity with $N/f_s$. In essence the FFT doesn't calculate the spectrum of the sine wave, but the spectrum of a signal that consists of 1000 sample long snippet of a 105Hz cosine that gets repeated infinite times. At the splice point there is a massive discontinuity and that causes the energy to spread over all bands
In order to minimize this spreading effect, one can apply a window, that tapers down the edges of the time domain frame and reduces the discontinuity. On the other hand, the window also smears out the energy in the frequency domain, so you can't get any sharp lines any more.
Many people think that the FFT will give you an "exact" representation of the signal's frequency content. That's not the case unless the time domain signal is either periodic with the FFT length or if it's finite in time. For all other signals the results will depend heavily on the exact process and parameters chosen.

EDIT: reworded item 7 per the comment


Object recognition using color information I want to use the L*a*b color model for recognition of a particular colored object by calculating the euclidean distance between the desired color in current image and the same color in next image frame.

As the LAB color space was designed in keeping with human vision, I think it can be used for the purpose.  Is this correct?

Yes Lab color is designed to approximate human vision. But there are also RGB and CMYK color models, and their transforms like HSL and HSV. Which color base can best make your object feature distinguishable needs your experiment to validate.


How to calculate the channel gain if the channel is a $N$ by $M$ matrix,$N,M>1$ How to calculate the channel gain if the channel is a $N$ by $M$ matrix,$N,M>1$

For example,i create a rayleight fading channel

H_AB  = 1/sqrt(2)*[randn(2,6) + j*randn(2,6)],and H_AB is a 2 by 6 matrix,if it is a $N$ by 1 matrix,i can just use $|H_{AB}|^2=H_{AB}^HH_{AB}$ to calculate a gain value,however,now $H_{AB}$ is a matrix,if i still use $|H_{AB}|^2=H_{AB}^HH_{AB}$,i will get a matrix,and power should be a value ,not a matrix.So how can i calculate the channel gain of $H_{AB} ?$ is it $trace(|H_{AB}|^2) ?$

It is important to understand the what the channel matrix $\mathbf{H}$ physically means. The channel matrix is a symmetric matrix that contains the fading coefficients, that is, in row $i$ column $j$ the value ${h}_{ij}$ is the coefficient between the $i^{th}$ transmit antenna and the $j^{th}$ receive antenna. The gain between any transmit/receive antenna pair can be found by $|{h}_{ij}|^2$.


Not able to reach minimum phase using Hilbert transform My problem is pretty simple, I've designed a magnitude response and I would like to find the corresponding minimum phase filter. I'm using the code below and unless there is a bug my eyes don't want to see, I'm pretty confident on the methodology.

nbTaps = 100;
fs = 96000;
desiredMagn_dB = [0 0 -100 -100];
f = [0 24000 26000 fs/2];

FIRFreq = (0:nbTaps/2)*fs / nbTaps;

FIRMagn = interp1(f, 10.^(desiredMagn_dB./20), FIRFreq);
FIRMagn = [FIRMagn fliplr(FIRMagn(2:end-1))];

minPhase = -imag(hilbert(log(FIRMagn)));

freqResp = FIRMagn.*exp(1i.*minPhase);

FIRCoefs = ifft(freqResp, 'symmetric');

And yet, my filter has not a minimal phase, since I have zeros outside the unit circle. I am a little bit lost, here, and I would like to understand why is so? or where I am wrong?

EDIT from r b-j:

nbTaps = 16384;
fs = 96000;
desiredMagn_dB = [0 0 -100 -100];
f = [0 24000 26000 fs/2];

FIRFreq = (0:nbTaps/2)*fs / nbTaps;

FIRMagn = interp1(f, 10.^(desiredMagn_dB./20), FIRFreq);
FIRMagn = [FIRMagn fliplr(FIRMagn(2:end-1))];

minPhase = -imag(hilbert(log(FIRMagn)));

freqResp = FIRMagn.*exp(1i.*minPhase);

FIRCoefs = ifft(freqResp, 'symmetric');

figure(1)
plot(20*log10(FIRMagn))
figure(2)
plot(minPhase)
figure(3)
plot(FIRCoefs)

i don't see anything unexpected.

More taps. You don't have anywhere near enough taps for a filter that steep. Start large with 8192 or so cut to desired accuracy, if needed
Due to the low number of tabs you are seeing the effect of "circular" hilbert transform. See for example: http://andrewduncan.net/air/
How do you know you have zeros outside the unit circle? Calculating the roots of a high order polynomial is a numerically difficult problem and using a Matlab function like "roots()" will given you often inaccurate or plainly wrong results.


How to know if the noise was removed in an image? I am looking for a technique for measure the quality of my denoise filter, is the SNR formula useful for that purpose?

My problem is, I have an image with noise, I apply the "adaptive median filter", how to know if the noise has decreased?

You can measure that with SNR, yes.  What you call 'signal' can get a little fluffy, but you should at least be able to compare two versions of the same image (before and after).

If your image has a nice homogeneous area in it, you can measure the variance of the pixel values in this region, and compare this before and after - so long as the image doesn't otherwise change.  Do make sure that you account for any scaling or shifting that a filter might be doing (an adaptive median filter shouldn't do either of these if it is implemented properly).

A common method for testing denoising filters is to start with a pure image (one that you define yourself), then add noise to it artificially then compare the denoised version to the original.  You will need to make sure to add the correct type of noise for your context (is it simple Gaussian, or Poisson, or ...?).


Are there DC removal techniques that do not introduce a filtering transient? I have an application in which I need to remove DC components from a sampled signal in real time (basically the same constraints of ECG, but it has nothing to do with the brain). My measured signal in in the single µV range and the undesired offset wanders around and can be as high as tens of mV.

In this application, fast electrical stimulation is introduced which can reach hundreds of volts and can introduce a measurement artifact of nearly a volt. This stimulus is basically an impulse from the point of view of the processing.

Currently I am using a 200Hz 1st order highpass, as this has an acceptable artifact (i.e., impulse response) after the stimulus. Higher orders are not needed, and these are undesirable as it would make the artifact worse.

My problem is that the actual application will require filters to go down to 1Hz or even lower (in some cases we use 0.05Hz, as this is above the frequency content of the "wandering"). A filter at these frequencies would have an unacceptably long impulse response.

Are there any other techniques that would remove the DC-wandering without having the impulse response artifact problem?

Addendum:

The way this problem is usually handled involves knowledge of the location and duration of the disturbance (whose description I simplified for this question). This enters the category of evoqued response and involves some non-linear processing. Some of the alternatives I have tried:

Blanking the interference before it hits the filter. This requires precise synchronization that is not always available.
Remove the real-time aspect and process the data non-causally. Long data records are still required due to filter transients, initial conditions are critical (and, of course it is not real time).
Model the transients and remove them with data fits. Precise synchronization is no longer needed, but this might introducers own artifacts.

Any filter with a DC gain of zero and, for most other frequencies, a gain of (or close to) 1, is a DC blocking filter.

So just make a high-pass filter (HPF) and put the corner frequency as close as you can to zero.

Of course, when the system is switched on, there will be a transient from blocking DC of 0 to blocking DC of whatever value it is.  Without prior knowledge of what the initial DC is, I do not see how you can assume the initial value is anything other than zero.  So all of the initial states of your HPF are zero.  But once the turn-on transient is over, there should be no other transients unless the DC component has a sudden change (like a step function) and is not just "trending" to another value.

If this is done in fixed-point, you will have some serious numerical issues.  You might even have serious numerical issues if it's floating-point and your HPF corner frequency is very very low. You might want to look at this old thing from me or this article from IEEE Signal Processing magazine.  I dunno where to find it free.


Is it possible to get the envelope of a signal via FIR filter? I was reading an article about ultrasonic signals, and I read in it that it is possible to get the envelope of a signal by using a band pass FIR filter.
My question: Is it possible to get the envelope of a signal by passing it through a band pass FIR filter?
Thank you in advance.

Getting the envelope of a signal requires more than a simple band-pass filter.

However, you could rectify the signal and bandpass-filter the rectified signal to get the enveloppe. That's how I used to do it with ultrasound signals.

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

You could also use more elaborate techniques, like the Hilbert transform (or Hilbert filter) to get a complex signal that would allow use to use the modulus (or magnitude) of the complex signal to get the envelope.


Convolution Kernels for Image Filtering As I understand it, in image processing, when it is desired to apply a filter of some sort to an 2D image, a kernel is applied to all the pixels of the 2D image in a process called convolution. Median, Gaussian, Bilateral, all are examples of kernel. These kernels are applied in a similar fashion to that shown in this example http://www.songho.ca/dsp/convolution/convolution2d_example.html (which I found here, searching Signal Processing at StackExchange).

Is there a good list of all/many/most kernels, there uses, and an example before/after image?

Also, is it correct to say that all kernels are just matrices which get convolved against a source?

Be careful, a median filter cannot be expressed as a convolution, and thus is not considered a kernel in this respect. This is because the median filter is based on order statistics of an image patch, and the resulting pixel at the output of a median filter is not a linear combination of other pixels within a patch.

Otherwise, you are right, kernels are generally thought of as just other 'small images', that one would convolve with the image to be analyzed. Those kernels can be chosen for a variety of reasons. Gaussian kernels are used for smoothing, laplacians are used for edge detections, etc.

There are as many kernels as you can imagine. However usually we start with a particular application in mind, "I need to blur the image", or "I need to detect edges in the x-direction", and then find kernels that magnify those features we are after.


FFT for bearing fault analysis Why is FFT a good method for detection of bearing faults? Why can't we do the analysis in time-domain?

Short answer: We can, but it is not always as robust as frequency domain techniques.

Longer answer: This is very broad topic, but let me try to throw some light. Time domain techniques tend to be naive. For example some are based on the RMS energy monitoring of the signal being recorded by an accelerometer. Main disadvantage is that on early stages of fatigue you can't distinguish between faulty and working bearing. Additionally setting a threshold is a real pain. Also some people are using methods based on crest factor (ratio of peak to RMS), knowing that crest factor for impulsive sounds is high. Although there are some 'fancy' methods for analysis in time domain, such as: New Time Domain Method for the Detection of Roller Bearing Defects

Frequency domain techniques (such as Cepstrum or Bispectrum) on the other hand are more robust and have properties that are desirable, such as you can very easily correlate detected fault frequency with a specific bearing and type of damage. For example when you use methods such as Bispectrum, which is (simply speaking) based on detecting presence of higher harmonics. Because it is also taking phase into consideration, in theory even if harmonic is almost covered with noise, you can still detect if its correlated with fundamental frequency - smart isn't it?


Determining pulse shape for Viterbi Equalization I believe I understand well enough how the Viterbi algorithm can reverse the effects of ISI when the combined impulse response (transmit filter, multipath, receive filter) is known. I've read How I Learned to Love the Trellis by Bernard Sklar, and I also asked about it in a previous question: Applying Viterbi algorithm to compensate for ISI in PSK31.

But how can this technique be applied when the channel distortion isn't known? From my superficial understanding of equalizing algorithms like CMA, they attempt to infer the channel distortion and compensate for it dynamically. My particular interest is HF channels, which would certainly seem to benefit from some dynamic adjustment.

However, in my study of Viterbi equalization so far, I've not found an explanation of how the unknown impulse response of the channel is found. Without that information, the best we can accomplish is to compensate for known elements like the receive and transmit filters. And in a case with no inherent ISI (like RRC filters on each end), then does Viterbi equalization have any value at all?

I'll answer my own question with what I've learned in the past 28 days. Hopefully I've got it right!

Viterbi equalization is just one way to equalize the channel, as an alternative to some other method such as convolving the signal with an equalizing filter. Estimating the channel is an entirely separate problem which must be solved independently of the equalization method.

The trouble with an equalizing filter is it compensates for dips in the channel frequency response by amplifying those frequencies. Unavoidably, the noise is also amplified.

There's a tradeoff to be made: the flatter the frequency response is made, the lower the error due to ISI, but the higher the error due to noise. There's an optimal amount of equalization which minimizes the error due to noise and ISI. An equalizer that tries to hit this sweet spot is a minimum mean square error (MMSE) equalizer.

Viterbi equalization offers an alternative solution. The channel estimate will produce an estimated channel impulse response. The coefficients in this impulse response can be used to build the trellis, and the Viterbi algorithm can then find the most likely sequence of transmitted bits, taking into account the estimated distortion from the channel. The receive filter remains a simple matched filter, with no equalization. The Viterbi equalizer is not linear (unlike the equalizing filter alternative), so it's able to compensate for deep spectral nulls in the channel response in a way the equalizing filter can't. The cost is increased complexity.


Fourier transform of function division in time domain $F_1(\omega)$ is the Fourier Transform of $f_1(t)$. $F_2(\omega)$ is the Fourier TRansform of $f_2(t)$. Can I obtain the Fourier Transform ($F_3(\omega)$) of

$$ f_3(t) = \frac{f_1(t)}{f_2(t)} $$

directly from $F_1(\omega)$ and $F_2(\omega)$? I mean, is there anything similar to the equivalence between Multiplication in the Time Domain to Convolution in the Frequency Domain but for the Division operation?

For division there is no equivalent to the duality between multiplication and  convolution. Note that from the existence of the Fourier transforms of $f_1(t)$ and $f_2(t)$, you cannot conclude anything about the existence of the Fourier transform of $f_3(t)=f_1(t)/f_2(t)$. So $f_3(t)$ might not even have a Fourier transform, and if it exists, it cannot be expressed in any simple way in terms of $F_1(\omega)$ and $F_2(\omega)$.


What is the significance of first zero upcrossing of a cross correlation? What is the significance of the first zero up-crossing of a cross-correlation? The signals being correlated are cosines signals with some phase difference.

EDIT 1 : Cross-correlation plot attached. The signals are cosines with phase difference and may be with some low frequency components.

Sampling frequency of the measured signal is taken as 50Hz

Cross correlation obtained using:

Cxx=fftshift(ifft(fft(x,N).*conj(fft(y,N))))/(norm(x) * norm(y));

Not sure what you are asking specifically. The cross correlation of two cosine waves is just another cosine. The triangular envelope in the picture is simply a consequence of the limited length of the signals.

The phase of the cross-correlation cosine is given by the phase difference of the two input cosines. This phase difference determines the position of the zero crossings. So in effect the position of that zero crossing is some metric of the phase difference of the input signals


Does an audio equalizer consist of a perfect-reconstruction filterbank? Within the theory of perfect reconstruction, I discovered that it is not obvious that the combination of all filters in a filterbank gives you the original signal. Therefore, I was wondering whether an audio equalizer (either analog or digital) distorts the spectrum of a signal, even when the equalizer has 0 dB gain in each frequency band. Is this ever considered a problem?

A filter bank system is generally composed of:

an analysis filter bank, that split an input signal into components with filters and may reduce their rate
a synthesis filter bank (FB), that takes the components, potentially increase their rate, and feed them into filter elements

so that the novel outputs can be combined into a signal. It is perfect reconstruction (PR) when, for every input, the output is a potentially delayed or scaled version of the input. Filters can be nonlinear, the rate reduction can be different for each component. Min-Max or morphological wavelets are examples of the above. There are even hybrid filter banks, where the analysis (A) or synthesis (S) FB are either analog or digital. Their design is quite involved, so often A and S are linear, with upsampling and downsampling operators, and an output obtained by summation, as explained by @Marcus Müller. However, the filters may have different bandwidths and the subsampling can be different for each filter.

Many signal processing systems, such as (most) equalizers, are made of an analysis FB, a processing unit (often nonlinear) and a synthesis FB (which could be a mere summation). The non-perfection arises, for PR FB, because of the middle filtering component.

Assuming an equalizer has linear A and S filters, the non PR aspect may arise from two aspects. First, what are the filters? It seems to me that the traditional filters (first or second order) do not allow easy PR. However, with sufficient frequency overlap, they together can form a near PR filter banks, with an error below the hearing threshold.
The second aspect is the processing. If you suspect saturation, or nonlinear phase alteration, then PR is not achieved in  general.

Phase distortion is sometimes considered harmless for audio in some   textbooks. However, picky sound engineers are cautious about audible phase distortion. The article Plugging into EQ discusses those issue, around FIR or IIR filters. The first ones can have linear phase, but have less sharp frequency transitions, and require a handful of computations. The second ones do not have linear phase in general. It also provides a comparative evaluation of several EQ.

So, as Markus, I'll say no in general, with a modest yes (in cases) if below the hearing sensitivity.


Imposing sparse structure in Gabor space I have a 2-D signal collected in the frequency vs time domain. There is a feature of interest in this signal (a chirp) appearing at some arbitrary times which I have successfully extracted using a properly orientated 2-D Gabor filter.

What I would like to do is exploit the detection of this feature to represent the data in a sparse domain, so that Compressed Sensing can be used to capture the data regions which contain the chirp trait without sampling the entire data set a-la Nyquist.

I have searched for a way to do this online but haven't found very much. The literature I've read says that the basis must be orthogonal for perfect reconstruction, and that Gabor domain is most times not, but given that I'm only interested in a faster detection (yes contains/no, does not contain chirp) and not really a perfect reconstruction, is such a thing possible? If the approach above is not the best way to go about doing this, recommendations would be welcome as well.

In my experience, you do not necessarily need an orthogonal dictionary, see for example Compressed Sensing with Coherent and Redundant Dictionaries (DOI: 10.1016/j.acha.2010.10.002).

The paper Quantifying the performance of compressive
sensing on scalp EEG signals (DOI: 10.1109/ISABEL.2010.5702814) seems to be an example of the use of Gabor dictionaries (in EEG) in compressed sensing. I would just go ahead and try it.

If the chirp you are looking for can have a fairly fixed location in the Gabor domain, it might be possible to use compressive detection or classification instead of actual reconstruction. This is explained in Signal Processing With Compressive Measurements (10.1109/JSTSP.2009.2039178. You can see an example of how we have used it for symbol detection in IEEE 802.15.4 communication.


Harmonic oscillation and Fourier I am back to college after 6 years of professional life and I need help with the following problem

  [...]
  The standard form of the function is $S_{E}(t) = e^{j2\pi f_k t}$, whose derivation results from the convolution theorem:
  $$
h(t)\star e^{j2\pi f_k t} = \int_{-\infty}^{\infty} h(\tau)e^{j2\pi f_k (t-\tau)}d\tau = e^{j2\pi f_k t} \int_{-\infty}^{\infty} h(\tau)e^{-j2\pi f_k \tau}d\tau = e^{j2\pi f_k t} \cdot H(f_k)
$$
  According to the equation above, the transfer function in the time domain can be determined when the channel is excited with an Eigenfunction $s_E(t)=e^{j2\pi f_k t} $. Is this relationship also valid for the harmonic oscillation $s_O(t)= \cos(2\pi f_kt)$? Please provide your derivation (Fourier correspondence table)."

I need some guidance for this problem or some hint for textbooks.

When you express the cosine as $\cos(2\pi f_kt)=\frac{1}{2}(\exp(j2\pi f_c t)+\exp(-j2\pi f_ct)$ and insert it into the convolution integral you get:

$$
h(t)*\cos(2\pi f_ct) = \frac{1}{2}(\exp(j2\pi f_ct)H(f_c)+\exp(-j2\pi f_ct)H(-f_c))
$$

In case your $h(t)$ is real-valued, you have $H(f_c)=H(-f_c)^*=|H(f_c)|\exp(j\phi)$, and you get
$$h(t)*\cos(2\pi f_ct)=\cos(2\pi f_ct+\phi)|H(f_c)|.$$

In case you have a complex impulse response, this does not hold.


How to recognize for small video sequence with a large sequence? I'm working on the idea of recognizing a small video sequence (a sequence of 18 seconds) within a long one (a sequence of 15 minutes), the first idea that come to mind stats that I have to characterize this small video with some kind of identification. so I will use the identification to look this small sequence inside the long sequence, I'm asking this question to the community, to share with me an approach that can help me in characterize this small video, then recognize it within the long sequence.

thanks.

The most intuitive approach here would simply be taking every pixel's color values, and finding peaks in the individual pixels' reference video/short video cross correlation function. That's basically a W×H-dimensional crosscorrelation function; you could speed that up by doing fast convolution.

You could algorithmically reduce the effort by first reducing the dimensionality of your problem, for example by reducing each video frame to a global brightness value and a mean motion vector of e.g. the frame center region – then you'd only have to calculate a threedimensional crosscorellation function.


Envelope of time domain signal // Hilbert transform to get the envelope of a time domain signal, s(t), I use the absolute of its analytical time domain signal. But there occur the problem that sometimes the envelope is not as desired a true envelope but it follows more or less the original time domain signal.

Just for illustration some python code where the plot looks good:

from pylab import *
import numpy
import scipy.signal.signaltools as sigtool

x = linspace(0,5,1e4)

s= numpy.sin(2*numpy.pi*x)*numpy.exp(-x**2/2)
env = numpy.abs(sigtool.hilbert(s)) # hilbert(s) actually returns the analytical signal

plot(x,s,label='Time Signal')
plot(x,env,label='Envelope')
legend()
show()

By changing the data for example to (taking abs):

s = numpy.abs(numpy.sin(2*numpy.pi*x)*numpy.exp(-x**2/2))

the resulting envelope is following the "ripples" instead of being an envelope.

--- WHY? ---

Hilbert envelope, also called Energy-Time Curve (ETC), only works well for narrow-band fluctuations. Producing an analytic signal, of which you later take the absolute value, is a linear operation, so it treats all frequencies of your signal equally. If you give it a pure sine wave, it will indeed return to you a straight line. When you give it white noise however, you will likely get noise back.

Example 1:

$$
\lvert\cos(\omega t)\rvert^2 + \lvert i\sin(\omega t)\rvert^2  = \cos^2(\omega t) + \sin^2(\omega t)= 1
$$

Here, you have $\cos(\omega t)$ as you input signal, and its analysis counterpart is $\cos(\omega t) + i\sin(\omega t)$. In it easy to see from expression above that its envelope is always a constant. Let's look at another example.

Example 2:

Let $x(t) = \cos(\omega_1 t) + \cos(\omega_2 t)$. Then its analytic counterpart is
$$
a(t) = \cos(\omega_1 t) + \cos(\omega_2 t) + i\sin(\omega_1 t) + i\sin(\omega_2 t).
$$
Let's take the absolute value of this signal.

$$\begin{align}
\lvert a(t)\rvert^2 &= \lvert\cos(\omega_1 t) + \cos(\omega_2 t)\rvert^2 + \lvert i\sin(\omega_1 t) + i\sin(\omega_2 t)\rvert^2\\&=\left(\cos(\omega_1 t) + \cos(\omega_2 t)\right)^2 + \left( \sin(\omega_1 t) + \sin(\omega_2 t) \right)^2 \\ &= \cos^2(\omega_1t) + \cos(\omega_1t)\cos(\omega_2t) + \cos^2(\omega_2t)\\& \quad + \sin^2(\omega_1t) + \sin(\omega_1t)\sin(\omega_2t) + \sin^2(\omega_2t) \\&= \cos^2(\omega_1t) + \frac{1}{2}\left[ \cos(\omega_1t-\omega_2t)+\cos(\omega_1t+\omega_2t)\right] + \cos^2(\omega_2t)\\&\quad + \sin^2(\omega_1t) + \frac{1}{2}\left[ \sin(\omega_1t-\omega_2t)-\sin(\omega_1t+\omega_2t)\right] + \sin^2(\omega_2t) \\ &= \left[ \cos^2(\omega_1t) + \sin^2(\omega_1t)\right] + \left[ \cos^2(\omega_2t) + \sin^2(\omega_2t)\right] \\ &\quad + \frac{1}{2}\left[\cos(\omega_1t-\omega_2t) + \cos(\omega_1t-\omega_2t) \right] + \frac{1}{2}\left[\sin(\omega_1t+\omega_2t) - \sin(\omega_1t+\omega_2t) \right] \\ &= 1 + 1 + \cos(\omega_1t - \omega_2t) + 0 \\ &= 2+\cos\left( \left[\omega_1-\omega_2\right]t\right)\end{align}$$

As you can see, this is not a constant as one would expect. Note that if $\omega_1$ and $\omega_2$ are close to each other, only low-frequency oscillations are produced as an artifact. hope this helps.


TV Norm - What Would Be the Formula? I have a 2-D discrete signal in which each point can be represented as $(x, y)$. These points are varying with time $t$. Can we represent the TV-norm using the following formula?

$$\sum_{t}|x_t - x_{t-1}| + |y_t - y_{t-1}|$$

Your definition looks fine to me, if your goal is to treat the signal at every fixed $(x,y)$ location as a function of time and then compute 1D total variation along time. Is that what you intend to do?

Here are a couple of other ways you can define it (again, I have no clue what your application is, and why you want to compute the TV norm in the first place, so I don't know which of these make sense in your application).

To fix notation, let's denote your signal of interest as $z(n_x,n_y,n_t)$ where $z$ is the quantity measured at discrete locations indexed by $(n_x,n_y)$ in space and at a time sample number $n_t$.

TV norm for each fixed spatial location:
Fix a location $(n_x^*,n_y^*)$.
$$
TV = \sum_{n_t} |z(n_x^*,n_y^*,n_t+1) - z(n_x^*,n_y^*,n_t)|
$$

TV norm for each 2D time snapshot:
Fix a time instant $n^*_t$.
$$
TV = \sum_{n_x,n_y} |z(n_x+1,n_y,n^*_t) - z(n_x,n_y,n^*_t)| + | z(n_x,n_y+1,n^*_t) - z(n_x,n_y,n^*_t)|
$$

TV norm in 3D space-time:

$$
TV = \sum_{n_x,n_y,n_t} |z(n_x+1,n_y,n_t) - z(n_x,n_y,n_t)| + | z(n_x,n_y+1,n_t) - z(n_x,n_y,n_t)| + |z(n_x,n_y,n_t+1) - z(n_x,n_y,n_t)|
$$

See also: https://en.wikipedia.org/wiki/Total_variation_denoising#2D_signal_images


Audio wave amplitude decreased after FIR filter I am doing sampling and filtering for a 17KHz sine audio wave signal, the Nyquist frequency is about 20.8889KHz. I found the fact that this audio signal's amplitude decreased dramatically after it went through an FIR bandpass filter( pass band 16.8khz-17.2khz), the amplitude also decreased after an FIR highpass filter(cutoff frequency 17khz). I think that's because the audio signal frequency is kinda close to Nyquist frequency, am I right?
Is there a filter can preserve the original amplitude of the audio signal within the frequency I want?
Shown below is a pic consists of 2 signals before and after filtering(highpass filter)
Thank you in advance

As per Peter K. recommendation, I will reformulate my comment as an answer, instead of starting by asking:

I think the fact that the cutoff frequency of your HPF is 17KHz (your tone frequency) can perfectly explain why your signal is being attenuated.

I would recommend you to try to move down the cutoff frequency of the High Pass Filter in order to avoid attenuating your tone.

BR


Image processing for eye movements I have series of images extracted from a eye tracker video for both eyes.
is there any image processing algorithm/method/technique available to find the following from those images:

Horizontal eye movement for both eyes.
Vertical eye movement for both eyes.
Eye blinking for both eyes.

Thanks in advance for the help.

You can see for algorithms in ieee sites there are tons of algorithm which might suit your requirement and

An Algorithm for Real-Time Eye-Movement Tracking with Free Head Movement

Robust algorithm for video based eye tracing.
Are some examples.


How to determine a suitable threshold when using amplitude thresholding method for spike detection? I have seen number of approaches that have been introduced in the literature for spike detection including, amplitude thresholding method, Nonlinear Energy Operator method, and template matching method. However, If I am to use the amplitude thresholding method, how can I identify a suitable threshold value for my dataset? I have seen several works using this threshold to be X times the standard deviation above the mean noise level. But this X seems to be different from one study to another. For example, commonly used values are 5, 5.5, and 6.
My question is how can I determine a suitable X for my dataset? My dataset is a set of MEA recordings taken from organoids in which I do not have any prior knowledge about neuron firing.

For threshold-based spike detection, the threshold defines your performance in the detection vs false alarm performance. That is to say, a lower threshold will improve your probability of detecting (smaller) spikes but eventually, some of the noise may be considered a spike, increasing the probability of a false alarm, detecing a spike where there is none. A common term for showing this trade-off is ROC.
Therefore, it is highly application-dependent how to choose the threshold. If you have a specific false alarm rate that is tolerable, you can use the noise distribution to determine the threshold that gives you the desired false alarm rate (refered to as the constant false alarm rate (CFAR) approach). Same is true if you have a tolerable missed detection rate (àla I need to find at least 99.99% of the spikes).
If not, you should consider the threshold a variable and do your evaluations as a function of the threshold. This allows you to pick a threshold that optimizes your final figure of merit.


How to filter out noise from non-stationary signal I have this non-stationary signal.

the mean is roughly constant but the second moment (autocorrelation) does not depend only on the time lag $tau$.
Correct me if I am wrong in the above statement.
Anyway as you see I have high frequency noise components and I would like to clean it but avoiding any sort of phase lag introduced by a filter.
What can I do to achive my purpose?
I thought of using wavelets to extract the coefficients and then reconstruct the signal, but maybe there is a more efficient and correct way of doing it.
Can you also provide with a link/guide to code this in python?

Break your signal into short overlapping frames over which it is stationary. Say, the length of the frame is $L$ samples. You can window it to smooth out discontinuities. For example, you can use a Hann window with an overlap of $L/2$ samples.
On each windowed frame, use the Wiener filter to find the optimal inverse filter coefficients to cancel your noise. See scipy.signal.wiener.
Overlap add your filtered signal frames to recover the original non-stationary signal.


What are the cosine functions in JPEG's DCTII table? As I understand it from this video, in a JPEG image an 8x8-pixel block is made up of weighted cosine waves, calculated using DCTII. There are lots of visualizations of these waves, such as this one from the Wikipedia article on JPEG compression:

I'd like to be able to visualize each of the waves for a given 8x8-pixel block in an image, weighted with values calculated by DCTII.

Where can I find (or how can I calculate) the 64 cosine functions used in JPEG compression?

Those 2D cosine functions are independent of your input image. They are just "cosine waves", or the 64 basis functions that yield 64 coefficients when transforming "$8\times 8$" blocks. Given a $D$ $8\times 8$ matrix for the DCT-II, and a $8\times 8$ image patch $I$, you'll get $64$ coefficients $C$ by:

$$C=DID^T$$

[EDIT] If you want to display the 2D DCT-II functions only, you have to create 2D arrays $D_{h,v}(m,n)$ ($0\le m< M$ and $0\le n< N$) for each couple of horizontal/vertical integer indices $(h,v)$, with $0\le h< H$ and $0\le v< V$, with $H=M$ and $V=N$. Classically for JPEG, $M=N=8$. A formula is (up to a transposition, and possibly a scale factor):

$$D_{h,v}(m,n) = 2\sum_{m=0}^{M-1}\sum_{n=0}^{N-1} \frac{1}{\sqrt{M}}\frac{1}{\sqrt{N}}\eta_h  \eta_v \cos\left(\frac{\pi h}{2M}(2m+1)\right)\cos\left(\frac{\pi v}{2N}(2n+1)\right)$$

with $\eta_x = \frac{1}{\sqrt{2}}$ if $x=0$, $\eta_x = 1$ if $x\neq0$. There are a lot of clever implementations, I'll provide the most straightforward. The goal is to get a picture like the following:

A pseudo-code version (from Matlab) could be (with the issue of 0-based or 1-based indices):

nRow = 8;nCol = 8;
nFreqHoriz = 8;
nFreqVerti = 8;
for iFreqH = 1:nFreqHoriz
    for iFreqV = 1:nFreqVerti
        iFreqHoriz = iFreqH-1;
        iFreqVerti = iFreqV-1;
        normFreqDC = 1/sqrt(2);
        matFreq2D = zeros(nRow,nCol);
        for iRow = 1:nRow
            for iCol = 1:nCol
                iRow0 = iRow-1;
                iCol0 = iCol-1;
                matFreq2D(iRow,iCol) = 2/(sqrt(nFreqHoriz)*sqrt(nFreqVerti))*cos(pi*iFreqHoriz*(2*iRow0+1)/(2*nRow))*cos(pi*iFreqVerti*(2*iCol0+1)/(2*nCol));
                if ~iFreqHoriz ,  matFreq2D(iRow,iCol) =  matFreq2D(iRow,iCol)*normFreqDC;end
                if ~iFreqVerti ,  matFreq2D(iRow,iCol) =  matFreq2D(iRow,iCol)*normFreqDC;end
            end
        end
        sum(matFreq2D(:).^2)
        subplot(nFreqHoriz,nFreqVerti,8*(iFreqV-1)+iFreqH)
        imagesc(matFreq2D);axis off;axis square;colormap gray
    end
end


What does poles in unit circles center mean? Suppose i have all my poles in unit circle center. What kind of information this gives me?
Can i determine if this filter is IIR or FIR?

If you have $N$ poles at $z=0$ this means that your transfer function has the form
$$H(z)=\frac{P(z)}{z^N}\tag{1}$$
where $P(z)$ is a polynomial. What does this tell you about the filter type? (Hint: write $H(z)$ as a sum of positive and negative powers of $z$).


What is the theoretical probability of error for OOK transmission? I am trying to simulate an optical wireless communication channel which uses OOK modulation.

Looking at equation 14 in this literature, I found that:

For electrical SNR at the receiver

$$\textrm{BER}=Q\left(\sqrt{\rm SNR}\right)$$

where (according to this page)

$$Q(z)=\frac{1}{2} \text{erfc}\left(\frac{z}{\sqrt{2}}\right)$$

Please note that:

  electrical signal-to-noise ratio is proportional to the square
  of the received optical power Pr.

Next, I implemented $Q(z)$ as follows:

getOokBer (double snr){
    double ber = 0.5 * erfc (std::sqrt (snr/2));
    return ber;
}

At this point I got stuck. I think the next step is to calculate $P_\rm e$ and/or $P_{\rm err}$ which has been discussed here too.

In this paper $P_\rm e$ is given as:

  $P_{\rm e}(k, l)$ -- function that represents an upper bound on the probability that an error is present in the chunk of bits located in interval $l$ for packet $k$.

I found the following implementation for $\rm BER$ calculations which takes $\rm SNR$ in mW. It is not clear to me whether the SNR given in the first equation above is for log scale or linear scale.

All in all, I would like to know how to calculate $P_\rm e$ and $P_{\rm err}$ for OOK?

It is not clear whether you want a derivation of the formula for the error probability of a wireless OOK system, or ideas about how the formula might apply to an optical communication system.

For a wireless OOK system, there are two main types of receivers called
coherent receivers and noncoherent receivers and these have different error probabilities. Neither error probability formula is suitable for application to an optical communication system.

Wireless OOK: The wireless OOK system model is that if a $1$ has to be be transmitted, the transmitted signal is an RF pulse of duration $T$ (the ON part of the keying) while if a $0$ has to be transmitted, then no signal is transmitted (the OFF part of the keying).  Assuming no signal distortion in the channel, then in the interval $0 \leq t \leq T$, the
receiver input can be expressed as
$$r(t) = \begin{cases} A \cos(2\pi f_c t+\theta) + n(t), & 1 ~~\text{transmitted},\\ n(t), & 0 ~~\text{transmitted}, \end{cases}
\tag{1}$$
where $n(t)$ is white Gaussian noise (with two-sided power spectral
density $\frac{N_0}{2}$) that is independent of the signal and
whether $1$ or $0$ is transmitted.  It is generally assumed that $1$ and $0$ are equally likely to be
transmitted.

Coherent receiver: A coherent receiver know the exact value of
the carrier frequency $f_c$, the received phase $\theta$ and the signal
timing, that is, it knows both the duration and location of the pulse
on the time axis. To understand how it knows all this will take us too
far away from this discussion. The receiver constructs the signal
$s(t) = \cos(2\pi f_ct+\theta)$ in a local oscillator and then correlates
this locally generated signal with $r(t)$, that is, it computes
$$X = \int_0^T r(t) s(t)\ \mathrm dt = \begin{cases}\frac{AT}{2} + N,
&1 ~~\text{transmitted},\\ N, & 0 ~~\text{transmitted},\end{cases}
\tag{2}$$
where $N$ is a zero-mean Gaussian random variable with variance
$\frac{N_0T}{4}$. The receiver compares $X$ with the threshold
$\Theta = \frac{AT}{4}$ and decides
at a $1$ or $0$ is transmitted according $X$ is larger than or smaller
than the threshold. The probability of an erroneous decision when a
$0$ is transmitted is $P\{N > \Theta\}$,
the probability that $N$ exceeds $\Theta$, while the probability of
an erroneous decision when a $1$ is transmitted is $P\{N < -\Theta\}$
which is the same as $P\{N > \Theta\}$
We readily get that
$$P_e = P\{N > \Theta\} = Q\left(\frac{AT}{4\sqrt{N_0T/4}}\right)
= Q\left(\sqrt{\frac{A^2T}{4N_0}}\right) \tag{3}$$
where $Q(z)$ is area under the right tail of the standard Gaussian
density to the right of $z$. It is, I believe, built in to MATLAB's
library; no need to fool around with erf or erfc.

It is possible to massage the right side of $(3)$ a little to give
a different perspective on the matter. Note
that $E_b$, the received signal energy per bit is $\frac{A^2T}{4}$
(because for half of the transmitted bits, there is no signal energy received).  Thus, we can re-write the right side of $(3)$ as
$$P_e = Q\left(\sqrt{\frac{E_b}{N_0}}\right)$$
which expression might be familiar to some readers.

Noncoherent receivers: Coherent receivers for OOK are difficult
to implement because the estimation of the RF phase $\theta$ is not
easy when the incoming signal (from which the phase is estimated)
keeps disappearing from time to time.  A noncoherent receiver is easier
to implement, and is sometimes referred to as a square law detector
because the earlier implementations merely measured the passband energy
and decided on $0$ or $1$ according as this energy is small or large.
More modern methods use I-Q signals, but regardless of the method, the
analysis is messy. The error probability expressions involve things like
Marcum's Q function which those knowledgeable about radar will likely
know. It is best to hide this stuff behind a veil.

Optical OOK communication:  This is completely different and the
system model is quite different. In particular, in contrast to
coherent wireless OOK communications, the probability of error
when a $1$ is transmitted is quite different from the error
probability when a $0$ is transmitted. Indeed, most of the errors are
mistaking $1$'s for $0$'s, and it is hardly ever the case that a
$0$ is mistaken for $1$. The reason is that the detector is essentially
a photon counter which can fail to detect a photon sometimes, but hardly
ever sees photons where none exist. As such, for the OP's simulation
of an optical communication system, the results from the analysis of
wireless RF communication systems are totally inapplicable.


How do I modulate complex time domain signal to a carrier frequency? I want to modulate a complex time domain signal to a carrier frequency 1880MHz but the sampling frequency is 30M which is smaller than Fs. I did upconversion using

exp(1i*2*pi*Fc*(0:length(signal)/Fs))

but the PSD isn't centered at Fc. I need to plot PSD and find the relation between the centered freq and the target Fc.

It is not clear why you would want to do this, and perhaps you are trying to achieve something at a bigger picture level where this is the wrong approach. However to clear up your question:

Since your sampling frequency is 30MHz, your unique digital spectrum is only from 0 to 30 MHz (or -15MHz to +15MHz depending on which portion of the spectrum you want to take). The rest of the spectrum, out to 1880MHz and beyond is replicas of this spectrum stated above.

1880 MHz/30MHz is 62.6666... so the boundary at 30*62 = 1860 MHz is identical to what exists at 0, and then from there up to 1860 + 30 MHz = 1890 MHz is identical to what exists at 30 MHz. So I suspect from your approach above that the result was a frequency at 1880-1860 = 20 MHz, which is identical to the frequency at 1880 MHz!

This is because the following frequencies are all aliases of the same when the sampling rate is 30 MHz:

20 MHz, 20+30*1=50 MHz, 20+30*2 = 80 MHz, 20+30*3 = 110 MHz .... all the way to 1880 MHz, which is 20 + 30*62.

The following post on the unique span for any digital spectrum may help to further clarify this result:

https://stackoverflow.com/questions/44429592/what-is-the-result-inside-fft-bin-at-different-index/44434979#44434979

UPDATE: From the comment you said that you "want to see multiple signals at the same frequencies on the same spectrum". And I assume your spectrum is centered at 1880 MHz, and you may be simulating a multicarrier system in the 1850 MHz to 1910 MHz band (my guess) such as CA_2-4 or PCS. Assuming that is the case, then your maximum bandwidth of the signal you want to simulate is 60 MHz, and to simulate that I would recommend simulating using complex signals at
a sampling rate of 100 MHz (greater than 2x the single-sided BW of 30 MHz to appease Nyquist with additional room for evaluating effects beyond the band such as spectral regrowth or simulating filtering there, etc).  In this case 1880 MHz can just as easily be represented at DC, and another signal at 1181.25 MHz (for example) can be placed at 1.25 MHz in your simulation. Why? Because the items of concern for us with regards to the carrier is the complex envelope of the carrier only; with regards to anything in our system simulation, we don't really care about the actual high frequency cycle rate of the carrier itself- only the magnitude and phase of the envelope versus time.

Alternatively, if you want to simulate with a real carrier, this is easily achievable too by centering the real signal at $F_s/4$ (where $F_s$ is your sampling rate) in your digital frequency span. The frequency span needed to simulate a real 60 MHz BW is at least 120 MHz (leaving no room at band edges), and for that case I would likely use 200 MHz which would be more than sufficient. This would result in the "1880 MHz carrier" being placed at 100 MHz.

Here is a close example I have that demonstrates this which was for a multi-carrier CDMA system (reverse link).

This is actually the radio receiver design but I draw your attention to the 35 MHz Digital IF, that is then represented with a 20 MHz sampling rate. In the end this is providing a real digital output (not a simulation) of what is at the 35 MHz carrier, which is also the same signal at the 1880 MHz carrier. The digital output is sampled at 20 MHz and centered at 5 MHz, which is $F_s/4$. The reason to center a real signal at or close to $F_s/4$ is that it maximizes the distance in frequency to the images, which in practical implementation such as what is shown, it simplifies the analog filtering required. In a simulation it will center your signal in the unique digital span.

This drawing may also help illustrate the point of not needing to simulate every cycle of the actual carrier:

A bandwidth modulated signal is shown in time and in frequency. For the case of time we see it as a vector moving around on the complex I Q plane versus time, changing in magnitude and phase (which is what we are interested in!). It moves relatively slowly due to being bandwidth limited (as opposed to white noise which would instantaneously be at any other arbitrary point at every next instant in time; we would just see a 2 dimensional Gaussian cloud in that case). In frequency we see the signal is centered at 0, and has some bandwidth as given by the rate of modulation and any filtering or similarly pulse shaping that was applied.

This next slide demonstrating frequency translation also helps to make the point I am trying to make here. On the left we see the same signal except at a carrier frequency of $F_o$.  The IQ time domain plot in this case is showing how the phase and magnitude of the signal is changing versus time relative to the carrier.

Consider the carrier alone with no modulation: When the carrier is at 0 as in our first case, the IQ plot would be a vector of constant magnitude and angle, not moving versus time (If the angle was 0 and magnitude 1 then it would be e vector on the real or I-axis centered on the origin pointing to the right). When the carrier is a positive frequency 1880 MHz signal (and to simplify this description consider the signal to still be complex), the IQ diagram would now be a vector rotating counter-clockwise at the rate of 1880 MHz.  If we could observe this with a strobe-light flashing at the same 1880 MHz rate synchronized to when the vector passes $1\angle 0$, we would then see the same case as when the carrier was at 0: a vector centered on the origin pointing to the right and not rotating.  Consider the same thought process with the modulated signal as in the picture below. When the signal is at the 1880 MHz carrier, the IQ diagram would be spinning at the 1880 MHz rate, but if we could observe it with the strobe light as previously described, we would see how the magnitude and phase deviates from 1880 MHz over time at each sample of the strobe light, and it would look identical to the same signal centered at 0. Since the signal is changing so slowly compared to 1880 MHz, there is no need for the strobe light to be going that high of a rate, we could equivalently run the strobe light at some sub-multiple of 1880 MHz (such as 100 MHz) and still observe the same result with the signal centered at 0; as at the 1880 MHz rate the signal has changed insignificantly in magnitude and phase from the previous sample.


Extended Kalman filter - why not use a 2nd order Taylor expansion? I'm implementing my first EKF and using a Taylor expansion to linearize the model. I'm somewhat new to this. Based on a paper I've read, I'm attempting to linearize the state transition matrix using a Taylor Expansion.

AFAIU, given a state vector X and state transition matrix F, we can write a linearized state transition matrix Φ (probably not the name) in the following manner:

$$
\textbf{F} = \frac{\partial \dot{\textbf{X}}}{\partial \textbf{X}}
$$
$$
\Phi = \textbf{I} + \textbf{F}\Delta t + \textbf{F}^2 \frac{\Delta t^2}{2!} + H.O.T
$$

Such that the state and covariance estimates are:

$$
\hat{\textbf{X}}_k = \Phi \hat{\textbf{X}}_{k-1}
$$

$$
\hat{\textbf{P}}_k = \Phi \hat{\textbf{P}}_{k-1} \Phi^T + \textbf{Q}
$$

Using a 2nd order expansion instead of a 1st, I've found that I get a much more accurate solution. It also allows me to use a larger timestep (minutes instead of seconds). I worry though that in doing so, I'm not actually linearizing, since the 2nd order TSE corresponds to a quadratic approximation.

Assuming the above steps are correct, is there anything wrong with taking the 2nd order TSE?

There is nothing wrong with taking a higher order approximation than first order.
An issue in doing so is that you no longer have Gaussian statistics, so some of the standard Kalman filter equations (which are derived using Gaussianity) no longer hold strictly true.
The usual approach to dealing with (otherwise ignored) higher order terms in the standard EKF is to increase either the process noise or the measurement noise variances from what the model tells you they should be. This is like saying,

I know what the noise power is, but I'm going to say the second and higher order terms in my Taylor series expansion add to it so I'm going to increase the noise variance I use in my filter to cater for it.


OFDM signal spectrum and characterisitics Is it correct to say that the spectrum of an OFDM signal is rectangular pulse?

What about the spectrum of each transmitted signal over subcarriers?

If the OFDM signal has a rectangular frequency domain spectrum, then the time domain after IDFT is much like a sinc function. Do we need pulse shaping in this case, since sinc function already looks like pulse shape?

Your statement isn't strictly correct; the spectrum of an OFDM signal isn't a rectangular pulse. Think about it this way: OFDM is really just a method to efficiently modulate multiple low-rate narrowband signals that are sent simultaneously. Therefore, the OFDM signal's spectrum is really just the sum of the spectra of each of its narrowband components.

For example, if your OFDM signal has 64 BPSK subcarriers, then its spectrum is equal to the spectra of 64 narrowband BPSK signals, spaced at the OFDM subcarrier spacing. If you were to estimate the OFDM signal's power spectrum with a wide enough resolution bandwidth, then the spectrum would appear to be rectangular. However, if you look closely enough at the signal's spectrum, you'll notice that it isn't flat.

Here is an example MATLAB script that demonstrates the phenomenon:

% generate random BPSK symbol vector; this is equal to the "low-resolution"
% spectrum, calculated with bin spacing synchronized with the OFDM
% subcarrier grid
X = 1 - 2 * (randn(64,1) > 0.5);
f = -0.5 + (0:length(X)-1) / 64;
% create time-domain OFDM symbol using the IFFT
x = ifft(X);
% calculate interpolated spectrum to see what it looks like between OFDM
% subcarrier centers
Sh = fftshift(fft(x,4096));
fh = -0.5 + (0:length(Sh)-1) / 4096;

figure;
plot(f,abs(X)); hold all; grid on;
plot(fh,abs(Sh));
legend('Regular', 'Interpolated');

As the modulated OFDM signal doesn't have a perfectly rectangular spectrum as you guessed, it does not necessarily have a sinc-like shape in the time domain either.


Scale invariant template matching For my application I am looking for a free, fast and scale-invariant template matching algorithm. The templates can only be scaled within the scene (no skew or rotation) and the scale is unknown. Scale can go up from 400% to 25%. All pictures are noise-free.

I have considered the following:

SIFT: Highest accuracy so far. However it is licensed and pretty slow (3sec)
ORB: Performance still not fast enough (1.5sec). Accuracy also poor on scaling.

What is the best algorithm in terms of accuracy for this particular problem? And which is the quickest? Are there any algorithms specially designed for this particular problem?

To increase the speed, I think the dimensionality of the feature descriptors should be reduced. SIFT feature descriptor is of 128 dimensions.

So you can try PCA-SIFT which reduces the descriptor number to 20. It uses Principal Component Analysis to the SIFT descriptors.

You can also try using Speeded Up Robust Features (SURF). The feature descriptors obtained from SURF are invariant to rotation and scale, and it can reduce the dimensions of feature descriptors significantly. it enables faster feature detection as well as template matching. Also, MATLAB has a specific command for SURF i.e. detectSURFFeatures, so if you are using MATLAB for your job you are in luck. I think OpenCV also has the SURF feature command so that you can use it on python too. I have used the SURF features which I think is pretty fast.


High performance DSP processor with better on-chip memory I am currently using a combination of 4 ADSP 21469 because of the complexity and huge chunk of code that I currently have. BOOT and application will execute from external Flash memory. The code for radar processing are equally distributed across 4 DSPs

Even after splitting the code across 4 DSPs, I still need to perform code optimization because of memory constraint. By doing a certain level of optimization, my variables are not accessible due to which I have to read from the registers, which vary all the time.

My thought was to use a single advanced processor that is robust enough to contain all the code so I can use common libraries wherever I can to make my code efficient on target.

To summarize, the bottleneck is that ADSP 21469 has 5 Mbits of on-chip RAM, 4 Mbits of on-chip ROM which is way less than my requirement. If I can use a multicore processor(Ex:quad code) that works at >=400 Mhz, which has an on-chip memory of >=20 Mb to hold all my data fr a faster access or supporting external memory without being slower in accessing the data.

That depends a bit on the constraints: are you out of code space or data space or both?

Anyway, the typical way of dealing with space problems is to use external memory in conjunction with a fast access internal cache. An alternative would be something like this http://www.ti.com/lsds/ti/processors/dsp/c6000_dsp/c66x/overview.page with a few gigs of DDR attached.

Obviously external memory is somewhat slower than internal memory, so the trick is to move data efficiently into internal memory, so that it's there when you need it. In many cases the processor will do a pretty decent job in managing this automatically, especially if it's mostly code (not data) and most MIPS consumption is in loops.

If that's not good enough, there are certainly strategies to manually optimize code for cache efficiency as well.


What makes an image feature stable? Feature detection is an essential task in low-level vision.

Good features are those that resist to different perturbations such as

noise addition,
blur,
geometric transforms (3D rotation with perspective, scaling),
radiometric transforms (monotonic or non-monotonic grayscale remapping).

By resistance, I mean that despite those alterations, the same features will continue to be detected.

My question is how can we mathematically express these types of resistance in general terms ? What would be the form of the equations, assuming grayscale images and features defined from the pixel values in a region of arbitrary shape ? (I am not asking for an assessment of the stability of known features described in the literature.)

Answering the question for a 1D signal would already be a good start.

Can anyone help ?

Without being an expert in the image processing field there is something that comes to mind that maybe could point you in the right direction. Usually we express a feature as a number, to which we assign an estimator, which is a random variable, with some probability distribution (for which we normally care about mean value and variance). If we take the noise example which i think is the easier. You could try to express this resistance in term of the variation of the probability distribution of the feature estimator given the perturbing noise source. Let the output of your estimator for the feature be the random variable $\mathbf{X}$, and let the noise be characterized by $\mathbf{Y}$. Some measure of the resistance could be $\frac{P(\mathbf{X})}{P(\mathbf{X}|\mathbf{Y})}$. If this value is close to one, then the feature would be resistent to perturbations by noise, if it is not then it won't be. Maybe you could define a function of this to adjust the values to be more significant. For the other transformations such as rotation, scaling etc, you could use the same but do $\frac{P(\mathbf{X})}{P(\mathbf{f(X)})}$, where $f(\mathbf{X})$ denotes the transformation of the random variable that characterizes the estimator when you apply the rotation to the image. Of course, I think this $f()$ would be very hard to find in general but perhaps if the feature is not hard to estimate you can do some math and some change of variable and arrive at closed expressions for it.

Hope this helps.


Calculating power in frequency band using power spectral density In a text I am using it states:

  The fourier transform $F(\omega)$ has both positive and negative
  frequencies and $F(\omega) = F^*(-\omega)$. The power spectral
  density $S(\omega)$ is similarly defined to have the positive and
  negative frequencies, and the power spectral density is similarly
  equally divided between the positive and negative frequencies. Thus
  power spectrum $S(\omega)$ is symmetric about $\omega = 0$. So the
  contribution to the total power that is contained in the band
  $[\omega_1, \omega_2]$ is $$P_{[\omega_1, \omega_2]} = \frac{1}{2 \pi} \int_{\omega_1}^{\omega_2} [S(\omega) + S(-\omega)]d \omega.$$

Question: If $S(\omega)$ is the power per frequency $\omega$, even though $\omega$ takes both positive and negative values, I still don't see why the power contained in the band $[\omega_1, \omega_2]$ is not simply given by $$\frac{1}{2 \pi}\int_{\omega_1}^{\omega_2} S(\omega) d \omega?$$ I am only beginning to learn about these things, hence there might be something trivial that I am missing. Thanks.

It's less about power spectrum and more about Fourier transform of real signals. Let's take an easy signal, $g(x) = \cos(10\pi x)$.

Then:
$$ F(\omega) = \frac{\delta(\omega -10) + \delta (\omega +10)}2.$$

When you ask "what is the energy in the band $ \omega = [9,11]$ you ask what is the energy of that cosine (of course this is an example).

If you look only at the positive side of the frequency you lose half the energy.


Why is the close looped model under control called plant? I am not sure if the question even belongs here. Please correct me if it is out of context.

Consider a closed feedback loop.

Usually, we say that $K^L(s)$ is a controller over the plat $G(s)$. Is there any reason it is called a plat?

The controlled system may be a car with cruise control, an actual plant, or even the famous F15. These are not necessarily plants. Was the theory originating from real plants and factories?

I am looking for something anchored in resources and not an intuitive explanation.

I have heard that the word comes from "plant" as in "steam plant".  But I don't know of any historical research that's been done to pin down the term.
I did just check my copy of "On Governors" by James Clerk Maxwell himself (that guy got around).  At least in 1868, Maxwell wasn't using the word "plant" for "thing being controlled".


Shelving and Peak filters name Why are Peak filters called "Peak" and why are Shelving filters called "Shelving". I am not a native English speaker. So excuse me if it is so dumb question.

This is probably the craziest, though important question on DSP SE ;) Answer is, because they:

look like peak:

and tend to have a shelf either in the beginning or the end of frequency range:


Image datasets used for salient objects? I would like to understand reasons for different image datasets used for salient object detection: their difficulties, their origins and likely future datasets. Does there exist any work explaining them such their progression and characteristics over time?

Datasets explicitly mentioning saliency

MSRA
SOD based on BSD Berkeley segmentation dataset

Apparently more general datasets (not researched in detail)

DUT-Omron Image dataset
THUS
ECSD
PASCAL
...

More challenging situations with more noise

Detection on magnetism such as MR: Medical image dataset and Extract voice from magnetic tape


Sampling and the creation of multiple images at integer multiples of the sampling frequency From what I thought I learned about sampling is that when you draw a sampled signal in the frequency domain, you can show copies or "images" of that signal at integer multiples of the sampling frequency. Is this true?

If so, I'm having a hard time understanding this conceptually. What does it really mean to make a copy at a different location on the frequency line?

For example, if I have a signal which when drawn in the frequency domain looks like a box centered at zero extending from $f_{low}$ to $f_{high}$, can I now say I can draw this box centered at integer multiples of the sampling frequency and so long as my box when redrawn does not touch another box I avoid aliasing? Why do I do this at integer multiples? For example if my sampling frequency is $20 \ \mathrm{Hz}$, what does the next integer multiple at $40 \ \mathrm{Hz}$ tell me? Does it mean I'll see something interesting if the signal were to run at $40 \ \mathrm{Hz}$?

I'm trying to get a better understanding of this so I can then understand bandpass sampling and what it means to select a location at $k$, $k-1$, etc.

Below is an example of two sine waves at 11 Hz and 1 Hz, both sampled with a 10 Hz sampling clock.  The sampled locations are drawn as circles. If you only had the samples, and had no knowledge of any filtering prior to the sampling process, you would not be able to know for certain if the samples came from a 1 Hz analog signal or an 11 Hz analog signal. In fact, any signal at N fs ± 1 Hz for any integer N would produce the same result. (9 Hz, 11 Hz, 19 Hz, 21 Hz, etc).

Thus, with appropriate analog bandpass filtering prior the sampling, we can directly sample ("undersample") a bandpass signal, as long as we follow Nyquist's criteria that the sampling rate be > 2x the signal bandwidth. Do not confuse bandwidth with highest frequency content!  The sine wave shown for example has 0 bandwidth in that it occupies an infinitely thin line in the spectrum. Actual bandpass waveforms of interest would have some bandwidth based on their modulation content.

What this picture also helps clarify is the increased sensitivity to any sampling clock jitter when performing bandpass sampling.  Note specifically the slope of the red lower frequency signal as compared to that of the blue higher frequency signal in vicinity of the sample locations, and imagine small variations due to clock jitter and how these variations would transform to amplitude error.

Please see this answer for specific details related to a conceptual understanding of having multiple frequency "images": Aliasing after downsampling