Confused on application of Anti-Aliasing Filter I have a question about the use of anti-aliasing filters prior to sampling a finite time signal and if they actually help anything.

The point of the anti-aliasing filter is to remove high-frequency components to reduce aliasing. It's know a signal cannot be both band-limited and time-limited thus time limited signals are first low-passed to "reducing aliasing",

This sounds good but it seems to create its own problems in so far as it will create infinite(or near infinite) duration output as a response to a finite input.

We can't sample forever, no matter what the signals samples will be finite in quantity. So aren't we creating as big of a problem as we are solving?

Math and Physics are not the same thing. They can be very close but infinity is hard to realize.  Going from analog to digital, a more accurate description is an average of the waveform over a small time window and the average is in turn quantized. Filters don’t have infinite stop bands but 80 dB is usually more than good enough.  There is always a small error, but one can bound it.


Python equivalent code for Matlab bandpass function I have used a bandpass filter function in Matlab as follows
Fs = 128;  sampling rate
Fd = [1 4]; passband freq
[data_1_bandpassed, digital_filter] = bandpass(signal, Fd, Fs);

How can we design this bandpass filter in Python to recreate the exact same output?
When using butter & filtfilt commands in Matlab & Python, I am able to get the exact same filtered signals between the two programming languages. However, I explicilty need to recreate this output of Matlabs bandpass function in Python.
I have also attached the output of the digital filter properties and values created using Matlab bandpass function.

A similar question is already asked by another user few months ago in different domain but with no solution, so I am adding that here for reference as well!

I know it is a bit late, but i had the SAME issue as you today. But from what you showed, i may got the solution similar to Matlab. I don't know if it is the right answer but you could tell me in a future!
Just to complement what you have already said, in the matlab documentation for the function, we have: stopband attenuation of 60 dB.
I did some changes from my variables, to fit your variables.

from scipy.signal import firwin, remez, kaiser_atten, kaiser_beta, kaiserord
Fs = 128.0
lowcut = 1 #In Hertz
highcut = 5 #In Hertz
stopbbanAtt = 60  #stopband attenuation of 60 dB.
width = .5 #This sets the cutoff width in Hertz
nyq = 0.5*Fs
ntaps, gb= kaiserord(stopbbanAtt, width/nyq)
atten = kaiser_atten(ntaps, width/nyq)
beta = kaiser_beta(atten)
a = 1.0
taps = firwin(ntaps, [lowcut, highcut], nyq=nyq, pass_zero=False,window=('kaiser', beta), scale=False)
filtered_signal = signal.filtfilt(taps, a, signal)

From that, i applied FFT in both: Matlab & Python and in spectral analisys i've got the same cut frequencies with my own 1-D data from EEG .


Best algorithm for time-scaling what is the best algorithm for time_scaling of a speech signal in which I've seen two of them and I don't know which one to use :

Psola:Pitch synchronous overlap and add.

time scaling modification using modified short Fourier transform magnitude .

what is the difference between them and which one is good for speech signal ?
In other words which one in more accurate and close  to the original signal.
If there are other good ones please tell me

The best depends exclusively on the type of signal you want to modify, if you just need apply time scale in monophonic (Speech), PSOLA and TDHS can give to you very nice results, but It is not as easy as you read in books, find Pitch Marks (Glottal Closure Instant) in PSOLA is really a big problem, in TDHS if your pitch track are consistent you can scale the signal to the hell, for Time scaling unsig FFT (commonly known as phase vocoder) this works better in polyphonic signals and to get best results here do you need use a locked phase vocoder.

Find the length is a easy step, for stretch factor between 0.5 and 2.0 you can use some time domain technics: ola, sola, wsola, psola, tdhs, you can try keith lent (I consider a bit different of PSOLA) in the Lent code you can try 50% of crossfade hann window that works nice for my ears.


Filter Percentage Uniform Noise from a DC-signal? I'm not good with signal processing but i've looked around and have got no clue how to approach this.

My Question is, If there is a Static value present - that is being corrupted by percentage uniform noise, is there a way to estimate the original static value given N corrupted measurements.

i.e

For a Constant X,
On each measurement, K = X + X*N
Where N is a noise signal with a Uniform PDF between bounds +0.50 to -0.50.

How to approach processing this noisy signal to estimate the original DC-value ?

The sample mean can be a choice for the estimator with the following pros and cons

Pros:

it is unbiased

Cons:
1. The variance of the estimator increases with increase in absolute value of $x$

The estimator will converge to true value of $x$ only when number of observations over which the mean is taken is large.

A different approach is illustrated below, though estimation in uniform noise of this nature is not trivial and many a times the best that can be done is the sample mean over large data samples

Let $$y[n] = x + xw[n]$$ where w[n] is uniform between $[-0.5 \: 0.5]$. Let us denote the equivalent noise $z[n] = xw[n]$. The correlation between these equivalent noise samples accross time assuming that $w[n]$ are IID is $$E(z[n]z[n+1]^*) = E(xw[n]xw[n+1]^*)$$ since $x$ is deterministic and w[n] are IID hence $$E(z[n]z[n+1]^*) = 0$$ so $y[n]$ are uniform  RVs between $[0.5x \:   1.5x]$. The mean of each of this RVs is $x$. Therefore the sample mean is an unbiased estimator of $x$ $$ \hat{x} = \frac{\sum_n y[n]}{N}$$. The variance of the estimator is $$var(\hat{x}) = \frac{var(y[n] -x)}{N}$$
$$var(\hat{x}) = \frac{x^2}{12N}$$

However consider this $$p(y[n];x) = \frac{u(y[n] + 0.5x) - u(y[n] - 0.5x)}{x}$$

Since, the PDF above is nothing but $z[n]$ which are IID, hence combined off of $Y$ is $$p(Y;x) = \prod_{n=1}^{n=N}\,p(y[n];x)$$ the PDF will remain non zero as long as $-0.5x<y[n]<0.5x$, if this is the case the PDF is given by $$p(Y;x) = \frac{u(0.5x - max(y[n]))u(min(y[n] + 0.5x))}{x^N}$$

So $max(y[n])$ and $min(y[n])$ are joint sufficient statistics. The next step would be to find unbiased estimator of these statistics. For the case of IID uniform RVs this approach is illustrated in fundamentals of estimation theory by Stephen M Kay in the chapter on sufficient statistics to be $\frac{N+1}{2N +1}max(y[n])$, you can modify that a bit to this problem. It is shown to have less variance than the sample mean.

I have illustrated this method for the case where the sample mean doesn't make sense due to the cons mentioned above. In general estimation over uniform noise of this percentage nature is tricky.


How can I improve the resolution of a spectrogram of a low-frequency signal with high sampling frequency (Python)? Imagine you have signals $x_1(t)$ and $x_2(t)$ that contain a limited number of periods of low frequency oscillations and are sampled at 44100Hz (coming from an external source) and you want to obtain their spectrograms.

I chose nperseg=8820 to capture one period of 5Hz (44100/5=8820) in an analysis window and noverlap=1 to move in steps of one sample.

As shown in the figure below, after zooming the frequency axis, the frequency resolution in this low-frequency region is poor. Are there ways to improve this resolution? By setting special parameters of the spectrogram function or by downsampling for instance?
And why does the spectrogram of $x_2(t)$ change color over time?

import numpy

import scipy.signal
import matplotlib.pyplot as pplot

fs = 44100
t = numpy.linspace(0, 100000, 100001) / fs
x1 = numpy.sin(2*numpy.pi*5*t)
x2 = numpy.sin(2*numpy.pi*5*t**2)

f1, t1, Sxx1 = scipy.signal.spectrogram(x1, fs, nperseg=8820, noverlap=1)
f2, t2, Sxx2 = scipy.signal.spectrogram(x2, fs, nperseg=8820, noverlap=1)

pplot.figure(1)
pplot.clf()
pplot.subplot(2,2,1)
pplot.plot(t, x1)
pplot.subplot(2,2,2)
pplot.plot(t, x2)
pplot.subplot(2,2,3)
pplot.pcolormesh(t1, f1, Sxx1)
pplot.ylim([0, 30])
pplot.subplot(2,2,4)
pplot.pcolormesh(t2, f2, Sxx2)
pplot.ylim([0, 30])
pplot.show()

I'm not sure but there's something called the time-frequency resolution limit, basically the shorter your windowing interval the broader your spectrum is going to be.

The frequency resolution in the low frequency area may be poor because of this. Things like wavlets attempt to resolve this problem but i dont know anything about them.


Amplitude decrease in FFT at higher frequency for multisine I am adding sine signals with different frequencies but same amplitude=100.
Sampling frequency is  50MHz.
I am not getting 100 in frequency domain for frequency 787KHz.
At 1K, 3K I am getting 100 amplitude,
but at 787KHz I am getting 75.9 .
All the signals have 100 V amplitude.
Any reason for this?

What you see is a result of spectral leakage, it can be mitigated by using multitone frequencies that are coherent with the sampling frequency.
There is a nice tutorial about coherent sampling from Renesas.


What is the value of the standard deviation $\sigma$? I am using the following formula to add the additive Gaussian noise.
noise= (x/255)*(randn(size(image))), where x is the variable and randn is random value generator. Intensity range is $0-255$. But, what is the value of the standard deviation $\sigma$ here? Is it x or (x/255)?

randn gives you Gaussian noise with std. dev. of 1.

(x/255)*(randn(size(image))) gives you a noisy image of std. dev. (x/255)


Relationship between spectral efficiency and PAR Low PAR means less fluctiations in the magnitude of a signal and hence analog signal processing (e.g., power amplifier) is tremendously reduced. For example, constant envelope modulation schemes like BPSK, GMSK can use a switching amplifier which are >90% efficient. The downside is low spectral efficiency.

On the other side, high-order QAM or OFDM have high spectral efficiencies but the PAR becomes very high (up to 12dB for OFDM) and hence power amplifiers only reach 10-20% effiency.

Specieal cases and tricks aside (like Crest Factor Reduction, Digital Predistortion, ...) and leaving SNR out of the picture, is there a first order relationship between spectral efficiency and PAR? Any theoretical bounds?

If not, are there at least published plots that plot spectral efficiency vs PAR for various standards/modulation formats? (I did not find any)

I'm with BlackMath here – a PAPR vs SE graph makes little sense, since you would be denying the differen BER at fixed SNRs you'd get.

Let me give you an example:

256-PSK has a PAPR of 1 – the "best" you can get. It's way lower than the PAPR of 256-QAM.

Both have the same spectral efficiency (which doesn't depend on the constellation shape at all, but only on the number of constellation points).

You would, however, never use 256-PSK, though, because it's terribly low on distance between constellation points.

Since you could define a PSK with any number > 1 of constellation points, and any PSK has PAPR of 1 ("doesn't get any better"), that graph would be mighty boring.

Aside from PSK, there's, for spectral efficiency-sensitive applications, only QAM (and things that are pretty much QAM, e.g. Amplitude shift-PSK). So, that graph, if only considering constellations that "make sense but aren't PSK" would in the end only compare different power-of-2 QAMs. There's little surprise in that – you simply get the PAPR of these QAM constellations – which is a bit ugly to calculate analytically, but holds little surprise: the higher the number of bits, the higher your PAPR in a QAM.

Instead of hoping for a constellation, I'd simply use e.g. Python to give me grids of $2^n$ points centered on 0+0j, and simply calculate PAPR on these over n. That's less than 50 lines of code :)


Time synchronous averaging - MATLAB I have a time domain signal that corresponds to a vibration signal of a machine. I have a second signal that corresponds to a tachometer signal (There is a pulse every one revolution of the shaft).
I'd like to use time synchronous averaging under MATLAB. I'd like to take all blocks of time from the time domain signal whose durations is the duration of one revolution of the signal and average them out.

The problem I have is that the speed of the shaft is slightly fluctuating. I've set it to 900RPM but it goes from 895RPM to 904RPM during the acquisition and so the duration of one revolution fluctuates as well (some miliseconds).
My question is would averaging blocks of signals of different durations be a problem? And if so, how should I deal with this?

The answer to your question "is it a problem" is really up to you to determine. Sometimes the simple things work well enough.

There is a technique called Dynamic Time Warping:

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

and Matlab has a function dtw

https://www.mathworks.com/help/signal/ref/dtw.html

which can use to figure out if it is a problem.

I suggest you pick a reference interval, that would be the closest to the average shaft interval and use that to warp the rest against.


Impulse Invariant mapping Can someone help me to better understand the Impulse Invariant mapping? :

$$h[n] = T \cdot h(nT)$$

I don't understand why the discrete-time impulse response, $h[n]$, is taken from the sampled impulse response in continuous time, $h(nT)$, multiplied by the sampling period $T$.

consider an linear time-invariant system (LTI) in continuous time.  the input-output relationship is:

$$ y(t) = \int\limits_{-\infty}^{\infty} h(u) x(t-u) \, \mathrm{d}u $$

and, in the frequency domain

$$ Y(f) = H(f) X(f) $$

where

$$ X(f) \triangleq \mathscr{F} \big\{x(t)\big\} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t  $$

and similarly for $H(f)$ and $Y(f)$.

now, suppose dimensionally $x(t)$ and $y(t)$ are the same species of animal, like they are both voltage or both position or both pure numbers.  then $X(f)$ and $Y(f)$ are also the same dimension of physical quantity.  then, dimensionally, what must $H(f)$ be?

$H(f)$, which is

$$ H(f) \triangleq \mathscr{F} \big\{h(t)\big\} = \int\limits_{-\infty}^{\infty} h(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t  $$

must be dimensionless.  $H(f)$ must be a pure number with no units attached.  that is the only way it can multiply $X(f)$ and result in $Y(f)$ and the latter two have the same units or dimension.

now, in the integral above, what dimension must $h(t)$ have for the resulting integral to be the dimensionless $H(f)$?  consider the relationship between velocity and position:

$$ v(t) = \frac{\mathrm{d}r}{\mathrm{d}t} $$

or

$$ r(t) = r(t_0) + \int\limits_{t_0}^{t} v(u) \, \mathrm{d}u $$

$u$ is a dummy variable of integration and is time, just like $t$ is.
dimensionally, what does the little differential element $\mathrm{d}u$ have?  $v(u)$ is velocity, length divided by time, and it is multiplied by tiny time increments, $\mathrm{d}u$, to become tiny length increments.  those tiny length increments are added up to become the length between positions $r(t_0)$ and $r(t)$.

so in the integral $\int_{-\infty}^{\infty} h(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t$ , what must $h(t)$ be, dimensionally, so that the result of the integral is the dimensionless $H(f)$?

and what is the dimension of $h[n]$?

so this is the first step, understanding what must happen dimensionally to $h(t)$ to scale it to $h[n]$.

perhaps the best way to see why the dimensional scaler should be the sampling period, $T$ is by viewing the convolution integral as a Riemann summation:

$$\begin{align}
 y(t) &= \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \, \mathrm{d}u \\
\\
  &\approx \sum\limits_{i=-\infty}^{\infty} h(i \Delta u) \, x(t-i \Delta u) \, \Delta u \\
\end{align}$$

$\Delta u$ is a vanishingly small period of time.  If the sample rate $\frac1T$ is sufficiently high or $T$ is sufficiently small that it can take the place of $\Delta u$, then the Riemann summation is exactly the same as the discrete-convolution summation:

$$\begin{align}
 y(t) &= \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \, \mathrm{d}u \\
\\
  &\approx \sum\limits_{i=-\infty}^{\infty} h(i \Delta u) \, x(t-i \Delta u) \, \Delta u \\
\\
 \to y(t) &= \sum\limits_{i=-\infty}^{\infty} h(i T) \, x(t - iT) \, T \\
\\
 \to y(nT) &= \sum\limits_{i=-\infty}^{\infty} h(i T) \, x(nT - iT) \, T \\
\\
  &= \sum\limits_{i=-\infty}^{\infty} T h(i T) \, x\big((n - i)T \big)  \\
\\
 \to y[n] &= \sum\limits_{i=-\infty}^{\infty} \Big( T h(i T) \Big) \, x[n-i]  \\
\\
  &= \sum\limits_{i=-\infty}^{\infty} h[i] \, x[n-i]  \\
\end{align}$$

there you can see (changing the dummy variable $i$ to $n$) that

$$h[n] = T \, h(nT)\ .$$


Interpolation filter gain If somebody could clarify why everywhere is written that the interpolation filter should have gain $L$ where $L$ is upsampling factor.

What this gain actually means?
Does it compensates for the amplitude scaling inherent to upsampling?
If I would hypothetically have an "interpolation" filter which is not attenuating spectral replicas (let's call it unity filter) would this filter still have gain $L$, or in other words is gain $L$ something what is by the definition within the interpolation filter and then the filter coefficients which define how well the replicas will be attenuated is something which depends on the design/quality of the filter?

What this gain actually means?

The factor between in- and output amplitudes for a signal in the passband

  Does it compensates for the amplitude scaling inherent to upsampling?

Exactly!

  If I would hypothetically have an "interpolation" filter which is not attenuating spectral replicas (let's call it unity filter) would this filter still have gain L, or in other words is gain L something what is by the definition within the interpolation filter and then the filter coefficients which define how well the replicas will be attenuated is something which depends on the design/quality of the filter?

You design a filter to have a certain gain; setting gain is easy, just multiply your filter with a scalar.


Mathematics / Signal theory behind billiard ball 'wave pendulum' effect This YouTube video shows a very interesting effect.

What is the underlying science?

I have recently started studying Fourier theory and DSP, and and trying to understand what is going on in terms of the material I've been learning.  I would like to understand this from as many angles as possible.

For example, at 0:39 the signal forms two lines moving against one another,  and shortly before that I can clearly see two sine waves in the signal.

This reminds me of the illustrations at http://www.dspguide.com/ch3/2.htm

PS  not quite sure how to tag this

As was pointed out in a comment above, this is a simple consequence of the dynamics of a pendulum. There's nothing particularly signal-processing-related about this problem, just some simple physics and trigonometry.

For a pendulum that is displaced from its angular equilibrium point by an angle of $\theta_0$ with an angular velocity of $\omega_0 = 0$ at $t = 0$, for small initial angles, the motion of the pendulum can be approximated by simple harmonic motion:

$$
\theta(t) = \theta_0 \cos\left(\sqrt{\frac{g}{l}}t\right)
$$

where $g$ is the magnitude of the force of gravity and $l$ is the length of the pendulum. In the video, a number of pendulums of varying lengths are given an approximately equal initial displacement $\theta_0$ and released simultaneously. Once they are released, each pendulum traces out a path whose displacement angle is a sinusoid with frequency related to the pendulum's length. The lengths were carefully chosen to give the interesting visual effect that you see in the video; the frequencies of the sinusoids are multiples of some fundamental frequency, so for particular values of $t$, they line up in interesting ways.

You'll notice toward the end, all of the balls "line up" again; this is because each of the sinusoids that describe their motion have moved through an integer number of periods, so the balls are back at a common angular displacement. If there were no effects of friction, air resistance, etc., the oscillation would continue, with the same patterns repeating indefinitely.


I have a noisy ECG Signal. want to use PSD to analyse the Signal to noise ratio I have a really noisy ECG Signal which contain, motion artifacts,muscle noise, a very high amplitude noise.

I am beginner to DSP. Can someone tell me in a simple language,how do i interpret the signals feature when i get this kind of image.
samping frequency is 1000hz. i used psd.welch in matlab.

Why it is plotting against only 500hz and what information i can get on this

Why it is plotting against only 500hz and what information I can get on this??

Plotting against $500Hz$: Because that is the frequency of the sampled signal. We usually sample a signal at least two times of its max frequency in order to avoid aliasing. In your case $fs=1000Hz$ (sampling frequency), so max frequency is (no higher than) $500Hz$.
The information you can get: Which frequency corresponds to high powers. You can avoid power below $0$ frequencies. High frequency components are usually our signal while low power is noise.


Control design: under what conditions can closed-loop poles be placed arbitrarily? Say we have a single-input linear system $\dot{\mathbf{x}} = A\mathbf{x}+Bu$. With full-state feedback ($u=-G\mathbf{x}$), it is straightforward to arbitrarily place the $n$ closed-loop poles (i.e., the eigenvalues of $A-BG$, or the roots of $sI-(A-BG))^{-1}=0$) by selecting the values of $G$.

I'll give a brief example:  $\dot{\mathbf{x}}= \begin{bmatrix}0&1\\-a&-b\end{bmatrix}\mathbf{x} + \begin{bmatrix}0\\1\end{bmatrix}u$.

Say we want the poles of the closed-loop system to be at $-1$, i.e., we want the characteristic equation to be $s^2 + 2s + 1 = 0$. We can achieve this by choosing $u = \begin{bmatrix}a+1&b+2\end{bmatrix}\mathbf{x}$.

As long as the system is controllable, the plant dynamics don't matter -- with full state feedback, we can always cancel them and place the closed-loop poles arbitrarily.

However, I don't fully understand why we aren't able to arbitrarily place the closed-loop poles in the classical transfer function/block diagram approach. Consider a unity feedback system with control transfer function $G_c(s)$ and plant transfer function $G_p(s)$.

The closed-loop transfer function from reference to output is

$$\frac{Y(s)}{R(s)} = T(s) = \frac{G_c(s)G_p(s)}{1+G_c(s)G_p(s)}$$

It seems like I can always pick a $G_c(s)$ to give me the desired closed-loop transfer function $T_d(s)$, i.e., by rearranging the above equation to solve for $G_c$:

$$G_c(s) \stackrel{?}{=} \frac{T_d(s)}{1-T_d(s)}\cdot\frac{1}{G_p(s)}$$

This seems like the Laplace domain analog of the state-space approach I described above. So why don't we do this? I can imagine several reasons this might go wrong:

Perhaps this approach always gives the desired response from reference to output, but gives undesirable disturbance or noise responses due to pole-zero cancellations between the controller and plant.
Perhaps the resulting control transfer function can't be implemented physically? Say the resulting $G_c(s)$ is improper - does that mean it is "unrealizable"? Would this imply that it can't be implemented?

Are there other reasons we don't use this approach?

--

Update. I've applied the proposed transfer function approach to the state-space system I described above. The plant transfer function is
$$G_p(s)= \frac{1}{s^2+bs+a}$$
My desired closed-loop transfer function is $T(s) = 1/(s+1)^2$ to place the poles at $-1$. So my required controller transfer function would be

$$G_c(s) = \frac{T_d(s)}{1-T_d(s)}\cdot\frac{1}{G_p(s)} = \frac{s^2+bs+a}{s(s+2)}$$

so I'm effectively just using the controller to cancel the plant poles in $GcGp(s)$. This makes me think that perhaps disturbance response or noise rejection would be a limitation with this design. However, I am still missing an insight into why this is a problem using the transfer function approach but not the state-space approach.

edit: added bounty.

Unnecessary additional information that might help

From a pragmatic point of view, there might be a problem when trying to control a given plant using only its transfer function $G_p(s)$. In general, you test the plant putting a series of inputs and seeing how it reacts, measuring the respective outputs. From that, you can model the plant as a transfer function in the Laplace domain. However, the intrinsic dynamics of the plant may have hidden unstable modes (for instance, poles that cancel with zeros, which means you are not able to notice their presence in the system) and thus the control system you design wouldn't work as expected.

But that's in real life. In Control Theory, we are often offered the transfer function of the plant with all the poles and all the zeros. So let's assume that $G_p(s)$ doesn't show this weird behaviour I described above. In that case, one must notice (at least) three advantages that state spaces (SS from now on) offer:

Working with SS instead of transfer functions can be way easier. If the order of the system is too high, operating with matrices can be the only way to deal with the system.
SS provide an easy way to deal with MIMO systems. Using transfer functions for each pair of input-output would be too tedious.
SS let us see those hidden unstable modes I referred to before (looking at the eigenvalues of the matrix $A$), while in the transfer function they cancel out.

The actual answer to your question

Now getting to your question, you've asked what happens when the system is of low order (in your example, the plant is of second order) and SISO.

You can definitely find your transfer function performing the variable isolation you've written for $G_c(s)$. The problem with that is that some things could go wrong. For example,

$G_c(s)$ could be impossible to implement.
$G_c(s)$ could be unstable, thus making the control system non-sense.
$G_c(s)$ could have zeros in the RHP that might cancel out with poles of $G_p(s)$.

If none of those things happen (maybe there are some other aspects that one should take into account, but I couldn't think of any additional one), then go ahead: you are free to find your control system, $G_c(s)$, without using full-state feedback at all.

The thing is that SS offer us a lot more information than transfer functions do. This is easy to notice, as when working with SS we are looking at not only the input and output variables, but the other state variables as well, which are totally left apart when working with transfer functions and block diagrams alone. When doing full-state feedback, we are using this additional information to control the system (i.e. you don't just use the output, but each one of the state variables). This means that there are other tools we can use to stabilize the plant and we were not using before. If the design using full-state feedback is good, then the closed-loop will work as a charm and none of the problems I've described above (unstability, unstable poles and zeros cancelling out, etc.) will appear.


How can I call the area under an impulse/spike? I am looking for an appropriate term to express the area under a spike or impulse for time series signals. The term "Briefness" is probably not appropriate since the amplitude should be taken into account

What you're describing "literally" is just the area under a curve, the integral. I don't think this helps you.
What you want, functionally, is not quite clear, but I'd bet it's either

the energy of the pulse or
the energy-bandwidth product.

The energy of the signal is the integral over the (magnitude)
squared signal in the continuous time case, or the sum over the
(magnitude) squared signal in continuous time.
The bandwidth is a frequency domain concept, and will tend to be larger the "shorter and sharper" the pulse is. The product of bandwidth and energy has important estimation quality effects for example in radar applications.


Philosophy of perfect inter-sample interpolation How would you define perfect (inter-sample) interpolation, and is it possible?

To quote Armifinn's prior answer:

  "I guess the most important result is that for signals with bandwidth limitation, you can have perfect reconstruction via sinc(⋅)sinc(⋅) convolution; the famous sampling theorem..."

If by perfect interpolation it is meant that an analog band-limited signal is recovered perfectly from digital samples, wouldn't a sinc interpolation lead to some problems? Considering that analog equipment is causal while sinc interpolation is anti-causal, shouldn't a minimum phase low pass filter be used instead?

On the other hand, I understand that the above method would not preserve the phase information of the system. Sampling theorem states that band-limited signals can be captured perfectly. So indeed that would be a contradiction.

Finally, a quick thought experiment: consider an analog dirac delta impulse that is low pass filtered using a perfect causal analog filter so that it satisfies the sampling criterion. The signal is then digitally captured by a DSP engineer. Now, the DSP engineer wants to reconstruct the analog signal in his lab, and after running the perfect sinc interpolation he discovers that the signal in fact started BEFORE the analog dirac delta even occurred, to be more precise an infinite time before the analog signal started. This casts some doubts on how the interpolation should in fact be done, and whether the band-limiting criterion is so easy to satisfy.

By the fundamental notions of signal processing, we can state that, a perfect intersample interpolation of a given perfectly bandlimited signal from its samples is not possible due to the fact that a perfectly bandlimited signal would be infinetely long whose infinite many samples are required to compute any one single interpolated intersample, hence impossible.

That being said, i.e. that the impossibility of obtaining infinite duration samples, a practical answer would be yes if you would allow an upperbound on the acceptable interpolation error, such as lower than the roundoff error within a digital system. This will be a truncation error whose spectral manisfestation is spectral smearing due to the convolution of the truncating window's frequency spectrum with that of the true spectrum of the bandlimited signal (i.e. multiplication in time, convolution in frequency)

Also a tradeoff can be sought that, instead of a perfectly bandlimited signal, you can work with bandlimited enough signal whose aliasing error is within acceptable bounds, just as in the previous paragraph, in this case you wont need infineteley many smples, but there will be spectral overlap.


Is there any simple algorithm to convert a sound source in 3d to volume and panning I'm looking for a simple algorithm which calculates an approximation of volume and pan at the position of the camera from the following 3d data:

Camera (position and look direction)
An item which emits mono-audio at a specific position.

Just a basic approximation is enough for this case.

inverse-square law for distance.  and Blumlein stereo for pan.


matlab impulse response i have a matlab file and i'm trying to understand what's going on.

when i run it, it shows this when i type the command :

this at

and this at

can someone explain to me what increasing the value in  does because i can't seem to understand some of what's going on in the graphs.

The two top figures show the response of an ideal low pass filter and (the center part of) its corresponding infinitely long impulse response. The bottom figures show the result of windowing the impulse response, i.e., cutting out the $N$ center samples and making the rest zero, which gives you a finite length impulse response that can be implemented. You can see that windowing obviously changes the frequency response, but the longer the window is made the better the approximation. However, note that the maximum approximation error (around the discontinuity) always has the same value, regardless of the window length. This is called the Gibbs phenomenon.


An effective algorithm for convolution of very many long signals I need to calculate the following linear convolution

$$ y[n] = h_1[n] \circledast h_2[n] \circledast h_3[n] \circledast \cdots \circledast h_{k-1}[n] \circledast h_k[n]$$

where $k$ exceeds $5000$ and the length of each $h_i[n]$ exceeds $100,000$.  I tried the method which is based on the circular convolution/DFT but it does not work because each $h_i[n]$ must be zero padded so DFTs then takes too much time and memory. Are there any other algorithm available?

So, I'm becoming a bit notorious with the whole "by using commonly available tools, how complex is your problem really?" business, but oh well:

I went ahead and simulated the "last" of your convolutions, i.e. I went ahead and picked up a FFT FIR (which is convolution by the overlap-save method Jason R mentioned) by dumping 500 Million random samples through a filter of 100,000 taps:

So, this flow graph produces an average total of 42 MS/s, or 420 (500x100,000) convolutions per second. The complete program, which includes initial FFT of the taps, calculating of the taps and the random pool, loading of libraries and set up of all, takes about 12 seconds. So, yes, doing this exact operation 5000 times would exceed the average coffee break's duration, and take nearly 17h.

I tried timing the initial FFT of the tap vector (which the FFT filter blocks automatically), but I failed - it took much shorter than producing the taps vector itself. So I made a separate flow graph that just does 100,000-transforms of input data, and includes all the data generation/copying-overhead in my speed estimate. With that flow graph, my PC does 1900 FFTs per second, so it should take roughly half a millisecond to do a single one, or 2.5 s to do 5,000 of those.

Which also means that if you're smart, and just used the overlap-add algorithm as contained in the FFT FIR block and hoped that performance scales roughly linear with length of the input vector you're transforming, then every convolution would take about $\frac{12}5\text{ ms}= 2.4\text{ ms}$ longer than the previous, giving us a rough estimate of $\sum\limits_{n=0}^{5000} n\cdot 2.4\text{ ms}\approx 8\text{ h},\,25\text{ min}$.

Obviously, that leaves a bit of headroom for optimization. Now, Operlap-Add and Overlap-Save should actually scale exceptionally well to multi-core/distributed solutions. Also, doing a set of long FFTs also sounds like something that GPUs are great at.


How to synchronize audio and video streaming for not fixed FPS streaming? It's a 3-stage process.

In the first stage, a video (with audio) is recorded at a fixed 25 FPS rate. In the second stage, it's demuxed in two different channels and the video is processed. In the third stage, the audio and video are again synchronized and muxed together.

The issue is for the third stage, the video after the processing is having varying FPS rate (like 16, 18, 30, 35 FPS). Is there any way to synchronize the audio with the video at original FPS rate. I'm trying to use FFMPEG for this but no success.

Container formats are usually storing video and sound information in chunks that are supposed to be triggered at the same time. In turn, each chunk might further be compressed according to the specifications of a given CODEC.

For example, consider a video at $FPS=30$ with its sound sampled at $Fs = 48kHz$.

In this example, one frame of the video sequence (one video "sample") should remain on the screen for $\frac{1}{30}$ seconds. But, $\frac{1}{30}$ seconds is equivalent to $\frac{48000}{30}=1600$ audio samples. And the idea here is that you put the frame on the screen and then trigger the sound and by the end of that sound clip you put the next frame up and trigger the next clip.

Or, to be more realistic, you can take the two streams and establish two functions that map audio samples to video samples and the reverse.

For instance, which video frame would be shown to the viewer at sound sample $n_{audio}=987600$?

That would be:

$$n_{frame} = \left \lfloor \frac{n_{audio}}{Fs} \cdot FPS \right \rfloor$$

Where $\left \lfloor \cdot \right \rfloor$ denotes the floor function.

And in this case:

$$n_{frame} = \left \lfloor \frac{987600}{48000} \cdot 30 \right \rfloor = 617$$

If either of your video frame rate or audio sampling frequency are variable, then the respective parts of the above equations become cumulative sums. So, to find where in time are you within either sequence, you have to accumulate the rates that brought you to that time (first 20 frames at 12FPS, next 40 frames at 20 FPS, ... and so on and similarly for the audio). CODECs will have that information in timestamps or even better, the decoder would be able to return "raw" streams and then you are back to the simple equations from above.

But in this case, unless you operate the CODEC yourself (even an uncompressed CODEC, just to specify which sound chunks match which video frames), it will be very difficult to "pipe" two variable streams to a third party software and let it figure out how to put things together.

The simplistic solution here would be to assume that your video runs at the highest FPS and handle other rates by repeating (or dropping if required) frames. The codec at the end of this will assume that the frames and sound it is passed are at a constant rate and "chop them up" accordingly. But I don't know how suitable this might be for your application (?)

Hope this helps.


DFT and Inverse DFT in Image Processing I have implemented DFT from this link.

1) Tried taking dft for the three Channels(R, G and B) and reconstructed the original image by taking inverse dft for all three channels and merged them together, but the image reconstructed doesnot have the same contrast as the original image.

But while altering the flags used in the link as such

  dft(input, Complex, DFT_SCALE);
  dft(Complex, InverseDFTImage, DFT_INVERSE + DFT_REALOUTPUT);
  InverseDFTImage.converTo(InverseDFTImage, CV_8U);

The output is same as the input, but while applying Gaussian filter, there are few colour differences near the edges.

 1. Original Input Image

 2. Output without any filters(followed the Procedure in the link)

 3. After applying Gaussian Low Pass Filter

 Can someone clarify this doubt.?
 P.S Thanks in advance

I saw same issue, and I found the answer here.
for short, you can set imaginary part of filtering kernel zeros before you call mulSpectrums.

Mat planes[] = {Mat::zeros(dftImage1.size(), CV_32F), Mat::zeros(dftImage1.size(), CV_32F)};
Mat kernel_spec;
planes[0] = your_filter; // real
// planes[1] = your_filter;
// imaginar should be zero! Don't touch it.
merge(planes, 2, kernel_spec);
mulSpectrums(...

Hope this helps.


PyWavelets CWT implementation I seek to understand PyWavelets' implementation of the Continuous Wavelet Transform, and how it compares to the more 'basic' version I've coded and provided here. In particular:

How is integrated wavelet, then differenced post-convolving, equivalent to convolving without int+diff?
How is resampling the integrated wavelet at increasing resolution equivalent to changing its scale?

Naive implementation:
def cwt(x, win_len=None, win='morlet', real=False):
    N = len(x)
    win_len = win_len or N // 8

    scales = _scales(N, nv=32)
    coef = np.zeros((len(scales), N), dtype='complex128')
    kernel = morlet_kernel(win_len)
    wl2 = win_len // 2

    for i, scale in enumerate(scales):
        coef[i, :] = np.convolve(x, kernel(scale)[::-1])[wl2:-(wl2 - 1)]
    return coef

PyWavelets (clipped to relevant parts):
def pywt_cwt(data, scales):
    out = np.empty((np.size(scales),) + data.shape)
    int_psi, x = integrate_wavelet(wavelet='morl', precision=10)

    for i, scale in enumerate(scales):
        step = x[1] - x[0]
        j = np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step)
        j = j.astype(int)  # floor
        if j[-1] >= int_psi.size:
            j = np.extract(j < int_psi.size, j)
        int_psi_scale = int_psi[j][::-1]

        conv = np.convolve(data, int_psi_scale)
        coef = - np.sqrt(scale) * np.diff(conv, axis=-1)

        d = (coef.shape[-1] - data.shape[-1]) / 2.
        coef = coef[..., floor(d):-ceil(d)]
        out[i, ...] = coef
    return out

Comparison on $f=1, 4$ sinusoids:

First see "Naive Breakdown" in a below section. Onto PyWavelets: the algorithm was found on Github to stem from an old MATLAB implementation, but it provides no details on coding the wavelet kernel. So I set on to explore; various visuals and code annotations here. Key findings follow. -- All code
See "PyWavelets Breakdown" below. It is implementing the following - namely Eq 4:

(1) How is integrated wavelet, then differenced post-convolving, equivalent to convolving without int+diff?
The formulation is in the right, deriving from the left, equations; Instead of continuously integrating over all $\mathcal{R}$, the integral is split into segments ssumming over $k$. The input is assumed constant over $[k, k+1]$ (though shouldn't it be $)$?), as it's discrete, whereas the wavelet is kept continuous. This enables to move the input outside the integral; visually:

So how do we code this? It's clear from Eq 4 that we need a diff somewhere. pywt does np.diff(conv), which is basically np.diff(coef); why does it seem to take difference of the would-be coefficient in Naive?
In Eq 4 we have the integrated wavelet up-to k + 1, minus up-to k, multiplied with the wavelet. This is equivalent to product with the wavelet at b - 1 and b, respectively (plug in t=k+1 and t=k, compare with plugging in b=b-1 and b=b instead`).
The answer's then in in conv; assume it was already trimmed so that it's only when the wavelet is at least halfway inside the signal. Then, every point in conv is the product of the wavelet with a portion of the input, at a different tau (or rather b). So two adjacent points in conv are products with wavelets at b and b + 1:
$$
\begin{align}
  \text{conv} &= [s \cdot \psi_{\text{int}}(t - 0),\ s \cdot \psi_{\text{int}}(t - 1),\ ...] \\
  \Rightarrow  \text{DIFF}(\text{conv}) &= [s \cdot (\psi_{\text{int}}(t - 1) - \psi_{\text{int}}(t - 0)),\ ...]
\end{align}
$$
However, we seek minus between t=k+1 and t, or a greater t and lower t, which is the negative of above, thus explaining step 7 in PW Breakdown.

(2) How is resampling the integrated wavelet at increasing resolution equivalent to changing its scale?
See figure under (5) in PW Breakdown. pywt appears to "scale" via the number of samples that define the wavelet, rather than recomputing the wavelet with a dilated argument. ... though, visually, this resampling clearly does dilate the wavelet, stretching it relative to the input.
So it's indeed a form of scaling; but how does it compare to naive? Naive is recomputed more intuitively per equations, using fixed wavelet length, and passing a scaled argument. Then, which is better? The opposite question is clearer to answer; consider what happens at higher and lower scales:

Naive higher: wavelet dilates too much, exiting the sampling frame far before its zero tails.
pywt higher: wavelet is just fine; the problem's with conv. Wavelet length is fixed at 1024, so if the input is any shorter, then higher scale wavelets can never fully multiply the signal. The greater the disparity, the more the wavelet is "seen" similar to "Naive higher" by the signal; this can be seen in the question's heatmaps differing by vertical shifts.

Maximum scale of the wavelet is also unlimited; len(j) grows proportionally with scale (but there's a normalization problem, see here).

Naive lower: all's good.
pywt lower: wavelet has small resolution.
Fixed vs variable window length: higher scales must correspond to greater dilation, or being-nonzero, of the wavelet relative to the input. Both accomplish this, except Naive discontinues to zero abruptly at higher scales.

The higher scale problems in both are solved by limiting the maximum scale; if both are in their completely safe zones, the Naive implementation seems preferable due to higher resolution. However, this "safe zone" for Naive is prohibitively small, if we are to consider "safe" as not jumping to zero at tails.
To contrary, pywt's problem is more easily solved by increasing maximum wavelet length, which was suggested to rid of the zipper-like artifacts (which can be seen in the question's plots). That'd, however, worsen higher-scale performance for short inputs - but in practice it's unlikely to be shorter than 1024.

PyWavelets Breakdown:

Wavelet, prior to integration, matches exactly with the shown code blob, which is an approximation of the complete real Morlet (used by Naive) assuming $\sigma > 5$ in the Wiki.
pywt integrates real Morlet via np.cumsum(psi) * step, accounting for the differential step size
The integrated wavelet, int_psi, is reused for all scales
For each scale, the same int_psi is resampled at increasing resolution, with j as an integer varying linearly between same min and max (almost) for all scales:

The resulting wavelet for increasing scale, side-by-side with naive (for later):

The convolution result first diff'd,
negated,
multiplied by sqrt(scale),
then trimmed on both sides to either (a) only include parts where the wavelet is at least halfway "inside the signal", or (b) force its length to len(x).

Naive Breakdown: if unfamiliar with CWT/STFT, I strongly recommend parts 1-3 here.
We need (1) wavelet type; (2) wavelet length; (3) wavelet increment/overlap. We'll use real Morlet, 96 samples, and increment by 20.
win_len = 96  # "win" == window == wavelet
win_inc = 20
n_wins = (len(x) - win_len) // win_inc + 1

coef = np.zeros((n_wins, len(scales)), dtype='complex128')
kernel = morlet_kernel(win_len)

for tau in range(n_wins):
    start = tau * win_inc
    end   = start + win_len
    coef[tau, :] = _transform(x[start:end], kernel, scales)

def _transform(x, kernel, scales):
    coef = np.zeros(len(scales), dtype='complex128')
    for i, scale in enumerate(scales):
        psi = np.conj(kernel(scale))
        coef[i] = np.sum(x * psi / np.sqrt(scale))
    return coef

Visualized for scale=2:

The difference between this and naive implementation in the question is, latter uses win_inc=1, equivalent to convolving with a flipped wavelet and trimming the result to only include parts where the wavelet was fully "inside the signal". Also, here we have tau ("increment") in the outer loop, there scale, but the two are equivalent.
Lastly, shapes along timeshift dim differ; question's implementation trims so that only parts of convolution where wavelet is at least halfway "inside the input" are included; this one trims further to where the full wavelet is included.

UNSOLVED:

Why * sqrt(scale) (step 8 pywt)? Eq 4 shows 1 / sqrt(scale). I considered * scale per chain rule, except we've integrated the wavelet numerically... Makes sense to have per-scale treatment; do we get away with not integrating a per-scale wavelet because the per-scale wavelet is... the same resampled wavelet?  Resolved.


Hilbert transform of unit step function How to calculate Hilbert transform, if it exists, of the signals like $u(t)$, $sgn(t)$.

What properties should a function satisfy for existence of Hilbert transform.

Absolute integrability of a signal $x(t)$ is necessary condition for HT to exist?

Absolute integrability is not a necessary condition for the Hilbert transform to exist. Simple counterexamples are the Hilbert transforms of $\cos(\omega t)$ and $\sin(\omega t)$, given by $\sin(\omega t)$ and $-\cos(\omega t)$, respectively, if $\omega>0$.

Since the functions $u(t)$ and $\textrm{sgn}(t)$ are related by scaling and an additive constant, and since the Hilbert transform of a constant is zero, the Hilbert transforms of the two functions - if they exist - must be scaled versions of each other:

$$\textrm{sgn}(t)=2u(t)-1\\\Rightarrow \mathcal{H}\{\textrm{sgn}(t)\}=2\mathcal{H}\{u(t)\}\tag{1}$$

I believe that the Hilbert transform of these two functions can be computed in the following way. [Edit: this turned out to be true. See below and the comments to RBJ's answer.] We know that the Hilbert transform can be written as a convolution:

$$\mathcal{H}\{x(t)\}=\frac{1}{\pi t}\star x(t)\tag{2}$$

This can be rewritten as

$$\mathcal{H}\{x(t)\}=\frac{1}{\pi}\left(\log(|t|)\right)'\star x(t)=\frac{1}{\pi}\log(|t|)\star x'(t)\tag{3}$$

where $'$ denotes differentiation w.r.t. $t$. The rightmost equality in $(3)$ is true because we generally have

$$(f\star g)'=f'\star g=f\star g'\tag{4}$$

With $x(t)=u(t)$ we get $x'(t)=\delta(t)$, and from $(3)$ we obtain

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi}\log(|t|)\star \delta(t)=\frac{1}{\pi}\log(|t|)\tag{5}$$

From $(1)$ and $(5)$ we get

$$\mathcal{H}\{\textrm{sgn}(t)\}=\frac{2}{\pi}\log(|t|)\tag{6}$$

The results given by Eqs $(5)$ and $(6)$ can also be derived in the frequency domain. The Fourier transform of the convolution

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi t}\star u(t)$$

is given by

$$\mathcal{F}\left\{\frac{1}{\pi t}\right\}\mathcal{F}\{u(t)\}=-j\;\textrm{sgn}\{\omega\}\left(\pi\delta(\omega)+\frac{1}{j\omega}\right)=-\frac{\textrm{sgn}(\omega)}{\omega}=-\frac{1}{|\omega|}$$

So we have

$$\mathcal{H}\{u(t)\}=\mathcal{F}^{-1}\left\{-\frac{1}{|\omega|}\right\}$$

From this document (p. 27, Ex. 19.4), we know that

$$\mathcal{F}^{-1}\left\{-\frac{1}{|\omega|}\right\}=\frac{1}{\pi}\log(|t|)$$

confirming Eq. $(5)$.

(Note the different scaling in the cited document due to the use of the unitary Fourier transform.)

EDIT:

The problem can also be solved using the definition of the Hilbert transform of bounded functions:

$$\mathcal{H}\{f(t)\}=\frac{1}{\pi}\left[\mathrm{p.v.}\int_{|\tau|<1}\frac{f(\tau)}{t-\tau}d\tau+\mathrm{p.v.}\int_{|\tau|>1}f(\tau)\left(\frac{1}{t-\tau}+\frac{1}{\tau}\right)d\tau\right]$$

where $\mathrm{p.v.}$ denotes the Cauchy principal value.

For $f(t)=u(t)$ we have

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi}\left[\underbrace{\mathrm{p.v.}\int_{0}^{1}\frac{d\tau}{t-\tau}}_{\displaystyle I_1}+\underbrace{\mathrm{p.v.}\int_{1}^{\infty}\left(\frac{1}{t-\tau}+\frac{1}{\tau}\right)d\tau}_{\displaystyle I_2}\right]$$

If $t$ is outside the interval $(0,1)$, the first integral is simply

$$I_1=\int_{0}^{1}\frac{d\tau}{t-\tau}=\log\left(\frac{t}{t-1}\right)$$

If $t\in (0,1)$ we get

$$I_1=\lim_{\epsilon\to 0^+}\int_0^{t-\epsilon}\frac{d\tau}{t-\tau}+\int_{t+\epsilon}^1\frac{d\tau}{t-\tau}=\log\left(\frac{t}{\epsilon}\right)+\log\left(\frac{-\epsilon}{t-1}\right)=\log\left(\frac{t}{1-t}\right)$$

Since in both cases the argument of the logarithm is positive, both results can be summarized as

$$I_1=\log\left|\frac{t}{1-t}\right|$$

If $t<1$ the second integral is

$$I_2=\lim_{\epsilon\to 0^+}\left[\int_1^{1/\epsilon}\frac{d\tau}{t-\tau}+\int_1^{1/\epsilon}\frac{d\tau}{\tau}\right]=\lim_{\epsilon\to 0^+}\left[\log\left(\frac{t-1}{t-1/\epsilon}\right)+\log\left(\frac{1}{\epsilon}\right)\right]=\\=\lim_{\epsilon\to 0^+}\log\left(\frac{t-1}{\epsilon t-1}\right)=\log\left(1-t\right)$$

If $t>1$ we have

$$I_2=\lim_{\epsilon\to 0^+}\left[\int_1^{t-\epsilon}\frac{d\tau}{t-\tau}+\int_{t+\epsilon}^{1/\epsilon}\frac{d\tau}{t-\tau}+\int_1^{1/\epsilon}\frac{d\tau}{\tau}\right]=\\=\lim_{\epsilon\to 0^+}\left[\log\left(\frac{t-1}{\epsilon}\right)+\log\left(\frac{-\epsilon}{t-1/\epsilon}\right)+\log\left(\frac{1}{\epsilon}\right)\right]=\\=\lim_{\epsilon\to 0^+}\log\left(\frac{1-t}{\epsilon t-1}\right)=\log(t-1)$$

Again, both results can be summarized as

$$I_2=\log|1-t|$$

Combining these results, we obtain the Hilbert transform of $u(t)$:

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi}(I_1+I_2)=\frac{1}{\pi}\left(\log\left|\frac{t}{1-t}\right|+\log|1-t|\right)=\frac{1}{\pi}\log|t|$$

which agrees with our previous result given by Eq. $(5)$.


Determining Causality From Discrete Impulse Response Trying to wrap my mind around the concepts of this one...

  Consider the following impulse response $h[n]$ for a linear, time-invariant system:
  $$
h[n] =\left\{\underline{1} , -2, 2, -1\right\}
$$
  where n=0 starts at 1 for h[n]

  (a) Is the system causal? Why or why not?
  (b) Is the system stable? Why or why not?
  (c) Is the system FIR or IIR? Give a reason for your answer.
  (d) Does the system have memory? why or why not?

Based on my understanding I've concluded...

(a):  ?
(b): Yes the system is stable because the summation of |h[n]| is less than infinity.
(c): The system is an FIR (Finite Impulse System) because the duration of h[x] is not infinite. It has a finite number of values.
(d): This system is not memoryless, therefore has memory. This is because in order to be memoryless (to my understanding) an impulse response should not have values at any times other than n=0. h[n] has multiple values

Updated the question to reflect how the values relate to time. The underlined 1 in h[n] denotes the origin n=0

As soon as you've figured out how the given values for $h[n]$ relate to the time index $n$ (which is unclear in the current problem formulation), it should be easy to decide whether the given system is causal or not by considering the input-output relation of a discrete-time LTI system:

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

where $x[n]$ is the input sequence, and $y[n]$ is the output sequence.

For the system to be causal, the output at time $n$ must only depend on the current value or on past values of the input $x[n]$. What does that imply for the impulse response $h[n]$?


MFCC extraction last step DCT or IDCT I'm confused with the last step of extraction MFCC where we have to take DCT or IDCT.

Cepstrum is defined as IDFT of log(power spectrum). Since MFCC are cepstral coefficients so I think we have to take IDCT at the last step.

It's confirmed by several articles. For example this one (point 2.4 we see F^(-1)) and this one (scheme on the second page). In fact in first link formula below F(^-1) is DCT-2. Why? As wiki says DCT-3 is an inverse transform.

And if we look at implementations they take DCT-2 too. I checked two of them: pyAudioAnalysis (stMFCC method) and implementation (mfcc method) by the guy who writes for practicalcryptography.com.
And this is also confirmed by several links: here and here (second page).

Well.. who is right?

I don't believe there is a right answer.  DCT or IDCT will achieve the same purpose: decorrelating to put most energy in fewer coefficients.  Whatever you do next with your MFCC (compression, feature extraction) will work with either approach.

Similarly, I think there are implementations that do an FFT instead of an IFFT for the first stage of the MFCC extraction.


What is really the Mel-filterbank? After applying STFT to a signal, you typically need to convert the frequencies to the Mel scale. Typically the frequency values range from $0$ to $N/2$, where $N$ is the width of each sliding window.
I have read so many articles about Mel Spectrogram. Basically they all say to convert frequencies to Mel Scale in the following way:

Partition the Hz scale (i.e. partition a sequence of numbers ranging from $0$ to $N/2$). In this way we will have a number of bins. This number is for example indicated in librosa by "n_mel". These bins are called frequency bands.

Covert each bins into Mel Scale using the formula $2595*log(1+f/700).$

Application of triangular filters for each bins to capture the energy at each frequency band and approximate the spectrum shape.

It seems that when this is buried in the articles, nobody really knows why these operations are carried out. I find it difficult to fully understand:

Why divide the Hz scale into bins?

What exactly is meant by a frequency band?

Why do we apply triangular filters?

What exactly is a triangular filter? How is it defined mathematically and from what concept?

Can anyone explain exactly these things to me or recommend a book where these steps are explained in detail? There seems to be no book that explains these concepts exactly.

This was my first google hit: https://medium.com/analytics-vidhya/understanding-the-mel-spectrogram-fca2afa2ce53
It seems to explain both how and why?
Mel frequency binning is afaik only used when human perception is involved. Either directly, or for signals that have probably been adapted to human hearing (such as human speech). It is a way to represent signal information in a way that mimics «low level auditory perception», or could be seen as removing perceptually redundant information.
I have not understood why it is that we generally have the same temporal reolution in narrow low frequency bands, as in wide high frequency bands. Historically, it may have been due to limited processing power and the efficiency of a single STFT windowing process. And/or it could reflect the limited neuron firing rate after the acoustic-mechanical filtering carried out by our ear.
But from a dsp perspective, I would expect a variable time/frequency resolution, where high & wide bands gets more temporal resolution.


Estimating the wavelength of a sinusoid I have a sinusoidal signal embedded in noise, and I would like to estimate it's amplitude and frequency, so I can subtract it out.

I realize that an FFT can do the job for me, and that's what I've been doing - I took an FFT of the input signal, and looked at the frequency bin with the highest amplitude, used that as an estimate for the frequency of the wave, and then subtract it out.

But the problem is that the frequency is accurate to only half a bin, and the bin size can be fairly large (since my sample size is typically ~ 100 points), and hence the error in wavelength is big. This leaves fairly large residuals once I subtract.

So as step two, I used the frequency from the FFT as an estimate, and I used the Levenberg-Marquardt algorithm to fit for the frequency, phase and amplitude. This solution worked wonderfully for toy data that I generated, but it fails quite frequently on real data, because of the noise. The amplitude of the noise can be as large (at most) as the amplitude of the signal.

I figured that I can estimate the wavelength using an autocorrelation, but I'm not sure how to interpret the results in the case of more than one frequency component.

Question: Are there any methods of estimating the wavelength of the wave that aren't FFT? I figure since fitting for the frequency is failing with noise, I could try fitting for the wavelength, but any other suggestions are welcome!

Note: I haven't really tried padding with zeros while performing the FFT because the noise seems to go up in the fourier plane if I do, but I'll try it and report the success here.

The first thing to note is that the LMA algorithm is an estimator, and tries to minimize error in the least squares sense, between your data, and a model function of variable parameters.

Typical of LSE based frequency estimators, you need to seed them at values close to the global minima of the cost function, otherwise they will find another local minimum to lock on to, and hence they need to be seeded properly, or 'in the ball park', ostensibly using some apriori information. How to attain this a-priori information? See below.

  This solution worked wonderfully for toy data that I generated, but it
  fails quite frequently on real data, because of the noise. The
  amplitude of the noise can be as large (at most) as the amplitude of
  the signal.

LMA will provide you with an amplitude, phase, and frequency estimate in the LS sense. Of course, the more noise there is in your data, the worse off this estimate is. This is why you will need to pre-process your data to remove as much noise as possible, before the application of the LMA.

You mention zero-padding your data, and while this will not increase your frequency resolution, this will certainly increase your frequency granularity. (Interpolation in the frequency domain). With those facts, I would attempt the following solution:

Pre-Process the data via an (FIR) Band-Pass Filter centered around my frequency of interest.
Perfrom a zero-padded FFT so as to ascertain the peak frequency of interest within the BPF'd data's band.
From this peak, ascertain the frequency, amplitude, and phase, with which to seed LMA with. (Note, make sure to remove the phase offset the BPF filter added for the seeding value).
Run the LMA.

  I'm guessing you mean throw away the highest frequency components etc?
  Would the statistics of my noise matter in that case? I can assume
  it's gaussian, but that may not be strictly true in all cases

The statistics of the noise does not matter in this case. All you are doing is removing frequency bands where noise exists, independent of their statistics. In this way, you are not utilizing DFT bands where your signal does not exist, and retaining the DFT bands where your signal exists.


FFT of a AM modulated signal I am generating an AM modulated sinusoidal wave. The carrier frequency is set at $1000 \ \rm Hz$; the modulation frequency is set to $40 \ \rm Hz$ and $100 \%$ of the amplitude is modulated.
I generated 2 signals with 2 slightly different equations.

Signal 1 follows:
$$\sin(2\pi f_m t) \cdot \sin(2\pi f_c t)$$
while signal 2 follows:
$$\big(1-\sin(2\pi f_m t)\big) \cdot \sin(2\pi f_c t)$$

I don't understand why the first signal doesn't have a $1 \ \rm kHz$ component on the FFT and why both signals do not have the same period.
Both $\sin(x)$ and $1-\sin(x)$ have the same period, and I can see that signal 1 has the amplitude sinus and the carrier sinus in phase, but I can't get my head around those plots and interprete them correctly. Thanks for all the information and explanation you can provide :)

Code snippet to produce the 2 signals:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
import os

#%% Signal 1
fs = 44100        # sampling rate, Hz, must be integer
duration = 1.0    # in seconds, may be float

# Define the time series
t = np.linspace(0, duration, int(duration*fs), endpoint=False) # time variable

# AM - Amplitude Modulation
fm = 40 # Modulation frequency
amplitude = np.sin(2*np.pi*fm*t)

# Output signal
fc = 1000 # Carrier frequency
signal1 = amplitude * np.sin(2*np.pi*fc*t).astype(np.float32)

# Apply FFT
fft_freq1 = np.fft.rfftfreq(signal1.shape[0], 1.0/44100)
fft1 = np.abs(np.fft.rfft(signal1))

#%% Signal 2
fs = 44100        # sampling rate, Hz, must be integer
duration = 1.0    # in seconds, may be float

# Define the time series
t = np.linspace(0, duration, int(duration*fs), endpoint=False) # time variable

# AM - Amplitude Modulation
fm = 40 # Modulation frequency
amplitude = np.sin(2*np.pi*fm*t)

# Output signal
fc = 1000 # Carrier frequency
signal2 = (1-amplitude) * np.sin(2*np.pi*fc*t).astype(np.float32)

# Apply FFT
fft_freq2 = np.fft.rfftfreq(signal2.shape[0], 1.0/44100)
fft2 = np.abs(np.fft.rfft(signal2))

#%% Plot
f, ax = plt.subplots(2, 3, sharex=False)
ax[0, 0].plot(t[:4411], signal1[:4411])
ax[0, 0].set_title('Signal 1')
ax[1, 0].plot(t[:4411], signal2[:4411])
ax[1, 0].set_title('Signal 2')

ax[0, 2].plot(fft_freq1[900:1101], fft1[900:1101])
ax[0, 2].set_title('Signal 1 FFT')
ax[1, 2].plot(fft_freq2[900:1101], fft2[900:1101])
ax[1, 2].set_title('Signal 2 FFT')

ax[0, 1].plot(t[:4411], amplitude[:4411])
ax[0, 1].set_title('Signal 1 AM')
ax[1, 1].plot(t[:4411], (1-amplitude)[:4411])
ax[1, 1].set_title('Signal 2 AM')

There two definitions for the AM modulated signals
The first one is called as the classical-AM (or conventional AM) and is given by
$$x_{AM}(t) = (A_c + m(t)) \cdot \cos(2\pi f_c t) \tag{1}$$
and the second one is called as the DSB-SC (double side-band suppressed carrier) and is given by:
$$x_{AM}(t) = m(t) \cdot \cos(2\pi f_c t) \tag{2}$$
In your code, you use Eq.2 (DSB-SC) for the first example, and Eq.1 (classical AM) for the second example. Their difference is that the first one includes the carrier spectrum (the middle spike at $f_c$) at the output in addition to the modulating signal spectra (side-band spikes). The latter DSB-SC, however, only includes the modulating signal side-band spectra
and not the carrier spectrum, hence the name suppressed carrier.
DSB-SC has the advantage that the tranmission energy is reduced as the permanently radiating carrier energy is saved yielding a more efficient transmission. Whereas the classical AM, which is less energy efficient to broadcast, has the advantage that its demodulation (by analog hardware) is extremely simple requiring what's known as the envelope detector (a diode, a capacitor and a resistance, in addition ot the antenna receiver and audio amplifier circuitry).


Trouble showing Time Invariance of recursive system The system is described with the following recursive differences equation:

$$y[n]-4y[n-1]+4y[n-2]=20x[n]+10x[n-1]$$

now lets say the input is delayed by k, then:

$$y[n]-4y[n-1]+4y[n-2]=20x[n-k]+10x[n-1-k]$$

and now the output by the same k:

$$y[n-k]-4y[n-1-k]+4y[n-2-k]=20x[n-k]+10x[n-1-k]$$

here is the problem , I cannot see how same expression will be acquired .
I've tried performing simple substitution of n-k=m , which leads to:

$$y[m]-4y[m-1]+4y[m-2]=20x[m]+10x[m-1]$$

which is still not quite the same.
Obviously the equation is linear differences equation with constant coefficients , therefore I suppose the system indeed TI but how to mathematically show this specific case.

You've actually proved time-invariance already. It's just a matter of clean notation to see this. Let's use $y_1[n]$ to denote the response to a delayed input $x[n-k]$. The sequence $y_1[n]$ satisfies the following difference equation:

$$y_1[n]-4y_1[n-1]+4y_1[n-2]=20x[n-k]+10x[n-k-1]\tag{1}$$

The delayed response to the input $x[n]$ satisfies

$$y[n-k]-4y[n-k-1]+4y[n-k-2]=20x[n-k]+10x[n-k-1]\tag{2}$$

Now we have to check if $y_1[n]$ and $y[n-k]$ satisfy the same difference equation. Comparing $(1)$ and $(2)$ we see that this is indeed the case.


How to calculate distance transfrom? What does the distance transform encode and how to calculate distance transfrom in the picture below?

My though is that there is a 2 in position (5,7) and the rest is 1 inside the black boundary. Outside would be 0.

The output of the distance transform depends on the applied vector norm (euclidean, maximum, manhattan, ...), see e.g. here or here.

From your example solution, I assume you take the Manhattan distance. Then, your result for the interior of the shape is correct. However, for the exterior of the region, why would it be zero? For the distance transform, the shape of the image does not count, it is only a pixel-based algorithm, where each pixel gets the value that matches its distance to the next non-zero pixel. Apply this to the outer pixels and you have the correct result.


Effect of bandwidth reduction in bit growth I just want to ask, whenever you reduce a bandwidth of signal by filtering, I wonder why do you need to improve the number of bits used to represent that signal.
The number of bit growth is
$$\frac{1}{2}\log_2\left(\frac{BW_{in}}{BW_{out}}\right)$$

Is it to represent to small noise due to filtering.
Thanks!

Consider a two-sample moving average of one-bit data with values {0,1}. The bandwidth is reduced by a factor of two. The possible outputs are {0,0.5,1}. You have gained half a bit, as predicted by the formula.

John


Which transformation in frequency domain equals a x-axis shift of a signal in time domain? I have discrete Fourier transformation results from measurements. Looking at the signal from the time domain perspective, I want to shift the signal on the $x$-axis to the left or right.

Which transformation in the frequency-domain results in a left/right shift in the time-domain?

If $s(t)$ is your signal in time domain, you want to make the operation $s(t+\Delta t)$, which according to the Fourier's transform properties, is equivalent to the operation :
$$
S(f)e^{j2\pi f\Delta t}
$$
$S(f)$ being $s(t)$ Fourier's transform.

In case of a discrete Fourier's transform, which is defined as
$$
S[k] = \sum^{N-1}_{n=0}{s[n]e^{\frac{-i2\pi kn}{N}}}
$$
with $s$ the signal and N the number of points of both signal and DFT.

The operation to be applied in order to induce a shift by a number c of samples to $s[n]$ and obtain $s[n+c]$ is :
$$
e^{\frac{-i2\pi kc}{N}}S[k]
$$

Which in matlab code would be done with :

exp(-j*2*pi*(1:(N-1))*c/N).*S_fdt;


Extended kalman filter for linear system I studied about Kalman filter(KF) for about a week. And I understand that Kalman filter(KF) is suitable for linear system and extended Kalman filter(EKF) can be used for nonlinear system. However, I wonder if extended Kalman filter(EKF) can be used for linear system. I understand that Kalman filter(KF) would be okay for the linear system. But I wonder if there's logical error of using the extended Kalman filter(EKF) for linear system.

the word "extended" might suggest you that EKF is a generalization of the result developed for the linear case. Hence the EKF applied to a linear system reduces to the standard KF.


Gaussian random generator I have quite a straight-forward question. What I aim for is the generation of a certain set of random numbers with a normal distribution (mu = 0, sigma = 1).

Now, the best way to approach the above mentioned gaussian bell is to generate quite a large number of samples.

Unfortunately, in my case, I can only generate a retained subset of samples: usually around 2048 samples, which lead to a poor fit.

I am actually wondering if there is any way to improve the fitness, by first generating a larger set of samples or more subsets of the same size (as 10 times 2048 samples) and then pick those 2048 which approximate the best for desired mean and standard deviation (and the gaussian of course).

What would you suggest or how would you proceed?

To elaborate on Dilip's answer (which is perfectly correct, though in practice the Ziggurat method is much more computationally efficient than Box-Mueller):

One key ingredient missing from your reasoning is whether you want your samples to be independent. It is not clear from your question, but this is the most common situation...

If you want your samples to be independent, then you will have to accept that their empirical mean and variance will not exactly be 0 and 1 - what is known is that the larger the sample, the closer the empirical values will be to 0, 1 - this is absolutely fine - this is a "feature" of the normal distribution rather than a "bug". It is a bad idea to try to select the samples that look more "random" or "well-behaved"; because this very selection process is making them less random! For example, if you come up with a process that always results in 2048-sample blocks with exactly $\mu = 0$ and $\sigma = 1$, it would mean that there are only 2046 possible degrees of freedom in your generation process, and thus that your samples are not independent (given the 2046 first samples, you can guess the values of the 2047th and 2048th using the $\mu = 0$ and $\sigma = 1$ equation... so these last two values are not random!).


inverse continuous wavelet transform and [Parm] in cwtft what is 'parm' means when you set the name of wavelet function in cwtft or icwtft. wave = {wname,[7.6]}. also can I change Fb and Fc when I use 'morl' function in cwtft or icwtft transform? and If not, then how can I reconstruct my signal with cwt transform? cause cwt let me to select optional value for fb and fc (cmorfb-fc). Matlab doesn't have direct function for inverse wavelet transform.

N = 1024;
t = linspace(0,1,N);
y = sin(2*pi*8*t).*(t<=0.5)+sin(2*pi*16*t).*(t>0.5);
dt = 0.05;s0 = 2*dt;ds = 0.4875;NbSc = 20;
wname = 'morl';sig = {y,dt};sca = {s0,ds,NbSc};
wave = {wname,[7.6]};
cwtsig = cwtft(sig,'scales',sca,'wavelet',wave);
sigrec = icwtft(cwtsig,'signal',sig,'plot');

Why won't cmorwavf work?

http://www.mathworks.co.uk/help/wavelet/ref/cmorwavf.html

Could you then construct a matrix of (complex) Morlet wavelet co-efficients (suitably dilated and translated) and simply invert it for the inverse wavelet transform?


How to model transfer functions using Python (without using higher level packages) and plot the magnitude and phase response? I'm trying to model a transfer function in Python and thought I could do it by simply plotting the transfer function at many frequencies. This seemed to work for a 2nd order LPF. See the below figure.

A bit of sample code would be like:
s = np.arange(10 * 2 * np.pi, 1000 * 2 * np.pi, 10) # Creates vector of frequencies in rad/s
lpf = np.zeros(len(s))  # Initializes lowpass filter vector

for k in range(len(s)):
   lpf[k] = lpfilt_func(s[k])  # Calls LPF function with TF = cutoff_freq ** 2 / (s + cutoff_freq) ** 2

To get the frequency response, instead of using s as a real number, I changed s to a complex number using the built-in function complex(0, s) and calculated the magnitude and phase as follows:
for k in range(len(s)):
   lpf[k] = lpfilt_func(complex(0, s[k]))
   gain_lpf[k] = np.absolute(lpfilt_func(complex(0, s[k])))
   phase_lpf[k] = np.angle(lpfilt_func(complex(0, s[k]))) * 180 / np.pi  # Phase in deg

I think this is essentially like setting $s = j\omega$ and then calculating the magnitude as $|LPF| = \sqrt{A^2 + B^2}$ and the phase response as $\angle{LPF} = \angle{(B/A)}$.
Is this a valid method to calculate the frequency response of a transfer function? Is this conversion of s from a real number to a complex number valid?

If the response is real you need not convert it to complex in order to compute the magnitude, but the process is correct for complex and real number since the abs function will use complex conjugate product prior to taking the square root.  If it was real, there would be a trivial phase response since the phase would always be 0 degrees. This would not be a causal system so not a realistic example to have a real frequency response. For complex cases the approach used is correct for magnitude and phase, just not the most efficient.
To confirm results (or to compute directly), use the freqs and freqz commands that are part of scipy.signal for continuous-domain (s) and discrete-domain (z) transfer functions.
For example the transfer function
$$H(s) = \frac{s+1}{s+4s+4}$$ could be plotted using:
w,h = scipy.signal.freqs([1, 1], [1, 4, 4])

Returning the frequency (w) and complex response (h)
Read the help for both of those functions (also available under same name in MATLAB and Octave) for further info.
For MIMO systems the control package has this all worked out as part of the bode function.


Downsampling and Then Upsampling Given this system:

I need to show the $\mathcal Z$-transform of $y[n]$ as a function of the $\mathcal Z$-transform of $x[n]$.

Now I know that for downsampling alone:

$$Y(z) = \frac1M\sum_{m=0}^{M-1} X\left(e^\frac{-j2\pi{m}}M\cdot z^\frac 1M\right)$$

And that for upsampling alone:

$$Y(z) = X\left(z^M\right)$$

My problem here is that I don't quite understand how upsampling would affect the downsampling equation. Say we call the point between downsampling and upsampling $w[n]$, we get:

\begin{align}Y(z) &= W\left(z^M\right)\\
W(z) &= \frac1M\sum_{m=0}^{M-1} X\left(e^\frac{-j2\pi{m}}M\cdot z^\frac 1M\right)\end{align}

I'm not sure if I will be getting:

$$Y(z) = \frac1M\sum_{m=0}^{M-1} X\left(e^\frac{-j2\pi{m}}M\cdot z\right)$$

Or:

$$Y(z) = \frac1M\sum_{m=0}^{M-1} X\left(e^{-j2\pi{m}}\cdot z\right)$$

Or even something else.

I'll be very happy for some clarification on this subject because I'm very baffled by it and I cannot find a solution anywhere, thank you.

You will get

$Y(z) = \frac1M\sum_{m=0}^{M-1} X(e^\frac{-j2\pi{m}}M*z)$

If you were to upsample first and then downsample you would get

$Y(z) = \frac1M\sum_{m=0}^{M-1} X(e^{-j2\pi{m}}*z)$

which just simplifies to (since $e^{-j2\pi{m}}=1$ for $m\in\mathbb{Z}$)

$\frac1M\sum_{m=0}^{M-1} X(e^{-j2\pi{m}}*z) = \frac1M\sum_{m=0}^{M-1} X(z) = X(z)$

Basically, if you upsample and downsample you can perfectly reconstruct the signal but if you downsample and then upsample, you are aliasing in the summation over rotations in the z-domain of the input.


An explanation of "phase" of a filter I always read the word "phase" (like linear phase, phase shift...) in DSP but I'm still not sure what it supposes to mean, in intuition and also in practice.

The phase of a sinusoid $s(t)=A_0\cos(2\pi f_0 t + \phi_0)$ is $\phi_0$ radians.

If this sinusoid goes through an LTI system with frequency response $H(f)$, then the ouptut is $y(t) = |H(f_0)| A_0 \cos(2\pi f_0 t + \phi_0 + \angle H(f))$. So, the phase of the output is different than the phase of the input -- the system introduced a phase shift in the signal.

Note that you can interpret the phase shift as a time delay. If $T_0=1/f_0$, then a phase shift of $2\pi$ corresponds to a delay of $T_0$. This allows you to calculate that the delay $\Delta$ for a phase shift $\theta$ is $$\Delta = \frac{T_0 \theta}{2\pi}.$$ Note that the delay is a function of the sinusoid's frequency. An LTI system that produces the same time delay for all its inputs is said to have linear phase.

Note that LTI systems without linear phase will in general distort their input, even if they have constant gain; this is the reason we want to design linear phase systems. In the case of filters, usually we only require linear phase in the filter's passband.

(Note that this explanation is for continuous time, but the ideas are also valid for discrete time).


Intersymbol interference due to limited channel bandwidth What is the difference between limited-bandwidth induced intersymbol interference (ISI) and ISI due to multipath propagation.
I noticed that the delay spread in indoor environment is in the range of nanoseconds, and is smaller than the symbol duration, yet ISI problem appears. Is this due to the limited bandwidth of the channel? Are there other causes for ISI?

Here is a graphic I have that demonstrates ISI due to multipath propogation:

Notice the eye diagram in the upper right hand corner showing what we should receive if there was no ISI, or other noise for that matter. (This eye diagram is showing the waveform vs time over 2 symbols again and again so that we can see the effects of ISI clearly. If you have any questions about eye diagrams in general, ask it as another question and we can go into further details on that).

This waveform reaches our receiver over multiple paths, each with their own delay and attenuation as shown by the multiple waveform vs time plots in the lower right hand corner of the figure. The eye diagram as received looks like the horrible mess due to this interference that is depicted in the lower left hand corner. (Luckily this can easily be undone through channel equalization in the receiver).

The variation in time of all the different propagation paths is known as the delay spread. Below is a plot I have shown in other posts showing the typical delay spreads for macro-cellular channels, as tabulated by Professor Paul Walter Baier of the University of Kaiserslautern, Germany (Vielen Dank, Herr Professor Baier!):

The significance of delay spread is whether we end up with Frequency Selective Fading (which means deep selective frequency nulls within the passband), or Flat Fading (which means the entire passband power level will vary, equally in frequency across the passband).  In the figure below this is demonstrated by two dominant paths with a Delay between the paths before summing together at the receiver. On the right is the resultant magnitude of the received signal versus frequency for two different delay conditions, the red plot is the result of a relatively long delay, while the blue plot is with a relatively short delay. The shaded region is our example channel at a particular frequency of interest. Notice for the red (long delay) given our channel bandwidth that we have the possibility of a deep null, selectively within our passband, while for the blue (short delay) the power level across the passband is generally much flatter. If you are familiar with the characteristics of an FIR comb filter then you know exactly what is going on; a fixed delay is a phase that varies linearly with frequency, so at some frequencies the phase will be $\pi$ radians so that the receive signals will cancel each other, while other frequencies the phase will $0$ (or equivalently multiples of $2\pi$) so that the signals will coherently add. Since the phase of a fixed delay varies linearly with frequency we see the repetition of this over frequency as the delayed path passes repeatably through $0$ and $\pi$ conditions as we vary the frequency of interest. A longer delay will increase more in phase for a given change in frequency, hence the tighter frequency response for that case.

Plot showing delay vs fading:
RED: Frequency response for a long delay.
BLUE: Frequency response for a relatively shorter delay.
SHADED AREA: Example frequency channel of interest.

Limited bandwidth induced ISI is easily understood in the time domain: A relatively tight filter in the frequency domain MUST have a long impulse response in the time domain (consider filtering as a low pass filter, and in order to really smooth the signal, you must average for a long time). So if you have a filter with a long impulse response, it goes to say that it's response to one symbol will extend into other symbols, causing Inter-Symbol-Interference. The trick is to provide relatively tight filtering while not causing ISI, and this is done with the class of Nyquist filters. They do this by having zero crossings in the impulse response at the symbol spacing. Thus the response from one symbol is always zero at the decision-sampling location of all other symbols. In fact the clean eye diagram in the first plot was created with such a Nyquist filter (the raised cosine filter specifically). All the deviations of the waveform are the result of trailing impulse responses of earlier symbols.

For further details on that see: Why root raised cosine filter can eliminate intersymbol interference (ISI) ?


Applying Circular Cross Correlation in MATLAB I have searched Google for circular cross correlation using Matlab, and I have found it only for one dimensional signals.
Could you please help me implementing circular cross correlation between two images using MATLAB?

Cross correlation can be implemented in the frequency domain using FFT by multiplying signals.
For a cyclic cross correlation you can use:
cc = ifft(fft(x).*conj(fft(y)));
You can also calculate the linear cross correlation using FFT by zero-padding the signals before the FFT.
Related: Cross correlation with FFT and fftshift
Edit:
In order to calculate the cross correlation of two images you have 2 options:

x and y should be the grayscale representation of the images
Repeat the process 3 times. Each time x and y are a single channel (r,g or b) from each image


How can I automatically find the boundary of a high intensity region inside an image? I have an image I(x,y):

Red indicates high values and blue low values. Everything in black is drawn onto the image afterwards.
The black dot (xp,yp) represents a point clicked by an user and the black curve a boundary around the high value region around the user clicked point.

How can I automatically find this boundary?

As can be seen from the picture; I do not want the boundary to propagate through narrow gaps.

One method that I have had lots of success with, and others have too, is Circular Shortest Path (CSP). Here is one paper, but there are quite a few more and possibly some code floating around the web. Another technique that are perhaps more popular, and for which code might be more common, are snakes (active contours).

These techniques work by finding a connected path that minimises a certain function. In your case, it would be to find the path that goes through many edge pixels while still remaining somewhat round. The roundness constraint is what keeps the path from going off down the narrow gaps. Snakes are different from CSP in that the boundary is 'grown' out from a seed area; CSP performs a polar transformation about the centre point.

Another possibility is some kind of morphological operator. These can be used to find regions of a particular size. One would convert the image to binary, the apply a morphological opening using a structuring element of a big enough size to fit inside the area of interest, but not anywhere else.


Why Does the Odd Multiple of $ \frac{\pi}{4} $ on Gaussian Cause Loss in Repeatability Under Image Rotations? I couldn't figure out below paragraph on SURF paper and hope that someone can help me to understand it. Why image rotations around odd multiples of $\frac{\pi}{4}$ lead to a loss of repeatability?

  Bay H., Ess A., Tuytelaars T. Van Gool L. - Speed-Up Robust Features
  (SURF), page 3, column 2

  Gaussians are optimal for scale-space analysis, but in practice they
  have to be discretized and cropped (figure 2 left half). This leads to
  a loss in repeatability under image rotations around odd multiples of
  $\frac{\pi}{4}$. This weakness holds for Hessian-based detectors in general.
  Figure 3 shows the repeatability rate of two detectors based on the
  Hessian matrix for pure image rotation. The repeatability attains a
  maximum around multiples of $\frac{\pi}{2}$. This is due to the square shape of
  the filter.

The answer boils down to 2 issues with the practical approximations of the Gaussian Kernel:

Though the Gaussian Kernel is radially symmetric its discrete approximation has a rectangle support. Unless this support will have infinite length a rotation by any angle different from a multiplication of 90 degrees will yield a shape which has to modified to fit a rectangle. Either choice (Cropping or extending) will yield different results as the weights will be different from the original approximation.
One must pay attention that the longer the support of the filter the less this issue is pronounced.
If I remember correctly, the SURF algorithm uses Box Blur based approximation of the Gaussian Kernel. This approximation, though can be very close visually, doesn't retain the radial symmetry of the Gaussian Blur hence even with large length ha different results when rotated.

I hope this clarifies the issue.


IQ Mismatch and Image I was reading the following paper about direct conversion receivers from TI: link.

In this paper they compare the performance of direct conversion and heterodyne receiver with IQ mismatch.

The paper says that heterodyne architecture performs better with IQ mismatch because the image is separated by 2 times Intermediate frequency (IF). But for direct conversion receivers the image and the desired signal lie on top of each other (figure 6 in the document).
The paper says that the image signal is created by the IQ mismatch. But I don't understand what is this image signal.

How can IQ mismatch create an image signal? And:
Why it is separated by IF in heterodyne systems?

Consider an IQ signal with amplitude or quadrature mismatch such as these QAM waveforms depicted below:

These can be decomposed into the sum of a perfectly balanced IQ signal and another signal that is modulating 0 to 180° with any fixed phase rotation. The perfect balanced IQ signal will be translated from the positive RF spectrum to baseband (as I depicted in the other post here Can Carrier Offset cause Image Problems) as a complex signal and thus if there was a carrier offset that would be maintained as an offset in the complex baseband signal but easy to correct for with a subsequent complex frequency translation. The signal that is on one axis only is a (potentially phase rotated) real signal and this will have alias distortion.

Consider an example with:

$$A\cos(\omega t)- jB\sin(\omega t))$$

If $A=B$ then this is simply $Ae^{-j\omega t}$. But if $B = (1+\Delta)A$ for example representing an amplitude imbalance of $\Delta$ then we get:

$$A\cos(\omega t)- j(1+\Delta)A\sin(\omega t))$$

$$=Ae^{-j\omega t} - j\Delta\sin(\omega t)$$

So the down-conversion that would have been image free with perfect IQ balance for a real input as follows:

$$r(t) = 2\cos(\omega_c t)$$
$$LO = \hat r(t) = e^{-j\omega_{LO} t}$$

$$r(t)\hat r(t) = 2\cos(\omega_c t)e^{-j\omega_{LO} t} = (e^{j\omega_c t}+e^{-j\omega_c t})e^{-j\omega_{LO} t}= e^{(\omega_c-\omega_{LO})t}+ e^{(-\omega_c-\omega_{LO})t}$$

Where the first component $e^{(\omega_c-\omega_{LO})t}$ is the only one that would appear at baseband.

With the amplitude impbalance as depicted above this example now becomes:

$$(e^{j\omega_c t}+e^{-j\omega_c t})(e^{-j\omega_{LO} t}- j\Delta\sin(\omega t)) $$

Consisting of the same result above with a new term given as:

$$(e^{j\omega_c t}+e^{-j\omega_c t})(-j\Delta\sin(\omega t))) = -(e^{j\omega_c t}+e^{-j\omega_c t})\frac{\Delta}{2}(e^{j\omega_{LO} t} - e^{-j\omega_{LO} t})$$

$$= -\frac{\Delta}{2}[e^{j(\omega_c+\omega_{LO})t}+ e^{j(\omega_c-\omega_{LO})t}+e^{j(-\omega_c+\omega_{LO})t}+e^{j(-\omega_c-\omega_{LO})t}]$$

Where we see we get two components at baseband, one $-\frac{\Delta}{2}e^{j(\omega_c-\omega_{LO})t}$ is constructive with our original signal (in this case decreasing the amplitude but no other distortion, while the other is an image $-\frac{\Delta}{2}e^{j(-\omega_c+\omega_{LO})t}$ which would cause serious distortion if the RF spectrums were not symmetric (as is typical for quadrature modulated signals) or if we had any frequency error between the carrier $\omega_c$ and the LO $\omega_{LO}$.

You can create a similar example with quadrature error, as well as the appearance of LO leakage with DC offset.

Consider this with the spectrums in a Zero-IF receiver as depicted below; the LO is depicted as a negative frequency only - with quadrature or amplitude error the positive frequency would start to appear starting with a very low level. You can follow the related spectrums when that occurs knowing the output is the convolution of the RF input spectrum with the LO spectrum to see the effect on the output.

With the super-heterodyne those overlapping images above would be separated appearing at their IF frequencies (and balanced since an IQ down-conversion would not be necessary from RF). But even a super-heterodyne receiver will have this same issue ultimately with quadrature and amplitude balance from the down-conversion at the IF frequency, since it does the down-conversion in two stages: first from RF to IF where we would then effectively have a Zero-IF down-converter from the IF frequency to baseband (either digitally or in the analog). However traditionally it has been a lot easier to achieve better amplitude and phase balance in the down-converter at the lower IF frequencies, as well as tighter filtering hence the superior performance of the super-heterodyne.

Thankfully with modern signal processing this is easy to detect and correct for, which has increased the popularity and use of Zero-IF receivers (and Direct RF transmitters). Below shows a typically correction approach for amplitude and phase imbalance:


Extracting top & bottom fields from a YUV 4:2:0 Frame I am wondering what would be the right approach to extract top and bottom fields from a YUV 4:2:0 frame.

In the YUV 4:2:0 representation, Luma and Chroma samples are located as shown below:

Y Y Y Y Y Y Y Y
c     c     c    c
Y Y Y Y Y Y Y Y

Y Y Y Y Y Y Y Y
c     c     c    c
Y Y Y Y Y Y Y Y

Hence, one needs to be careful while extracting chroma top and bottom fields (due to the chroma offset problem). I have tried a few methods (using asymmetric filters and such), but none have been satisfactory.

Are there some good resources (papers) on this problem?

Assuming the encoder properly low-pass filtered the chroma before sampling (which may not always be true), chroma samples from the adjacent macroblocks might need to be used to properly reconstruct (interpolate) any per pixel color.  Those samples may exist due to motion compensation from previous fields, or might need to be inferred by reflecting, duplicating or extrapolating edge data.


What are i-vectors and x-vectors in the context of Speech Recognition? I have read that i-vectors and x-vectors are widely used in speaker recognition tasks but I don't get the difference between them and how exactly they work. Can someone explain it starting from the ground to a bit technical?
I came across following papers talking about these terms:
https://www.danielpovey.com/files/2018_icassp_xvectors.pdf
https://arxiv.org/pdf/1611.00514.pdf
https://dl.acm.org/citation.cfm?id=2846832

The i-vectors and x-vectors share the ability to represent speech utterance in a compact way (as a vector of fixed size, regardless of length of the utterance).

The extraction algorithms of i-vectors and x-vector are quite different. The x-vector concept is newer and the name of the method is similar to "i-vector" to suggests that this representation can be used instead of i-vectors in state-of-the-art speaker (or language) recognition systems.

There are many good sources explaining how i-vectors work, I recommend this presentation from Najim Dehak (who is the inventor of the i-vector concept).

The x-vectors are embeddings extracted with DNN, but more detailed explanation can be found here (there is also a youtube lecture explaining the concepts presented in this paper).


For 2D signals can it be said that the frequency response is the same as the Fourier transform? Say I have a signal: [1, 1 , 1; 1 100 1; 1 1 1];,

I wanted its frequency response, then can I take its FFT?
Why does MATLAB have two functions to do the same thing, fft2 and freqz2?

Yes for 2D signals you can take a 2D fft, and if the 2D signal is represented in the time domain, then its fft is represented in the frequency domain.  2D FFT's have many other interesting applications, for example image creation in synthetic aperture radar (SAR), where an inverse 2D FFT of radar reflections results in the creation of an image.

If your interest is in the frequency response, then you most likely want a continuous function of frequency, which would be the Discrete Time Fourier Transform, as opposed to the Discrete Fourier Transform (which the FFT computes) which is a discrete function of frequency. freqz returns samples of the DTFT, while fft returns samples of the DFT.

The fft (and fft2) is an algorithm that implements the Discrete Fourier Transform (DFT), while freqz (and freqz2) is an algorithm that returns samples of the Discrete Time Fourier Tranform (DTFT). The DFT is a discrete function in frequency while the DTFT is a continuous function in frequency.

Below are some graphics that will help illustrate the difference between the DFT and DTFT, as well as show for reference the Fourier Transform from the continuous time domain (The CTFT).

A key take-away to see from these and waveforms in general whether they be time or frequency domain waveforms, is that if a waveform is repeating in one domain, it must be discrete in the other domain. (Or stated another way if a waveform is continuous in one domain, it must not be repeating in the other). Pay attention to this relationship as you review the plots below. And this relationship really helps see the difference between the DFT and DTFT.  Another example some may be familiar with is sampling in the time domain ---- this is a repeating function in the frequency domain (called aliases), which is why we only need to concern ourselves with the frequency response from $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate; to bridge the analog and frequency domain worlds I like to still view the digital frequency domain as extending to infinity, but knowing that everything beyond $f_s/2$ keeps repeating. This has been a helpful intuitive view for me when dealing with multi-rate processing, or when going from Analog to Digital with undersampling etc...  I have other posts that dive into that path in more detail, but specific to CTFT, DTFT and DFT:

First the CTFT. Notice that it extends in time from minus infinity to plus infinity (even if the waveform itself is zero). There is no implied repetition in time, therefore the frequency domain is a continuous function. The frequency domain also goes from minus infinity to infinity, with no implied repetition, so the time domain is also continuous.

With the DTFT, the only difference from above is we now sample in time the non-repeating time domain function. This one change causes the frequency domain to repeat. But notice that the frequency domain is a continuous function. (Because the time domain is not repeating).

Now finally the DFT, in this case we limit the time domain over a finite duration (similar to the Fourier Series Expansion), which I argue is identical (mathematically and intuitively) to repeating in time, as evidenced by the frequency domain becoming discrete (and the frequency domain is repeating as well for the same reason). This detail may require further explanation than i am giving here; for further details on the implied repetition please see my other post at this link:

Intuition for sidelobes in FFT

Regarding practical uses for the DTFT, which the freqz command implements: The freqz command is used to plot the frequency response for digital transfer functions. In this case we are interested in a continuous function in frequency, as the system can take any value frequency (even though it is a sampled system) at its input within its allowed Nyquist range. If you didn't use the freqz command, you can emulate it's effect with zero padding. Reviewing the plot of the DTFT above should illuminate how zero padding can approximate the DTFT (for the "real" DTFT you would have to pad to infinity!). This also shows clearly why zero padding results in interpolating more samples in the frequency waveform (without changing the underlying waveform itself).  If we add twice as many samples by adding zeros in the time domain equal to the number of original time domain samples, we end up with additional samples in the frequency domain, except each new sample will be in between the original frequency domain samples (interpolated)--- as we add more zeros we get more samples in between, approaching a continuous function in frequency.