2-d circularly symmetric low-pass filter For a square pixel grid, the ideal 2-d low-pass filter with a horizontal and a vertical cut-off angular frequency $\omega_c$ in radians has an impulse response (kernel) $h_{\small\square}(x, y)$ that is the product of a horizontal and a vertical stretched and scaled sinc function, with $x$ and $y$ the integer horizontal and vertical pixel coordinates:

$$h_{\small\square}[x, y] = \frac{\omega_c\operatorname{sinc}\left(\frac{\omega_c x}{\pi}\right)}{\pi}\times\frac{\omega_c\operatorname{sinc}\left(\frac{\omega_c y}{\pi}\right)}{\pi} = \begin{cases}\frac{\omega_c^2}{\pi^2}&\text{if }x = y = 0,\\\frac{\sin(\omega_c x)\sin(\omega_c y)}{\pi^2 x y}&\text{otherwise.}\end{cases}\tag{1}$$

If $\omega_c = \pi$, the kernel of Eq. 1 is simply:

$$h_{\small\square} = [1]\quad \text{if}\quad\omega_c = \pi.\tag{2}$$

For real-valued $x$ and $y$, Eq. 1 would not describe a circularly symmetric kernel, and neither is the square-like frequency response $H_{\small\square}$ circularly symmetric. Sometimes an isotropic kernel may be desired, for example if we want to remove the effect of the choice of the direction of the image axes.

What is a kernel $h_\circ[x, y]$ such that the frequency response $H_{\circ}$ is $1$ inside a circle of radius $\omega_c$ (angular frequency in radians) and $0$ outside it?

Excerpted from Jae S.Lim 2D signal and image processing ch.1, as an example of $2$-D circularly symmetric lowpass filter with a cutoff frequency of $\omega_c$ radians per sample, whose impulse response is given by:
$$h[n_1,n_2] = \frac{\omega_c}{2\pi \sqrt{n_1^2 + n_2^2} } J_1 \big( \omega_c \sqrt{n_1^2 + n_2^2} \big) $$

where $J_1$ is the Bessel function of the first kind and the first order...

Interested readers may consult the book for a derivation which is not untitive but nevertheless tractable; familiarity with Bessel functions is required but also provided as is.[Derivation is added.]

[Olli below]

At $n_1 = n_2 = 0$ the limiting value must be used:

$$h[0, 0] = \frac{\omega_c^2}{4\pi}$$

A slice through the middle of $h[n_1,n_2]$ with $\omega_c = \pi$:

Python source for 2-d filter kernel (you may want to apply a 2-d window function):

from scipy import special
import numpy as np

def circularLowpassKernel(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 x, y: omega_c*special.j1(omega_c*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2))/(2*np.pi*np.sqrt((x - (N - 1)/2)**2 + (y - (N - 1)/2)**2)), [N, N])
  kernel[(N - 1)//2, (N - 1)//2] = omega_c**2/(4*np.pi)
  return kernel

Example with $\omega_c = \pi$:

import matplotlib.pyplot as plt

kernelN = 11  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-1, vmax=1, cmap='bwr')
plt.colorbar()
plt.show()

Example with $\omega_c = \pi/4$:

kernelN = 41  # Horizontal size of the kernel, also its vertical size. Must be odd.
omega_c = np.pi/4  # Cutoff frequency in radians <= pi
kernel = circularLowpassKernel(omega_c, kernelN)
plt.imshow(kernel, vmin=-np.max(kernel), vmax=np.max(kernel), cmap='bwr')
plt.colorbar()
plt.show()

[fat32 below]

Ok, for those who would benefit from a (not so intuitive) derivation, here I make a (n almost) verbatim copy of it from the same book.

First, lets simply write the 2D discrete-time inverse Fourier transform to define the impulse response as:

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} {\int \int}_{\omega_1^2+\omega_2^2< w_c^2} 1 \cdot e^{j(\omega_1 n_1 + \omega_2 n_2)} d\omega_1 d\omega_2 \tag{1} $$

Let's make a change of variables $\omega_1 = r \cos(\theta)$ and $\omega_2 = r \sin(\theta)$ (effectively write the integral in (1) in the cyclindrical (or polar) coordinates) :

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} \int_{\theta=a}^{a+2\pi} e^{j r (\cos(\theta) n_1 + \sin(\theta) n_2)} r ~ dr d\theta \tag{2} $$

Now, make a further change of variables $n_1 = n \cos(\phi)$ and $n_2 = n \sin(\phi)$, with $n = \sqrt{ n_1^2 + n_2^2 }$ and obtain

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r ~dr  \int_{\theta=a}^{a+2\pi} e^{j r n \big( \cos(\theta) \cos(\phi) + \sin(\theta) \sin(\phi) \big) } d\theta \tag{3} \\$$

which is now:
$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r ~dr  \int_{\theta=a}^{a+2\pi} e^{j r n \cos(\theta -\phi) } d\theta \tag{4} \\$$

defining $f(r) = \int_{\theta=a}^{a+2\pi} e^{ j r n \cos(\theta-\phi) } d\theta $ , then we get:

$$ h[n_1, n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r f(r) dr  \tag{5} $$

Now, exploring $f(r)$ we can apply the Euler's identity :

$$ f(r) = \int_{\theta=a}^{a+2\pi} \cos(r.n.\cos(\theta-\phi)) d\theta + j \int_{\theta=a}^{a+2\pi} \sin(r.n.\cos(\theta-\phi) ) d\theta \tag{6} $$

And we can notice that the imaginary part is zero (can check an integral table for that) and with $a =  \phi$  then $f(r)$ becomes

$$f(r) = \int_{\theta=0}^{2\pi} \cos(r.n.\cos(\theta) ) d\theta  \tag{7}$$

Now the integral in (7) is recognized as the the Bessel function of order zero, kind one $J_0(x)$ which is given as:

$$ J_0(x) = \frac{1}{2\pi} \int_{\theta =0}^{2\pi} \cos( x \cos( \theta) ) d\theta \tag{8} $$

from (7) and (8) we see that $f(r) = 2\pi J_0(r n) $...

And the last identity is given as is:
$$ x J_1(x) |_a^b = \int_a^b x J_0(x) dx \tag{9}$$

Now putting the impulse response into the form

$$h[n_1,n_2] = \frac{1}{(2\pi)^2} \int_{r=0}^{\omega_c} r 2\pi J_0( r n) dr \tag{10}$$

applying (9) onto (10) with the recognition that $x = r n$ and $dr = dx/n$ and $n= \sqrt{n_1^2+n_2^2}$ yields the result:

$$\boxed{ h[n_1,n_2] = \frac{\omega_c}{2\pi \sqrt{n_1^2+n_2^2}} J_1( \omega_c \sqrt{n_1^2+n_2^2}) }  \tag{11}$$


Rectangular or trapezoidal signal filtering. Beginner I have the following signal:

I know, that there are 2 states ('High' and 'Low') in my real signal. I need to convert my signal to the sequence of these states (Low from 2 to 5 second, High from 5.5 to 7 and so on).
Could you give me a piece of advice about classic and well established method to extract such signal from mine?
Median filter was as an example. Is it the most classic way to do it? Or you may provide some information for a better approach?
I don't like the median filter approach, because I should rather voluntary define kernel and 2 thresholds. I am searching for a method, with more natural parameters.
UPDATE:
Some knowledge about signal:
Amplitude and mean duration of the state may vary from approx. 1 to 3 fold. But I have some knowledge about CV and skewness/CV. CV~0.5-0.9, skewness/CV~2. Distribution of States duration is gamma-like.

Depending on your noise model, a median filter might be quite good or not. Have a look at this example:

import scipy.signal

# from http://stackoverflow.com/questions/23289976/how-to-find-zero-crossings-with-hysteresis
def hyst(x, th_lo, th_hi, initial = False):
    hi = x >= th_hi
    lo_or_hi = (x <= th_lo) | hi
    ind = np.nonzero(lo_or_hi)[0]
    if not ind.size: # prevent index error if ind is empty
        return np.zeros_like(x, dtype=bool) | initial
    cnt = np.cumsum(lo_or_hi) # from 0 to len(x)
    return np.where(cnt, hi[ind[cnt-1]], initial)

Fs = 300
T = 9
a = np.zeros(T*Fs)
# setup the state changes
a[0*Fs] = 1
a[0.5*Fs] = -2
a[2*Fs] = 2
a[3*Fs] = -2
a[4.5*Fs] = 2
a[7*Fs] = -2
a[8*Fs] = 2
a[8.5*Fs] = -2
# create the rect-waveform
data = np.cumsum(a)

# add some noise to the signal
rx = data + 2*np.random.randn(len(data))

plt.subplot(121)
plt.plot(rx)
plt.plot(data, 'r', lw=2)

plt.subplot(122)
# plot the filtered signal (if you have access to the purple signal in your plot, just use this signal)
filtered = scipy.signal.medfilt(rx, Fs//2+1)
plt.plot(filtered, lw=3)

# Perform thresholding with Hysteresis
detected = 2*hyst(filtered, -0.3, 0.3)-1
plt.plot(detected)

# From the thresholds, get the state changes via the derivative
diffed = np.diff(detected)
print ("estimated times of state switch:", ["%.2f" % (x/Fs) for x in np.nonzero(diffed)[0]])
print ("actual times of state switch:", ["%.2f" % (x/Fs) for x in np.nonzero(a)[0]])

estimated times of state switch: ['0.18', '0.60', '2.05', '3.10', '4.56', '7.07', '8.07', '8.59']
actual times of state switch: ['0.00', '0.50', '2.00', '3.00', '4.50', '7.00', '8.00', '8.50']

It creates an artificial signal, adds some noise, does a recovery of the original signal via median filtering. Eventually, it detects the state changes in the signal.


Power or Energy Signal? Can you tell me whether the following signal is a energy or power signal?
$$
x(n) = e^{j (n \pi/2 + \pi/8) }
$$
I've solved it and found it as it was neither power nor energy signal was that right?

The energy of a discrete-time signal is defined as

$$E_x=\sum_{n=-\infty}^{\infty}|x(n)|^2\tag{1}$$

and its power is given by

$$P_x=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{n=-N}^{N}|x(n)|^2$$

With $x(n) = e^{j (n \pi/2 + \pi/8) }$ we have $|x(n)|^2=1$ which implies that the sum in (1) does not converge, i.e. $x(n)$ has infinite energy. For the power we get

$$P_x=\lim_{N\rightarrow\infty}\frac{2N+1}{2N+1}=1$$

The signal $x(n)$ has finite power and is consequently a "power signal".


Kalman Filter - Updating the Covariance Matrix Step I am trying to simulate the Kalman Filter. I have the covariance matrix P_{0|0}. Tell me please, how can I get the predicted (a priori) estimate covariance matrix on the (k-1) step? P_{k-1|k-1}?

Given the Covariance Matrix $ {P}_{k - 1 \mid k - 1} $ then:

$$ {P}_{k \mid k} = {F}_{k} {P}_{k - 1 \mid k - 1} {F}_{k}^{T} + {Q}_{k} $$

Where $ {F}_{k} $ is the Model Matrix at iteration $ k $ and $ {Q}_{k} $ is the Process Noise Covariance Matrix at iteration $ k $.


Can one generate a signal with an asymmetric power spectral density? My question is as follows. I know that classically speaking, if you have some time varying signal $X(t)$ its power spectral density $S_X(f)$ will be symmetric around $f=0$ as the signal is real. This is why your spectrum analyzer only shows the positive part, the negative part is redundant.

But if one has a complex signal (such as often encountered in quantum mechanics) the PSD can be asymmetric around 0. This for example leads to the asymmetry between absorption and emission in atoms: atoms have spontaneous emission (emission in the absence of photons in the environment) but not spontaneous absorption.

What I am interested in is if there is a way to physically make such a complex signal with asymmetric PSD. Maybe this can be done by using two channels that output the signal with some phase difference? I've tried looking for some source material on this but I have not been successful.

@user129412, there are many ways to generate both continuous signals and discrete sequence having asymmetric PSDs. One way, for discrete sequences, is shown in the following diagram.

$x_r(n)$ is a narrow bandpass signal having spectral energy in the vicinity of $±f_c$ Hz. The complex-valued $x_c(n)$ output sequence is your desired asymmetrical sequence having spectral energy only in the vicinity of $+f_c$ Hz. (If your Hilbert transformer has M taps, where M is odd, then the upper path delay must be (M-1)/2 samples to achieve time synchronization between the two signal paths.)  If you desire sequences whose asymmetrical spectra are centered at zero Hz, many ways to obtain such sequences can be found at:

https://www.dsprelated.com/showarticle/153.php


Absolute Integrable Sinc function How do we prove that $$\int_{-\infty}^{\infty}\bigg|\dfrac{\sin t}{\pi t}\bigg|dt\to \infty$$
This comes in the context of stability of LTI system with impulse response $h(t) = \dfrac{\sin t}{\pi t}$.

Thanks in advance.

The question has an answer on Math SE. The argument goes as follows:

\begin{eqnarray}
\int_{-\infty}^\infty \left| \frac{\sin t}{\pi t} \right|dt &=& \frac{1}{\pi}\sum_{n=-\infty}^\infty \int_{n\pi}^{(n+1)\pi} \left| \frac{\sin t}{t} \right| dt \\
&>& \frac{1}{\pi} \sum_{n=-\infty}^\infty \frac{1}{|n+1|\pi}\int_{n\pi}^{(n+1)\pi} |\sin t|dt \\
&=& \frac{2}{\pi^2} \sum_{n=-\infty}^\infty \frac{1}{|n|}\\
&\rightarrow& \infty
\end{eqnarray}
where in each interval $[n\pi, (n+1)\pi]$ we lower bound $1/t$ with $1/(n+1) \pi$.


Kalman filtering for position using GPS,accelerometer and speed sensors I am working on tracking a vehicle under tunnel when GPS is lost.
Whenever the vehicle in on the road, the GPS works fine and gives good accuracy but when the vehicle is under tunnel, the GPS is lost and its difficult to track vehicle. This is where i have decided to use kalman filtering.
I have a 9 axis IMU sensor(accelerometer,gyro,magnetometer) and speed value from Candata and would like to predict the location using kalman
So far, this is wat i have done

1. Px(t+1) = Px + delta_t * vx + 0.5 * ax * delta_t 2
2. Py(t+1) = Py + delta_t * vy + 0.5 * ay * delta_t 2
3. Vx(t+1) = Vx + ax * delta_t
4. Vy(t+1) = Vy + ay * delta_t

where px,py are my positions and vx and vy are my velocities this is my statemodel
my input is the abs acceleration in x and y direction calculated from 9 axis IMU using all 3 sensors.
and in my update step i use the gps value i received.
So my question here is, where do i use my speed value? Is accelerometer enough?

In my experience at the very least you need a nine or ten state Kalman (velocity, position, and angle -- that's nine states if you attempt to use Euler angles, ten if you use quaternions).  Practically, you need to use IMU offset, so add six more states.  If you use the IMU's compass function, add three more.  It may not be a bad idea to monitor the acceleration due to gravity -- one more.  If you use the vehicle's odometer function (from CAN), then add a state for proportionality error, and possibly offset error.
You need vehicle orientation.  That's because the ground is accelerating upwards at $9.81 \mathrm{m/s^2}$ (if Einstein is to be believed), or the accelerometer's proof mass is being pulled down with a force equal to $9.81 \mathrm{m/s^2/kg}$ (if Sir Isaac Newton is to be believed).  Either way, to get the actual acceleration with respect to the ground you need to null out the acceleration with respect to gravity, and you can only do that if you have the vehicle's orientation in 3D.
It's been over a decade since I worked on this last, but there has to be canned packages available open source -- have you done a search?


LTE questions (subcarriers and cyclic prefix) Assume we have transmission bandwidth of 10 MHz, in this case the total number of subcarriers (FFT size) present is 1024. According to LTE standards, 600 out of 1024 subcarriers are data subcarriers.

data subcarriers= total number of subcarriers- pilot/null subcarriers

Moreover, 0.07 % of 1024 are used as Cyclic Prefix extension

1) Is it correct that we have 1024-600 are therefore unused subcarriers?
(where 600 = total number of Resource Blocks (RB) * number of subcarriers/RB = 50*12)

2) An OFDM symbol is made up of data part and CP part. Is it correct to say that the number of QAM symbols for each OFDM symbol is then 600 (equal to number of data subcarriers) and the number of CP symbols are 72? How is the actual transmission over the subcarriers taking place. For example in downlink, does the transmitter allocate a certain number of subcarriers over which it transmits?

I hope I can get some feedback! Thanks

The FFT size in LTE is defined as 2^(nextpow2(Nsubcarriers/2)+1) or the nearest power of two that is also higher than the number of subcarriers for a given bandwidth. For the 10 Mhz case, for example, this means that 424 FFT bins(1024-600) are zero padded apart including the DC bin. The usual practice is to zero pad the central 423+DC bins and use the rest for filling the data so that you get a continous spectrum of 10 MHz around DC. You can also fill the central 600 FFT bins with data and do a fftshift. See table 1 in this document.
Yes. Each Resource Element(i.e.1 subcarrier and 1 OFDM system in the Time- Frequency Grid) has one QAM Symbol in shared channel.

Check out this visualization to get a better clarity about the LTE resource grid. The resource is allocated both in time and frequency domain and in the case of Multi User MIMO(MU-MIMO) transmission mode, spatial resource allocation of resource elements to different user equipment in the downlink also take place. To keep things simple, let's just focus on a single-antennae transmission. The following diagram will give you a better idea of the whole transmission process.


MATLAB: Generating and Visualizing Mother Wavelet of different scale In MATLAB, there is a function wavefun that outputs both the scaling and wavelet functions. To generate a Daubechies, the following command can be used:

[phi,psi,xval] = wavefun('db40',10);

The question I have is: How to generate a Daubechies at different scales and visualizing them? For instance, generating a similar Daubechies with frequency twice larger.

Subsequent wavelets, spanning a specific range of scales are generated by shifting and scaling the Mother Wavelet.

Therefore, to produce "...a similar Daubechies with frequency twice larger..." you would have to "take every other sample" of the Mother Wavelet. For more information, please see this link and the Wavelet Toolbox User Guide on the section about the discrete wavelet transform.

It may also be useful to have a quick look at decimation and upsampling concepts in DSP.

Hope this helps.


Generating spectrogram for band-pass filtered wav file python I have some .wav files of heart sounds. I want to apply band-pass filter [20Hx 100 Hz]. Then generating  PCG images size of 560x420 pixels and RGB format. How could i do this?
Any help will be highly appreciated.

Using scipy you can easily:

bandpass signals, for instance with scipy.signal.butter
obtain a spectrogram, using for instance scipy.signal.stft (or LTFAT, The Large Time-Frequency Toolbox (LTFAT) in Python)
scale the values of this 2D array so that each channel has values between $0$ and $1$
reshape the array into a $560\times 420$ grid, eg with scipy.ndimage.affine_transform or scipy.ndimage.geometric_transform
save the image, but I don't know the PCG format


What is the DFT of $[x_1 -x_2 x_3 -x_4....x_n]$ If DFT of $[x_1 x_2 x_3... x_N]$ is $Y(k)$, what is the DFT of $[x_1 -x_2 x_3 -x_4,....x_N]$ in terms of $Y(k)$?
I have tried to formulate it but I cannot get a simplified expression for DFT of modified sequence in terms of $Y(k)$.

Hint: Consider a factor 2 down-sampled sequence : $\tilde{x} = \{x_1, x_3, x_5, ...\}$ and consider another sequence which is circular shifted by $-1$ sample and downsampled by a factor of 2 : $\hat{x} = \{ x_2, x_4, x_6, ...\}$. Observe that the required sequence is difference of $\tilde{x}$ and $\hat{x}$. Each of the DFT can be computed by observation using properties.


Converting poles/zeros to differential/difference equation solutions Does anyone have a reference handy on how to convert poles and zeros of a system to differential/difference equations.

Here is a quick draft of math, but I am not sure if it's at all correct. First, solving the characteristic equation gives the eigen values (equal to poles). Applying rudimentary knowledge of differential equations, the solution regarding only the poles should be:

$$\text {Poles Diffrential}: p(t)= \sum_{i=1}^{n_1} c_ie^{t\times \text{p}_i} $$

$$\text {Poles Difference}:p[n]= \sum_{i=1}^{n_1} c_i\text{p}_i^n $$

Where the $p_i$ is the value of ith pole and C is some constant, possibly scaled to the inverse value of the number of poles,  the index runs from 1 to $n_1$, the number of poles. I denote the solution that considers only poles/characteristic equation as $p(t)$. (To be exact, c is the input/initial value at time zero so we still have the denominator of characteristic equation as cz^0=c).

Next, zeros should be solved in respect to time/samples. This should be an easy task, considering that the input equation which zeros are solved from, directly outputs the signal. First solve the Y(z) (equivalent to the numerator of the transfer function):

$$Y[z]=\prod_{i=1}^{n_2}(z-q_i)$$

Where $q_i$ is the value of ith zero and $n_2$ the number of zeros. The function needs to still be converted to time doimain by converting each individual elements $az^{x}$ into $z[n(-x)]=a$. Things get more complicated when both elements are considered simultaneously. Unless my interpretation is incorrect, the input equation (used for zeros) acts as a signal that is fed into the recursive filter(poles portion). Thus the final interpretation can be had by simply convolving the equations (with $c_i=1$):

$$ \text {Zeros, Poles}: y[n] = p[n] * z[n] $$

Where the $p[n]$ and $z[n]$ from the preceding equations and $*$ is the convolution operator. To solve the differential equation variation, simply substitute $t$ in place of $n$ for zeros and the output.

Okay, I'm gonna attempt an answer to this while I am still convinced the question is ill-formed.  I am answering the title of the question without really referring to the text of the question, which I don't understand.

First, for discrete-time:  A general IIR (or even FIR in which all of the poles are at the origin, where $z=0$) transfer function is:

$$\begin{align}
 H(z) \triangleq \frac{Y(z)}{X(z)} & = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_M)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \\
\\
 & = H_0 \frac{z^M (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{z^N (1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{H_0 \ (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})}{(1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = z^{M-N} \ \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} }
\end{align}$$

$N$ is the number of poles, $p_i$, and $M$ is the number of zeros, $q_i$.

The difference equation comes from moving a few of the above things around.

$$\begin{align} Y(z) & \left(  1 + a_1 z^{-1} + a_2 z^{-2} + ... + a_N z^{-N}  \right) \\
& \quad = X(z) z^{M-N} \left( b_0 + b_1 z^{-1} + b_2 z^{-2} + ... + b_M z^{-M} \right) \end{align}$$

$$\begin{align} Y(z) = \ & b_0 z^{M-N} X(z) + b_1 z^{M-N-1} X(z) + b_2 z^{M-N-2} X(z) + ...  \\
 & \quad ... + b_M z^{-N} X(z) - a_1 z^{-1}Y(z) - a_2 z^{-2}Y(z) - ... - a_N z^{-N}Y(z) \end{align}$$

Inverting the Z Transform, this becomes

$$\begin{align} y[n] = \ & b_0 \ x[n+M-N] + b_1 \ x[n+M-N-1] + b_2 \ x[n+M-N-2] + ...  \\
 & \quad ... + b_M \ x[n-N] - a_1 \ y[n-1] - a_2 \ y[n-2] - ... - a_N \ y[n-N] \end{align}$$

For a causal system or filter, we cannot look into the future for any input sample $x[n]$.  That means $M \le N$ so that there is no negative delays on the $x[n+M-N]$ term.

Often $M=N$ (the numbers of poles and zeros are the same) and you get the more familiar

$$\begin{align} y[n] = \ & b_0 \ x[n] + b_1 \ x[n-1] + b_2 \ x[n-2] + ... + b_N \ x[n-N] \\
 & \quad - a_1 \ y[n-1] - a_2 \ y[n-2] - ... - a_N \ y[n-N] \end{align}$$

and the more familiar transfer function

$$\begin{align}
 H(z) & = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_N)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \\
\\
 & = \frac{H_0 \ (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_N z^{-1})}{(1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})}  \\
\\
 & = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_N z^{-N} }{ 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... +  + a_N z^{-N} }
\end{align}$$

That demonstrates where the difference equation (and the feedforward, $b_i$, and feedback, $-a_i$, coefficients) come from the poles, $p_i$, and zeros, $q_i$.

The problem remaining, which is the hard problem, is that I didn't bother to multiply it out.  So you have to blast this out:

$$\begin{align}
 H_0 \ & (1-q_1 z^{-1})(1-q_2 z^{-1})(1-q_3 z^{-1})...(1-q_M z^{-1})  \\
 & \quad = b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + ... + b_M z^{-M}
\end{align}$$

and

$$\begin{align}
  & (1-p_1 z^{-1})(1-p_2 z^{-1})(1-p_3 z^{-1})...(1-p_N z^{-1})  \\
 & \quad = 1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + ... + a_N z^{-N}
\end{align}$$

You can see right away that $b_0 = H_0$ and that $a_0=1$ (I never really had an "$a_0$" term anyway), but you have to blast out the multiplication of the binomials to get the polynomials.  You can also see that $b_M = (-1)^M H_0 q_1 q_2 q_3 ... q_M$ and $a_N = (-1)^N p_1 p_2 p_3 ... p_N$   That's not a simple expression, but I hope you understand what you must do here.


Complex Data type not accepted in Simulink I am having an issue with the interaction of the Matlab function block in the Simulink with my model block. Below you can also find the picture.
I have created a Matlab function with three inputs and 1 output. Since in the MatLab function block, I have solved a differential equation, it has complex expressions in the form of for example.

(exp(sqrt(-a/b)))

This generates the output of complex data types from the function block. When passing this to the integrator which is the input to my plant, I can not run the Simulink due to the error

    "Cannot pass the complex value to non-complex location"

I tried changing the data types of the integrator, then the same thing happens to other blocks in the Simulink.
I also tried using absolute values, which of course gives real numbers, but turns out to be a very large number and generates high magnitude. I really need help in solving this issue.
Can anyone help me in this/

Thanks in advance.

The data type of the signal must be fixed in Simulink. Since you use a square root of a real signal, the output will be real if the data is positive, imaginary if negative. Simulink cannot allow that.
The workaround is to use force your real signal to be a complex signal by setting its imaginary part to 0. Then feed the signal to the square root block. The output will always be complex.


Why do we use sinc interpolation rather than just multiply the frequency domain signal by rect? I know if we have a continuous-time bandlimited signal $x_{ct}(t)$ that is sampled frequently enough we get a signal that we can call $x_s(t)$.

If we take the Fourier transform of $x_s(t)$ to get $X_s(j\Omega)$, we end up with a Fourier transform of the sampled input that resembles $X_{ct}(j\Omega)$ at $\Omega = 0$ and repeating at $2\pi/T$ and $-2\pi/T$. If we want to get $X_c(j\Omega)$, we can multiply by a $\mbox{rect}()$ function.  It seems to me that we could then reconstruct the original signal by taking the inverse Fourier transform to get $x_{ct}(t)$

However, when reading on signal processing theory, we instead use a $\mbox{sinc}()$ function to perform reconstruction in the time domain since it is the "ideal interpolator". Why can't we instead use the method described in the frequency domain? It seems like it would be easier?

The "two" operations you talk about are one and the same:  low pass filtering using the $\mbox{rect}()$ function in the frequency domain and $\mbox{sinc}()$ interpolation in the time domain are the same.


What is the instantaneous frequency of this signal Suppose we have the signal

$x(t) = \cos(\pi t) + \cos(\pi t^2) +  \cos(\pi \cos(\pi t))$

How can I find the formula which gives the instantaneous frequency $f(t)$?

Concerning the comment from Juancho, we can assume that the time is normalized to e.g. 1s, meaning that $t=1$ corresponds to $1s$. Then, we dont have a problem with the units (this a very common operation and not specific to your problem).

Now, regarding the "instantaneous frequency" (also check Wikipedia and another question on math.se), it is not rigorously defined. Due to the uncertainty principle you normally cannot assign an accurate specific frequency to an infinitely short interval of time (which would the term "instantaneous" correspond to). So, in general, the question for an "instantaneous frequency" needs to be asked with care and you need to understand that its interpretation might not be as straight-forward as it seems.

However, in case of a simple sine wave, consider the simple expression

$$x(t)=\cos(2\pi \cdot f\cdot t),$$

where we know that its frequency is $f$. The argument to the $\cos$ is $\phi(t)=2\pi f t$. If we calculate $\frac{d\phi(t)}{dt}=2\pi f$, which tells us the phase of the cosine changes with rate $2\pi f$. This term can be considered as the instantaneous angular frequency (i.e. frequency is $f$, as we would expect) of the sine wave.

Let's extend this expression to the general one

$$x(t) = \cos(\phi(t)) $$

where $\phi(t)$ is the phase of the cosine. We can now define the rate of change of the phase as the instantanous angular frequency: $\omega(t)=\frac{d\phi(t)}{dt}$.

With this idea, you can calculate the instantaneous frequencies (there are multiple) of your signal

$$\phi_1(t)=\pi t, \quad\phi_2(t)=\pi t^2,\quad \phi_3(t)=\pi\cos(\pi t).$$

So, the instantaneous angular frequencies are given by

$$\omega_1(t)=\pi,\quad\omega_2(t)=2\pi t,\quad\omega_3(t)=-\pi^2\sin(\pi t).$$

I.e. your signal consists of three different instantaneuos frequencies.


Variance of white noise and the STFT I hope someone can help me understand the meaning of variance in this case. I am reading an old research paper by Ephraim and Malah (for my DSP project) and it says:

Knowing that $D_k$ is the FFT of white noise, what I am trying to find out is $\lambda_d(k)$. Can anyone explain to me the meaning of $E\lbrace|D_k|^2\rbrace$? I think it is the variance of a particular $D_k$ with respect to the mean, but since this is white noise, shouldn't it be zero?

Please help, I am stuck!!!!

Thanks.

$E\{|D_k|^2\}$ is the expected value of the squared magnitude of the noise in the $k^{th}$ frequency bin. Since they assume that the noise is zero mean, this equals the noise variance. Remember that the variance of a random variable $X$ is defined as $\sigma_X^2=E\{|X-\mu_X|^2\}$, where $\mu_X=E\{X\}$ is the mean of $X$. So if $\mu_X=0$ (as assumed here for the noise) then $\sigma^2_X=E\{|X|^2\}$. Note that $E\{|D_k|^2\}$ is also the noise power in the $k^{th}$ bin, so it can't be zero unless there is no noise at all in that bin.


Extracting DCT coefficients from H.264/AVC bitstream Is it possible to extract DCT coefficients from H.264/AVC bitstream without doing a full-decoding of the video?

There are a lot of papers on compressed domain video processing which are based on the assumption that DCT coefficients can be extracted through partial decoding and thus a lower computational complexity and faster than real time decoding can be achieved by using the DCT coefficients for various computer vision tasks (e.g. tracking, action recognition). However, nobody really describes the exact process of such partial decoding, and I seriously doubt it can be achieved with the recent video codecs such as H.264/AVC or H.265/HEVC. Does anybody have an insight on this one?

It is possible to extract the transform coefficients without fully decoding the entire picture, but you would still have to find your way through the bit stream, parsing whatever is necessary in order to arrive at every transform unit, which might also not be present for every coding unit (if some macroblock matches exactly one in a reference picture then you don't need to encode any residual, for example).
How you would use only the transform information to do your computer vision tasks, I have no idea, but it is entirely possible to only extract the transform coefficients. What i suggest you do is grab whatever standard you want (H.264 or H.265 or some other) and go through all of the tables with the syntax elements, and the parsing process for arriving at the transform coefficients would be to ignore every fixed width syntax element not useful for decoding the transform (anything related to motion vectors or reference pictures) and parse whatever is necesary to obtain all the transform units in the coded picture. I don't know if there is any tool out there that does this, perhaps you can go into ffmpeg or one of the libraries like x264 or x265 and modify them to extract what you want. It is a difficult task, but definitely possible.

Hope it helps


2nd Order Type 2 PLL for Tracking Frequency Error This is a "DSP Puzzle", Please preface your answer with spoiler notation by typing the following two characters first ">!" as the intention is for readers to think through what the solution might be first (the ones that want to participate, others can hover over immediately to see the solution).

A second order "type 2" Phase Lock Loop is required in order to track a phase ramp with time with zero error. A phase ramp is a frequency step. Similarly a third order "type 3" PLL would be required to track a frequency ramp (phase acceleration), which I have summarized in the graphic below. The "Type" refers to the number of poles in the open loop transfer function at the $s=0$ in continuous time or $z=1$ in discrete time loops, since these represent pure integrators. An integrator has infinite gain at DC and hence infinite suppression of DC error.

A second order type 2 PLL is often implemented with a proportional-integral loop filter, adding the second pole in addition to the oscillator which has an inherent pole since it integrates phase (controlling the output frequency of the oscillator directly is controlling the integral of the phase of the output).

The question is: There is a way to implement a similar tracking system to track a phase ramp with zero error using two first order PLL's. How do we do this?

Notation: Incoming carrier frequency : $f_c$ and Local Oscillators synthesized frequency : $f_o$.   A First Order PLL to track a constant phase offset $(\phi)$ in the received carrier sinusoidal is based on the following PLL daigram. In working of the single PLL Loop it is essential that $f_c = f_o$:  (image courtesy : Software Receiver Design by Johnson, Sethares & Klein) Here, $r_p(kT_s) = BPF_{2f_c} \{ r(kT_s)^2 \}$, meaning received signal is squared and then bandpass filtered with BPF's center frequency at twice carrier frequency $2f_c$. This preprocessing is done to emphasize the carrier, since this will boost carrier SNR. So, input to the PLL is basically, $$r_p(kT_s) = cos(4\pi f_c k T_s + 2\phi)$$Now if instead of constant phase offset of $\phi$ if we are having a constant frequency offset of $\Delta f = f_c - f_{o}$ and a constant phase offset of $\phi$ in carrier frequency then, we can extent this single PLL Loop to a combination of 2 PLL loops as follows: (image courtesy : Software Receiver Design by Johnson, Sethares & Klein)  Here, input to this frequency Tracking PLL is as follows: $$r_p(kT_s) = cos(4\pi (f_c - \Delta f)k T_s + 2\phi)$$  Now, if the first Loop processing is worked out, it can be shown that it tracks a ramp:$$\theta_1[k] \rightarrow 2\pi (f_c - f_{o})k + \beta$$, where $\beta$ is the y-intercept of the ramp being tracked by first PLL.This is added to $\theta_2[k]$, which makes it possible for the $2^{nd}$ PLL to track the frequency accurately :$$\theta_2[k] \rightarrow \phi - \beta$$ , which mean, $$\theta_1[k] + \theta_2[k] \rightarrow 2\pi (f_c - f_{o}) + \beta + \phi -\beta$$ $$\theta_1[k] + \theta_2[k] \rightarrow 2\pi (f_c - f_{o}) + \phi$$, And this sum is being used to drive the sinusoid of Local Oscillator, hence the local oscillator is now generating a sinusoid which is synchronized with $r_p(kT_s)$ in both frequency and phase. Basically, top PLL Loop is tracking a ramp and hence generating the correct frequency for the bottom loop, and bottom loop is further correcting the residual phase error.(image courtesy : Software Receiver Design by Johnson, Sethares & Klein)


Frequency array feed real values FFT EDIT

The first version of my question has not worked as I expected, so that I will try to be a little bit more specific.

The final goal I am trying to achieve is the generation of a ten minutes time series: to achieve this I have to perform an FFT operation, and it's the point I have been stumbling upon.

Generally the aimed time series will be assigned as the sum of two terms: a steady component $U(t)$ and a fluctuating component $u^{'}(t)$. That is

$$u(t) = U(t) + u^{'}(t);$$

So generally, my code follows this procedure:

1) Given data

$time = 600 [s];$

$Nfft = 4096;$

$L = 340.2 [m];$

$U = 10 [m/s]$

$df = 1/600 = 0.00167 Hz;$

$f_{n} = Nfft/(2*time) = 3.4133 Hz;$

This means that my frequency array should be laid out as follows:

$$ f = (-f_{n}+df):df:f_{n} $$

But, instead of using the whole $f$ array, I am only making use of the positive half:

$$ f_{+} = df:f_{n} = 0.00167:3.4133 Hz; $$

2) Spectrum Definition

I define a certain spectrum shape, applying the following relationship

$$ S_{u} = \frac{6L/U}{(1 + 6f_{+}L/U)^{5/3}}; $$

3) Random phase generation

I, then, have to generate a set of complex samples with a determined distribution: in my case, the random phase will approach a standard Gaussian distribution $(\mu = 0, \sigma = 1)$.

In MATLAB I call

nn = complex(normrnd(0,1,Nfft/2),normrnd(0,1,Nfft/2));

4) Apply random phase

To apply the random phase, I just do this

$$ H_{u} = S_{u}*nn; $$

At this point start my pains!

So far, I only generated $Nfft/2 = 2048$ complex samples accounting for the $f_{+}$ content. Therefore, the content accounting for the negative half of $f$ is still missing. To overcome this issue, I was thinking to merge the $real$ and $imaginary$ part of $H_{u}$, in order to get a signal $H_{uu}$ with $Nfft = 4096$ samples and with all real values.

But, by using this merging process, the $0-th$ frequency order would not be represented, since the $complex$ part of $H_{u}$ is defined for $f_{+}$.

Thus, how to account for the $0-th$ order by keeping a procedure as the one I have been proposing so far?

To get a real signal valued time domain signal: create the negative frequencies as the complex  conjugate of the positive frequencies

To get a complex time domain signal: create negative frequencies by copying the amplitudes of the positive frequencies and applying a random phase.

If you construct the signal in frames you will been apply windowing and overlap to avoid discontinuities at the frame boundaries.


Relation between SER and the channel capacity Is there a relation between Symbol Error Rate (SER) and the channel capacity?

What is {Modulation Order x (1-SER) / symbol duration [bps]} called?

Yes, there is a relationship between SER and channel capacity.  The channel capacity equation is-

$$
C = B\log_2(1 + \frac{S}{N})
$$
where $C$ is the channel capacity in bits/s, $B$ is the bandwidth in Hz, $S$ is the average received signal power, and $N$ is the average noise power across the bandwidth $B$.

Looking at the equation, it is clear that increasing the SNR will increase the channel capacity.  Increasing SNR does, of course, also decrease the symbol error rate (SER).  Thus, there is an inverse relationship between SER and channel capacity.

  What is {Modulation Order x (1-SER) / symbol duration [bps]} called?

As far as I am aware it does not have a name.  Having a bit stream where some of the bits are good and some of the bits are bad is generally not very useful.  Shannon's theorem assumes that error correction codes are used to reduce the number of errors to an arbitrarily low rate.  This reduces the number of generated bits/s, but increases the reliability of those bits.  The "channel capacity" of a channel is the number of bits it can produce with an arbitrarily low error rate.


How to insert \ell in text box on figure in matlab? I have a plot in matlab and want to insert a text box write inside ell symbol. Till now i couldn't do it. I wrote '\ell','interpreter','latex' but it doesn't accept. I appreciate If anyone could help. Thanks.

Here is an example how you may do it;

plot([0 2 1]);
text(2,1,'$\ell$','interpreter','latex','fontsize',30);


Which is correct representation for function of variable? I have a function as
$$E=\int_\Omega -\log\big( p_i(x)\big) dx$$
where $p_i(x)$ is density distribution which estimated by Parzen window method.
$$p_i(x)=\frac{1}{\Omega_i} \int_{\Omega_i}K_\sigma\big(I(x)-I(y)\big) dy$$
where $\Omega_i$ is the image domain, with
$$K_\sigma=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(\frac{-z^2}{2\sigma^2}\right)$$
and $\sigma$ is scalar parameter. $I$ is an image $I(x):\Omega \to \mathbb{R}$, $i$ is given.

So, which ones is correct for left side: $E(\sigma),E(x),E(x,\sigma),E(p_i)$ or other form?

Update: I upload the original equation for better visualization

Based on your new explanations, I think the answer is just
$$
E_{u_i}
$$
given that what you're really trying to find a shorthand for is the expression inside the summation sign in (5) of the screenshot. That is:

$$
\min_{U} \left\{ \lambda \sum_{i=1}^N \int_\Omega r_i(\chi, p_i(\chi)) u_i(\chi) d\chi \right\} = \min_{U} \left\{ \lambda \sum_{i=1}^N E_{u_i} \right\}
$$

so that
$$
E_{u_i} = \int_\Omega r_i(\chi, p_i(\chi)) u_i(\chi) d\chi
$$

Original answer below for reference

As per @Fat32's comment, any of:

$E(\sigma, z, p_i)$
$E_{\sigma,p_i}(z)$
$E(z; \sigma,p_i)$

is acceptable.

But you need to tighten up your notation. What is $x$ vs $z$? How does $\Omega$ relate to them? What is $y$?


How many bins to include when calculating SNR from FFT? I'm a scientist conducting an experiment that requires some signal processing. My expertise is not in signal processing, thus here I am. We've basically re-created an experiment conducted by other scientists, attempting to check their results. Here is a link to their paper: Ultrasensitive Inverse Weak-Value Tilt Meter
In short, a laser bounces off of some mirrors, one of which is oscillating at a controlled sinusoidal frequency, onto a quadrant detector, which outputs an electrical signal to an oscilloscope where we record it. So, you end up with a noisy record that has a tiny, known sine wave hiding in it.
My question has to do with calculating the SNR from the FFTs of our records. How many bins do I include for noise? Do I include all of the non-signal bins (besides DC and Nyquist)? Or is there some standard for this type of thing?
As a follow on question, when determining the noise floor from our spectral densities, is there a certain standard calculation, or is it more of an eyeball the plot type thing? Is the floor determined for bins near the signal, or should you look at the whole record?
Thanks in advance for any help you can provide. As a side note, I have done quite a bit of due diligence trying to figure this stuff out. But, anytime I find some source that seems authoritative that says one thing, I'll find another that says something different. I can't tell if I'm misunderstanding what they're saying, or if there are just lots of different definitions for this stuff.
Here are a couple of our PSD plots.
Full Spectrum:

Zoomed in on the reference signal:

The OP's question requires further details to provide a definitive answer, but the following will give the considerations involved. If the noise is white and stationary then the answer is clear in that we can simply use a power spectral density and work with the SNR as an SNR/Hz quantity. If the noise is not white (meaning the average power across all bins is not the same) then the concern will ultimately be with the measurement bandwidth of your final system that is providing the estimate of the quantity you are looking for. You would include all noise that is within this bandwidth as that is what will give you the statistical error in your result (which is why we are concerned with SNR in the first place). Usually the choice of bandwidth comes down to a time /frequency decision and the coherence time of your signal, and for how long it can be assumed to be stationary. For bandwidth limited signals that can be assumed to be stationary for a long time, then we can use a relatively longer averaging time (reducing the bandwidth down  to that of the signal, and therefore the amount of total noise in our result, as long as the noise is a spectral density and not a "spur" that exists at one frequency).
Ultimately to answer your question, determine accurately what the bandwidth is of your measurement. Use that bandwidth in the analysis where you predict the total noise and you will get a useful SNR. (This makes sense, consider the result you will get from your measurement system with a noise-free signal, and then consider the result you will get with noise alone-- the ratio will be your SNR that matters as it is based on the final result-- any noise that got filtered out by your measurement system does not apply).  This is analogous to a communication system where we will use the bandwidth of the signal, since the ideal receiver will have filtered out anything outside of that bandwidth but not reject any of the valuable signal. Your ideal measurement system should do the same, which means you will have interest in the spectral occupancy of your signal (and for how long it can be assumed to be stationary, which leads to concern with phase noise since that is typically not stationary... read on below).
You are likely doing an experiment where you are trying to estimate the presence of a specific signal, in the presence of noise, and ultimately will need to determine the best observation time in which to maximize your SNR (and therefore the accuracy of your estimate). Phase noise as your first measurement gave is a great indication of this as it is typically (if not always) NOT stationary over long durations, (and therefore for very long duration the typical power spectral density measurement is not valid). A much better measurement you can make that will allow you to optimize your estimation duration is the Allan Variance (also known as a two-sample Variance) since that will indicate the duration of time in which you can use an assumption of stationarity, and will then give you the averaging time that will maximize your SNR.
I go into more details of the Allan Variance and Allan Deviation for this particular application here: What determines the accuracy of the phase result in a DFT bin?


Stft - Alternative way of visualising In theory I have 3 axis: time, frequency and amplitude (which is scaled in by $20\log$). I have an audio signal. I noticed that I have very high amplitudes in the range of $1$-$1000$ $\mathrm{Hz}$. So I made 3 frequency bands: $y_{low}$, $y_{medium}$ and $y_{high}$, where $y$ refers to the signal and the subindex to the frequency band.

Now I defined a frequency vector for my bands:

LOW = 1:200
MEDIUM = 201:600
HIGH = 601:1000

Afterwards I wanted to know the values of my bin, depending on the frequency band. Basically y_LOW/4096   --> (fft points). Now I did an arithmetic mean with the values of my bins, which is saved in a single variable.
Before that happens, I am doing a STFT. So, a spectrum is shown every $50 \ \mathrm{ms}$. I want to plot my 3 values (since I have 3 frequency bands) every $50 \ \mathrm{ms}$. The values should be shown as a point or something similar.

If the code is required, I can post it.

Instead of calculating the STFT, you could just implement a 3 band filter bank as shown in the image below.

Your input signal x[n] could then be broken into 3 subbands: x0[n], x1[n], x2[n] which contain the energy associated with each frequency range. You could then feed each individual signal into an accumulator block in which the output would be the sum of the last x number of points. Or in you case, the x number of points that occurred in 50 ms. This would give you a single number in which you could map to some visualization output. You could implement normalization factors however you would like.

All you would have to do is define the filters H0, H1, and H2 as bandpass filters corresponding to your desired cutoff frequencies.


How to find phase shift and process I/Q channel I am currently working of 24GHz FMCW radar signal processing. The Module I'm using has 2 receiving antennas and each has 2 outputs (I and Q channel). So overall I've 4 channels. I'm using combination of harware and software for the processing in frequency domain. Range and velocity is fine, if I sample any one channel out of those 4 (either I or Q) I can find the R and V values.

Now I'm interested in finding the direction of a object, when it is moving (away from sensor or towards sensor) and the position of the object (right of the sensor or left or the senor, at the best the horizontal angle).

In the datasheet of the radar sensor I read 2 lines in the features, that gives me some idea in finding the direction and location. Below are the lines:

"two receiving antenna for phase comperision operation"
"I/Q channels for direction of motion discrimination"

Now I'm thinking how to find out the position or the horizontal angle of the target using phase comperision of the 2 antennas. Below is the diagram I made for analysis.

Surely the 2 antennas receives reflection from a single target with a phase shift. Now the question is how to I find the phase shift processing the I channel (let's say) from the 2 antennas. I'm not aware of any DSP technique, so want suggestions.

Also the direction of motion, what technique can be ised with the I and Q channel of a antenna to deduce it. I'm again on poor dsp knowledge, how do I process the I and Q channel.

I appreciate any suggestion of improve my knowledge and know about more dsp techniques in signal analysis.

The way to do it is to correlate the received signal with the signal you sent out- likely a chirp or something.  You correlate it with the signal from antenna 1, which will give you a time (location of the cross correlation peak) and phase (phase of the complex value of the peak).  You do the same with the signal from antenna 2, getting a time and phase.  If all is right the time should be the same (or, depending on your time resolution, very, very close).  The phase difference is simply the difference of the phases of the cross correlation peaks from the two antennas.


Finding the phase and frequency of a tone with Doppler. What is wrong in the code? I am trying to find the phase and frequency of a single tone having a small Doppler (10 Hz). I have the following code. The integration time for the operation (the length of the time vector) is defined in the code. It works if the integration time is a correct multiple otherwise not. Why is this failing?

Here is the code:

clearvars; close all; clc;

 Fc = 10.005e3;
 Fd = 10; % random doppler
 Fc2 = Fc+Fd;
 phi = deg2rad(37.8);
 disp(['True Freq: ',num2str(Fc2),' Phi: ',num2str(rad2deg(phi))]);
 Fs = 1e6;      % Max. Sampling Freq allowed
 Ts = 1/Fs;     % Sampling interval
 T_intg = 1.2347e-3; % Max. Integration time allowed

 t  = 0:Ts:T_intg;                        % Time Vector
 s  = exp(1j*(t*2*pi*Fc2+phi));           % Signal Vector

 Fs = 1/Ts;                               % Sampling Frequency
 Fn = Fs/2;                               % Nyquist Frequency
 L  = length(s);

 fts = fft(s, L)/L;                       % Normalised Fourier Transform
 Fv = linspace(0, 1, fix(L/2)+1)*Fn;      % Frequency Vector
 Iv = 1:length(Fv);                       % Index Vector

 amp_fts = abs(fts(Iv))*2;                % Spectrum Amplitude
 phs_fts = angle(fts(Iv));                % Spectrum Phase

 [val, index] = max(amp_fts);             % locating the peak in the spectrum
 rad2deg(angle(fts(index)))                % this gives 104 instead of 37 degrees

Your FFT method of finding the phase of a single tone is only valid if your 's' signal contains an exact integer number of cycles. I.E., no spectral leakage. And that is NOT the case for your 's' signal.


Transposed Direct Form II VS Direct Form II IIR filters? I read in "Discrete-Time Signal Processing" Reference, in transposed Direct Form II filters section, that it implements zeros first then poles(unlike Direct form II which implements poles first), Which is important in the presence of quantization of finite precision or noise without adding further explanation.

Why implementing Zeros first is better than implementing poles first?

Compared to the Direct Form I, Direct Form-II has pros and cons.

The advantage of DF-II is its more efficient usage of the delay lines. Although both use separate all-pole and all-zero sections, the DF-II can share the delay lines between them and reduce the number of required delays to implement the same transfer function.
The disadvantage of DF-II is that it precedes the poles and imposes a higher dynamic range on the intersecting delay line in some frequencies. Therefore, unlike the DF-I which first implements the zeros and is almost immune to overflow, DF-II is prone to such effects.

In the transposed DF-II the zeros are implemented first AND the delay lines are shared between the all-pole and the all-zero section. Hence, the transposed DF-II has the advantage of DF-II and also offers a better robustness similar to the behavior of DF-I.

Also, it is stated here that they are more robust when filters with a sharp-transition are implemented.


What distribution is the easiest to compress? I'm currently playing around with some compression algorithms and I'm asking myself if there is a type of data distribution / noise distribution that is easier to target with quantization (meaning less distortions at same rate). To my understanding i.d.d. Gaussian is the upper bound on "compression difficulty".

Are there distributions that especially efficient to compress, maybe via quantizers designed for those (I know about LLoyd-Max but I'm looking for something that works generally well with a known distribution and doesn't need to be optimized samplewise)? Does it make sense to transform the data into such a distribution prior to quantizing it?

It's very easy to make fundamental mistakes on the topic of information theory, unless you have some PhD or alike in the field. That being said as far as I remember the following could roughly be stated about your question.

Considering a discrete probability and source model of iid, the most efficient compression of a source happens when the source has minimum Entropy (information) in it; i.e, its average codeword length (in bits) per symbol is minimum. Such a variable length code would have the highest compression and the most efficient coding to represent the source in a lossless manner.

Such a case would happen when one of the symbols of the source alphabet has a probability close to 1, and all other probabilities are close to zero. Such a distribution is said to be highly skewed. But I'm not sure if a unique minimum can be shown to exist unlike the case of maximum Entropy (minimum compression) which happens when all the symbols have equal probability.

On the other hand, you are talking about lossy Quantizers and probability distributions. PDF optimized lossy quantizer design is about finding the best input decision and output reconstruction levels for a given soure PDF such that the the expected quantization error (MMSQE) per given source distribution a minimum. And that's the Max-Lloyd quantizer for the scalar case. Rate-Distortion theory shall be invoked for various considerations here.

Merging these two concepts you would say that, for a given source power (or dynamic range) the one which required minimum number of quantization levels while still providing least possible quantization error (distortion) would be given by a source whose PDF is dense at a certain region. In that case you would put a few fine quantizer levels at those dense signal positions and then use larger step size levels outside of that region, which would still provide minimum expected distortion.


YV12 to RGB - What is wrong with my algorithm? I am attempting to convert YV12 to RGB using the following fragment shader:

precision mediump float;
uniform sampler2D tex0, tex1, tex2;
varying vec2 v_texCoord;

void main(void) {
    float r, g, b, y, u, v;
    y = texture2D(tex0, v_texCoord).x;
    u = texture2D(tex1, v_texCoord).x;
    v = texture2D(tex2, v_texCoord).x;

    y = 1.1643 * (y - 0.0625);
    u = u - 0.5;
    v = v - 0.5;

    r = y + 1.5958 * v;
    g = y - 0.39173 * u - 0.81290 * v;
    b = y + 2.017 * u;

    gl_FragColor = vec4(r, g, b, 1.0);
}

Below is the original image followed by the converted image.  Based on the output, can you determine where I might look to determine what is causing the improper decoding?

The output image has lines shifted one against another, which is typically due to incorrect [byte] distance between the lines (stride). That is, the color conversion math is perhaps correct, however you are either taking source pixel from wrong position and/or putting output pixel into wrong place.

For performance reasons YV12 images might be using extended strides in hardware to facilitate performance. For example, 320x240 YV12 image might have a stride of 512 and this means that the data buffer is a 512x240 data space (512x240 bytes for Y, followed by two 256x120 arrays for V and U) where the payload image on the left, and there is invisible padding on the right.


Is it possible to observe more area in azimuth than physically traveled in stripmap SAR? I have been examining stripmap SAR algorithms, mainly Range Migration Algorithm. When I check the textbook examples so far I always saw the data trajectory is larger than scene that is observed. An example from Melvin, Scheer et. al page 309.

As you can see, the platform trajectory roughly travels -20 to 20 and scene is half of it. Similarly, if you check this RMA implementation, this RMA Matlab implementation from the book "Spotlight SAR" by Carrara et. al., you can see the along-track axis of the reconstructed image is actually some subset of trajectory points. Just like the case above.

My question is, what happens in the opposite scenario? For instance in image above, we travel a small distance [-5 to 5] and try to produce, maybe not so precise, the image of the scene [-10 10]. If our antenna does not have high directivity, we should also able to get reflections from targets that are not directly in front of the antenna, but how we should map the frequencies in RMA?

I know that in Stolt mapping of RMA, we shift and strech each azimuth column, but in that case I should retrieve more azimuth columns than I have.. I don't know, if there is a such imaging type in SAR literature, please refer me to correct terminology, I am new to that area so I might be missing correct nomenclature.

It's so that the cross-range (along-track) resolution is consistent for all targets in the scene. There is no reason why you could not use the flight path between [-5,5] to form the same size image. However, in that case, the targets on the edge of the scene will have been illuminated for a smaller angle than the targets in the middle of the scene, which will given them a smaller region of support in the frequency domain and, hence, a coarser resolution.

When choosing the first and last sensor locations that are used to form the image, it's common to select the first point to be that where the beam first enters the imaged scene and the last point to be when it leaves it. Picture the footprint of the beam on the ground as the sensor flies by; it's usually an oval or cone or something like that. At -20 is where the edge of the oval (or whatever shape it is) begins to intersect the imaged scene. At +20 is where it stops intersecting it.

I should add that, in general, the size of the imaged scene is limited only by the PRF of the sensor.


How many DFT bins must I use to get best precision when analyzing melody scale I need to analyze some signals containing chords.

How big DFT spectrum is needed to identify each possible musical note in such a signal sampled at 44100 Hz? I can see that 8192 bins gives some frequencies close to some notes like C3 A4 etc., but how big should it be to contain the whole scale, let's say to C5 and with good precision, like for a very good tuner?

You can always improve the accuracy of frequency estimates by quadratically interpolating the peaks in the DFT. This page should explain how to do this: https://ccrma.stanford.edu/~jos/parshl/Peak_Detection_Steps_3.html

You still end up with the issue of converting peaks in the spectrum to pitch, i.e. which peaks correspond to the fundamentals. The fundamental may not even be present as well!


Demodulation of OFDM symbol in LTE for multi-user case I was trying to understand what exactly happens when you try to modulate an ofdm symbol in the downlink in LTE for multi-user case.

When we modulate, usually we allocate different sets of RBs for different users. Since an RB is 12 subcarriers in frequency, a given ofdm symbol (say 2048 SCs in a 20 MHz case) can contain data from different users. And we do take IFFT of this 2048 sized data vector comprising of different user data. Is that right?. OR do we append 2048-12 zeros to a specific user data and then take IFFT?

Thanks in advance.

The zero padding (or the cyclic prefix) takes care of inter-symbol interference that occurs due to delay spread. It is hence something you append to your signal in time domain. Hence the order is this:

Assemble all the data symbols along the subcarriers (can contain data for multiple users, as you said).
Compute DFT (via FFT, typically)
Add cyclic prefix (or zero padding)

Then transmit.


Why Is the Total Time Equal to $ N \cdot {T}_{s} $ and Not $ \left( N - 1 \right) \cdot {T}_{s} $ In the Context of DFT? In the definitions of the DFT
DFT
$$
X(j)=\sum_{k=0}^{N-1} x(k) \exp \left(-i 2 \pi\left(\frac{j}{N}\right) k\right)
$$
Let us say, if we have $10$ points, $N=10$, each sampled at $0.2$ seconds, why is the total time used for calculating the frequency resolution is equal to
$$
\frac{1}{N\Delta t }
$$
where $k$ will run from $0, 1, 2, \ldots , 9$.
If the first point was at zero, the sampled time will be at $0.2$, and the last sampled point will be at
$$0.2\cdot (N-1)= 0.2\cdot 9 = \mathbf{1.8 \ \rm s}$$
Rather the total time is equal to $0.2\cdot N= 0.2\cdot 10=\mathbf{2.0 \ \rm s}$ in the frequency step.
P.S.: I have seen the query and the discussion How do I measure the time duration?

Here $\Delta t = 0.1 \ \mathrm{s}, N= 11 (\text{Eleven data points}),
k= N-1$; So
$$\text{total signal duration} = k\cdot \Delta t = (11-1)\cdot \Delta t= 1 \ \text{second}$$
This agrees with the $10 \rm \ Hz$ sampling rate, i.e., $10$ points were collected in $1 \ \text{second}$ and the $11^{th}$ point belonged to the next cycle.

You are right, the duration associated with taking $N$ uniform samples of a signal is
$$ D = (N-1) \cdot T_s$$
where $T_s$ is the sampling period.
A concrete example is sufficient; assume your sampling period $T_s$ is 1 hour long, and you want to take 3 samples of a slowly changing process, such as the height of an iceberg tip while it's melting.
Your first sample is taken at $t=0$ (the electronic sampling process itself takes about a micro second or less, so ignore it compared to an hour of sampling-period!). Then your second sample arrives at $1$ hour later and your third (and final) sample arrives at $2$ hours later.
Hence your $3$ samples long observation takes $D = (3-1) \cdot 1 = 2$ hours long. As soon as you take your last (third) sample, you shut down the sampling system. You do not wait one more hour (one more sampling interval) after taking your last sample.
And this methodology of calculation is exactly on par with calculating distances within crystal lattice structures. What's the distance between N atoms? What's the total length of N atoms (regularly placed on x-dimension) ?
Nevertheless, in the literature, you can find expressions involving $D = N \cdot T_s$. Some applications may require that; i.e., block based signal processing, DFT, sample rate conversion make use of such a viewpoint, which is justified in their processing of data blocks one after other.
To understand why $D = N \cdot T_s$ may be used in DFT analysis, consider the following example. Assume you have long data set, such as $4 \cdot N$ samples, divided into 4 blocks of $N$ samples; i.e, you will have 4 blocks of $N$ samples each. The blocks are adjacent, their sample orders are (1,N)  ,(N+1,2N) , (2N+1,3N) , (3N+1,4N). The sample $N+1$ belongs to second block, but the duration of first block is measured beginning from the sample 1 up to sample N+1. Because, the duration between the samples N and N+1 belongs to the first block, and this explains why the duration of that block is taken as $D = N \cdot T_s$. However, for the last block of samples (3N+1,4N), the duration will be $(N-1)\cdot Ts$, as there are no more adjacent blocks.
Last but not least, this is a topic of debate. :-)


Time-frequency representation of sound signal using Matlab I am currently doing a project on sound signal processing by using MATLAB. I have a problem regarding creating the command line so that I can get a graph containing frequency by its time. But this graph is the result of an audio input.
example:

x = audioread('samson.wav')

But I still can't get this kind of a result. How can I get this colour graph? Thank you

I believe that this "color graph" you are looking for is a spectrogram (although it looks to me more like a scalogram, but you did not mentioned wavelets). Let me give you an example in MATLAB of obtaining such plot:

load handel
nfft = 512;
noverlap = 128;
win = hamming(nfft);
spectrogram(y, win, noverlap, nfft, Fs, 'yaxis')
colormap('jet')

So first line is loading of some standard test recording, it gives you two variables: Fs (sampling frequency) and y (your signal). Next you define length of your analysis: nfft in your Short Time Fourier Transform and by how many samples you want to shift your analysis window. When defining window you can use variety of them (simply check in MATLAB help), I've used most common one: Hamming window. Last line is for setting colour representation.

Computation of spectrogram is very straightforward, please refer to help for more details: click. The last argument 'yaxis' tells MATLAB to use horizontal time axis and vertical for frequency.

Having that you can play with length of your analysis window and the overlap. You must though understand one, most important thing: better resolution in time (short window) yields poorer resolution in frequency domain and vice versa - that's a trade-off.

In the end you should get something like:

EDIT:

Because you would like to obtain the scalogram, here's how to do it in MATLAB. Most important thing - you must have Wavelet Toolbox installed!

% Definition of signal consisting of two sinusoids
f1 = 10;
f2 = 40;
Fs = 1000;
t  = 0:1/Fs:1;
y =  sin(2*pi*t*f1) + sin(2*pi*t*f2);

wname = 'morl'; % Choosing the wavelet
scales = 1:1:128; % Defining scales
coefs = cwt(y, scales, wname); % Get wavelet coefficients

wscalogram('image', coefs, 'scales', scales, 'ydata', y); % Get the scalogram, together with time domain signal overlayed
colormap(jet) % Set the colormap

As a result you will get something like:


Basic question on numerical precision I appreciate this opportunity to submit a query on this forum.

When studying the continuous-time & discrete-time distinction, specifically with reference to discrete signals being identical when separated by 2*pi, it has struck me that a basic premise, per google calc, doesn't hold:

exp(2 * pi * sqrt(-1)) = 1
exp(4 * pi * sqrt(-1)) = 1
exp(6 * pi * sqrt(-1)) = 1
exp(8 * pi * sqrt(-1)) = 1

but...

exp(10 * pi * sqrt(-1)) = 1 - 1.2246468 × 10^-15 i

Any thoughts on the above (discrepancy?) would be appreciated.

Regards,

wirefree

Writing another way

$$1 - 1.2246468 × 10^{-15} i = 1 - 0.0000000000000012246468i$$

It is just rounding errors adding up due to having a limited number of bits in whatever computer is doing the calculation.


A General Question regarding OFDM symbol I have cross across an OFDM modulator that says that

A frequency domain sequence is passed into a 64 point IFFT. A cyclis prefix is added (0.8 $\mu$sec) and creates as a result a 4 $\mu$sec time domain sequence (OFDM symbol). I am wondering The input of IFFT is discrete the output is also discrete of length 64. Then what do the authors mean by 3.2 $\mu$sec OFDM symbol? Is it a continuous time domain signal?

Your question is a little confusing, but I'll give it a shot.

The 64 sample sequence output from the IFFT is a discreet sampled time domain sequence of complex values.  (A signal at complex base band)  The 3.2usec length comes from the sample rate (Fs) of the system.

In this case it sounds like Fs is 20MHz.  That Fs comes from the fact that there are 64 sub-carriers (FFT = 64)  therefore the frequency spacing to maintain orthogonality between sub-carriers is 1/Tu = 1/3.2usec = 312.5KHz.  Tu is the "used time", which in this case is that 3.2usec.  The total bandwidth is then 64*(sub carrier frequency spacing) = 64*312.5kHz = 20MHz.  Since the signal is complex Fs = BW.

Fs = BW = NFFT / Tu = 64 / 3.2usec = 20MHz

In my experience an OFDM symbol (or FFT symbol) refers to the data without the cyclic prefix added with its length equal to the FFT size.  I call the OFDM symbol plus the cyclic prefix and OFDM frame. Other people may have a different definitions.  It looks like this author uses the same definition as me not including the 0.8usec cyclic prefix in what they call an OFDM symbol.  That 0.8usec is also called the guard time and acts like a spacing in between successive OFDM symbols offering multipath protection.

Not sure if that will help or not.


Calculate amplitude in dB and dBA from audio signal Progress
I'm a product developer so kind of light on DSP inner workings. I have captured an audio signal - 1024 samples at a time. What I want to know is the amplitude of the wave in dB, and in dBA.
dB
For dB, I've created a 10 deep circular buffer recording data every 100 ms, and averaging the maximum reading over the last 10 reading. Therefore my output looks less jumpy and is more human readable. To get the amplitude I just record the max and min of the signal and take the difference.

Q1: Is this the correct way (or a valid way) to get a dB amplitude from a sample? Analternative is to work out the running average and then do average RMS to get energy in the signal. But I don't think that really corresponds to the maximum audio sound level?

dBA
For dBA, I have windowed 512 captured samples with a flat-top window and done an FFT on the data. It all looks fine so far. ( I am resource constrained and can only fit 512 samples into my FFT).
Sampling at 44.1 kHz for standard audio, with 512 samples, I end up with 256 usable frequency bins of 86 Hz, 172 Hz, etc... up to 22 kHz. So far any corrections, please jump in.

Q2: How do I interpret this bin? As representing the band from 0 Hz to 80 Hz? Or as 80 Hz? Since the next step is to apply A-Weighting to the frequency measurements, I need to know what frequency to use.

A-Weighting: I've generated an A weighting table in code based on the values of the bins. In this case, 256 values for 86, 172 ... up to 22 kHz. Clearly the point is to add the weighing on to compensate for the human frequency response - which translates PSL to phons.
What is confusing me is that I have seen some formulae that suggest that the dBA measurement is done by taking only single measurements for each octave starting at a suitable low frequency i.e. 33.25 Hz. Then you create a table with 10 entries in it taking you up to 16 kHz and apply the A-Weighting to those. You basically do an antilog, average, log of the 10 values and bob's your uncle.

Q3: I have 256 bins, and I'd like to make best use of them for improved accuracy. The approach above would suggest that each entire octave contains an equivalent amount of energy for the averaging process. On that basis, if I were to average all 256 bins, the result would be heavily (and incorrectly) weighted towards the highest octave in my sampled set, since it contains half of the bins! Which approach is correct? Is there some magical reason that energy isn't distributed linearly across the frequency domain, albeit that we like to look at frequency in logarithmic scales.

Accepted Answer And Additional Detail
I have accepted @ZR Han 's answer as it succinctly provided the required information to solve the stated problem, deftly sidestepping the misdirection in the errant assumptions I originally made in the question.
Using the formula provided, the time domain approach works perfectly, yielding the answer in dB. I wanted to add that the 2nd order IIR cascaded 3 times with the coefficients taken from  jenshee.dk/signalprocessing/aweighting.pdf also works perfectly. I have tested this using a signal generator at various frequencies (old school I know) and the ratio is exactly 1 at 1 kHz, as expected for A Weighting. Problem solved.

Q1: Is this the correct way (or a valid way) to get a dB amplitude from a sample?

Assuming you want to build a sound pressure level meter the answer would be a resounding "no". The correct way is to build a running energy detector with a proper time management as defined in the IEC standard IEC61672. https://webstore.iec.ch/publication/5708. This is actually pretty simple. By far the most difficult part is the calibration of the microphone. As ZR Han in their answer already pointed out: you get to a sound pressure level in dB (dBSPL) by using a reference pressure of $20 \mu Pa$, but that means you have to convert your input signal from Volts or "number" to Pascals as well and this is very device and settings dependent.

Q2: How do I interpret this bin?

That's a complicated question. Each digital bin contains energy from ALL analog frequencies EXCEPT the center frequencies of the other digital bins and the exact relationship isn't trivial. In your case you should just ignore this since there is no reason to do go into the frequency domain in the first place.

Q3: I have 256 bins, ...

Again, that's a very complicated question and I won't get into the answer since its completely unnecessary for what you want to do.
Building a sound level meter is relatively straight forward

Calibrate your input signal to represent actual sound pressure in Pa (Pascal)
Apply time domain weighting filter (A,C, Z, none)
Square the signal
Apply exponential (first order lowpass) averaging with proper time constant
Convert to dB using $p_{ref} = 20 \mu Pa$

Step 1 is by far the most complicated.


Does flat surface detection and orientation require special hardware to compute? Modern augmented reality platforms such as Google's ARCore and Apple's ARKit seem to only operate on mobile devices, I'm guessing, because their underlying algorithms require specialized hardware that is typically available on these devices (accelerometers, gyroscopes, etc.). Since this hardware isn't available on laptops/desktops, I'm guessing these libraries would never be able to work off of mobile platforms.

Having said that I'm wondering if the following AR/image processing "capabilities" can be achieved via algorithms that do not require such specialized hardware:

Flat surface detection inside an image or video; and
Orientation/angle of the camera view inside an image or video

Meaning, if I am given an image or video (set of sequential images), can I detect flat surfaces in that media without specialized hardware (as in, on a Linux box running on server or PC)? Same goes for orientation/angle determination (meaning, determining the (x,y,z) coordinates of the camera within the "domain" of the image/video and determining which angles its pointed towards)?

Modern augmented reality platforms such as Google's ARCore and Apple's ARKit seem to only operate on mobile devices,...

They only make sense on such devices.

  ... I'm guessing these libraries would never be able to work off of mobile platforms.

They do work outside of a mobile device... (and here).

  Flat surface detection inside an image or video

For anything to be detected by a computer vision algorithm it has to be in contrast to its environment. If not, then extraneous information would have to be used that does make contrast, e.g. obtain images in different wavelengths.

Provided that surfaces are in contrast to their surroundings it is possible to reconstruct 3D geometry (e.g. multiple surfaces coming together forming one object), by obtaining more than one images around the object of interest.

If a surface is flat, its orientation can also be be detected easily with a marker anyway.

  Orientation/angle of the camera view inside an image or video

This is called "Simultaneous Localisation and Mapping (SLAM)", it can be done with or without a marker (usually, when the camera is used as a sensor in robotics, it is done without a marker) and there are lots of libraries already on places like github that you can look for inspiration or further details of their inner workings.

Hope this helps.


What does an ADC DNL specification mean, exactly? I'm looking at the differential non-linearity specification for an Analog-Digital converter. The spec sheet claims that the DNL ranges from -1LSB to +1LSB, with a typical value of $\pm$0.6LSB. What exactly does that mean? Is this some value that's taken statistically? Does this mean the 3$\sigma$ value is 0.6LSB? I.E., for any given code, the range of voltages over which it has jurisdiction is 99.7% likely to be no greater than 1.6LSB or less than 0.4LSB?

Suppose I wanted to model a 16-bit ADC with these specifications. How would I go about injecting DNL into the output codes? The spec sheet is completely unclear as to what 0.6LSB means. Is there some statistical convention? Or would it just be some random number which caps at $\pm$1LSB for every single code?

Differential non-linearity specifically is the deviation from ideal between successive bits. For an ADC, the ideal would be the converted signal will not transition to another bit value until the input deviates by 1 bit. If the DNL is greater than 1 bit, non-monotonic behavior could result (an linearly increasing input could result in some transitions actually going down).

If the spec claims ±1 bit, this is a absolute maximum limit and is therefore guaranteeing monotonic behaviour. The "typical" limits are not so stringent in what that means and you would need to ask the manufacturer directly, but in essence it is what you could typically expect, on average (so therefore could be the the mean of the absolute value of the deviation).

Similar to if I said I had a 9V battery that had a spec of +/-0.5V and provided 9V typical. The average is 9V, the peak limit that the manufacturer will guarantee is 8.5V to 9.5V; now as far as the statistics that is completely up to them in balancing yield at final test in what they will ship.

If I was to model this, I was first tempted to generate a Gaussian random sequence with a standard deviation of $\sigma = (0.6) 4/\pi$ and clamp it at $\pm 1$ since the mean of the absolute value of a Gaussian R.V. Is $(4/\pi)\sigma$... but the clamping reduces the mean, so increasing that to $\sigma = (0.7) 4/\pi$ results in a random (clamped Gaussian) distribution that never exceeds ±1 but also has an average magnitude (therefore typical) of 0.6.


Why the angle of a real signal in frequency is continuous? This original signal is shown below (time domain):

It is a real signal. According to the signal processing theory, the form part of the angle is conjugate to that of the latter part. That is, $\sphericalangle X(e^{j\omega}) =\sphericalangle X^*(e^{-j\omega})$ when $x[n] \in \mathbb{R}$.

However, I implement X = angle(fft(x)); and plot(X) in MATLAB. The figure is shown below:

The phase is continuous! Could anyone explain the reason? Thanks in advance!

X = angle(fft(x))

is numerically very questionable. If your signal is periodic in exactly the length of vector x you're observing, then the (theoretical) DFT would be zero in all but two bins. You can't really take the angle of $0$.

If your signal has a period that doesn't exactly fit in x, then you're observing the phase of two superimposed sinc pulses, as that's what a rectangular windowing leads to.

A sinc is for all points (but the maximum) defined as $\frac{\sin ax}{bx}$ (with $a$ and $b$ depending on who you ask and how you scale), and if you look at that, it definitely does have at least a strongly linear component.


MATLAB phase of 2D rectangular pulse's Fourier transform I'm trying to plot the graph of the phase of the Fourier transform of a 2D rectangular pulse. I've been able to evaluate the FFT but I'm not sure if the phase is correct because there are some tilts I can't explain.

What are those diagonal lines that appear in the phase?
How can I eliminate them?

Here there is the code:

clc
close all
clear
npoints=512;
perc=1;
dt=6*1E-7/(npoints*perc);  % Tempo di campionamento
df=1/(npoints*dt); % Frequenza di campionamento

t(1)=0;
f(1)=0;
for k=2:npoints/2
    t(k)=(k-1)*dt;
    t(npoints-k+2)=-t(k);
    f(k)=(k-1)*df;
    f(npoints-k+2)=-f(k);
end
t(npoints/2+1)=t(npoints/2+2)-dt;
f(npoints/2+1)=f(npoints/2+2)-df;

ts=ifftshift(t);
fs=fftshift(f);

figure
[X,Y]=meshgrid(ts,ts);

D = npoints/2;        % to indicate origin at the center of the function
a = 5;          % change it to enlarge or reduce the pulse
y = repmat(1:npoints,npoints,1);
x = y';
rect = zeros(npoints);
rect(D-a:D+a-1,D-a:D+a-1) = ones(2*a);
rect=(rect);
surf(X,Y,rect);
shading interp
axis tight
title ('Rect 3D');
rect=ifftshift(rect);
figure, surf(X,Y,rect);
shading interp
axis tight
title ('Rect 3D shifted');

R = fft2((rect));
R = fftshift(R);
[X,Y]=meshgrid(fs,fs);
figure;
surf(X,Y,abs(R));
shading interp
axis tight
title('Fourier Transform of Rectangular function');
%plot real part
figure;
surf(X,Y,real(R));
shading interp
axis tight
title('Real part');

Rm=abs(R);
imm=imag(R);
re=real(R);
re(abs(re) < 1e-8) = 0;
imm(abs(imm) < 1e-8) = 0;
R=re+imm;
phase=angle(R);

%plot phase
[X,Y]=meshgrid(fs,fs);
figure
surf(X,Y,phase);
shading flat
axis tight
title ('Phase of the rect');

The issue is that the phase value found is between $-\pi$ and $\pi$ (or $0$ and $2\pi$) but it needs to be "unwrapped" to be continuous.

In 1D, the unwrap function will help. Your mileage may vary when applying it to a 2D Fourier transform.

Also, there is a bug in your code that forces the R variable to be real valued:

R=re+imm;
phase=angle(R);

is not going to get you the right phase at all.


Downsampling by block averaging Many spectroscopic detectors used to study light absorbance as a function of time. These detectors cab sample data at 160 to 250 Hz. There is an option of choosing lower sampling rates such 80, 40, 20, 10 Hz to decrease noise. The manual uses an interesting term called "data bunching" when lower sampling rates are chosen. I asked the engineers and he clarified that it is akin to block averaging. The data is still being sampled at 160 Hz, however digitally it is being downsampled to 80, 40 Hz etc.

Just for the sake of simulating data bunching I would like to downsample data simulated at 200 Hz to 50 Hz by averaging blocks of 4 and plot it versus time. For example, if y=[1,2,3,4,5,6,7,8,9,10,11,......], I would like to do averages of (1,2,3,4), then (5,6,7,8), and so and plot it versus t=[0:1/50:60], so that the downsampled data appears to be sampled at 50 Hz."

Is there a simple syntax for Matlab for that?

Thank you.

Averaging by N but keeping the same sample rate is equivalent to filtering input "y" with an FIR filter "B". Where B = 1/N * ones(1,N);

yFiltered = filter(B,1,y);

Now you simply need to downsample the data by N. I'm not sure exactly how the engineer implemented the original algorith, but I would do this.

yDownsampled = yFiltered(N:N:end);

Basically start at index "N" and keep 1 out of N samples until the end.

Your starting index can be anything between 1 and N, you can toy around with it.

You'd need to do the same with your time vector.


Using MFCC in spoken words recognition We're trying to implement a "simple" speech recognition application in MATLAB (isolated words from a very limited dictionary). We've been trying the following methods:

Extract MFCC coefficients for each frame of the word, and then compare to templates that we have recorded previously using DTW, and take the template that gave the minimal distance as the recognized word. (note: the "templates" are 10 recordings of each of the 6 dictionary-words).
Extract MFCC coefficients for each frame, and run an SVM on the coefficients (for each dictionary word we had a different SVM classifier).
Extract MFCC coefficients for each frame, and then define the feature vector as the vector of all distances from all the templates, then run an SVM on these features.

All three of these methods didn't work well (and gave almost everything wrong results).
We don't know what the problem is.. A few questions that arose, but we don't know the answers, are:

Do we really need 10 samples of each word? or can we "combine" them into a single template? and if so, how?
Should we run the SVM directly on the MFCCs vectors? On the MFCC vector of vectors (for each word)? On the DTW values? or should we combine all the MFCC vectors into a single one?
If we should look for a minimal match for the DTW between the input and the templates, should we give different weights to the templates depending on length, or even account for the differences between the templates (the templates of the word "1" are more distant from one another than the templates of the word "2").

We would appreciate any help, or references to good sources. We thought about using HMMs, but it seems much more difficult and not really possible in our time-frame..

Thanks.

Do we really need 10 samples of each word? or can we "combine" them into a single template? and if so, how?

You can use single template, 10 templates will make result more accurate.

  Should we run the SVM directly on the MFCCs vectors? On the MFCC vector of vectors (for each word)? On the DTW values? or should we combine all the MFCC vectors into a single one?

SVM has no meaning for your application because it is not invariant to time stretch so it's not suitable for classification of time series. There are SVM modifications to deal with time series, for example

Hidden Markov Support Vector Machines

but it's completely different story

  If we should look for a minimal match for the DTW between the input and the templates, should we give different weights to the templates depending on length, or even account for the differences between the templates (the templates of the word "1" are more distant from one another than the templates of the word "2").

No, it's not necessary. A simple comparison of DTW values should work.

  All three of these methods didn't work well (and gave almost everything wrong results).

Probably you made a mistake in DTW implementation first of all. You need to verify if you implemented everything correctly. You might make a mistake in DTW and go to incorrect conclusion

  or references to good sources.

Read a good text book. At least Rabiner "Fundamentals of Speech Recognition"


Relation between power of the baseband signal and the low pass equivalent Let $x(t)$ be a real bandpass signal and $x_l(t) = (x(t) + j\hat{x}(t))e^{-j2\pi f_0t}$ be the lowpass equivalent. Is there any relation between $P_x$ and $P_{x_l}$ where $P_x$ refers to the average power of $x(t)$ and $P_{x_l}$ refers to the average power of $x_l(t)$? I tried to use the definition $$P_x=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}|x_l(t)|^2dt = \lim_{T\rightarrow\infty}\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}(x^2(t) + \hat{x}^2(t))dt \tag{1}$$
And I got stuck here. Maybe we should add some assumptions on $x(t)$ in order to be able to simplify $(1)$.

Since the Hilbert transform is an allpass transformation, the powers of $x(t)$ and its Hilbert transform $\hat{x}(t)$ are the same. Consequently, the complex baseband signal as defined in your question has twice the power of the real-valued bandpass signal. That's why in some textbooks the bandpass signal is defined as
$$x(t)=\sqrt{2}\textrm{Re}\big\{s(t)e^{j\omega_ct}\big\}\tag{1}$$
where $s(t)$ is the complex baseband signal.
The scaling in $(1)$ ensures that $x(t)$ and $s(t)$ have the same power.


Upsampling question- Aliasing/imaging problem I'm trying to study vibrations in a non-stationary shaft by upsampling the original signal.

The original sample rate is at 25.6kHz, and it is being upsampled to 36kHz. In the FFT of the upsampled signal I am getting freq. components above the Nyquist frequency of 12.8 kHz. I'm using linear interpolation when upsampling.

Is this imaging/aliasing? What are common techniques for dealing with such behaviour? I'm working in python/matlab.

Would scipy.signal.resample_poly.html be sufficient?

Yes I would suggest to use resample_poly in scipy. When doing upsampling, you would get artefacts outside 12.8kHz, which you would remove via Low Pass Filtering. This is what is done by scipy.signal.resample_poly. You can enter the upsampling factor value as 36k/25.8k = 1.39534, and downsampling factor = 1.

In the above method while doing low pass filtering, the effect in time domain is to do sinc-interpolation. But in Linear Interpolation, you are not doing interpolation using sinc functions. You are just using neighboring points to compute the interpolated value. This is like moving average filtering, whose effect in frequency domain is sinc. It will not totally filter out all the artefacts outside 12.8kHz. But you still can improve it by taking more number of points for linear interpolation.

Here is the figure showing remains of original signals copy even after doing linear interpolation beyond 12.8kHz (not an exact mathematically correct shape).


Fundemental scintific publication in up conversion down conversion Prior to performing DSP on a signal and afterwards, up and down conversion of the same in order to transmit and receive the signal is important.

I would like to know if anyone knows a fundamental scientific publication of good clarity (preferably from a reputed publisher e.g. IEEE, ACM, Science Direct, SIGCOMM, etc) which introduces this concept. I understand that this publication will be an old one.

Certain Topics in Telegraph Transmission Theory by Nyquist (1928) might be what you're looking for (emphasis mine):

  The following topics are discussed:

  The required frequency band is directly proportional to the signaling speed.
  A repeated telegraph signal (of any length) may be considered as being made up of sinusoidal components. When the amplitude and phase,
  or real and imaginary parts, of these components are plotted as
  ordinates with their frequencies as abscissas, and when the frequency
  axis is divided into parts each being a frequency band of width
  numerically equal to the speed of signaling, it is found that the
  information conveyed in any band is substantially identical with that
  conveyed in any other; and the bands may be said to be mutually
  redundant.
  The minimum band width required for unambiguous interpretation is substantially equal, numerically, to the speed of signaling and is
  substantially independent of the number of current values employed.
  A criterion of perfect transmission is selected; and a discussion is given of the characteristics which the received wave must have to
  be nondistorting with the requirement that the frequency range should
  not be greater than necessary.
  Directions are sketched for specifying systems to meet this requirement.
  Several alternative criteria of distortionless transmission are considered and a method for computing the corresponding transmission
  characteristics of the circuit is explained and illustrated.
  An analysis is given of the carrier wave, and it is shown that the usual carrier telegraph requires twice as much frequency range as the
  corresponding d-c. telegraph, other things being equal.
  A discussion is given of two alternative methods for overcoming this inefficiency of carrier telegraphy, namely, the use of phase
  discrimination and of a single sideband.
  After the d-c. and carrier waves have thus been analyzed a corresponding analysis is given of an arbitrary wave shape, including
  these two as special cases. Calculations are given on the shaping of
  the transmitted wave so as to make the received wave perfect.
  A discussion is given of the dual aspect of the telegraph wave. The wave may be looked on either as a function of , requiring the
  so-called steady-state method of treatment, or as a function of
  requiring the so-called method of transients. It is shown that the
  steady-state theory can be made to yield the information necessary to
  specify the characteristics of an ideal system.
  A discussion is given of the effect of interference and departures from ideal conditions.


Meaning of lines on x- and y-axis of a centered FFT-spectrum of an image I am currently learning about the Fourier transform of digital images, but I am wondering about one aspect.
In most visualisations of fourier transforms of images, a horizontal and a vertical line in the middle seem to be rather dense compared to the rest of the image. Below is an image of two pattern which, from my understanding, should only have three entries in the fourier spectrum each.

That there are more than three lines should be due to the fact that there are horizontal and vertical sinusoidal waves of most frequencies used to represent the original image.
These two lines in the FFT-spectrum seem to appear in most images, not only in such with regular patterns, so my guess is it has something to do with the discrete pixels, but i can't quite wrap my head around why it would use all frequencies and only in those two directions.
Even if it is about "sampling" the pixels, wouldn't diagonal lines be as well suited to represent quadratic pixels (for less regular images than above)?

For an image containing an oriented 1D sine pattern looking like a corrugated sheet like the ones displayed, the 2D Fourier transform magnitude spectrum generally contains information along a line (aligned dots). One indeed may expect (at least) three dots: one in the center for the average or DC component (which could be zero if the mean is subtracted from the image), and two symmetric ones whose distance is in inverse proportionality with the frequency of the sine.
Since the image is limited in size (hence windowed) and the spectrum is computed using an FFT (which hypothetizes periodicity of the input), one may observe "ripples" (corresponding with a cardinal sine kernel in the frequency domain, dual to the 2D rectangular window) and border effects (with a centered cross).
You can find additional descriptions in Fourier Transform or
slides 8 to 14 of 2D Fourier transforms.


How to unify the color and brightness measure of two given cameras? I have two cameras, and I'm up to perform some matching stuff on videos captured by these two cameras since I'm a CV guy.
The problem is that two cameras have very different performance in terms of color and brightness. The same object in the same illumination appears different in two cameras.
Is there any academic work on measuring something like the "rgb response curve" of a camera? I want to unify the rgb response of these two cameras to guarantee the same object "looks" the same in two cameras.
Note Below is just an intuitive illustration of my problem. Same object appears different in different cameras.

Here you have two options.

1) Color calibrate your cameras (radiometric calibration) to the same reference. Then try matching. There are many academic works on this topic, one of them being:
http://research.microsoft.com/en-us/um/people/yasumat/papers/rankcalib_PAMI12_preprint.pdf

2) Use a vignetting invariant (or color invariant) feature descriptor. For example you could try to work with gradients instead of intensities and so on. This is harder, but if you could achieve it you get rid of the cumbersome calibration procedure. There are also valuable works done in this field, such as:

http://www.cvip.louisville.edu/wwwcvip/research/publications/Pub_Pdf/2006_2/CSIFTA%20SIFT%20Descriptor%20with%20Color%20Invariant%20Characteristics.pdf

Good luck,


Implementing Wavelet Transform using Equations I want to implement Wavelet Transform from the scratch, that mean breaking the wavelet transform into its equations to implement in any Programming language. Matlab Comes with built-in functions to implement Wavelet Transform but It is really hard to understand which processes are exactly involved in the implementation if one wants to develop their own functions.

I know there are Low and High Pass filter involved and the other step called , Down-Sampling but I still have so many doubts on how exactly to combine these filters and samplers to design one's own Wavelet Transform function.

The wavelet tranform block diagram looks like this,

but there are so many Wavelet transform like Haar, db1, db2 etc. Which wavelet tranform does this blog diagram defines anyway.

I found the Wavelets for Computer Graphics: A Primer a good introduction to the Haar wavelet and its role in image processing.