Estimation of accelerating target using position measurements only I am currently thinking about approaches to estimating the position and velocity of an accelerating target.

At this time, I have tried a few approaches that work alright. I have tried two variations of Kalman filters, the first where it estimates position and velocity given the following:

Dynamics:
\begin{align}
\dot{x} &= v + \eta_1\\
\dot{v} &= \eta_2
\end{align}
Measurement:
$$y = x + \eta_3$$
where $\eta_i$ are noise from a Gaussian distribution.

The other Kalman Filter approach was based on the following:

Dynamics:
\begin{align}
\dot{a_1} &= \eta_1\\
\dot{a_2} &= \eta_2
\end{align}
Measurement:
$$y_{i+1} = y_{i}+a_1 \Delta t + a_2 \frac{\Delta t^2}{2}$$

I also implemented a filter using local regression based on the last $N$ position samples in the time series.

All of these filters do fairly well, but I still question if there might be something better. I have tried looking for any papers or lectures on the subject but struggle to find anything worthwhile.

One other thing that coworker and myself discussed was potentially using Euler's Equation for Rigid Body Dynamics to add some constraints into a target's ability to turn too rapidly based on using a ballpark inertia matrix (assuming we're going after a small car).

This obviously makes the dynamics nonlinear and would push us to using something like an Unscented Kalman filter. Does anyone know whether this would be worth it? We suspect it could help our estimate not make any unrealistic jumps that it might otherwise make with the purely linear model, but don't know if this would make a big enough difference.

I don't have a really good answer for you, but I can see a possible problem with your second attempt:  the value of $a_3$ has the potential to completely blow up the dynamics because it is multiplying $t^3$.

In cases like this, you are better off defining an interval over which you wish to do the estimation, $[T_1,T_2]$ and then choosing well-behaved polynomials on this interval (e.g. the Legendre polynomials, $P_k(x)$)  and multiplying them.

That way your measurement becomes
$$
y(t) =  \alpha_1 P_0(x;T_1, T_2) + \alpha_2 P_1(x;T_1, T_2) + \alpha_3 P_2(x;T_1, T_2)
$$
where $P_k(\cdot;T_1, T_2)$ is the (scaled) $k^{\rm th}$ order Legendre polynomial:
$$
P_k(\cdot;T_1, T_2) = \frac{1}{T_2 - T_1} P_k\left( \frac{2t -T_1 - T_2}{T_2 - T_1} \right)
$$
where $P_k(\cdot)$ is the standard Legendre polynomial.

I've used the technique in a formulation of the EKF for frequency tracking.


Undecimated DWT vs. CWT. In what cases is one preferable over another? Talking about undecimated DWT let's focus here on SWT algorithm.

SWT is said to be:

numerically more stable due to Mallat algorithm. And the inverse operation is capable of returning the original signal which is not the case with CWT
SWT is computationally faster than CWT

On the other hand:

SWT works on dyadic scale while with CWT you can choose any scale sampling

Are there any other factors supporting the choise of one or another alternative?

Are there any other factors supporting the choice of one or another alternative?

possibility to use almost any wavelet shape (and associated sampling) or signal-matched: in favor of the CWT
possible parallelization in cycle-spinning form: in favor of SWT
tunable redundancy factor: in favor of the CWT
possibility to derive sharp signal/noise estimators for denoising: in favor of SWT
geometrical extensions (direction, shear, anisotropy) in higher dimension: in favor of the CWT


Obtaining power spectral density without having knowledge about frequency bandwidth or sampling frequency of the signal in Matlab? Is there a way to plot PSD of a signal which show true frequency amounts ( not normalized one) without having knowledge about frequency bandwidth or sampling frequency of the signal in matlab ?

Without more information about the signal, it is not possible.

The set of samples is just a set of numbers with no further information associated with them. If you know the frequency at which the samples were taken, then you can associate "real-world" frequencies to their FFT bins. Otherwise, it may be impossible.

To give an example, consider the samples of $\cos(2\pi 100t)$ taken at $f_s=500\,\text{Hz}$, and the samples of $\cos(2\pi 200t)$ taken at $f_s=1000\,\text{Hz}$. They are exactly the same! There is no way to tell the two signals apart just from their samples.

It may be, however, that you have some extra piece of information that may help deduce $f_s$. For example, you may know that the signal's spectrum has a peak, or a particular phase, or some other property at a certain frequency. If you find the normalized FFT of the samples and identify this property, you can extrapolate from there to estimate the actual signal frequencies.


why total energy of a finite duration continuous signal becomes infinite after sampling If I have a continuous-time energy signal and multiply it with an impulse train, the frequency domain representation has infinite replicas. The total energy if we sum them up becomes infinite. What causes it to have infinite energy now?

I understand the derivation provided in many texts but I am trying to get intuitive understanding.

If you multiply a continuous-time finite energy signal $f(t)$ with an impulse train you get
$$\tilde{f}(t)=\sum_{n=-\infty}^{\infty}f(nT)\delta(t-nT)\tag{1}$$
where $T$ is the sampling interval and $\delta(t)$ is the Dirac delta impulse. Note that the "signal" $\tilde{f}(t)$ is a mathematical fiction, it cannot exist in practice, and it cannot even be evaluated for any value of $t$. Of course it is pretty safe to say that $\tilde{f}(t)$ is zero for $t\neq nT$, but you can't assign any value at $t=nT$. This comes from the fact that the Dirac delta impulse is no ordinary function but a generalized function or distribution, which only makes sense under an integral. So we can do calculations with $(1)$, we can take its Fourier transform, etc., but we cannot evaluate $(1)$ directly.
As for the energy of $\tilde{f}(t)$, since the square of $(1)$ is undefined (because the square of a distribution is undefined), we cannot compute the energy in the time domain. As you know, the spectrum of $(1)$ is the sum of shifted spectra of $f(t)$. Consequently, if you compute the energy in the frequency domain you get a diverging integral. So you could state that the energy of $(1)$ is definitely not finite. But the question remains if that has any meaning, given that $(1)$ is not even a function or signal in the ordinary sense.
Of course, the actual discrete-time signal $f[n]=f(nT)$, which is a sequence of finite numbers (the samples of $f(t)$), does have finite energy. If $f(t)$ is band-limited with maximum frequency $f_{max}<\frac{1}{2T}$ the energy of the sampled signal is equal to the energy of the continuous-time signal $f(t)$ (up to a constant):
$$E_f=\int_{-\infty}^{\infty}|f(t)|^2dt=\frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|^2d\omega\tag{2}$$
$$\begin{align}\sum_{n=-\infty}^{\infty}|f(nT)|^2&=\frac{T}{2\pi}\int_{-\pi/T}^{\pi/T}\left|\frac{1}{T}\sum_{k=-\infty}^{\infty}F\left(\omega-\frac{2\pi k}{T}\right)\right|^2d\omega\\&=\frac{1}{2\pi T}\int_{-\pi/T}^{\pi/T}|F(\omega)|^2d\omega=\frac{E_f}{T},\qquad\textrm{if } F(\omega)=0\textrm{ for }|\omega|>\pi/T\end{align}$$
A different way of showing this can be found in this answer.


Human ear - freq log scale I am not sure if it's DSP question or physic.

But I wonder is there any math equation (method) that describe that graph? (It could be anything approximately like that).

I am working on audio freq analyser. And for typical pop song I get graph that Low End is very big but Middle and High End are very low. That's the mathematical reality, but I want to transform it to human ear reality. Could anyone help me?

(Of course I know how to make freq log scale on horizontal axis and dB scale on vertical axis, just have problem to modify it to human ear characteristic)

So Robert Wannamaker did a real bang-up job parameterizing the 0 dB curve and came up with, what he called F-weighting.  It's in this AES paper from the 90s:

you can see that this can be turned into an $H(s)$, and with the bilinear transform, into an $H(z)$.

here is the curve on a dB vs. log/linear freq scale.  the log frequency plot is the upside-down version of your 0 dB curve above but is more updated from Fletcher-Munson to the ISO standard which cam from i dunno where.

here's a table of poles and zeros:


Actual measurements of the influence of fading and multipath on HF radio communication I'm looking for papers or books that contain actual measurements of the influence one has to expect in HF radio communication. I'm interested in this, because i've simulated the influence of a physical channel using the Watterson Model on the detection performance of a feature-based approach. The approach is not robust for the parameters I've picked. Now I'm interested in how realistic the Watterson parametrzation I've picked is.
For you information and the sake of completeness, the parameters I picked:

The information I'm gathering serves the purpose of anwering the question:
Is it just the average case of distortion that is to be expected in HF communication or is it a more challening scenario.
EDIT:
PSK,QAM and FSK signals were tested. The bandwidth was between 100 and 1000 Hz.The parametrization was chosen for the frequency range of 3kHZ to 30 MHz.

Matlab have in his comm toolbox some scripts dedicated to HF ionospheric channel simulation, check it out


shifting and reversing the signal Say I have a signal $x[n]$.

What is $x[3-n]$?

My solution #1.

$x[-n]$, reverse first.   Then, $x[-(n-3)]$, shift to the right by 3.

My solution #2.

$x[n-3]$, shift right by 3. Then $x[-(n-3)]$, reverse the entire shifted signal.

Both of my solution is incorrect. Why?

Independent variable transformations (such as exemplified in your question) play an important role in mathematical formulation of signal processing concepts and therefore should be understood well and practiced without difficulties before any real mastery of subsequent concepts. Fortunately they're easy to master! as long as a few tricks has been followed carefully. The most important rule is to break a composite transform into basic transform steps, one at a time.

Method-1:  given $x[n]$ to find $x[3-n]$ proceed as follows:

Define $x_1[n] = x[-n]$ , (reversal of $x[n]$ about origin)
Define $x_2[n] = x_1[n-3] = x[-n]|_{n=n-3} = x[-(n-3)] = x[3-n]$ , (right shift of $x_1[n]$ by 3 samples)

Hence $x[3-n]$ can be obtained from $x[n]$ by first reversing $x[n]$ about origin and then shifting the reversed signal $x_1[n]$ to right by 3 samples.

NOTE-1: While decomposing a composite transform such as $n=3-n$ into its simple steps, you can define a new signal at each step such as $x_1[.]$ and $x_2[.]$, which results in a more cumbersome notation but eventually pays off so that you can clearly see the result of every simple transformation independently.

Method-2: The following set of steps could also be used to get the same result: Given $x[n]$ to find $x[3-n]$ proceed as follows:

Define $x_1[n] = x[n+3]$ , meaning a left sihft of $x[n]$ by 3 samples.
Define $x_2[n] = x_1[-n] \rightarrow x_2[n]= x[3-n]$ , meaning a reversal of $x_1[n]$ by about origin.

Hence $x[3-n]$ can be obtained from $x[n]$ by first shifting $x[n]$ left by 3 samples and then reversing the resulting signal $x_1[n]$ about the origin.

NOTE-2: Be careful to observe the following during step-2: $$x_1[n]|_{n=-n}=x_1[-n] = x[n+3] |_{n=-n} = x[(-n) + 3] = x[-n+3] = x[3-n]$$ and not $$x_1[-n] = x[-(n+3) ],$$ a commonly practiced mistake...


Theory behind Linear Predictive Coding (LPC) What is the theory behind LPC?
Why are(were) certain implementations of LPC said to be more tolerant of transmission or encoding errors quantization than other compressed voice encoding schemes?
Can LPC methods also be used for smoothing or short term "prediction" similar to the use of Kalman filter methods?
Under what conditions or constraints is the use of LPC valid?

First, to say linear predictive coding (LPC) is "more tolerant of transmission or encoding errors" isn't entirely true. The form in which the coefficients is transmitted makes a big difference. For example, if the linear prediction coefficients are solved for, they can be very sensitive to quantization, much like high order IIR filter coefficients (this is because the synthesis filter will be IIR, but more on that later). However, if they are transmitted in some other form, this problem can be mitigated easily.

One way is to transfer the reflection coefficients. If you recursively solve for a k-th order linear prediction filter, the highest order coefficient at each stage is called the reflection coefficient. These can be used together to completely characterize the system (which can easily be seen from the Levinson recursion). In fact, you can use all of them together to form a lattice filter. These filters are often used when quantization is a concern, as they're much more robust to low bit counts. In addition, if the magnitude of these reflection coefficients are bounded by unity, you're guaranteed a BIBO stable filter which is critical for LPC where the filter is used to synthesize your signal. There are other methods such as line spectral pairs which are frequently used, but aren't as intuitively defined as the reflection coefficients in relation to AR modeling for LPC.

Now, to address the first question, the theory of LPC revolves around vocal tract modeling. Essentially, we're modeling speech as air vibrating as an input to a tube of some structure. You can look for some resources that go into much more detail to flesh out this model (length of tubes, intensity of air, structure etc). These resources relate these structures directly to IIR filters responding to various stimuli, white noise for example.

So when we solve for the linear prediction coefficients, we're looking for the coefficients such that if we input our signal (voice for example) into an FIR filter created from the coefficients, we get white noise as an output. So think about what that means. We're inputting a highly correlated signal, and outputting a white noise sequence. So in effect, we're removing all linear dependence of that signal. Another way to look at this, is that all of the meaningful information is contained in the coefficients that remove this linear dependence. Therefore, we can transfer these coefficients (or some form of them as above), and the receiving end can recreate the signal. This is done by inverting the linear predictive FIR filter to create an IIR filter, and inputting white noise. So the compression comes in from removing this linear dependence, and transferring the coefficients. This is why the Burg method is also sometimes referred to as the maximum entropy method, as it aims to maximize the "randomness" or whiteness of the output noise in the linear prediction filter. Another way to look at this, outside of physical modeling, is that vocal spectrum (with exception to unvoiced sounds such as hissing) tend to be very impulsive, making them excellent candidates for AR modeling.

To answer your final question, I'm not sure what you're asking entirely. LPC, or linear predictive coding is meant to "compress" the signal assuming that it can be efficiently modeled as previously discussed. You can certainly use linear prediction to do "shortterm prediction" as you've mentioned. This is the implicit basis behind the high resolution AR methods used for power spectral density estimation. The autocorrelation sequence may be recursively extended from it's finite form from the limited data record to infinity as the theoretical autocorrelation sequence of the unwindowed sequence. This is also why AR methods of PSD estimation don't exhibit sidelobe phenomena.


How do I prove that delta - sinc function is the same as an (-1)^n times the sinc $$\delta(n) - \frac{1}{2} \mbox{sinc} \left(\frac{n}{2}\right) = (-1)^n \frac{1}{2} \mbox{sinc} \left( \frac{n}{2} \right)$$

The picture shows what I've tried

Hint : Use the fact that, $\sin(-\frac{n\pi}{2}) = (-1)^n \sin(\frac{n\pi}{2}), \forall n \in \mathbb Z$


Periodogram frequencies not matching with signals I am new to signal processing world.  I am working on a project where we have to find change points using spectral densities.  I am running some simulations where in I created a wave which is a linear combination of a few sinusoidals of frequency ranges (3.5-7.5, 7.5-9, 9-12, 12-14).  I find the periodogram of these sinusoidals using the usual formula which goes from 0 to pi.

My understanding is, I will be able to observe peaks in the periodograms which corresponds to the frequencies of the signal.  But since the x axis in a periodogram ranges only from -pi to pi, how can I compare my input frequencies and the frequencies I got in periodogram?

The output you get from the periodogram are the power spectral density for the angular frequencies W in the range from 0 to pi. With the sampling rate Fs you can compute the frequencies f from W with: f = W * Fs / (2pi).

Be aware that you can only display frequencies from 0 to Fs/2 due to the Nyquist criterium which states that you need at least two points of a signal to be able to detect it.


Formula of frequency modulated sine wave from data points I'm trying to determine the formula of a variable frequency sine wave based on the values of max1, max2, min1 and min2. y(0) isn't necessarily 0 and I can assume the frequency increases or decreases linearly (I only want to plot the first couple oscillations).

I've figured the formula for the non-modulated signal (based only on max1 and min1) as the following:

$y(x) = sin(2π*f*x + C)$
where
$f=\frac{1}{2\left(\max _1-\min _1\right)}$
and
$C = -\left(max_1-\frac{\left(max_1-min_1\right)}{2}\right)\cdot \frac{\pi }{\left(max_1-min_1\right)}$

I've also figured that plotting:

$y(x) = sin(2π*(f+M(x))*x)$ where
$M(x)= a*x + b$

results in an increasing frequency signal but I can't figure out the exact formula to have the wave match the data points.

Any help would be much appreciated, my math skills are rusty...

If your function is $y(x)=\sin(2\pi fx+C)$, and you have data points $\{x_i,y_i\}$, then finding $C$ is easy, if you suppose $x=0$ then $C=\arcsin(y(0))$.

Now, according to your plot your frequency varies with $x$. The first thing you need to do is to know how does your frequency vary with $x$. Supposedly this variation is linear, then you write: $f(x)=ax+b$ and your function model becomes $y(x)=\sin(2\pi f(x)x+C)$.

The problem now is finding (an estimate for) $f(x)$. You have two unknowns to determine $a$ and $b$, you need two data points ($\{x_1,y_1\}, \{x_2,y_2\}$) and you solve a system of two (nonlinear) equations with two unknowns.

$$\begin{cases}
y_1=\sin(2\pi (ax_1+b) x_1+C)\\
y_2=\sin(2\pi (ax_2+b) x_2+C)
\end{cases}$$


2D-otsu algorithm in Digital image processing I'm working on image project using Matlab & I heard about 2D-otsu algorithm as effective segmentation with low SNR and low contrast images, but I can't find enough information about it or how to implement it. I'll be thankful if anyone can help me in such case.

The Otsu algorithm most responders are referring to, is the 1D algorithm indeed, exactly for the reason mentioned above. The main drawback is that it does not work very well in images with low s/N ratio or low contrast. Also, when you only have a small object with a big background, Otsu will incorrectly select background pixels.

The 2D Otsu provides a solution to this. I don't know the maths involved, but it performs '2D histogram projection', I found it was well documented in:
J. Zhang, J. Hu, 2008 International Conference on Computer Science and Software Engineering
From what I understand it projects the histograms in such way and then applies some filters to it so that the peaks relating to the 'object' will appear more profound in the resulting histogram, enabling better segmentation.

Both algorithms can be applied in 2D and 3D ( to add some more confusion).


What is about the circular convolution in OFDM In an OFDM system, serial-to-parallel conversion for data is done, then the DFT is performed and then adding the cyclic prefix (CP).

My question is related to that step of adding a CP. As I know, adding the CP is
My question, if $H$ is the vector of $y$ in frequency domain, "$H = Fy$", so Can we diagonalize $H$ by left-right multiplying by proper DFT?  Which mean $A = F \, H \, F^H$ and $A$ is a diagonal matrix of $N \times N$ ??   Does that right in OFDM? And what's about after removing the CP when performing DFT operating in receiver side, will this propriety of circular convolution will be kept?

thank you

Assuming you have $N$ symbols to transmit encoded in block $k$,
\begin{equation}
 s(k)
 =
 \begin{bmatrix}
  s_1(k) \\
  \vdots \\
  s_N(k)
 \end{bmatrix}
\end{equation}
Performing $N-$FFT at the transmitter, you now have
\begin{equation}
 x(k) = F^Hs(k)
\end{equation}
It is here that you insert a cyclic prefix of size $D > L - 1$, where $L$ is the number of taps in your FIR channel, i.e. at the receiver you have
\begin{equation}
 r(k)
 =H_{ISI}\begin{bmatrix}
  \overbrace{x_{N-D+1}(k)} \\
  \vdots \\
  \underbrace{x_N(k)}_{\text{CP}} \\
  x_1(k) \\
  \vdots \\
  x_N(k)
 \end{bmatrix}
 +
 H_{IBI}
\begin{bmatrix}
  x_{N-D+1}(k-1) \\
  \vdots \\
  x_N(k-1) \\
  x_1(k-1) \\
  \vdots \\
  x_N(k-1)
 \end{bmatrix}
\end{equation}
where $\lbrace c_k \rbrace_{k=0}^{L-1}$ are the channel taps. Note that the channel matrix is decomposed into two parts, ISI (intersymbol interference) and IBI (interblock interference)
\begin{equation}
 H = H_{ISI} + H_{IBI}
\end{equation}
i.e.
\begin{align}
 H_{ISI}
 &=
 \begin{bmatrix}
  c_0 & 0 & \ldots & 0\\
  \vdots & \ddots & \ldots & \vdots \\
  c_{L-1} & \ddots &  & \vdots \\
  \vdots & \ddots &  & \vdots \\
  0 & \ldots & c_{L-1} & c_0
 \end{bmatrix} \\
 H_{IBI}
 &=
 \begin{bmatrix}
  0 & \ldots & 0 & c_{L-1} & \ldots & \ldots & c_1 \\
  \vdots & \ddots & &  &  \ddots &  & \vdots \\
  \vdots & \ddots & &  &  \ddots &  & c_{L-1}\\
  0 & \ldots & &  &  \ldots &  & 0\\
  \vdots & \ddots & &  &  \ddots &  & \vdots \\
  0 & \ldots & &  &  \ldots &  & 0\\
 \end{bmatrix}
\end{align}
If we remove the first $D$ samples of $r(k)$, we have
\begin{equation}
\tilde{r}(k)
=
 \begin{bmatrix}
  r_{D+1}(k) \\
  \vdots \\
  r_{N+D}(k)
 \end{bmatrix}
 =
 \begin{bmatrix}
  0 & \ldots & 0 & c_{L-1} \ldots & c_0 & 0 & \ldots & 0\\
  \vdots &  &  &  \ddots &  & \ddots &  & \vdots \\
   0 & \ldots &  &  \ldots & 0 & c_{L-1} & \ldots & c_0\\
 \end{bmatrix}
\begin{bmatrix}
  x_{N-D+1}(k) \\
  \vdots \\
  x_N(k) \\
  x_1(k) \\
  \vdots \\
  x_N(k)
 \end{bmatrix}
\end{equation}
After a little bit of manipulation, you can see that
\begin{equation}
 \tilde{r}(k)
 =
 C F s(k)
\end{equation}
where $C$ is a circular matrix that contains
\begin{equation}
 \begin{bmatrix}
  c_0 & 0 & \ldots & c_{L-1} & \ldots c_1
 \end{bmatrix}
\end{equation}
in its first row. Now, applying a DFT, we get
\begin{equation}
 y(k)
 =
 F
 \tilde{r}(k)
 =
 F
 C F^H s(k)
\end{equation}
Since $C$ is circular, then $F
 C F^H$ is diagonalizable because any circular matrix is diagonal in a Fourier basis, hence you get
    \begin{equation}
  y(k) = D s(k)
 \end{equation}
    where $D$ is diagonal containing the eigenvalues of $C$.


Properties of a custom wavelet family for continuous wavelet transform First of, I'm new to signal processing. I have a signal which is a linear composition of several basis signals, whereas the same basis signal can occur several times, that is translated, but not scaled. These basis signals look very similar to wavelets. I'm thinking of using a continuous wavelet transform to get the decomposition of the signal with respect the basis signals. However, I assume that a set of basis signals needs to fulfill certain properties to be considered a wavelet family. Could you please point out to me what these may be?

It sounds like you don't need a wavelet. A wavelet's translation and dilation properties allow you to get both time and frequency resolution.

A matched filter takes an inner product of a signal with another signal, allowing you to detect a known signal in noise (by correlating the signal you have with the known signal). You may want to use this idea instead.

Another idea is matching pursuit type algorithms -- if you have an overcomplete set of signals, and you want to find an approximation of a signal with respect to these overcomplete set of signals, you can do this through a matching pursuit algorithm. Basically, iteratively form a residual by finding the largest inner product signal in the set and then subtracting a scaled verison of the signal.


IIR Adaptive Filter in MATLAB Suppose I have a IIR filter represented by $$G_0\left(z\right)=\frac{1}{1-0.2z^{-1}-0.1z^{-2}}$$
I would like to use the LMS algorithm to model an FIR filter $G\left(z\right)$ of order $N = 15$ such that it would adaptively reach $ G_0\left(z\right)$ coefficients values
So we know that $G_0$ can be represented by:
num = 1;
den = [1; -0.2; -0.1];

my issue is that I don't know how to initiate $G(z)$ given its order is $N=15$.
Is it valid to say:
G(z) = zeros(N,1)??

because my confusion is that if $G_0(z)$ has only 3 coefficients, how would $G(z)$ converge to 3 coefficients given that it is 15 taps long?
Please give me some clear insights into this. Thank you in advance.

An adaptive FIR filter is a FIR filter, that uses some kind of an adaptive algorithm to change the filter weights and reach a desired state. In case of using an LMS algorithm the general update equation is the following:
$$
\mathbf{w}(n+1) = \mathbf{w}(n) + \mu\cdot e(n) \cdot \mathbf{x}(n),
$$
where $\mathbf{w}(n)$, $e(n)$ and $\mathbf{x}(n)$ are respectively filter weights, an error and tap values at $n^{th}$ sample and $\mu$ is an adaptation step size.
That means, that the result of adaptation of a $N$-tap filter will be $N$ new values of its weights.
In your case, you can make a FIR filter $G(z)$ with $N=15$ try to adapt, such that its impulse response will be close to the impulse response of the IIR filter $G_0(z)$. And it should be able to do so starting from zero weights. However, after the convergence it will certainly still have an error.
Useful references on the LMS algorithms are Chapter 9 "Adaptive Filtering" in Statistical Digital Signal Processing and Modeling by M.H. Hayes and a Least_mean_squares_filter entry on Wikipedia.


Spectrogram from scipy.signal with Python To test the python spectrogram (from scipy.signal) , I've created a signal with 2 harmonics: 2 Hz and 8 Hz. Then I've added 50Hz noise and a trend (https://i.imgur.com/t1SnQck.png). I wanted to retrieve those frequencies with spectrogram. I only find 10 Hz and about 50 Hz with the noised signal (https://i.imgur.com/I0hMtAC.png). I tried on the original signal without trend and without noise. I only find about 10 Hz (https://i.imgur.com/MVokSSS.png). Do you have any explanation of this ?
The code is only 3 lines ("new_sa" is the noised signal with the trend and "PerEch" period of sampling).

from scipy.signal import spectrogram
f,t1,Sxx = spectrogram(new_sa,1/PerEch)
plt.pcolormesh(t1,f,Sxx)
plt.colorbar()

Thank you

The input to your spectrogram is too short in terms of time.  To get 2 Hz plot resolution, you need an input vector long enough to represent at least 0.5 seconds (preferably 3X more) worth of samples at your sample rate.


FFT has unexpected DC component I have a mixture of Gaussians and I want to look at their power power spectrum. The spatial distribution looks like this:

It's already been convolved with a Gaussian window function. I subtract the mean over all values and run an FFT in matlab.
For the power spectrum I used $|\mathcal{F}|^2 $ and the image under is a section through the middle of the 2D phase space.

The zero phase component is orders of magnitude bigger than the rest of the spectra; given the nature of the problem I expected the power spectral density to look like a combination of Gaussians.

Edit:
The data is initially given as a bivariate poisson process, defined on $ \mathbb{R}^2$. The weights are assumed to be equal and the Gaussians are equally spaced on a hexagonal grid. Now, I'm supposed to estimate the variance of the Gaussians, and wanted to do it using the empirical characteristic function,because the modified EM doesn't work well with data (I guess it's because sample is relatively small ~300 points in total). Following supervisor's suggestion, I estimated a density using Gaussian kernels, and then ran the FFT.
Concerning the Matlab implementation:
The 2D spatial density matrix is 'x'

nfft= 2^(nextpow2(max(size(x))));
fftx = fft2(x,nfft,nfft);
mx = abs(fftx)/length(x);
mx = mx.^2;

I selcted arbitrarily a section and plotted plot(ffshift(mx(2048,:)))

Edit2: I added a link to the data.
https://www.dropbox.com/s/y5sgz0ataxe9z3s/FFT_problem_DataSet.mat

The first and the obvious algorithm for the solution of this problem is the EM algorithm. This algoritm is getting worse in the number of parameters. It is also too much dependent on the initial conditions.

The drawbacks of this algorithm can be circumvented by parallel computing. This means you have $N$ EM algorithms run simultaneously which have different initialization points. The final algoritm takes all parameters estimated by N processors and based on the wellness of the fit gives the final decicion.

Another method is to use the central moments until Mth order. When written in a parametric form each xth order moment gives a linear/non-linear combination of parameters which are to be estimated. From here one can construct a test. Of course the moments should also be estimated from the data.

Another method would be to employ a better search algorithn. That is called particle swarm optimization. It will run similar to EM. However it is a stochastic algorithm have different characteristics than EM such as to be able to deal with more number of parameters. It is a kind of team optimization algorithm.

Thats all in my mind at the moment. Good luck with the experiments.


fliter doesn't filter In this code I try to decompose music from singer sound with wfilters but when I hear it it is still mixed where Iam I wrong ?

 function s = scale3(q,n)
  p = length(q);
   if p==n
   s = q;
    else
    s = resample(q,n,p);
   end

the program

s = wavread('K:\Khaleel_3aaza (1).wav'); % whatever here

 [l,h] = wfilters('sym2','d');
 [z.x]=wavedec(s,2,l,h);

p = scale3(z.x(2245001:22350000),2245000);

Background music and human voice occupy the same frequency ranges, so they cannot be separated by filters. See this Q&A for a few suggestions of techniques to separate human voice from background music.


Signal equation using signal waveform and Fourier series? Part A: Obtain signal waveform from mathematical equation of the signal

Let a sinusoidal periodical signal is represented by an equation

$$y=f(t)=10+10\cos\left(\frac{2\pi f_1t}{T} +\frac{\pi}{6}\right)+10\sin\left(\frac{2\pi f_2t}{T}+ \frac{\pi}{3}\right)\tag{1}$$.

Here, let us take $f_1=10 ,f_2=5 ,T=100$

Now, with this we can get the waveform of the signal $y$ on Cartesian axis with magnitude as $Y$-axis and time as $X$-axis.

(Sorry but I don't have any software tool now so that I could show the waveform diagram)

Part B: Obtain mathematical equation of the signal from its waveform and Fourier series formula.

i.e. part B is reverse of part A.

Here, my question is that how can we obtain (probably the same) mathematical equation $(1)$ for signal $y$ from waveform of signal $y$ and Fourier series formula
?

Fourier series formula for any sinusoidal periodical signal $y$ is given as

$$ y(t) = \sum\limits_{m=0}^{+\infty} a_m \cos(w_m t) + \sum\limits_{m=0}^{+\infty}b_m \sin(w_m t). \tag{2}$$

Did you check the usual equations?

Responding to:

  The steps I should follow are

  Find $f(x)$ from the signal waveform
  Compute $a_0$, $a_n$, $b_n$
  Put above values into the (2) and see whether I get (1) or not.

Essentially, yes.

For your equation to relate to the $f(x)$ in the image above you may want to think about
$$
w_m = \frac{2\pi 5 m}{100} = \frac{\pi m}{10}
$$

Also note that there needs to be an integer (or rational) ratio between $f_1$ and $f_2$ in your equation (1).


DFT provides coherent integration? I read from material that coherent integration can be provided by DFT. But it doesn't give further explanation on this remark. As far as I know, coherent integration is the summation of discrete samples with phase considered and it is usually used in time domain. How to understand that DFT's capability of providing coherent integration?

Thanks.

@ecook: Your time-domain notion of coherent integration is correct. (Where multiple periodic sequences must be time-synchronized before the averaging operation.)

The only way I could justify the remark that "coherent integration can be provided by DFT" is in the following scenario: let's say your DFT's input is an $x1(n)$ sinusoidal sequence and your DFT result is the complex valued $X1(m)$ sequence. Then some moments later you collect a new $x2(n)$ time sequence of the same sinusoidal signal and produce a new complex-valued DFT result of $X2(m)$. To average (integrate) those $X1(m)$ and $X2(m)$ sequences you would add their magnitudes. That's a form of incoherent averaging (incoherent integration), also called "scalar averaging." (Averaging $X1(m)$'s and $X2(m)$'s real and imaginary parts individually may lead to spectral energy cancellation and you don't want that to happen.)

If, on the other hand, your $x2(n)$ time sequence was time-synchronized with the previous $x1(n)$ time sequence you can add $X1(m)$'s and $X2(m)$'s real and imaginary parts individually. That's a form of coherent averaging (coherent integration), also called "time-synchronous averaging."

Know this however, as Peter Kootsookos points out, averaging phase angles is tricky. Have a look at Peter's explanation at:
http://www.dsprelated.com/showarticle/57.php


Applying FIR filter to data with different sample rates I have a filter design, and it filters over a 1-2 kHz range.
What should I do if I want to apply it to data with a different sample rate than the one for which it was designed?

Let's say it consists of Bessel and Chebyshev filters. How do I find a function that defines each filter's coefficients at an arbitrary sample rate? Or should I do this by hand?

Since you mention sampling, you are presumably talking of a digital filter.

The cut-off frequency or half-power frequency of a digital filter is
actually relative to the sampling frequency $f_s$.  If your digital
filter is passing signals in the $1$ kHz to $2$ kHz range when you
feed it a signal sampled at $f_s = 20$ kHz, then the pass band is from
$5\%$ to $10\%$ of $f_s$.  These ratios do not change if
$f_s$ changes to some other value, say $40$ kHz.  The same digital
filter will become a filter with passband $2$ kHz to $4$ kHz without
your having to do anything.


What Kind of Features Can I Extract from a Signal I want to monitor (automatic-)gearbox failures on some vehicles. For each vehicle I have a captured signal representing the selected gear at each one millisecond. An example of two signals are shown on the figure bellow.

I would like to compute extract some features from each signal (and represent it as a feature vector) in order to compare a specific signal to other signals (compute its distance to other signals) and see if it is deviating from most of the other signals.

Do you have any idea on how to represent this type of signals as feature vectors ? Which kind of features can we extract from those ?

Some Features:

Mean.
Variance.
Skewness.
Kurtosis.
Dominant 3 frequencies in the DFT.
Energy of the 3 dominant frequencies.
Max Value.
Min Value.
Median.
Total Variation.

Usually I'd compute them in running windows.
Another great information is the Histogram of the Derivative.
Or just all the above of the Derivative.


Is linear convolution same thing as aperiodic convolution? I am reading an online article and what i understood from there is that linear convolution is aperiodic
Is it right?
I have also attached snapshot and higlighted two instances

What they call noncircular or aperiodic convolution is more commonly called linear convolution. Linear convolution describes the input-output relation of linear time-invariant (LTI) systems. Linear convolution can be computed via circular convolution if we use sufficient zero-padding. Circular convolution is important because it can be computed using fast algorithms (FFT).


Second order section digital Bandpass Does anyone have a ressouce on how I can calculate second order section bandpass coefficients?

What I tried allready:

I used following formulas from the Book Electronic Circuits (U.Tietze, Ch. Schenk) for the Lowpass:
My coefficients match exaclty (difference is max 2e-14) with the second order sequions of MATLAB for the same lowpass.

For a bandpass filter I used following formulas:

The results differ quite alot from MATLAB, expecially if the order is higher than 2.
Also I dont think I can apply the formulas direcly for second order sections because the coefficients only depend on constant values.

So it seems that the formulas from the textbook for the bandpass are not accurate.

I used following approach to calculate second order sections given an analogue filter lowpass. As an example following code displays the coefficients for a 4th order bandpass, calculated from a second order lowpass protopye.

Paramters

fs = 1000;
fs1 = 200;
fs2 = 300;
w1 = 2*pi*fs1/fs;
w2 = 2*pi*fs2/fs;
Omega1 = 2*tan(w1/2);
Omega2 = 2*tan(w2/2);
syms s z Om2 Om1 a1 b1 A0
% apply LP to BP transformation
G_s_BP = subs(H_s, s, (s^2 + Om1*Om2)/(s*(Om2 - Om1)));

Start with analogue transfer function H_s

H_s = A0 / (1 + a1*s + b1*s^2);
% apply LP to BP transformation
G_s_BP = subs(H_s, s, (s^2 + Om1*Om2)/(s*(Om2 - Om1)));

Digital conversion with bilinear transformation

G_z_BP = subs(G_s_BP, s, 2*(( z - 1)/(z + 1)));
% order polynomial such that coefficients can be read out
coeffs_z = collect(G_z_BP, z)

Those coeffs_z are very similar to the ones given by the iowahills software
The band-pass filter coefficients then can be calculated with any analogue low-pass prototype, just enter numerical values on read out coefficients from coeffs_z


How can I design a very narrow filter? Suppose I've got an audio signal sampled at $48000$ Hz, and I'd like to design a low-pass filter that isolates everything below ~$60$Hz. In the digital world, this is a low-pass filter with the passband at $[-\frac{\pi}{400} , \frac{\pi}{400} ] $. Also, the transition band should be reasonable as well. Building a FIR filter for this can have a lot of taps which in the long run affects precision. An IIR filter is too not ideal because audio suffers for non-linear-phase response in filters, so unless the signal is filtered, then reversed and filtered again, it is not really an option.

Could a wavelet transform be better at this than in-one-go regular filtering?

If you are optimizing for engineering time and are on a platform that supports large FFTs well (i..e not fixed point), then take hotpaw2's advice and use fast convolution.  It will perform much better than a naive FIR implementation and should be relatively easy to implement.

On the other hand, if you have some time to spend on this to get the best implementation or are on a fixed-point platform, you should use a multirate down-filter-up-subtract structure.  But it's a bit trickier to get everything right.

I've got access to trusted and highly optimized implementations of both fast convolution and multirate filtering tools. The fast convolution takes about 3x longer to get equivalent signal performance compared to the multirate structure.  Furthermore, that is even on a floating point platform. The gap would widen considerably on a fixed point dsp.

In general terms:

Down-conversion:

Use 8 stages of halfband,decimate-by-2 filters to transform your 48kHz signal into a 187.5 Hz signal.  The first stage of this downsampling can have a very wide transition band, allowing energy to alias as long as it doesn't alias back into the sub 60 Hz range. As the stages progress, the number of taps needs to increase, but they will be applied at a progressively lower sampling rates, so the overall cost of each stage remains small.

Filtering:

Perform your tight filtering around the 60 Hz bw to keep the energy you will eventually want to subtract.  There is a double advantage to doing your tight filtering at the low rate:

1Hz of transition bandwidth is 256 times larger in terms of digital frequency at the low rate vs. the original rate. So every tap of your filter is 256 times as powerful.
The signal itself is at a lower rate, so the filter only needs to process 1/256 the data.

Up-conversion:

Essentially, this is the reverse of the decimation stages.  Each of the 8 interpolator stages doubles the rate by estimating the sample that goes between consecutive input samples.  The transition band gets wider as the sample rate gets higher.

Subtract:

Subtract your full-rate low-pass filtered signal from the original signal. If you've adjusted for all the group delays properly, the overall structure will be a highpass filter with a narrow transition bandwidth.


Taking audio frequency and picking out peaks (Peak Detection) So I have a unique challenge. I want to take an audio file, and pick out the peaks in the samples and plot that onto a graph. Each point based on the amplitude/pitch(?) will be represented as either a small dot (lower amplitude) or a bigger dot (higher amplitude).

So if I had to look at a drum track, I should notably see big dots where the drum is hit, but smaller dots when a cymbal is hit.

From a technology stack perspective, I am using Python - and using a Library called Aubio (although I am not sure if there is a better library out there).

EDIT:

Forgive me on the lingo - very new at this audio stuff.  Okay so the issue I have currently is that I have a waveform graph that shows the samples over time using a standard sampling technique.

But I found an article on Peak detection, which is exactly what I am looking for (I think).

The challenge I have now is pulling those peak values out via Python and Aubio and plotting them on a Graph.

I found the answer finally. I found a great article that explains many different libraries that can be utilized for peak detection.  I now have the peaks I am really interested in, and can now create the output I am requiring.

Finding Peaks in Python


Scipy.signal noise with rfft compared to fft I'm trying to get the fourier transform of a signal with real values, however the results I get with rfft are noiser than those with fft.

I wrote the following code:

import numpy as np
from scipy.fftpack import rfft, fft
import matplotlib.pyplot as plt
# Number of sample points

N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
f1 = 50
f2 = 80
y = np.sin( f1* 2.0*np.pi*x) + 0.5*np.sin(f2 * 2.0*np.pi*x)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
plt.plot(xf, 2.0/N * np.abs(yf)[0:N/2])
#plt.savefig('sin.png')
plt.show()

yf2 = rfft(y)
xf2 = np.linspace(0.0, 1.0/(2.0*T), N)
plt.plot(xf2, 2.0/N * np.abs(yf2))
plt.savefig('sin2.png')
plt.show()

and I get the following results:

I thought that for real values I would get the same result with fft and rfft, do you know why there are some differences ?

The DFT of a real sequence is complex-valued. The array output by scipy.fftpack.rfft consists of the real part of the 0th entry, followed by the imaginary part of the 0th entry, followed by the real part of thre 1st entry, followed by the imaginary part of the 1st entry,...

Convert the output of rfft into the appropriate complex array and plot the absolute values of that array, and you will have what you seek:

yf3 = yf2[0:-2:2] + 1j*yf2[1:-1:2];
plt.plot(np.abs(yf3));
plt.show()


Transfer function of a PLL Loop Filter that can support a linearly increasing (chirping) frequency What is a transfer function for a the loop filter of a PLL that can track a sinusoid of the form:

\begin{equation}
x(t)=\cos\left(2\pi\left(\frac{1}{2}c_0t^2 + f_0t + \phi_0\right)\right)
\end{equation}
where $c_0$ is the chirp (units of Hz/s), $f_0$ is the center frequency (units of Hz) and $\phi_0$ (units of cycles) is the initial phase?

I am currently using a loop filter followed by a direct digital synthesize (DDS). The loop filter I am currently using is a proportional plus integrator loop filter with transfer function
\begin{equation}
F_{\text{loop}}(z) = k_1 + \frac{k_2}{1-z^{-1}}.
\end{equation}
The DDS has a transfer function of
\begin{equation}
F_{\text{DDS}}(z) = k_0\frac{z^{-1}}{1-z^{-1}}.
\end{equation}

I would like to replace the loop filter with one that, combined with the DDS, can track the aforementioned sinusoid. The feedback loop is configured in the usual manner, namely, conjugate-multiplying the output of the DDS with the input sinusoid and then the output of that is fed into the loop filter, and, finally, the output of the loop filter is fed into the DDS. (Sorry I don't have a block diagram.)

This PLL works well with no noise present and $f_0 < 10$ and so, is clearly not very useful. I need something that can handle a phase offset, a phase ramp (nonzero frequency), and a frequency ramp (nonzero chirp) simultaneously, in addition to some noise.

NOTE: I have simplified this down to tracking a sinusoid. I'm actually tracking the phase of a modulated signal and so, rather than conjugate multiplying, I feed the match filter and the DDS outputs into a phase error detection function. However, for my purposes here, this simplified approach should suffice.

(I thought about posting this to the Electrical Engineering stack exchange, but since I'm working with a completely digital implementation, my question felt more natural here.)

To track a frequency ramp with a Phase lock loop, with zero steady state error requires a type 3 PLL Loop; which means three integrations (DC Poles) in the open loop gain (your NCO would be one of the integrators and your loop filter needs to provide the other two). Stabilizing such a system becomes more challenging but here is one reference paper detailing an approach to the design:

http://vbn.aau.dk/files/72548181/Final_PLL.pdf

This picture helps clarify tracking and loop order; given a phase lock loop, phase is the parameter being tracked. A frequency ramp as you are dealing with is an accelerating phase, hence the requirement for a type 3 loop.  The detailed loop filter design would depend on your overall system and related gain coefficients for your specific design, but I hope the above reference offers some guidance to help you.

Also note that many people often confuse "Type" with "Order", so to be very clear a Type Three Loop is not the same as a Third Order Loop. Specifically a Type Three Loop must be Third Order (or higher), but a Third Order Loop need not be Type 3. Type specifically is number of integrators in the loop, which are poles at z=1 in the open loop gain, while order is the total number of poles, regardless of location (specifically the order of the denominator polynomial of the open loop gain, which ends up being the order of the closed loop gain given proper transfer functions).


Doppler Shift of Phase Coded Pulse Compression Waveform Suppose I want to transmit the signal:

$$ x(t) = g\Big(\frac{t-t_{0}}{T}\Big) \cdot \cos(2 \pi f_{0} t + \theta_{0}) + g\Big(\frac{t-t_{1}}{T}\Big) \cdot \cos(2 \pi f_{1}t + \theta_{1}) $$

Where $g(\cdot)$ is some window function of duration $T$ and $t_{1} = t_{0} + T$. Suppose that it encounters some target at range $R$ traveling with some velocity $v$, giving rise to a Doppler shift $d$.

My question:

Can I then represent the return signal as (with $c$ the speed of light, and $a,b \in \mathbb{R}_{++}$ constants representing compressing the waveform in time):

$$ x(t) = g\Big(\frac{t-t_{0}-\frac{2R}{c}}{aT}\Big) \cdot \cos(2 \pi df_{0} t + \theta_{0}) + g\Big(\frac{t-t_{1}-\frac{2R}{c}}{bT}\Big) \cdot \cos(2 \pi df_{1}t + \theta_{1}) $$

If so, what are these constants a and b?

I suppose what I really want to find out is:

How exactly does the Doppler shift compress these pulse compression waveforms? Will the chips still be connected? ie It just compresses the entire waveform in time?

Or will there be a gap in the middle because the Doppler shift compresses each chip from the chip's midpoint? (thereby leaving a gap between chips that did not exist before)

So this is something we run into with the radar community on a pretty regular basis.  For LFMs, the solution is a bit trivial, but for coded waveforms it becomes a real issue.  For us radar folk, we're dealing with targets that are at escape velocity (right around 11 km/s typically). I actually have some pending intellectual property being released via my company on this subject, so I'll try to discuss the solution without going too far into it.

Depending on the direction of the velocity vector, the Doppler distortion on an LFM will result in either an elongation of the waveform, or a "squishing" of the signal.  It's important to note here that for the LFM case, this pretty much just results in an effective chirp slope change, so I like to think of it as the start/stop frequencies have remained the same, but the effective pulse duration it took to get there has been distorted.

The solution for an LFM is simply just to account for the effective slope change that you expect. For example, if you're expecting an elongation by a factor of 2, "squish" your waveform by a factor of 2 when you transmit the waveform; we call this process pre-distortion, and the goal is to equalize the distortion effect from the target such that at your pulse compression filter, you get maximum pulse compression gain.  This of course relies on some a priori knowledge of the target velocity, which isn't totally unreasonable given how most of these systems perform tracking operations.

Coded waveforms present an issue though, we're now non-linear, so the solution isn't as easy on the surface.  If you chunk through the algebra a bit using generic phasor notation, you'll start to see some terms pop out.  It turns out the result is fairly similar: the chip widths change.  How the chip widths change, I'll leave that for you to figure out with your algebra.

The important thing to note here though is that, at least in my experience, you do not punch holes/zeros in your signal like you've asked.  The coded waveform chips will now either be shorter or longer, so you may have added or lost some phase information, but if you can account for that on transmit just like we did with the LFM, you can negate the distortion via pre-distortion of the waveform.  In short: the chips will change their duration, but there will not be gaps since all of the chips change.

Now as for any Doppler shift that's on the phase coded waveform, that can easily be taken care of if you have some expectation of target velocity: you simply just shift the pulse compression filter in frequency (or the received waveform itself, whichever is easiest for you to implement in your system) so that the Doppler shift is "matched" to your pulse compression filter.  Again, this requires some knowledge of target velocity.  If you didn't have knowledge of target velocity, you could do a multi-hypothesis approach.

The reason we typically don't have to worry about this for LFMs is evident in their ambiguity function: they are incredibly Doppler tolerant, i.e. Doppler shifts just cause slight range migration and decrease in pulse compression gain, but not signal loss.  With coded waveforms, their ambiguity function typically looks like a spike (often called "thumb tack ambiguity"), which results in complete signal loss with even the slightest Doppler shift.  Therefore, it's pretty important that we get that Doppler hypothesis correct!

This method is in open literature, and I'm pretty sure some of the Mark Richards signal processing books touch on it. This paper from MITRE touches upon a multi-hypothesis technique.  A lot of that paper is poorly written and sometimes incorrect in my opinion, but they have a good explanation of the multi-hypothesis approach, and a nice block diagram for it towards the end.


Missing Data in Machine Learning Classification I have a training set of 4000 for ML with 36 features, two-class classifier.

Consider that due to some discrepancy feature no. 25 is missing for half of the data. But during data acquisition I know what class it belongs to.
I have defined a specified value (-1000) for that missing data. Other feature values are relatively smaller than -1000.
If I train the data with the same data set, it might not train the
data with high accuracy. Can I somehow extrapolate the data to create
the data set that may be present (near-by values) or is it fine to
train the model with that data-set? (or any other suggestion).

Thanks

Labelling a missing value as -1000 (assuming a non-missing value is much smaller in absolute value) will not cause problems with anything based on decision trees/stumps (bagged/boosted). You might also have some success with mixture models, for example a GMM - the learning algorithm will allocate some components to cover the missing data. Same for non-parametric density estimators, nearest neighbours, or SVM with RBF kernels.

Note that in all these cases, even if the model works, there will be a waste of resources ("model space budget") if feature #25 is always present in your testing data. In other words, your model will have learnt how to deal with tricky examples that will never occur when classifying new data! But if feature #25 will sometimes be missing in some of the data you will encounter "in the field", this is the way to go.

Do not expect things to work with algorithms such as neural networks or support vector machines with linear or low order polynomial kernels.

If feature #25 will never be missing "in the field", it might be simpler to just extrapolate feature 25 in your training data. To do so, you can proceed as follows:

Option 1: Nearest neighbor

Normalize all features (except #25).
For each training example in which feature #25 is missing, search for the nearest training sample, according to features (1 .. 24 and 26 .. 36), for which feature #25 is present.
Assign to feature #25 the value taken from the nearest neighbour.

Works best when you have lots of data (but not too much so as to make the search feasible) with kinky distributions.

Option 2: Regression

Use any regression method of your choice to find a relation ship between feature #25 and the other features, on the subset of your dataset in which feature #25 is not missing. For example, a linear regression.

Use the result to predict feature #25 on the remaining training instances.

Works best when your features are "well behaved" for the regression method you plan to use (for example, see the prerequisites for linear regression).


Marginal Stability based on Poles We know that a discrete-time system with a (Z-transform) transfer function that has a pole of magnitude 1 (i.e. $|z|=1$ is a pole of the transfer function) is marginally stable if the pole at $z=1$ is a single pole and unstable on every RoC otherwise.

My question is: Suppose a transfer function that has two simple poles at $z=1, z=-1$. Is it marginally stable or not? And, therefore, what is the most general form of the stability theorem based on poles of the (Z-)transfer function?

A causal discrete-time LTI system is marginally stable if none of its poles has a radius greater than $1$, and if it has one or more distinct poles with radius $1$. So a system with poles at $z=1$ and $z=-1$ is marginally stable (if there are no other poles outside the unit circle).

A causal discrete-time system with all its poles strictly inside the unit circle is called asymptotically stable. If at least one of the poles is outside the unit circle, the system is unstable.

For anti-causal systems just replace 'outside' by 'inside' and 'greater than $1$' by 'smaller than $1$' in the above definitions.

For non-causal systems with two-sided impulse responses the region of convergence (ROC) is an annulus centered at $z=0$. Such a system is asymptotically stable if the ROC contains the unit circle. If the ROC is limited by the unit circle and if one or more distinct poles are on the unit circle, then the system is marginally stable.

In sum, we get the following ROCs for marginally stable systems:

causal: $1<|z|<\infty$
anti-causal: $0\le |z|<1$
non-causal (two-sided): $1<|z|<R,\; (R>1)\quad\text{or}\quad R<|z|<1,\; (R<1)$

where in all cases it is assumed that all poles on the unit circle are distinct.


Selecting the Nyquist rate for a combination of signals If several sine waves were combined but each had a separate frequency, for example 10 Hz, 20 Hz and and 30 Hz, what would the required Nyquist rate be for analysis?
From my research, it suggests that Nyquist rate should be double the maximum frequency and for this case it should be 60 Hz. Is this correct or does the fact that it is a combined signal change this concept?

According to the lowpass sampling theorem, the required minimum sampling-rate (a.k.a Nyquist rate) is twice the highest frequency present in the frequeqncy spectrum of the bandlimited signal.
Note however that, the lowpass sampling theorem assumes there is no impulse (no pure sinusoidal component) at the highest frequency (a.k.a Nyquist frequency) of the signal being sampled. If there's an impulse, then you should choose a rate arbitrarily larger than the Nyquist rate.
Therefore, for your combined (linearly superposed) signal, the highest frequqency is 30 Hz; and thus the Nyquist rate is 60 Hz, however, since there's a pure sinusoidal at 30 Hz, then your sampling rate should be arbitrarily larger than 60 Hz.
In practice a 10% tolerance is already applied for various reasons. Hence you may use about 66 Hz for sampling your combined signal.


Weird results when using MATLAB's $\tt butter$ function? I'm trying to learn some of the basic MATLAB filter design functions and have been experimenting with the butter function. This function calculates the transfer function coefficients of a Butterworth filter given a desired filter order and normalized cutoff frequency.

I have previous experiences with designing prototype filters and then modifying them, however, I'm having a difficult time understanding MATLAB's definition of normalized frequency.

I'm used to designing Butterworth filters with a normalized frequency of $1 \textrm{ rad/s}$, yet MATLAB insists that the normalized frequency has something to do with a sample rate? I assumed the butter function was used for designing analog filters?

Furthermore, when I plot the frequency response of the function, the cutoff frequency I specify (Wn parameter) never seems to correspond to the $-3\textrm{ dB}$ mark in the plots.

Any thoughts on these?

The answer to your question can mainly be found in the Matlab documentation. The butter function can design filters either in the analog domain as well as in the discrete domain. In order to make an analog filter, you can do so by giving the string 's' as a third parameters. The transfer function can be defined by the zeros/poles/gain or by numerator/denominator coeefficients. The two available syntaxes are

[z,p,k] = butter(n,Wn,'s')
[b,a] = butter(n,Wn,'s')

When designing a filter in the discrete domain, the normalized frequency is a ratio of the cutoff frequency over the biggest frequency that you can represent according to Nyquist theorem. For instance : if you have a sample time of 1ms, hence a sampling frequency of 1Khz, the biggest frequency you can represent is 500Hz.  Designing a filter with a normalized frequency of 0.5 will put the cutoff frequency at 250Hz, since 250/500 = 0.5

In your question, you specify that the cutoff frequency didn't seems to match the -3dB point. Can you give us more detail about your issue. Maybe show us the work you have done until now ?


Detect lanes from a filtered image I am working on lane detection system and I would like some suggestions from you as to how to detect lanes from a filtered image which I got after performing certain operations. Firstly, I should mention a few things:

The lane detection would be in a real road environment unlike here.
As of now, I am running it in real-time (O(n)) to be precise, which I would like to maintain.

The following image shows the original image and my filtered image.

As you can see, there are false positives which I want to remove and that is my main problem. Now, these are the few approaches I have thought of:

Check the intensity on the left and right hand side for desired areas. It should be high for desired area and low on both sides.
**Problem: I am finding it difficult to define 'high' and 'low' since the method has to be independent of time of the day.
One approach I thought was to check shapes but I can't think of any algorithm since those are not rectangles, those are essentially distorted parallelograms(?).

One more, I have seen people generate lines on these type of images using hough, this same question does that. If I run hough, I am not getting lines. Is the intensity of the image a problem or should I run Hough on binary mask image? I am currently doing this with MATLAB but I will be finally implementing with OpenCV to keep it real-time.

Finding shapes / lines is a good idea. Threshold tweaking is almost certainly a bad idea.

Hough transform typically expects your lines to be (1) in a binary image and (2) of one pixel width. Try to apply an edge detector (Matlab: edge) on your filtered image in order to get the boundaries of your white regions as one-pixel lines. Then Hough transform should find something.

I would also consider doing edge detection on your original image rather than your feature image. Then you can try to tweak the edge detector, rather than first running your "certain operations" and then detecting edges.

Once you have some lines you will face the problem of finding the good ones. You could look at (i.e. google for) "vanishing point detection", and read about RANSAC.


Rician factor effect The Rician factor is the ratio between the mean power of the LOS path to the mean power of all other (non-LOS) paths. Increasing/decreasing this factor means increasing/decreasing the mean power of the LOS path. Does increasing the Rician factor improve or degrade transmission quality?

A strong LOS component, i.e. a high Rician factor, will always improve receiver performance. If $K$ denotes the Rician factor, i.e. the ratio between the power of the LOS component and the power of the other (scattered) components, then for $K\rightarrow\infty$ you have no multipath propagation at all and you get a simple AWGN channel. For $K=0$ you have a Rayleigh fading channel (with no LOS component).


Find the impulse response from the frequency response So I'm having a problem here which gives me the frequency response and asks for the impulse response:

$H(\Omega ) = e^{-j\frac{\pi }{2}}$ for $\Omega>0 $ and
$H(\Omega ) = e^{j\frac{\pi }{2}}$ for $\Omega<0$

I need to find the impulse response. My first thought was to (obviously) apply Inverse Fourier Transform but I'm thinking that since it's a branch function maybe that's not the correct approach. Also, I see that the two forms look a lot like each other so maybe I can somehow combine them into one and THEN apply the Transform. Any thoughts and guidance?

HINT: (because this is a homework type problem)

Using the sign function the frequency response can be expressed as

$$H(\Omega)=-j\,\textrm{sgn}(\Omega)\tag{1}$$

Direct computation of the inverse Fourier transform is tricky because of convergence issues of the corresponding integral.

If you can use Fourier transform tables, find the inverse Fourier transform of the step function $u(\Omega)$ in the frequency domain. If you can't find it, use the Fourier transform of the step function in the time domain $u(t)$, and figure out a simple transformation to obtain from it the inverse Fourier transform of $u(\Omega)$. Then express $\textrm{sgn}(\Omega)$ in terms of the step function $u(\Omega)$.

As a final hint, that system is called a Hilbert transformer.


Book suggestion for upsampling and downsampling I want to ask for a  good textbook or online material  which explains
in detail:

-the time and frequency domain mathematics
behind  upsampling and downsampling.
-easy to follow

Note: Oppenheim's discrete-time signal processing book is hard to follow.

The answer to my question is:
Essentials of Digital Signal Processing 1st Edition
by B. P. Lathi
https://www.amazon.com/gp/product/1107059321/ref=ppx_yo_dt_b_asin_title_o00_s00?ie=UTF8&psc=1
where in Section 6.6 you can find detailed discussions and good examples.


Power Spectral Density Reconstruction I am from an ML background and I am working on MP3 audio reconstruction and I have managed some amount of spectrogram reconstruction but I have noticed that the higher frequency end of the generated Power Spectral Density plot(Wave2) isn't matching with the original audio's plot (Wave1).

There is a similar distribution in the magnitude spectrum which is obvious since the two are linked. I wanted to know if there is any way I could regenerate the higher frequency power (or magnitude) because I cannot map the input PSD data to output PSD data as there is no way to convert the PSD data back to an audio clip even if I am successful in mapping them. The real and imaginary part re-constructions for the spectrogram reconstruction are visually pretty accurate but I am still unsure about how to approach this issue. Any guidance would be great!

A few observations

Your pictures looks odd. "Wave 1" doesn't look like an original audio file but something that has already been encoded/decoded and aggressively lowpass filtered in the process. Your reconstructed audio (Wave 2) actually as MORE high frequency energy, which I suspect is simply just quantization noise from the other bands. The expected spectrum of music is roughly pink (3dB/octave) up to about 10 kHz and then a steeper decline (maybe 10dB octave) to 20 kHz and above. There shouldn't be a step function in the spectrum.
The standard method of recovering for lost high frequency spectrum is Spectral Replication. That's been already implemented in some of the more modern codec standards (e.g. AAC) although not in MP3. See for example http://www.stockis.se/AES112th_SBR.pdf
If you want to do any commercial work in this area, I highly recommend studying the patent situation carefully.  A massive amount of work has already been done in this field and since the license revenue is quite substantial, the work has been well protected by patents and there is incentive to enforce them.


Laplace transform of $f\left(\frac{t - b}{a}\right)$ Consider the function $f\left(\frac{t - b}{a}\right)$. We want want to calculate its Laplace transform. There are two approaches:

Firstly,

let $g(t) = f\left(\frac ta\right)$.
Then $\mathcal{L}\left\{f\left(\frac{t-b}{a}\right)\right\} = \mathcal{L}\left\{g(t - b)\right\} = e^{-bs}G(s)$ and $G(s) = \mathcal{L}\{g(t)\} = \mathcal{L}\left\{f\left(\frac ta\right)\right\} = |a|F(as)$.
Therefore $\mathcal{L}\left\{f\left(\frac{t - b}{a}\right)\right\} = |a|e^{-bs}F(as).$

Secondly,

let $h(t) = f\left(\frac{t - b}{a}\right)$.
Then $\mathcal{L}\left\{f\left(\frac{t-b}{a}\right)\right\} = \mathcal{L}\left\{h\left(\frac ta\right)\right\} = |a|H(as)$ and $H(s) = \mathcal{L}\{h(t)\} = \mathcal{L}\left\{f\left(\frac{t - b}{a}\right)\right\} = e^{-\left(\frac ba\right) s}F(s)$.
Therefore $\mathcal{L}\left\{f\left(\frac{t - b}{a}\right)\right\} = |a|e^{-bs}F(as).$

So far so good. However, let $f(t) = e^t$. Then
$$F(s) = \frac{1}{s - 1}$$
and
$$\mathcal{L}\left\{f\left(\frac{t - b}{a}\right)\right\} = \frac{ae^{-\left(\frac ba\right)s}}{as - 1}$$
from Wolfram Alpha.

The problem comes from the fact that we introduce new data by shifting, i.e., the function is not zero for $t < 0$. What can we do then? Can we not use the time delay property?
What about the Fourier transform? The arguments above are valid for the Fourier transform for $s = j\omega$ and you do not have the requirement of $f(t < 0) = 0$ when applying the shift property. But when you set $f(t) = \cos(t)$ you get:
$$\mathcal{F}\left\{f\left(\frac{t - b}{a}\right)\right\}\bigg\vert_{\omega \ge 0,\ a > 0} = \frac{1}{2}e^{\left(\frac ba\right)j\omega}\delta(\omega - 1/a)$$
So it looks like it should be $e^{-\left(\frac ba\right)s}$ instead of $e^{-bs}$. But where have I made a mistake?

Your two alternate derivations of the Laplace transform (specifically the bilateral I'm referring to) of the signal $f((t - b)/a)$ seems right, both resulting in, the same , as a check.

$$
\mathcal{L}\{f((t - b)/a)\} = |a|e^{-bs}F(as)
$$

However, when you assume the signal $f(t) = e^t u(t)$ whose Laplace transform is
$$
F(s) = \frac{1}{s - 1}
$$

your conclusion is wrong
$$
\mathcal{L}\{f((t - b)/a)\} = |a|e^{-bs}F(as) \neq \frac{ae^{-(b/a)s}}{as - 1}
$$

which should be corrected as (according to your own derivation)
$$
\mathcal{L}\{f((t - b)/a)\} = |a|e^{-bs}\frac{1}{as - 1}
$$

So is this the problem you are referring to ?


create a Calibration Board How could I know the best size of the board for my camera?
My camera is Canon D-100 It's 18 megapix
I have created a board consist of 14*10 squares
the length of each square is 7*7cm
I need to know what the step of choosing the number of squares ?
and How could I know the best dimension for each square ?
Thank you very much

You should position the board at roughly the same distance from the camera as the distance between the camera and your objects of interest. The board should be big enough to cover a good portion of the field of view at that distance.

Another point: ideally, your board should have an even number of squares along one side, and an odd number of squares along the other.  Then the orientation of the board can be detected unambiguously.

Finally, since you have a "matlab" tag on your question, please check out the Camera Calibrator app, which is a part of the latest release of the Computer Vision System Toolbox for MATLAB.


Help find the bug in my SPLL code In an attempt to get myself up to speed with PLLs I'm trying to implement Example 1 from this tutorial, in Python:
https://wirelesspi.com/phase-locked-loop-pll-in-a-software-defined-radio-sdr/
(search for "Example 1")

I have output that is close to what the article shows, except that my cos_out[n]  seems to be "leaning" to the left:

My phase-difference signale_D[n] is flipped vertically from what the tutorial shows:

I must be missing a - sign or have a sin() where I need cos() but I've checked and triple checked my code.

My code:

import numpy as np
k = 1
N = 15
K_p = 0.2667
K_i = 0.0178
K_0 = 1

samples = np.linspace(0, 100, 100, endpoint=False)
input_signal = np.cos(2*np.pi*(k/N)*samples + np.pi)

integrator_out = 0
phase_estimate = []
e_D = [] #phase-error output
e_F = [] #loop filter output
sin_out = []
cos_out = []

for n, sample_in in enumerate(input_signal):

    # phase detector
    try:
        e_D.append(sample_in * sin_out[n-1])
    except IndexError:
        e_D.append(0)

    #loop filter
    integrator_out += K_i * e_D[n]
    e_F.append(K_p * e_D[n] + integrator_out)

    #NCO
    try:
        phase_estimate.append(phase_estimate[n-1] + K_0 * e_F[n])
    except IndexError:
        phase_estimate.append(K_0 * e_F[n])

    sin_out.append(-np.sin(2*np.pi*(k/N)*n + phase_estimate[n]))
    cos_out.append(np.cos(2*np.pi*(k/N)*n + phase_estimate[n]))

I've written articles like that, and it's incredibly hard -- and sometimes not the best exposition -- to get the graphs to match exactly with the text.  From the looks of things, the author was intentionally showing you idealized graphics rather than the actual output of code.

  My phase-difference signale_D[n] is flipped vertically from what the tutorial shows...

Eh, don't worry.  It's a control loop, there's always a stray minus sign in there someplace.  As long as you match it up with another stray minus sign so that it's stable, you're OK.

  I have output that is close to what the article shows, except that my cos_out[n] seems to be "leaning" to the left: cos_out

That's because your K_p is so high.  You can get a simple loop, fast settling, or a nice sinusoidal response -- not all three.  If you're using the author's numbers, or if you've selected gains to make your loop settle as fast as his idealized example (which was probably chosen more to make things fit in the page rather than because it's the best tuning) then what's happening is that strong signal at the sum frequency is bleeding into the NCO phase, and making the "lean".

Typically a PLL will be set up to take hundreds or thousands of cycles to settle, which implies a much smaller K_p (certainly less than 0.01, and almost arbitrarily small depending on what you're doing).  How fast a PLL settles depends on what you're doing, and what you're doing it with, so if your next question is "how fast should it settle" -- I can't answer.  But I would say try it for settling in around 100 cycles, 1000, and 10,000.  Keep in mind that the capture range (the frequency error between signal and NCO initial frequency over which it can actually acquire a lock) and its tracking performance (its ability to follow a varying input frequency) both go down with slower settling.


Relationship between Nyquist and Shannon channel capacity Nyquist channel capacity(Cn) used for theoretical noiseless channel $C_n=2B\log_2(M)$

Shannon channel capacity(Cs) for noisy channel $C_s=B\log_2(1+\text{SNR})$

So is there a relationship between Nyquist and Shannon channel capacity?
is Nyquist channel capacity($C_n$) the upper limit? so $C_s \le C_n$?

Well, I think both are unrelated unless you have some constraints.

The Nyquist channel capacity says you can transmit $2B\log_2(M)$ bits per second for given channel bandwidth $B$. It does not say if these bits would be received reliably at the receiver in the presence of noise or channel distortion. I am assuming here that bandwidth limit is imposed by channel. Your transmitter can afford to have higher bandwidth but channel does not.

The Shannon Channel capacity theorem says, in order to receive the transmitted bits reliably (that is Probability of Error $\text{P}_e \rightarrow 0$) in the presence of additive noise, the maximum rate at which you can transmit is $B\log_2(1+\text{SNR})$ bits per second.

So if your requirement is reliable communication and the channel has a bandwidth limit of $B$, then $C_n \le C_s$.


STFT output frequency does not match sin wave frequency I have the following function and pre definition for producing a sin and its' stft:
fs=48000
nfft=2048
hop_len=nfft/2
win_len=2048
nfft = next_greater_power_of_2(win_len)

def sinus_gen(f):
   np.expand_dims(np.sin(2*np.pi*f*np.arange(fs*20)/fs).astype(np.float64)...
    ,axis=1)

I am trying to generate and plot an 8KHz using the following stft command:
audio=sinus_gen(8000)
s=librosa.core.stft(audio[:, 0], n_fft=nfft, hop_length=hop_len,
                                        win_length=win_len, window='hann')
librosa.display.specshow(librosa.amplitude_to_db(abs(s)),x_axis='time',y_axis='log',cmap='jet')
plt.title('Power spectogram before beamforming')
plt.colorbar(format='%+2.0f dB')
plt.clim(-60,30)
plt.tight_layout()
plt.show()

I am getting the following spectrogram image, which shows 4000KHz for some reason:

Why am I not seeing the expected 8000KHz frequency? what am I doing wrong?

You have to pass sr=fs as well as hop_length=hop_length to the specshow function, as its default values are different from yours.


$e^{j\omega}$ on unit circle I am new to DSP, and I'm self studying. I am confused about the magnitude of $e^{j\omega}$ - where $\omega$ is the normalized angular frequency - when we are on the unit circle.

According to the text book, it's $|e^{j\omega}|=1$. But why is that? I'd really appreciate it is someone could please kindly explain this to me.

Warning: $|e^{j\omega}|$ is equal to $1$ if and only if $\omega$ is a real number. More generally, for $z$ complex,  $|e^{z}|=e^{\Re{z}}$.

This really depends in your prior knowledge, because functions like exponentials, sines or cosines, can be developed in different ways. Here, I'll assume you know about sines  and cosines, as you refer to angular frequency. Then, an historical way to define the exponential is via Euler formula, for instance:

$$e^{j\omega} = \cos \omega+j\sin \omega\,.$$

Then, through Argand's complex plane interpretation, you can interpret a complex number as a point with  2D coordinates $(x,y) \leftrightarrow  x+jy$: the real part refers to the $X$-axis, the imaginary one to the $Y$-axis. Thus, obviously, $e^{j\omega}$ has the coordinates of a point on the $(0,0)$-center circle with radius $1$, hence the magnitude  $|e^{j\omega}|$ is equal to the radius, or  $(\cos \omega)^2+(\sin \omega)^2=1^2$.

In more modern versions, $e^z$ is defined by a series (see Proof of complex conjugate symmetry property of DFT), and $\cos$ and $\sin$ are derived for it, so they don't come first.


Basic Algorithm for Image Classification I have been trying to implement image classification by extracting facial features. I tried starting with some already established methods like Multi Kernel Approach, using Independent Component Analysis etc. Is there a very basic algorithm for image classification with which I can proceed the implementation ?

The data set I am currently considering is CMU Face Images data set. So aim is to separate images with human faces to those who do not contain human faces. I think I need to learn parameters by supplying  positive data set and negative data set as well. I am aware that SIFT, SURF, VJ algorithm etc. are helpful fir extracting features from an image. But I am stuck on how to start i.e. what features are required for face detection, and how to proceed towards implementing them.

Have a look at this report by Marques (PDF). In particular, have a read through section 1.4. It describes several different aspects to the face detection problem. You might want to decide what scenario you're after. Fully unconstrained (uncontrolled environment) face detection is hard to do accurately.

As a first approach, try just doing color-based segmentation.  See this paper by Singh, Chauhan, Vatsa and Singh for some details about how to use color to segment your images for face detection.


Reconstructing Tree Limbs that are Occluded by Leaves How would would one approximately reconstruct the branches in a photo of a tree, where portions of the branches are covered by leaves?

I can extract the visible portions of the branches by threshholding then opening with line segments at a few different angles.  But how would I then (approximately) fill in the spaces that were covered by leaves?

Suppose we start with this drawing of a tree:

load the image as grayscale

uu = mmreadgray('tree2.jpg')

threshold the darker colors; convert to binary; remove connected regions under 200px area

vv = mmareaopen(mmbinary(mmthreshad(mmneg(uu), 200)), 20)

open vv with 10px line segments at 90 degrees, 180 degrees, 135 degrees and 45 degrees; add these images together to form a marker image zz

zz = mmaddm(mmaddm(mmaddm(mmopen(vv, mmseline(10, 90)), mmopen(vv, mmseline(10, 180))), mmopen(vv, mmseline(10, 135))), mmopen(vv, mmseline(10, 45)))

reconstruct the visible limb portions in vv from marker image zz

aa = mminfrec(zz, vv)

display result aa on top of original image uu; a reasonable extraction of limbs is achieved

mmshow(uu, aa)

I would appreciate any thoughts on how to now approximately reconstruct the entire limb.  (Or improvements to the steps given above).  Thank you!

The masked area hidden by leaves is very large and while the border and shading seems to imply what is hidden to a human being it is unsuited to "inpainting" which relies on pretty consistent boundary, such as where a foreground object temporarily obscures something continuing behind.

The branch orientation is a start but not enough data on its own to imply how the twigs extend without having some external concept of an oak tree versus a willow, poplar or pine.  That is more about 3d fractal modeling than about image restoration.

The least obviously visible algorithms for in-painting involve texture synthesis called "Healing" in Photoshop.  G'MIC and Resynth are popular open-source tools for that available as plug-ins for GIMP.  While visual continuity is quite good they have no concept of a tree or depth, just 2-d patch statistics so would need input from an artist as well in the example given.  The mask is white background on the outer edges and thick branches on the inside .. so texture interpolation gives a weird unnatural result.


Which time-frequency coefficients does the Wavelet transform compute? The Fast Fourier Transform takes $\mathcal O(N \log N)$ operations, while the Fast Wavelet Transform takes $\mathcal O(N)$. But what, specifically, does the FWT compute?

Although they are often compared, it seems like the FFT and FWT are apples and oranges.  As I understand it, it would be more appropriate to compare the STFT (FFTs of small chunks over time) with the complex Morlet WT, since they're both time-frequency representations based on complex sinusoids (please correct me if I'm wrong). This is often shown with a diagram like this:

(Another example)

The left shows how the STFT is a bunch of FFTs stacked on top of each other as time passes (this representation is the origin of the spectrogram), while the right shows the dyadic WT, which has better time resolution at high frequencies and better frequency resolution at low frequencies (this representation is called a scalogram). In this example, $N$ for the STFT is the number of vertical columns (6), and a single $\mathcal O(N \log N)$ FFT operation calculates a single row of $N$ coefficients from $N$ samples. The total is 8 FFTs of 6 points each, or 48 samples in the time domain.

What I don't understand:

How many coefficients does a single $\mathcal O(N)$ FWT operation compute, and where are they located on the time-frequency chart above?
Which rectangles get filled in by a single computation?
If we calculate an equal-area block of time-frequency coefficients using both, do we get the same amount of data out?
Is the FWT still more efficient than the FFT?

Concrete example using PyWavelets:

In [2]: dwt([1, 0, 0, 0, 0, 0, 0, 0], 'haar')
Out[2]:
(array([ 0.70710678,  0.        ,  0.        ,  0.        ]),
 array([ 0.70710678,  0.        ,  0.        ,  0.        ]))

It creates two sets of 4 coefficients, so it's the same as the number of samples in the original signal. But what's the relationship between these 8 coefficients and the tiles in the diagram?

Update:

Actually, I was probably doing this wrong, and should be using wavedec(), which does a multi-level DWT decomposition:

In [4]: wavedec([1, 0, 0, 0, 0, 0, 0, 0], 'haar')
Out[4]:
[array([ 0.35355339]),
 array([ 0.35355339]),
 array([ 0.5,  0. ]),
 array([ 0.70710678,  0.        ,  0.        ,  0.        ])]

You are correct that the FWT is better thought of as a "cousin" of the STFT, rather than the FT.  In fact, the FWT is just a discrete sampling of the CWT (continuous wavelet transform), as the FFT/DFT is a discrete sampling of the Fourier transform.  This may seem like a subtle point, but it is relevant when choosing how you discretize the transform.

The CWT and STFT are both redundant analyses of a signal.  In other words, you have more "coefficients" (in the discrete case) than you need to fully represent a signal.  However, a Fourier transform (or say a wavelet transform using only one scale) integrate a signal from -infinity to +infinity.  This is not very useful on real world signals, so we truncate (i.e. window) the transforms to shorter lengths.  Windowing of a signal changes the transform -- you multiply by the window in time/space, so in transform space you have the convolution of the transform of the window with the transform of the signal.

In the case of the STFT, the windows are (usually) the same length (non-zero extent) at all time, and are frequency agnostic (you window a 10 Hz signal the same width as a 10 kHz signal).  So you get the rectangular grid spectrogram like you have drawn.

The CWT has this windowing built in by the fact that the wavelets get shorter (in time or space) as the scale decreases (like higher frequency).  Thus for higher frequencies, the effective window is shorter in duration, and you end up with a scaleogram that looks like what you have drawn for the FWT.

How you discretize the CWT is somewhat up to you, though I think there are minimum samplings in both shift and scale to fully represent a signal.  Typically (at least how I've used them), for lowest scale (highest frequency), you will sample at all shift locations (time/space).  As you get higher in scale (lower in frequency), you can sample less often.  The rationale is that low frequencies don't change that rapidly (think of a cymbal crash vs. a bass guitar -- the cymbal crash has very short transients, whereas the bass guitar would take longer to change).  In fact, at the shortest scale (assuming you sample at all shift locations), you have the full representation of a signal (you can reconstruct it using only the coefficients at this scale).  I'm not so sure about the rationale of sampling the scale.  I've seen this suggested as logarithmic, with (I think) closer spacing between shorter scales.  I think this is because the wavelets at longer scales have a broader Fourier transform (therefore they "pick up" more frequencies).

I admit I do not fully understand the FWT.  My hunch is that it is actually the minimum sampling in shift/scale, and is not a redundant representation.  But then I think you lose the ability to analyze (and mess with) a signal in short time without introducing unwanted artifacts.  I will read more about it and, if I learn anything useful, report back.  Hopefully others will like to comment.


WSS Ergodic Process with Power Spectrum I was given a WSS ergodic process $x(t)$ with power spectrum :

$$
\begin{array}{rcl}
       G_x(f) &=&1−\left|\frac{f}{B}\right| &\mbox{for } |f|<B\\
       G_x(f) &=& 0 & \mbox{ elsewhere}
\end{array}
$$
and first order probability density function is uniform in the range $[−A, A]$.
How do I evaluate $A$?

It's not a home work but it was a question asked in my exam which I failed.
Any help will be appreciated.

You can compute the power of the process from its power spectrum as well as from its PDF. Equating the two gives you a relation between the constants $A$ and $B$. More specifically you get

$$\int_{-B}^BG_x(f)df=\frac{1}{2A}\int_{-A}^Ax^2dx$$

If I'm not mistaken this should give $A=\sqrt{3B}$.


what does frequencies in non periodic signals mean? What do the frequencies in the a Fourier transform of a non-periodic signal mean physically?
Are there another definition of frequency that doesn't include the FT?

The physical siginificance of the Fourier transform for nonperiodic (finite energy) signals is not something clear from any point of view.
The limiting case of a Fourier transform of a periodic signal, as the period goes to infinity, consists of line components (dirac impulses at a countable set of discrete frequencies) which becomes more and more crowded, eventually approaching a continuum of frequencies, at which are no more dirac impulses.
This infinitesimal frequency component at any $\omega$, has no physical significance, unlike the dirac impulse at $\omega$, which has a direct implication of a sinusoidal signal component, extending from minus infinity to infinity (in time for example).
Non-periodic signals are also known as transients, or non-stationary signals, and for their analysis, Wavelet transform family can provide a better insight than the Fourier family.
Nevertheless, the Fourier family of transforms can still be used to analyse periodic as well as nonperiodic signals with great success.


minimum requirement of SNR for fft algorithm i will like to ask that does fft algorithm requires any minimum snr for frequency estimation

please also refer some technical documentation in this regard

Some assumptions are made in this paragraph. The signal is sinusoidal. It lands exactly at a frequency domain bin, and has magnitude $a$. No frequency domain interpolation is needed. Frequency estimation will look for and report the bin with the largest magnitude. This can be the signal bin or a noise-only bin. The noise is additive, but its contribution to the signal bin will be ignored as an approximation. The real $\operatorname{Re}(Z)$ and imaginary $\operatorname{Im}(Z)$ components of the noise $Z$ in the noise-only bins are independent and identically distributed (IID) Gaussian random variables of zero expected value.

The probability that the magnitude of the noise, at a given noise-only bin, is lower than $a$ is given by the cumulative distribution function of the Rayleigh distribution:

$$P(|Z| < a) = 1-\operatorname{exp}\left(\frac{-a^{2}}{2\sigma^2}\right),$$

where $\sigma^2 = E\left(\operatorname{Re}(Z)^2\right) = E\left(\operatorname{Im}(Z)^2\right)$, the expected value of both the squared real and imaginary components of $Z$. If there are $N$ bins to look in, then the probability that none of the noise-only bins has higher magnitude than the signal bin is given by:

$$P(|Z| < a)^{N-1}.$$

A conclusion is that the correct bin will be chosen at a higher probability if you can narrow down the frequency range where to look for the signal, making $N$ smaller. It is up to you/someone to decide how much risk of reporting the wrong bin is acceptable. There is more to say about this I'm sure.


Bandpass butterworth filter in python is not working I've implemented this bandpass filter in matlab for the audio signals and it's working as expected. Now, I tried to implement the same thing in Python but it doesn't produced expected results, it just sort of discard the whole signal.

Here is the Matlab code:

[f,fs] = audioread('abc.wav');
n=6; %order
beginFreq = (300/(fs/2));
endFreq = (2300/(fs/2));
[b,a] = butter(n,[beginFreq,endFreq],'bandpass');
%filter the signal
filtered = filter(b,a,signal);

Here is the Python code (using scipy.signal)

[originalSignal, sampleRate] = sf.read('abc.wav')
lower = (300/(sampleRate/2))
higher = (2300/(sampleRate/2))
n=6 #order
[b, a] = signal.butter(n, [lower, higher], 'bandpass')
filtered = signal.lfilter(b, a, originalSignal)
plt.plot(originalSignal, 'b', filtered,'r')

There is one more weird thing with the output! I tried to change the values of the lower and higher frequencies. As the frequency range increases; only the amplitude of output signal increases from 'almost zero' to higher value. For example, when I give the lower and higher frequency values as '100.0' and '13000' respectively, the output signal is "same" as the input signal but with almost half the amplitude. So, apparently the filter is not affecting the frequencies at all but rather affecting the amplitude.....
I've no idea what's going on... :(

Filter results:


What is state of the art in multivariate 2D image display? If I have a 2D scalar image, it's trivial to display it with MATLAB or matplotlib type tools with the pixel values mapped via a color map, and there is any amount of literature (one example, another one) describing the perils of those tools' default "rainbow" colour-maps and the advantages of various other schemes (color ramps constructed to be more suitable for the significant proportion of the population which is colorblind, for example).

However, if I have two 2D scalar images and want to fuse/overlay them in some useful way (rather than just displaying them side by side), it's less clear what the best approach is.  The naive approach of using one channel as red and one channel as green in a combined image kind of works, but is pretty ugly looking.  Mapping one channel to luminance and one to hue in an HSV color space also kind of works, but it becomes impossible to identify the hue of the dimmest pixels (unless luminance is kept quite high, but that destroys contrast).

So my question is: is there some recognized and generally considered good/best way of visualizing such fused image pairs?  Or does it all depend on the details of the data and the specifics of the "visualization problem" to be solved?

NB I'm not considering the case of visualizing images of vector quantities or complex numbers; this is more about independent but spatially registered scalar image data from different sensors/"modalities" or GIS/mapping systems where the viewer wants to get some insight into spatial correlations between the imaged quantities.

If "kind of works" but slightly less ugly is acceptable, than any two not-too-near-to-co-linear RGB color space vectors will do, and you or the user can just play with them to get a preferred combination of differentiation and non-ugliness (say, from a suitable artistic palate coordinated with other artwork color schemes in or around the data presentation, etc.) for the particular data sets.  Perhaps let the user choose one color vector, and auto rotate the other color vector to some roughly orthogonal point in the color space if it's too close to co-linear.

If you invert the data vector (or use a subtractive color vector) and thus fade to white, the lighter pastel version of some color vectors may blend better, but this depends on the mean and distribution of the data.


Periodogram: an inconsistent estimator What is the reason that a periodogram is not a consistent method for obtaining information about frequency response of a system?

To be more precise, why even after increasing number of samples doesn't the variance of the estimator decrease?

The Periodogram implicitly uses a window to get an idea of the spectrum. Implicitly this is either a rectangular window, with frequency response which is like $\dfrac{\sin(\omega n)}{N\sin(\omega)}$, or a triangular window $(N-k)/N$.

Roughly, the variance of the periodogram is the square of the Fourier transform multiplied by the square of the window. Note that neither windows decrease as N goes to infinity, so the periodogram isn't a consistent estimator.


the position of image plane in optical imaging In the book Learning OpenCV, there is one figure that show the image plane always is focal image.

The description of the figure from the book is as follow.

We begin by looking at the simplest model of a camera, the pinhole camera model. In this simple model, light is envisioned as entering from the scene or a distant object, but only a single ray enters from any particular point. In a physical pinhole camera, this point is then “projected” onto an imaging surface. As a result, the image on this image plane(also called the projective plane) is always in focus, and the size of the image relative to the distant object is given by a single parameter of the camera: its focal length. For our idealized pinhole camera, the distance from the pinhole aperture to the screen is precisely the focal length. This is shown in Figure 11-1, where f is the focal length of the camera, Z is the distance from the camera to the object, X is the length of the object, and x is the object’s image on the imaging plane. In the figure, we can see by similar triangles that –x/f = X/Z,

From this link https://en.wikipedia.org/wiki/Lens_(optics)

There is thin lens equation. Only object in infinity can be imaged in focal image. Which one is correct?

Edit:
From my previous knowledge, camera calibration is performed for cameras using fixed focal length lenses using methods from the opencv library. In calibration process, the pinhole camera model is assumed and the image plane is f away from optical center. In the book Learning OpenCV, there is one chapter which gives detailed description of calibration process.

Nowadays I learn that the same procedure is used for cameras on mobile phones.
Such cameras have the autofocus functionality. AF is implemented by changing the position of the lens by using VCM(voice coil motor). The link have one figure which shows the process. http://electroiq.com/blog/2011/11/mems-alternatives-for-miniature-auto-focus-cameras/. The focal length of the lens is fixed. If you capture the chess board each time, the VCM will modify the lens to make objects at the center of the image in focus.
Can the calibration procedure apply to such system?

I agree with user lxg that it is a notational difference, as for the pinhole camera we always have:

$$S_2 = f$$

That is not necessarily true for a camera with a lens. A pinhole camera is not the same as a camera with a lens, because a pinhole camera is sharp at all $S_1$. What makes focal length a descriptive quality of both camera types is that the magnification $M$ of a distant object imaged onto the image plane when the object is in focus can be estimated just from the focal length and the distance to the object. For a single-lens camera:

$$M = -\frac{S_2}{S_1}$$

When an object is in focus:

$$\frac{1}{S_1}+\frac{1}{S_2} = \frac{1}{f}\\
\Rightarrow M = \frac{f}{f-S_1}$$

If $f$ is very small compared to $S_1$, as it is for distant objects, it follows that:

$$M \approx -\frac{f}{S_1}$$

For the pinhole camera we have similarly:

$$M = -\frac{S_2}{S_1} = -\frac{f}{S_1}$$

The magnification is negative because the image is upside-down compared to the object.

In mobile phones the focal length is in the order of 1.5 mm to 7 mm. If you are not doing close-up, then the pinhole model approximates the magnification well, with an error under 0.15 % to 0.7 % for objects at a distance of 1 m or greater.


How to set the frequency cut-off for a low pass RC filter? How do I set the frequency for a low pass filter illustrated below?

You probably know that the frequency response of this first order low-pass filter is

$$H(\omega)=\frac{1}{1+j\omega RC}$$

Now compute its magnitude

$$|H(\omega)|=\frac{1}{\sqrt{1+\omega^2R^2C^2}}$$

In order to find the 3-dB cut-off point $\omega_0$, you have to solve $|H(\omega_0)|=1/\sqrt{2}$:

$$\frac{1}{\sqrt{1+\omega_0^2R^2C^2}}=\frac{1}{\sqrt{2}}\Rightarrow
\omega_0=\frac{1}{RC}$$

And since the frequency in Hertz is $f=\omega/2\pi$, you get

$$f_0=\frac{1}{2\pi RC}$$