Hilbert FIR Filter for Q - Matching Filter for I For an ultrasound application, I am collecting data from an ADC using high spec Altera FPGA which grabs data and can then do some DSP processing that would be otherwise too slow for a desktop to do.

Part of the DSP required is the precursor to envelope detection. The incoming signals from the ADC contain only the in-phase component of the signal. As a result I also need to generate the quadrature component to be able to then do the sum of squares of I/Q to get the envelope. To get the Q component a Hilbert FIR filter is called for.

Using MATLAB I've been able to design a suitable Hilbert filter for the application. However there are limited resources, so the filter is limited in terms of taps. Using the "Altera FIR Compiler II" tool I've then generated the filter for the FPGA.

I've been doing some simulation and the filter does correctly produce the quadrature component - for example if I feed in a sine I get a cosine. The trouble is the filter causes an attenuation of the signal - i.e. it isn't unity gain. Part of this will be due to coefficient rounding, and part of it due to the calculations being done in fixed point.

Using this filter if I try to do envelope detection using sum of squares, the resulting envelope is not brilliant. in the sine/cosine example I would expect the result to be a DC level, however it is in fact a DC level with a sine wave floating on it. This is presumably due to the Q signal being slightly attenuated. In fact if I multiply all the I values by a constant scalar found by trial and error, the envelope improves in accuracy.

The trouble is the pass band of the Hilbert filter isn't perfectly flat, there are <0.1dB ripples in it - there is a limitation on the number of taps I can implement due to resource restrictions. As a result the envelope of signals that are not pure tones get distorted again because my scalar for the I data is perfectly flat (its a simple multiplication).

What I am wondering then is can I improve the envelope by using a FIR filter on the I data as well. Rather than just delaying the I data and multiplying it by a constant scale factor, is it possible to design a FIR filter whose pass band response is matched to the Hilbert filter? If so, is there a simple way to do it using MATLAB?

Clay Turner has an interesting couple of papers.

what you want to do is compute (using MATLAB firpm() or firls()) two filters that shift the phase, one at -45° and the other at +45° relative to some linear phase angle that corresponds to a constant delay.  you design one and then flip the coefficients around for the other.  they will both have the same magnitude frequency response and be 90° out of phase of each other for every frequency.


Transfer function of double exponential smoothing? Double exponential smoothing
a.k.a.
Holt-Winters smoothing
tracks level and trend of a time series in coupled IIRs:
$\qquad$ In:  $Y_t$, t = 0 1 2 ...
$\qquad$ State: $L_t, T_t \quad$ -- level and trend
$\qquad$ Out: estimate $\hat{Y}_{t+1}$
$\qquad$ Parameters: a, b  (a.k.a. $\alpha, \beta$)
$\qquad$ Step equations:
$\qquad$ $\qquad L_t = a Y_t + (1 - a) (L_{t-1} + T_{t-1})  \quad$ -- level
$\qquad$ $\qquad T_t = b (L_t - L_{t-1}) + (1 - b) T_{t-1}  \quad$ -- trend
$\qquad$ $\qquad \hat{Y}_{t+1} = L_t + T_t$

How can I calculate the transfer function of this smoother, for given $a$ and $b$ ?

(Two possible approaches:

manipulate the step equations into ratio-of-polynomials form
brute force iterate with input $e^{2 \pi i f t}$: converges slowly for small $a, b$.

The first of these is beyond me, the second unsatisfactory.)

I would use the $\mathcal{Z}$-transform. Let me change the notation slightly by using $n$ as the discrete time index, lower case letters for time domain quantities, $x[n]$ as the observed input data sequence, and $y[n]$ as the output sequence, which is an estimate of $x[n+1]$. With this notation, the time domain equations are

$$\tag{1}\begin{align}
l[n]&=ax[n]+(1-a)(l[n-1]+t[n-1])\\
t[n]&=b(l[n]-l[n-1])+(1-b)t[n-1]\\
y[n]&=l[n]+t[n]
\end{align}$$

Taking the $\mathcal{Z}$-transform of these 3 equations gives

$$\begin{align}\tag{2}
L(z)&=aX(z)+(1-a)z^{-1}(L(z)+T(z))\\
T(z)&=b(1-z^{-1})L(z)+(1-b)z^{-1}T(z)\\
Y(z)&=L(z)+T(z)
\end{align}$$

From the second equation in (2) we can express $T(z)$ in terms of $L(z)$:

$$T(z)=\frac{b(1-z^{-1})}{1-(1-b)z^{-1}}L(z)\tag{3}$$

Plugging (3) into the first equation of (2), we can express $L(z)$ in terms of $X(z)$:

$$L(z)\left[1-(1-a)z^{-1}\left(1+\frac{b(1-z^{-1})}{1-(1-b)z^{-1}}\right)\right]=aX(z)$$

from which, after some algebra, you get

$$L(z)=\frac{a\left[1-(1-b)z^{-1}\right]}{1+[a(1+b)-2]z^{-1}+(1-a)z^{-2}}X(z)\tag{4}$$

Finally, plugging (3) into the last equation of (2) gives

$$Y(z)=\left[1+\frac{b(1-z^{-1})}{1-(1-b)z^{-1}}\right]L(z)=\frac{1+b-z^{-1}}{1-(1-b)z^{-1}}L(z)$$

and combining with (4) results in a relation between the output $Y(z)$ and the input $X(z)$:

$$Y(z)=\frac{a(1+b-z^{-1})}{1+[a(1+b)-2]z^{-1}+(1-a)z^{-2}}X(z)=H(z)X(z)\tag{5}$$

where $H(z)$ is the desired transfer function, which is of course a second order IIR filter.


Rendering a spectrogram of audio data - resulting magnitudes are all too close together I'm working on creating an application which will render a spectrogram of the audio in an MP3 file. I'm having some trouble, though - after running the audio samples through my program, the resulting magnitudes are all very close together, and don't seem to be very patterned like one would expect. I'm coloring my spectrogram using a pretty typical black/blue/purple/red/yellow/white gradient. The current result I get is this:

The steps I'm performing to get this are:

Decode the MP3 file using libmad, getting raw signed 16-bit PCM samples back. The samples are stereo, but at this point (for debugging, mostly) I'm just ignoring the right channel. I've verified that this step works by writing the raw samples back out to a file, importing them into Audacity, and then playing the audio - and it sounds correct.
Compute the STFT of these samples (each sample is converted to a complex number of sample + 0i - I realize this isn't the most efficient, but it's very simple and my understanding is that it should produce the same result). I'm not scaling the samples, or converting them from signed to unsigned, or anything else. I'm using a window size of 8192 samples, which is about 185ms of audio data (sampled at 44.1KHz). I'm also applying the Hann function to each of these windows, and overlapping them by about 800 samples. I'm using my own implementation of the FFT, however I've verified it by comparing its outputs to those from Numpy over a very large list of inputs.
Plot the result. Each column in the graph is a DFT from list of STFT results (time), and each row is a particular value in the DFT result (a frequency). I'm only plotting the results in DFT bins [1, N/2 - 1], since I'm using real inputs. I'm computing the magnitude as log10(sqrt(r*r + i*i)).

The color range is from [minimum magnitude, maximum magnitude]. I've tried shrinking the range, thinking that outliers may be screwing it up, but even with a range of [5.6, 5.8] the data is still not a real spectrogram. The magnitudes seem to be more or less randomly distributed in that [5.6, 5.8] range.

What issue might be causing this type of output?

I'd love to post code or go into more detail, but as I'm not sure what area the problem might lie in I'm not sure what to expand upon. For what it's worth (although obviously I don't expect anyone to look at all of my code and fix my bug), this project is on GitHub here.

EDIT: I generated a wave file with a single 1000Hz sine wave as suggested, and for the first time ever my application actually produces some (maybe?) useful output. Here's what I get:

Perhaps this provides some more insight as to what's going wrong?

A periodic looking output from an FFT could mean a bug such that the FFT did nothing (output == input, unchanged), or that the input to the FFT only had a single large magnitude or non-zero data point (input == garbage).


causality of the system $y[n] = x(2n)$ Can somebody please tell me why the system $y[n] = x(2n)$ is non-causal ?

I know that causal systems depend on the past and present values of input and
this system satisfies the condition. So why is the answer that the system is non-causal ?

Thank you.

No it does not satisfy the condition. Simply take an example:

$$n = 1 \implies y[1] = x[2]$$

Hence the output value at the present time $n=1$ depends on a future value of the input at time $n=2$. This violates the causality principle.


Kalman filter - question about matrices In order to familiarize myself with Kalman filter I decided to see how it work step by step. I would gladly use your help.
Let's say I have a falling object with air resistance proportional to its velocity. Since my measurements are not perfect let's say that position measurement is gaussian with standard deviation of $15m$, so is velocity with deviation $5 \frac{m}{s}$. Equation for this system is $\ddot{x} = g - \frac{b}{m}\dot{x}$, where $b,m$ are constant. Now how do I find matrices $A$ and $B$. Do I simply convert to space-state form? Also, how do I determine the processing noise $\mathbf{v}$? And finally, how do I find covariance matrices of of processing($\mathbf{Q}$) and measurement ($\mathbf{R}$)? If you are kind enough to help me I would gladly ask you some more questions tomorrow.

According to the scenario you formulated, you only have measurement noise and no input noise. You are correct that you might need a state space model. For example,

$$
\begin{bmatrix}\dot{x} \\ \ddot{x}\end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & -\frac{b}{m}\end{bmatrix} \begin{bmatrix}x \\ \dot{x}\end{bmatrix} + \begin{bmatrix} 0 \\ 1 \end{bmatrix} u + \begin{bmatrix} 0 \\ 1 \end{bmatrix} v,
$$

$$
\begin{bmatrix}y_1 \\ y_2 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix}x \\ \dot{x}\end{bmatrix} + \begin{bmatrix}w_1 \\ w_2\end{bmatrix},
$$

with $u$ always equal to $g$, $v$ the input noise (which assumed to be zero), $y_1$ the measured position, $y_2$ the measured velocity and $w_1$ and $w_2$ zero mean Gaussian noise.

Since we did not take into account input noise, therefore you can also view it as a zero mean Gaussian noise with also zero (standard) deviation. When calculating the covariance matrices you can take the variance (standard deviation squared) on the diagonal, if all random variables are independent (which I assume $w_1$ and $w_2$ are). This would yield,

$$
Q = \mathbb{E}(v^2) = \begin{bmatrix}0 \end{bmatrix},
$$

$$
R = \mathbb{E}\left([w_1\ w_2]^T [w_1\ w_2]\right)= \begin{bmatrix}225 & 0 \\ 0 & 25 \end{bmatrix}.
$$

But depending on whether you want to implement a Kalman filter in continues or discrete time, you might first have to discretize the state space model.


Is real-time Hilbert envelope feasible? Is the Hilbert envelope (analytic signal) feasible to be calculated in audio frames/blocks in real-time, when using FFT with windowing.

What I'm specifically concerned with is that how does it not create or not be able to create discontinuities (in the end points) between successive frames?

As with most things in Engineering the answer is "depends on your requirements". In general you can't do a perfect Hilbert transform unless the signal is strictly periodic (see http://andrewduncan.net/air/ ). Since you can't do perfect you need to define what's good enough for your particular application.

For most typically audio applications, that it's pretty straight forward. The Hilbert transform can be represented as a simple linear time invariant filter. Unfortunately the impulse response infinite and non-causal. For an FIR implementation you need to truncate and delay it. The larger your truncation window and your delay the more precise it will be, however it will also become more expensive computationally and the latency goes up.

If your application requires a really long filter (typically when you need high precision at very low frequencies) than FFT can be useful in an overlap add configuration.

If you just care about relative phase and not absolute phase you can run the signal through a set of differential allpass filters, which is a very efficient way of implementing a band limited 90 degree phase shift.

You can also use advanced filter design methods (least square, Parks McLellan etc.) to optimize any Hilbert approximation for your specific requirements.


What is the noise whitening process? I came across this process when working on OFDM in wireless communication. It's mentioned if we have two copies of a signal resulted from SIMO system (Single Rx Multiple Tx), we can process the noise whitening, after the FFT, in order to improve the quality of the signal.

What does mean the noise whitening, and what's its process?

Here is an example,  in the paper.

After taking the DFT in OFDM, we supposed to get frequency response, so how w[k] is obtained?

Consider a WSS (wide sense stationary) discrete-time random process (RP) with i.i.d. samples $X[n,s]$ whose auto-correlation sequence (ACF function) is $\phi_{XX}[m]$. We call a random process as white-noise if its ACF is the following:

$$ \phi_{XX}[m] = \begin{cases}
{\sigma_X^2 ~~~~~, \text{ for } m=0 \\ 0 ~~~~~~~~, \text{ otherwise} }\end{cases} $$

This also means that the RP is uncorrelated. If a RP does not have this property then it's a colored noise. And it means its samples are correlated.

When a RP is colored, it may be possible to find an LTI filter which will try to remove the correlation between the samples of the noise and hence make it whitened. Such a filter is called either as a decorrelator or equivalently as whitening filter.

It's best described in the context of AR-p (Auto-Regressive) random processes generation through all pole filtering, where the whitening filter will be the inverse of that all-pole AR filter which generates the colored noise from an input white noise.


Dynamically changing cut-off and sampling frequency of a digital filter I have designed a low pass filter that smooths the output coming from an accelerometer attached to a vibrating machine. I designed this assuming cut-off frequency $f_{c1}$, and sampling frequency $f_s$. Now the issue is when I change the machine , I have to change the cut-off frequency and $f_s$ manually, else the output is not as smooth as expected. That is every machine is actually vibrating with slightly changed frequency.

Is there a way the filter learns the appropriate Cut-off frequency $f_{c2}$ and $f_{s2}$, itself, without me to change the code  and perform filtering based on this $f_{c2}$ and $f_{s2}$?

I shall appreciate if some one can provide some web-links so I can explore this area further.

There are a few ways that you could do run-time dynamic filters.

With a program like Matlab or Octave, but I'm guessing that they aren't available to your application.  There might be DSP libraries that have functions for generating filters though.
Calculating from scratch a windowed sinc function.
You could create somewhat crude low-pass filters based off of your starting filter.  For instance, if you wanted to halve the cutoff frequency you would halve the filter sample time.  You would thus create a new filter that had all of the original filter's samples and a sample in between those samples.  You could create the extra samples through splines.  This same approach could be used to create any fractional cutoff frequency change.  A more complicated example would be if you wanted the cutoff frequency to be $\frac{9}{8}f_c$.  If we say that your original filter's samples are at time 0, 1, 2, ..., $n$, then the new filter's samples would be at time 0, $\frac{9}{8}$, $\frac{18}{8}$, ...


different between MVG and joint MVG? Distribution for "joint multi-variate gaussian distribution" (joint MVG):
$$f_{X}(x) = \frac{1}{(2\pi)^{n/2}\prod \limits_{i=0}^{n}\sigma_i} ~~\text{exp}\bigg[-\frac{1}{2} \sum \limits_{i=1}^{n}\Big(\frac{x_i\mu_i}{\sigma_i}\Big)^2 \bigg]\tag{2}$$

Distribution for "multi-variate gaussian distribution" (MVG):
$$f_{\mathbf{x}}(x) = \frac{1}{(2\pi)^{n/2}~|\text{det }C|^{1/2}  }~~\text{exp}\bigg[ -\frac{1}{2} (\mathbf{x}-\mathbf{\mu})^T C^{-1} (\mathbf{x}-\mu)\bigg]\tag{1}$$
Where:
$[n\times1] ~~~\mathbf{x} \Rightarrow \text{sample vector}$
$[n\times1]~~~\mathbf{\mu} \Rightarrow \text{mean vector}$
$[n\times n]~~~ C \Rightarrow \text{covariance matrix (symmetric positively defined matix)}$

What's the difference between the two?

There's no difference in formula:
Try setting the $C=\text{diag}(\sigma_1^2, \sigma_2^2,\ldots,\sigma_n^2)$ in $(1)$ to see how $(1)$ simplifies to $(2)$. So, your formula $(2)$ is just a special case of $(1)$.
In the diagonal case, the elements of the multivariate distribution are jointly Gaussian, which is much stronger than them all just being Gaussian distributed.
Joint Gaussian distribution has a name of its own (unlike other multivariate distributions with identically distributed elements) because "joint gaussian multivariate random varaible" has the very interesting feature that simple MVG doesn't have: that uncorrelatedness implies independence. This is actually something that makes it useful in communications engineering (because rotation and normalization doesn't change joint gaussianness and doesn't introduce interdependencies between components - so, in an IQ receiver, the noise on I and Q are still independent, even after phase correction, for example).


Sample Dataset for Kalman Filter I'm a newbie to Kalman filter. I have found the code online but I was wondering if there is any sample dataset available online to get hands-on with it (for example: CIFAR-10 for classification etc. ).

You can easily create measurement for Kalman Filter model.
Just run the model itself on some noise with starting point and you'll have endless number of data points.

Namely, if your model is given by:

$$ x \left[ k \right] = F x \left[ k - 1 \right] + v \left[ k \right] $$

Then all you need is to generate $ x \left[ 0 \right] $ randomly according to $ P \left[ 0 \right] $ and start moving it according to the above model while generating $ v \left[ k \right] $.


Bit error rate curves, theoretical and simulation The BER is the average number of errors that would occur in a sequence of
$n$ bits.  When $n = 1$, we can
think of the BER as the probability that any given bit will be received in error. Basically, it lies between 0 and 1 --

when BER = 0, we can say no error in the received bits.
If at least one error is seen, we can say with 100
percent certainty that the BER is definitely not zero.
as long as at least one bit was successfully received, we can state that the BER is most definitely less than 1.

When we measure the performance of a digital communication system with and without equalization, with that of the theoretical BER for a particular kind of modulation, we prefer the simulation to be exact match with the theoretical.
My Questions are :

(1) Is theoretical BER the upper or lower performance bound?

(2) Is it desired that the BER tends to zero? Is it reasonable to say that  since this may not happen in real communication, so there is always a quest to come up with an equalizer designed with an objective that the BER tends to zero?

(3) What if the simulation curve with and without equalizer falls below the theoretical? What can be said about the equalization method playing a role and if it is desired that the simulation curve - with and without equalizer is below the theoretical.

(4) What can be said about performance of the system with and without equalizer if the simulated BER is above the theoretical?

(5) Does an equalizer essentially help to improve the BER? By improving -how much improvement is desired? Tradeoff? Do we wish the simulated curve to be an exact match to the theoretical and then say that the equalizer is functioning properly in removing intersymbol interference and noise?

What is the interpretation for the bit error performance?

Note: usually, the "theoretical BER" is called $P(e)$, the probability of bit error (sometimes with other symbols, like $P_b$). A rate implies experiment; the probability implies theory.

1) $P(e)$ is a lower bound for the BER: since practice is messier than theory, the BER will be higher than $P(e)$ (assuming the theory is correct; that is, that it models reality properly).

2) I'm not sure I understand the question. We always want $P(e)$ to be as low as possible, but it is ultimately bounded by the channel capacity.

3) If by "below the theoretical" you mean "better than $P(e)$", I'd say something is wrong with your simulation or your assumptions.

4 and 5) You use an equalizer when you have a wideband channel with intersymbol interference. Such a system will always have worse BER than a narrowband channel. The equalizer tries to eliminate the ISI; if the equalizer is perfect, the BER of the wideband system will be the same as that of a narrowband channel, assuming equal $E_b/N_0$. In practice, the equalizer will not be perfect. In this case, theory can help predict the equalizer's performance, and how much BER improvement to expect.


Will binarization / thresholding of image yield better results when performed in a linear color space? All thresholding algorithms I have read about make no mention about color space, except sometimes using pixel values between 0 and 255 as an example.

But since the intensity output by a monitor of a greyscale is $greyscale^{gamma}$, the perceived difference between 1 and 2 is less than the difference between 254 and 255. Provided this fact, should it not improve accuracy to perform binarization on intensities rather than greyscale values?

Disclaimer: What I am about to say is a general rule of the thumb, there are exceptions.

In short: Linear space is more useful for computer vision while non-linear space is more useful for showing the final image to humans, without outputing a quantitative analysis.

This highly depends on your problem.

Linear space is more related to physical measurements. Thus, thresholding in it is more natural when you are dealing with measurements such as: counting instances, location of objects, etc. Noise behaves in a predictable way - variance of the signal is linearly dependent on the mean value of the signal. There are many other properties that make your algorithms more straightforward and easier to implement. Among other things, when you are selecting a threshold, you can use the linearity of the signal.

On the other hand, the non-linear space is more related to human perception of light intensity. Thus, it is more suited for image processing problems in which there is no measurement, but rather showing a good image to the observer, or measuring something in a way that humans do.

Examples:

Linear space -
Counting cells in microscopic images, finding defects in a wafer, tracking cars on a road. (Some image processing operations must be done in linear space due to the physics of the process involved, but the final image is almost always in non-linear domain)

Non-Linear space - Performing image processing pixel crunching operations to show a nice image on your tablet, attempting to find defects on a carpet texture, performing image segmentation of medical images in order to assist a human interpreting a scene.


What are imaginary signals? Reading about complex signals, I learn a phasor is composed of real and imaginary components as $e^{2\pi j \cdot \omega}$ and $e^{2\pi ⁻j \cdot \omega}$. Why is this? How can a signal be imaginary?

How can a signal be imaginary?

It's not imaginary in the everyday sense of meaning something that doesn't exist.
"Real" and "imaginary" have technical meanings in signal processing, and more widely in mathematics, that are different from their everyday meanings.
In signal processing "real" denotes a signal component that is in phase with (i.e. has a 0° phase shift from) some reference signal.  Imaginary denotes a signal component that is in quadrature with (i.e. has a 90° phase shift from) the same reference signal. The reference signal can come from a local oscillator. (In DSP or SDR equipment the local oscillator might be a mathematical representation of one, rather than an analog circuit.) In a receiver we might seek to synchronise the reference signal to the received signal.
In signal processing equipment we also refer to real and imaginary as I and Q for in-phase and quadrature, respectively. Some equipment has sections which are divided into separate I and Q channels. You can measure the signals flowing through a Q (imaginary) channel, implemented in hardware, and find they are just as real (in the everyday sense) as the signals flowing through the I (real) channel.
The term "imaginary" was originally coined in the 17th century by René Descartes as a derogatory term and has been confusing students ever since.
For why $e^{j \omega t} = \cos  \omega t + j \sin \omega t $ see Euler's formula. Mathematicians use $i$ to represent the imaginary unit. Electrical engineers use $j$ to represent the imaginary unit because we use $i$ to represent (instantaneous) current.  You can use either representation but please choose one and stick with it.
(This question was originally posted in the amateur radio SE group. This answer was drafted as a reply there and leans towards radio usage.)


Signal to Quantization Noise Ratio (SQNR) I have this problem in my homework (following picture).

How do find the signal to quantization noise ratio at the output?

I don't want the solution. I just want to know a better approach to tackle it.

Because the way I am trying to do is -  take the inverse fourier transform of H(z) and convolve it with x(n) and find the power. Then divide this recently found power by the noise power. As easy it sounds (probably correct) it is computationally intense.

I end up with 4 convolutions just to find the y(n).

No need to do the convolution to find the output. Just use the fact that $z^n$ is an eigenfunction of the system with $z=e^{j\omega_0}$. This will lead to the following rule

$$x[n]=A\sin(\omega_0 n+\phi) \Rightarrow y[n]=|H(e^{j\omega_0})|A \sin(\omega_0 n+\phi+\angle H(e^{j\omega_0}))$$

That is, the output to a sinusoidal signal with angular frequency $\omega_0$ is a similar sinusoidal signal with the same angular frequency, but scaled by $|H(e^{j\omega_0})|$ and phase-shifted by $\angle H(e^{j\omega_0})$.

For SQNR calculation, we only need the scaling $|H(e^{j\omega_0})|$. You probably know how to find the power of a sinusoidal signal, given its amplitude. Right?

Then you need to know how to calculate SQNR when a sinusoidal signal with amplitude $|H(e^{j\omega_0})|A$ gets quantized.


why only impulse signal for convolution? Why do we only use impulse response for convolution for getting output of LTI systems. why we don't use unit step or any other signal for this purpose?

Because the impulse response completely characterizes an LTI system. The reason is an arbitrary input $x(t)$ can be written as an infinite sum of time shifted and scaled Delta functions:

$$x(t) = \int_{-\infty}^{+\infty}x(\tau)\delta(t-\tau)d\tau$$

or alternatively, using the convolution operator, $$x(t) = x(t)*\delta(t)$$

Properties of linearity and time invariance lead to the following sequence of input-output pairs:

$\delta(t) \longrightarrow h(t)$
Time invariance: $\delta(t-\tau) \longrightarrow h(t-\tau)$
Linearity: $x(\tau)\delta(t-\tau) \longrightarrow x(\tau)h(t-\tau)$
Linearity: $\int_{-\infty}^{+\infty}x(\tau)\delta(t-\tau)d\tau \longrightarrow \int_{-\infty}^{+\infty}x(\tau)h(t-\tau)d\tau$

The last input-output pair can be written as $$x(t) \longrightarrow y(t) = x(t)*h(t)$$

This way, any output can be computed using the input $x(t)$, the impulse response $h(t)$, and the convolution integral (or sum, in discrete time).

@a concerned citizen gave you some other hints regarding your question.

On the contrary, if you use the unit step function, you get

$$\int_{-\infty}^{+\infty}x(\tau)u(t-\tau)d\tau = \int_{-\infty}^t x(\tau)d\tau$$

or alternatively, using the convolution operator, $$\int_{-\infty}^t x(t) = x(t)*u(t)$$

Once again, just for comparison, you get the following sequence of input-output pairs:

$u(t) = \int_{-\infty}^t \delta(\tau)d\tau \longrightarrow s(t) = \int_{-\infty}^t h(\tau)d\tau$
Time invariance: $u(t-\tau) \longrightarrow s(t-\tau)$
Linearity: $x(\tau)u(t-\tau) \longrightarrow x(\tau)s(t-\tau)$
Linearity: $\int_{-\infty}^{+\infty}x(\tau)u(t-\tau)d\tau \longrightarrow \int_{-\infty}^{+\infty}x(\tau)s(t-\tau)d\tau$

where $s(t)$ is the step response, the system output when the unit step is presented at the input. The last input-output pair can be written as $$\int_{-\infty}^t x(\tau)d\tau \longrightarrow x(t)*s(t) = x(t)*u(t)*h(t)$$

The latter can be written as $$(x(t)*h(t))*u(t) = y(t)*u(t) = \int_{-\infty}^t y(\tau)d\tau$$

So after all, what you get at the output is the integral of the output you would get using the delta decomposition of an input signal. That is, your output must be differentiated to obtain the actual system output for a given $x(t)$.


how can i call the input data in wavelet transform and generate code for gabor wavelet manually? I am an ECE student, and i am doing my project

I want to generate a code for gabor  wavelet transform. I have no idea  could you help me with a sample program

I have DNA as input sequences and how can use it as input data

"sample sequence  s={ACGTACGTACCCCAGGGATTT} "

thanks and regards!

You're probably way out of your league if you have these problems, and just getting code won't help you with your project. your teacher will spot rather quickly that you didn't write it.

To start, you should read up on wavelet transforms in general. Gabor wavelets are just a specific kind.

You'll learn that wavelet transforms work on signals. That means you'll need to represent your DNA as a signal, i.e numbers.


How to make filters? I am trying to make filters like sepia filter in my app like other photo editing applications , till now I don't have any idea about it , Can any one help me out for it about its algorithms and how to do it ? and Can I implement in opencv using native c++ or we have any of its build in functions ?

You can look at implementations for open source image editing programs (e.g. GIMP). This link from MSDN Magazine, January 2005 by Stephen Toub, entitled "Sepia Tone, StringLogicalComparer, and More" in the .NET matters section specifies a linear transformation of the color space to get a Sepia effect: "New Red = R*.393 + G*.769 + B*.189 New Green = R*.349 + G*.686 + B*.168 New Blue = R*.272 + G*.534 + B*.131" where each channel is normalized to lie in $[0,1]$.

Its just applying the formula on each pixel, so you don't need any special libraries.


Integrating over an image? I am currently trying to implement the method described in this paper. In short we have a system of the form $a=B\times c$. Where

$$a_i = \int d^3r \space w(r)f_i(r)t(r)\quad\text{and}\quad B_{ij} = \int d^3r \space w(r)f_i(r)f_j(r).$$

The solution is given by $c = B^{-1}a$. In this case $w,t,f_i$ are volumes, specifically a sequence of images obtained from an x-ray simulation. When I solve these integrals using MATLAB's trapz function I get a poorly-conditioned $B$ matrix.

Is there any other way to integrate this type of data (images)?
Or a different way to calculate my $a$ and $B$ terms?

You can try to improve one both sides. One issue is that smooth separable kernels $w$ or integration scheme like trapz tend to reduce rank, and pose conditioning problems. Did you use pinv or X = (A)\(b) already?

You can invest on non-separable or higher order integration rules, see for instance n-dimensional simplex quadrature;
You can invest in regularization, for instance with the simple diagonal update: $B_c = B + c*\operatorname{eye}(\operatorname{size}(B))$. If this is not sufficient, there are quantities of constrained or penalized methods (LASSO for instance), implemented for instance in the CVX toolbox (e.g. Premade solvers for specific problems (vector variables)).


why human voice is band limited while it's time limited? According to Wikipedia, time-limited signals can not be band-limited and vice versa. Also we know that human voice is band-limited and has frequencies from $50\textrm{ Hz}$ to $4\textrm{ kHz}$.

So how is the human voice band-limited while it's time-limited ? (consider a man speaking in a period of time to have a time-limited voice signal)

It's true that band-limited signals can't be time-limited (and vice-versa). However, band-limited and time-limited in that context mean "having finite support", i.e. the respective function is identically zero outside a given interval.

In practice, we almost always deal with functions that are time-limited and band-limited, because we never have infinite time and infinite bandwidth available. However, in that context we use the terms "band-limited" or "time-limited" in an approximative sense. I.e., the contribution of the function outside a given interval becomes negligibly small compared to the noise in that interval. We don't require it to be exactly zero.

Using your example of human speech: if you utter a sentence, no matter how long it is, it will be finished at a certain moment. However, the acoustic signal you've produced will be reflected and reverberated, and it will decay exponentially. Soon the reflections will be inaudible and will be hidden in noise, but in theory they go on forever.

When you say that human speech is band-limited to $4$ kHz, then that's actually not true. We often use that frequency as a cut-off frequency because we know that speech remains intelligible if we filter out frequencies above $4$ kHz. Yet, there is a significant contribution above that frequency, and the harmonics continue indefinitely. But again, their contribution becomes negligible and they disappear in noise.


How to get the melody from a signal I'm trying to retrieve a melody from a signal (simple sine with 3 notes).

I've windowed the signal and performed autocorrelation to get the fundamental frequency of each window and then converted frequency to note.

So for a signal with notes "G A E" I am getting something like "G G G A A B E E E E F". It makes sense that notes are repeated and there are some mistakes... But how should I deal with this? If I use a bigger window then I might miss some notes with shorter duration. I can think of a really naive algorithm to get the melody for this example but what are some real techniques that are usually applied?

Some common approaches:

Simple post-processing

Apply a mode filter on your resulting sequence of pitches. This will increase temporal inaccuracies though!

Pre-segmentation

Instead of analyzing your signal according to a grid of fixed-size windows, try to detect first the boundaries between notes using an onset detection algorithm. Once you have performed the onset detection, your signal is divided into segments in which you can make the assumption that the pitch is constant. You then apply your pitch detection algorithm to each segment. This works best for instruments with a sharp attack - like guitar or piano (playing monophonic melodies) - but not very well in which the pitch of a held note is modulated without sharp attacks between notes (flute, whistling, singing voice).

Dynamic programming

You define a function $C_i(f)$ which indicates how unlikely it is that frequency $f$ is the fundamental frequency of the $i$-th signal window. You can use, for example, the reciprocal of the autocorrelation function; or sort your pitch hypotheses. For example, if the 4 highest peaks in the autocorrelation function of the 6th signal window are at 220 Hz, 440 Hz, 500 Hz and 250 Hz, $C_6(220) = 0$, $C_6(440) = 1$, $C_6(500) = 2$ and so on. The higher the function, the most unlikely the argument is to be true fundamental frequency.
You define a transition cost $T(f_1, f_2)$ which indicates how unlikely it is that a window with fundamental frequency $f_1$ will be followed by a window with fundamental frequency $f_2$. This transition cost can embed basic musicological knowledge (which notes are more likely to follow other notes) but will mostly serve as temporal smoothing.
You use dynamic programming to find the sequence $f^*$ that minimizes:

$$\sum_i C_i(f^*_i) + T(f^*_{i - 1}, f^*_i)$$

Intuitively, this searches for a sequence of pitches that is both compatible with your observations, but also doesn't change too much. Adjusting the relative weight between $C$ and $T$ will allow you to control the amount of smoothing.

HMM

Manually annotate a large collection of files. Train a HMM using autocorrelation vectors as your feature vectors, and musical notes as your states.

To transcribe a melody, find the most likely sequence of states that explain your autocorrelation vectors. This is done by the Viterbi algorithm, which will take exactly the same form as the dynamic programming method given above. The only difference is that there is this time a rigorous statistical framework for defining and learning the scoring and transition functions.


Working backwards from PSD to possible signal I have been trying to reconstruct a random signal from its PSD and am running into trouble. I know that many different signals in the time or spatial domains can result in the same PSD-- I am interested in any of these signals. I am not the greatest at signal processing, so I've gotten a bit confused. I am reading this paper. They state:

  Therefore, PSD is approximated as
  $$
S(k)=\frac{2\sigma^2L_c}{(1+k^2L_c^2)^{0.5+\alpha}}
$$
  where
  $$
k=2\pi f, f=i/(N\Delta Z), 0\lt i \lt N/2
$$
  Lastly, N is the number of points along the line.

  ...

  With the magnitude information provided by S(k), we can reconstruct random line edges by applying a random phase to each frequency component of the PSD to form a unique signal in the frequency domain. A line edge with roughness can be simulated by doing an inverse Fourier transform of this signal. Random lines are distinguished through applied random phases.

I have computed values for S(k), per the equation given. I then tried taking the square root of those values, multiplied them by the number of values, and subtracted their mean. Finally, I generated random phase information before doing an inverse FFT. Even after subtracting the mean I still see a very large spike at the beginning and end of the signal.

Does anybody have any suggestions as to what I'm doing wrong?

For those that speak Python, here is the code I've written so far:

import numpy as np
from matplotlib import pyplot as plt

if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.linspace(0,N/2,N)
    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    window = 25.0/46-21.0/46*np.cos((2*np.pi)/(N-1)*np.array(range(N)))
    magnitude = N*np.sqrt(PSD(k))
    phase = 1j*np.random.randint(0,6,N)+np.random.randn(N)
    FFT = magnitude-np.mean(magnitude)+phase
    FFT = window * FFT
    FFT = np.concatenate((FFT[::-1][0:-1],FFT))

    fig = plt.figure()
    ax = fig.add_subplot(111)

    ax.plot(np.fft.ifft(FFT))

    plt.show()

Edit: I forgot about Gibbs phenomenon! As per Bruce Zenone's suggestion, I applied a Hamming window and mirrored my data. I added plots and updated my Python code. I am still seeing some effects of Gibbs phenomenon. Truncating the affected parts, the signal does not seem to be varying with sigma, alpha, and Lc.

Here is my resulting signal:

And a figure from the paper showing the effect of the different parameters:

There are several issues with your code:

The "Gibbs" issue is a non-issue, so you shouldn't need to window.
You are better off generating the complete $-\pi$ to $+\pi$ spectrum that you want, rather than relying on getting the symmetry correct.
The code

FFT = magnitude-np.mean(magnitude)+phase

is changing the phase of the PSD magnitude by addition. It should be using complex multiplication:

FFT = magnitude * np.exp(1j*phase)

The code below is a modified version that seems to do the trick. The top graph shows the PSD and the (absolute value of the) FFT, showing that they overlap exactly.  The bottom graph shows the generated time domain signal for this particular realization of the phase noise.

UPDATE: I just realized that the original code probably shifted the PSD by $\pi$ because the peak was in the center, rather than at the extremities.

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

    if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.concatenate([np.linspace(0,N/2,N/2+1), np.linspace(-N/2+1,-1,N/2-1)])

    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    magnitude = N*np.sqrt(PSD(k))

    phase = 2*np.pi*np.random.randn(N)
    FFT = magnitude * np.exp(1j*phase)

    fig = plt.figure()

    ax1 = fig.add_subplot(211)
    ax1.plot(np.abs(FFT))
    ax1.plot(magnitude)

    ax2 = fig.add_subplot(212)
    ax2.plot(np.fft.ifft(FFT))

    plt.show()


Laplace Transform: zeros and corresponding impulse response $h(t)$ Poles and the impulse response
If our impulse response is in the form :
$$h(t) = e^{-\sigma_0 t}\cos(\omega_0 t) \, u(t)$$
(where $u(t)$ is the unit step function)
And its Laplace transform is :
$$H(s) = \frac{N(s)}{D(s)} = \int_{0}^{+\infty} h(t)e^{-st}dt$$
$$s = \sigma + j\omega$$
Poles are values of $s$ so that $$D(s) = 0 \rightarrow H(s) = +\infty $$
But to understand this, I prefer to look at the integral : it will go to infinity (poles) when $s$ reflects components of $h(t)$. In a way, $e^{-st}$ "probes" $h(t)$. Indeed :

A single real pole ($s = -\sigma_0$) means $h(t) = e^{-\sigma_0t}u(t)$ because : $$\int_{0}^{+\infty} e^{-\sigma_0t}e^{-(-\sigma_0)t}dt = \int_{0}^{+\infty} 1dt = +\infty $$.

Complex conjugate poles ($s = -\sigma_0 \pm j\omega_0$) mean $h(t)$ is an exponentially decaying sinusoid (say $h(t) = e^{-\sigma_0t}\cos(\omega_0t)$) because : $$\int_{0}^{+\infty} e^{-\sigma_0t}\cos(\omega_0t)e^{-(-\sigma_0)t}e^{-j\omega t}dt = \int_{0}^{+\infty}\cos(\omega_0t)e^{-j\omega t}dt $$ which is infinite at $\omega = \pm\omega_0$ (Fourier transform of $h(t)$ without its exponential component, which is a sinusoid).

Complex conjugate poles with $\sigma = 0$ ($s = \pm j\omega_0$) mean $h(t)$ has no decaying component (say $h(t) = \cos(\omega_0t) u(t)$) because : $$\int_{0}^{+\infty} \cos(\omega_0t)e^{-j\omega t}dt$$ which is infinite at $\omega = \pm\omega_0$ (Fourier transform of $h(t)$ which is a sinusoid).

Zeros : a dirac in the impulse response ?
Now, let's look at $H(s)$ for a Notch filter, as shown in ch.32,p.17 of "The Scientist and Engineer's Guide to DSP" and see if similar reasoning on the integrals can be done.

Let's use the following filter (figure above for illustration only, I use different poles and zeros here) :
$$H(s) = \frac{s^2+1}{(s-(-1+i))(s-(-1-i))}$$
This filter has 2 poles and 2 zeros :

Zeros : $z_1,z_2 =\pm i$
Poles : $p_1,p_2 =-1 \pm i$

Let's find $h(t)$ and see why the integral would indeed go to 0 or $+\infty$ for these values of zeros and poles, respectively.
If it makes sense, this tool gives the following inverse Laplace transform for $H(s)$ :
$$h(t) = \delta(t) - 2e^{-t}\cos(t) u(t) + e^{-t}\sin(t) u(t)$$

Poles : for $s=p_1$ or $p_2$ in the Laplace transform, the exponentials of h(t) get canceled and remain the Fourier transform of some sinusoid which is indeed infinite at $\omega = \pm 1$ (I'm not discussing the $\delta(t)$ but I suppose it won't change this result).

Zeros : for $s=z_1$ or $z_2$ in the Laplace transform, the result is 0 if real part and imaginary of the Laplace transform are 0. Real part is :

$$\int_{0}^{+\infty} (\delta(t) - 2e^{-t}\cos(t)+e^{-t}\sin(t))\cos(t)dt$$
$$=\int_{0}^{+\infty} \delta(t)\cos(t)dt +  \int_{0}^{+\infty} (- 2e^{-t}\cos(t)+e^{-t}\sin(t))\cos(t)dt$$
with
$$\int_{0}^{+\infty} (- 2e^{-t}\cos(t)+e^{-t}\sin(t))\cos(t)dt = -1$$
Imaginary part is :
$$\int_{0}^{+\infty} \delta(t)\sin(t)dt +  \int_{0}^{+\infty} (- 2e^{-t}\cos(t)+e^{-t}\sin(t))\sin(t)dt$$
with
$$\int_{0}^{+\infty} (- 2e^{-t}\cos(t)+e^{-t}\sin(t))\sin(t)dt = 0$$
Questions

If the inverse Laplace transform is correct, how to handle $\int_{0}^{+\infty} \delta(t)\cos(t)dt$ and $\int_{0}^{+\infty} \delta(t)\sin(t)dt$ to show that $H(s)$ is indeed 0 at $z_1$ and $z_2$ ?
If all of this is correct, what does it (physically) mean for an impulse response to have a dirac in its expression ? I thought impulse response of most physical systems was only a combination of decaying exponentials and sinusoids ?

For you first question you can use the following
$$
\int_{-\infty}^{\infty} \delta (t-a)\,f(t)\,dt = f(a),
$$
with $f(t)$ any function. In your case those integrals would thus yield the values one and zero respectively.
For your second question I will only considering linear time invariant systems. In that case the impulse response of such a system can only contains a Dirac delta function if the transfer function of that system has a numerator of the same order as the denominator. Namely, any transfer function of the form
$$
G(s) = \frac{b_n\,s^n + b_{n-1}\,s^{n-1} + \cdots + b_1\,s + b_0}{s^n + a_{n-1}\,s^{n-1} + \cdots + a_1\,s + a_0},
$$
with $b_n \neq 0$ can also be written as
$$
G(s) = b_n + \frac{b'_{n-1}\,s^{n-1} + \cdots + b'_1\,s + b'_0}{s^n + a_{n-1}\,s^{n-1} + \cdots + a_1\,s + a_0},
$$
with $b'_k = b_k - b_n\,a_k$. The inverse Laplace transform of the constant $b_n$ would contribute a Dirac delta term. For the remaining part of the transfer function one could use partial fraction expansion to show that it can't contribute a Dirac delta term.
If a physical system would have a numerator of the same order as the denominator then it would require that the output of the system is directly affected by the input. An example of such physical system might be some electrical motor where you input a voltage and measure the angular position with some voltage leakage from the input signal to the output. However, most physical systems have numerator of a lower order as the denominator. It is more likely that you might encounter equal order numerators and denominators in digital filters (though, those would be z-domain and not s-domain, but roughly the same argument holds) such as notch filters. Those filters are however often used in series with physical systems, so their combined transfer function would also have a lower order numerator.


How to Subtract One Signal from Another? I am capturing some body internal signals.

I am using 2 identical mics. One captures the in body signal and one open air.

I want to remove/delete/subtract the open air signal (noise for me) from  the main signal.

Please suggest an efficient algorithm.

Under the assumption that your background noise is statistically independent of the signal you want to recover and that your second microphone picks up virtually nothing from the body noise, you have the following model for your measurements:
$$ \begin{eqnarray} M_1(t) &=& L\{N(t)\} + B(t) \\ M_2(t) &=& N(t) \end{eqnarray}$$

Here, $M_1,M_2$ are the signals picked up by the microphone, $B(t)$ is the signal you want to recover, $N(t)$ is the environmental noise signal and $L$ is a linear time invariant system that maps the recorded environmental noise to what is picked up by the second microphone.

With $\langle\langle \cdot, \cdot \rangle\rangle$ as cross correlation, we also have the non-correlation constraint $\langle\langle N,B \rangle\rangle=0$ and $\langle\langle L\{N\},B \rangle\rangle=0$ as a consequence of statistical independence of $N$ and $B$.

These prerequisites guarantee a cross correlation matching pursuit to converge against $B(t)$:

Normalise $M_2$ with the standard inner product: $M_2(t)\leftarrow \frac{M_2}{\sqrt{\langle M_2(t),M_2(t) \rangle}}$
Find $c(t):=\langle\langle M_2,M_1\rangle\rangle$
Find the time associated with the maximum of the squared cross correlation: $\tau=\mathrm{argmax}\left( c^2(t) \right)$
Update $M_1$ to remove the contribution of the environmental noise: $M_1(t)\leftarrow M_1(t)-\langle M_2(t-\tau),M_2 \rangle\cdot M_2(t-\tau)$
If result not good enough, go back to 2.

In step 4 you may also use $c(\tau)$ instead of the inner product if the cross correlation uses the same inner product. Stopping the matching pursuit is a question of experience. Some of the assumptions you need for this algorithm are only approximately true, like the perfect uncorrelatedness of the environmental noise and the body signal. A good way to check if your matching pursuit runs into a non-productive state of signal modifications is to observe the natural square-norm of $M_1$. If it stops decreasing in step 4, terminate the loop.

After the algorithm ran, $M_1(t)$ will contain only your body signal $B(t)$.


Different results for different orders of estimating AR model using Yule-Walker equations I'm trying to use MATLAB to estimate the AR parameters to the following filter:

$$H(z) = \frac{1}{1-0.5z^{-1}+0.25z^{-2} -0.25z^{-4}}$$
As I can see, the process at the output of this filter depends on the previous 4 samples. So, i figured that i would need 4th order model of Yule-Walker equations.

I run a MATLAB code solving Yule-Walker equations for 4,8,16,64 order and when I compared the 4 options, I got that the higher the order, the higher the MSE error. How can i explain this fact?

Here is a plot I got when compared $\lvert S(w)-S_{\rm estimated}(w)\rvert^{2}$.

Your true AR model is order 4. If an AR(16) model did a good job, then 12 of its parameters should be zero. Because of noise effects / overfitting, those 12 parameters are NOT zero and, more than likely, the other matching parameters will be more inaccurate because of this. Hence the error in an AR(16) will be larger than that for an AR(4) fit.


Contradiction while using the convolution sum for a non-LTI system In a recent quiz, we were given the following problem:

  The cascaded LTI systems $\mathcal{T}_1$ and $\mathcal{T}_2$ respectively have impulse responses $h_1[n]=\delta[n+3]$ and $h_2[n]=\delta[5−n]$. What is the output when the input is $x[n]=n$, i.e., find $y[n]=\mathcal{T}_2\{\mathcal{T}_1\{x[n]\}\}$.

I am primarily confused about the "LTI-ness" of the system and at a contradiction; given the second stage to be LTI (assuming that for a system to be LTI as a whole, all its stages/subsystems must also be LTI). Thus, we can say that a sub-system with  $h \lbrack n \rbrack = \delta \lbrack 5 - n\rbrack$, i.e. $y \lbrack n \rbrack = x \lbrack 5 - n\rbrack$ should be LTI, which already seems to be false. To prove that, I considered the following convolution sum: $$\tilde{y} \lbrack n \rbrack = \displaystyle{\sum_{k = -\infty}^{\infty} x \lbrack k \rbrack \delta \lbrack 5 - n + k \rbrack} = x \lbrack n-5 \rbrack \ne x \lbrack 5-n \rbrack = y \lbrack n \rbrack$$

And this is a contradiction since, after the convolution, the output doesn't match the original output with which we started.

Also, I realized that interestingly, the convolution always gives a result that corresponds to an LTI system (as in the case above too: $y[n]=x[n−5]$ is LTI whereas $y[n] = x[5-n]$ is not).

Also, since $\delta[n]$ is even, any non-LTI system's impulse response, for example, $\delta[1−n]$, will equal $\delta[n−1]$ which corresponds the impulse response of an LTI system. This explains why I am getting LTI characteristics after the convolution. This is interesting too since in either way the impulse response implies delaying the signal by 1 (in the current example).

And so did the impulse response $\delta[5-n]$ stated in the quiz question just qualitatively imply a delay of 5, and is, technically, not the exact description of the underlying system?

So, in summary, I have the following two doubts:

Is the convolution sum only true for the output relation of an LTI system? If so, then can it be used to prove non-LTI-ness of a system in the same manner as above?
Is the quiz question incorrect?

Any help will be greatly appreciated!

You've some troubles about the LTI property and the impulse response. First of all let me give you the answer to your quiz question.

Given two LTI systems $\mathcal{T}_1$ and $\mathcal{T}_2$ with impulse responses :

$$h_1[n] = \delta[n+3]$$
$$h_2[n] = \delta[5-n] = \delta[n-5]$$

you can compute the output of the cascade $$y[n] = \mathcal{T}_2 \{ \mathcal{T}_1 \{ x[n] \} \} $$ for the input signal $x[n] = n$ as :

$$ y_1[n] = \mathcal{T}_1 \{ x[n] \} \implies y_1[n] = n \star \delta[n+3] = (n+3) $$

and

$$ y_2[n] = \mathcal{T}_2 \{ y_1[n] \} \implies y_2[n] = (n+3) \star \delta[n-5] = (n-2) $$

so your output is $$y[n] = n-2$$

Note that $h_2[n] = \delta[5-n]$ and $h_2[n] = \delta[n-5]$ are the same systems, and yield the same output.

Note again that both systems are declared to be LTI and are given their respective impulse responses $h_1[n]$ and $h_2[n]$ , therefore their cascade will also be LTI and furthermore the output does not depend on the ordering of the systems; i.e.

$$y[n] = \mathcal{T}_2 \{ \mathcal{T}_1 \{ x[n] \} \} = \mathcal{T}_1 \{ \mathcal{T}_2 \{ x[n] \} \}  $$

Now coming to your own questions:

1- You cannot use convolution sum to deduce LTI property. You cannot form a convolution sum if the system is not already LTI.

2- The question seems ok:

Note: your claim that $ h[n] = \delta[5-n] $ , signifies a non-LTI system is wrong...


Why do I have frequency leakage in DFT after zero padding if frequency resolution is fine? Let's consider this example:

Fs=1000;
Ns=500;
t=0:1/Fs:(Ns-1)*1/Fs;
f1=10;
f2=400;
x=5+5*sin(2*pi*f1*t)+2*sin(2*pi*f2*t);
X=fft(x);

In this scenario, frequency resolution is 2, and all frequency components are captured correctly.
However, if I do this:

  X=fft(x,1000);

frequency resolution is 1, but there is a spectral leakage. Similar effect is seen here.
It seems to me, that Fourier transforms of both windows (one with length 500, and one with length 1000) have zeros at the frequencies that are presented in the signal, so I don't see why leakage will happen?

This phenomenon has nothing to do with spectral leakage. What you are observing is the effect of zero padding. Given a number of samples $N$, there is a maximum possible frequency resolution $\Delta f$ that can be achieved:

$$\Delta f=\frac{f_s}{N} $$

In your case $\Delta f$ is exactly $2\;\mathrm{Hz}$. If you zero-pad your signal, there is no extra information to retrieve - you will only decrease the frequency spacing.

In the example above, when you increase $N$ to $1000$, you get a frequency spacing of $1\;\mathrm{Hz}$.  All extra observed samples are merely an interpolation, done by the window function ($\mathrm{sinc}$ in your case). You will start observing side-lobes of the window spectrum.  Since you implicitly multiplied your signal by a rectangular window, this will result in the convolution of the spectrum of your signal (two Dirac's + DC) with the $\mathrm{sinc}$ function.

Another way to look at it is to imagine that DFT is basically a filter bank, consisting of shifted $\mathrm{sinc}$ functions. Those are aligned in such a way, that peak of one is where zeroes of all the remaining ones are present. If you start looking in between those zeros, you will start taking those samples. Here is an example plot of such $\mathrm{sinc}$ filter bank.

Let's imagine that frequency corresponding to the blue filter is present. That will yield the amplitude in a corresponding bin. All the remaining frequencies are not present (orange and yellow), thus you multiply those $\mathrm{sinc}$'s by $0$ and get nothing in the bins. In the case of zero padding, that will be no longer the case. Samples of that blue $\mathrm{sinc}$ will fall in the intermediate bins and will be sinc-interpolated.

Here is what happens for $N=1000$ and $N=10000$:

And a zoomed part:

Things to notice:

For $N=500$, there is no leakage whatsoever. There are perfect spikes, representing each of your frequencies and DC offset.
We can also observe the FFT noise at the very bottom.
For $N=10000$ the shape of $\mathrm{sinc}$ function is clearly visible.

And obviously the code to reproduce the results:

Fs=1000;
Ns=500;
Ns2=1000;
Ns3=10000;
t=0:1/Fs:(Ns-1)*1/Fs;
f1=10;
f2=400;
x=5+5*sin(2*pi*f1*t)+2*sin(2*pi*f2*t);

X1 = abs(fft(x))/length(x);
X2 = abs(fft(x, Ns2))/Ns;
X3 = abs(fft(x, Ns3))/Ns;

F1 = 0:Fs/Ns:Fs-Fs/Ns;
F2 = 0:Fs/Ns2:Fs-Fs/Ns2;
F3 = 0:Fs/Ns3:Fs-Fs/Ns3;

plot(F1, 20*log10(X1))
hold on
plot(F2, 20*log10(X2))
plot(F3, 20*log10(X3))
xlim([0, Fs/2])
grid on
legend({'N=500', 'N=1000', 'N=10000'})


OFDM training symbol format For OFDM synchronisation, we are using two training symbols according to the Schmidl and Cox method of frequency synchronization and timing synchronization. The paper says:

The first OFDM training symbol has only even numbered subcarriers applying a PN sequence.The result is two identical half symbols in time domain each consisting of Nc/2 samples each.
The second training symbol consists of even numbered subcarriers that are differentially modulated with even numbered subcarriers of first training symbol using a PN sequence. The odd numbered subcarriers of second training symbol can be used for data, pilot or reference symbols.

I have two questions on this.

In OFDM, every subcarrier corresponds to one symbol which is a complex value. So what does it mean by saying  even numbered subcarriers are modulated in a symbol using a PN sequence? How does that ensure that it results in two identical halves of a time domain symbol.

what is differential modulation using a PN sequence. How is it achieved?

To answer your first question, what they mean is that the first training symbol only encodes data on the even-numbered subcarriers. The other subcarriers are set to zero. That is, the frequency-domain,

$$
X[k] =
\begin{cases}
s_k, &k \text{ mod } 2 = 0 \\
0, &\text{otherwise}
\end{cases}
$$

The symbols to encode on the even-numbered subcarriers $s_k$ are chosen from a pseudorandom noise (PN) sequence. They assert that this results in a time-domain symbol that has two identical half-symbols in it (i.e. the first half of the symbol period is equal to the second half). Why is that?

Recall one of the properties of the discrete-time Fourier transform: if you zero-stuff a signal in the time-domain (by inserting zeros between each pair of samples), then in the frequency domain, the spectrum is repeated periodically. You might have seen this before in a description of discrete-time interpolation: insert zeros, then add a lowpass filter to eliminate the spectrum duplicates.

The property is dual in the sense that it works the other way also. When assigning symbols to subcarriers, you're working in the frequency domain. An inverse DFT is used to generate the time-domain signal for transmission. If you set all of the odd-numbered subcarriers to zero, then you're generating a frequency-domain signal that is zero-stuffed by a factor of two. When you inverse transform that to the time domain, you get the same effect: the resulting symbol consists of two periodic half-symbol waveforms.

I haven't read the paper in detail, but I'm assuming their synchronization technique takes this redundancy into account in some way.


Is Gabor uncertainty a feature of the Fourier transform, or of nature? By Gabor uncertainty, I mean the principle of uncertainty as applied to signals — with the result that you can't have arbitrary time and frequency localization.

By way of background to my question, I wrote about this recently in a blog post. In the last paragraph, I stated:

  Signals do not have arbitrarily precise time and frequency
  localization. It doesn’t matter how you compute a spectrum, if you
  want time information, you must pay for it with frequency information.

The point here is that, as I understand it, this uncertainty is an irreducible property of signals. It is not a 'feature' (or bug) in the Fourier transform or mathematics generally. And you can't get around it by using non-Fourier methods like matching pursuit.

I'm now unsure this is really true. Someone wrote to me recently:

  ...the uncertainty is there because of the math, not nature. A
  'signal' is a mathematical construct... The timing information of a
  signal is in the phase, and for a short-lived pure tone you can
  exactly recover the timing and frequency from the temporal phase.

I guess I buy this. But when we're measuring unknown, real, broadband, noisy signals does it still hold?

Time/frequency uncertainty is a property of the meanings or definitions attached to the terms "time" localization and "frequency" estimation.  It also alludes to some apriori knowledge, or lack of it, about the characteristics of the signal in question.


Detection of lines in a point cloud What are the best solutions to detect line in a point cloud? Comparison being made with and between Hough Transform, Radon Transform, RANSAC (see wikipedia) and Brute-Force Search (see wikipedia).
Which one is the most robust against the dispersion of points in point cloud?
Note:
1- The question is about 3D point cloud not image.
2- Points in point cloud are randomly dispersed (completely sparse locations).
3- There is no information about the object (line) being explored in terms of preferred orientation, size, etc.
4- A tolerance needs to be considered around the candidate line.
Updates:
According to my experiments: RANSAC could miss easily some lines. It is good for quick detection of edges however the complexity of point dispersion could produce undesired outputs. Hough and Radon are very similar and I had not chance to try for 3D point cloud however they work well on 2D cases. There is a difficulty in extraction of segments of found lines. BFS is simply impractical for large data set.

It really depends, how will you measure the quality of the solutions? what are your requirements, real time, high accuracy? how large is the point cloud?

You mentioned valid, yet fancy signal processing based methods to address the problem.

Let me add three methods you have not mentioned that are classic and more statistical in nature: least squares, ridge regression and lasso.


Generic text recognition OCR is a quite an old subject in Image Processing & Vision. However, the most literature that I had come across 1, 2, 3 shows various matured techniques which I am sure will give some good results in document processing.

Unfortunately, most OCR begins with document and that's when its utility gets limited. If you look at more generic scenario when a person goes to airport figures out where the sign board is and "reads" the signs board; now if we want to get a robot to achieve something like this, how does one go about achieving this?

Is there a more generic "text recognition" systems and algorithms?

A short answer to your question is "NO". To my best knowledge, there is no such a smart end-to-end OCR system that can handle various OCR tasks. However, there are OCR products could reasonably handle some tasks under controlled environments.

You seem to be confused of different OCR related tasks and you erroneously believe these tasks are quite similar ( and thus something works well on document could also work well on reading signs ).

Here lists several levels of OCR tasks:

level 0: single character recognition
level 1: word recognition level
level 2: sentence recognition
level 3: paragraph recognition
level 4: document recognition
level 5: OCR in the wild

Depending on the source of characters, you may classify texts into 1) machine printed texts ; 2) handwritten texts; and 3) scene text. All these are OCR related tasks, and thus can all be called "OCR".

These tasks are clearly of different difficulties. Basically, the higher the level the more difficult it is. More importantly, all these tasks are NOT just about recognition. Sometimes pre- and post- processing are both important. For example, it is well known that the line segmentation problem in handwritten document recognition is very challenging when text lines are cursive, or touched, or mixed of texts from different directions, etc. If you cannot first segment text line correctly, how to recognize them? For another example, it is also well known that deciding the reading order of multi-column document mixed with pictures is challenging. Even if you recognize every sentence correctly, you may still make mistakes in organizing different paragraphs. In your interested scene text recognition tasks, text localization itself is challenging in the sense that texts can be any font in anywhere with any direction and perspective change.

For these reasons, a good end-to-end OCR engine is at least a mix of image processing, pattern recognition, and natural language processing, and thus this is a quite big system requiring knowledge from all these parts.

I put my comments of 'prior knowledge' below.
To discuss 'prior knowledge', it is better to clarify its definition.
- For recognition, it often means there is prior distribution to tell you how likely a character is expected to see. If you said you have no prior knowledge, then you believe it is equally likely to see all possible characters. Think about how many valid characters in all languages. This number is big, and the number of valid words in all languages is even bigger.
- For image processing, it often means you have (explicit or implicit) assumptions about a given image. Such assumptions could be 'one-column document', 'car plate', 'taken by iphone', etc. If you say I have no prior knowledge about where a car plate is, then this simply means you treat the possibility of seeing a car plate is equally likely for any image patch size and location. Think about how many possible image patches you will have.

It is clear that doing OCR without prior knowledge is the worst thing we can try. In other words, we could do better when we have more prior knowledge for our task. In practice, people are using various prior knowledge, but often in an implicit way, and thus they even donot notice it.

One typical example is training data. Because your training data is not the entire universe, any rules, distributions, characteristics etc you learnt from this dataset may not necessarily apply to the rest of the world. However, you are still use your trained models to do various tasks. Your prediction is based on the prior knowledge that you learnt from data.

Back to your 'car plate' OCR problem. We know a car plate is of a rectangular shape under some perspective transformation. We know it is very likely to locate at the lower part of a car. We know it cannot be to too thin nor too small. Such prior knowledge helps quickly locate a candidate plate while avoiding a lot false alarms. When you recognize a plate, there is more noticeable prior knowledge. If you know only digits are allowed in your plate, then you cut your solution space to the digit space. Even if you see something more similar to letter 'g' than digit '9', you can safely say "this must be '9'" because of your prior knowledge. For the same reason, if you know plate are valid, you can safely thrown away any candidate solution that is invalid. Even 'not that useful prior knowledge' helps. For example, if you know no plate in our world has more than N characters, you know a plate candidate you are recognizing is NOT a valid plate because its length is above N.

If you donot have any prior knowledge, what happens? For me, it is similar to ask me find some object named '%#$$' on an image, which I completely have no idea what it is. And this is an impossible task. However, if I have the knowledge that '%#$#' is the ciphertext of 'ball', then I can easily do this task.

I hope that I convince you that prior knowledge is useful and closely associated with each processing step, and we cannot do OCR without prior knowledge. Indeed, this is the exact reason why we need to combine models from different aspects for OCR.

Donot be fooled by the name of OCR "optical character recognition". It is not just character recognition: we are actually reading words, sentences, paragraphs from an image that is taken by someone using some particular device... and thus we need prior knowledge at many different levels.


Cross-correlation of maximum length sequence with noise Background

I am using a maximum length sequence to find the impulse response of a linear system. I want to calculate the expected SNR. I know that my signal amplitude is scaled by the length of the sequence after cross-correlation, but I am not sure what happens to noise.

Idea

Since the sign of noise is random, the sign of an element in a maximum length sequence doesn't matter: [1,1,-1] may as well be [1,1,1]. Cross-correlation with a maximum length sequence then is like using a moving average filter except the coefficients are all 1's instead of 1/n where n is length of the filter. A moving average filter reduces noise power by the length of filter, and scaling by the length of the filter increases noise power by the length of the filter squared. Therefore the noise power increases by the length of the filter.

Questions

What is the cross-correlation of a maximum length sequence and thermal noise?

What is the cross-correlation of a maximum length sequence with a discrete sine (e.g. EMI from a clock signal)?

I assume that this is a discrete-time problem where the maximum-length sequence is a pseudorandom sequence $x[k]$ of $\pm 1$ values and the noise is a sequence $n[k]$ of independent identically distributed (iid) zero-mean random variables with variance $\sigma^2$.  Then, $\sum_{k=0}^{N-1} x[k]n[k]$ is also a sum of $N$ iid random variables and its variance is $N\sigma^2$, that is, the noise power increases in proportion to the length of the filter, exactly as you have conjectured.  So, the crosscorrelation of the maximum length sequence and the noise is a zero-mean random variable with variance $N\sigma^2$ where $N$ is the length of the sequence being correlated against. If the noise is modeled as Gaussian, then crosscorrelation value is also a Gaussian random variable.


How to interpret the signal in frequency domain I'm having difficulties understanding some basics regarding FFT. I was hoping you could answer a few questions I have. These questions will probably seem silly to you, sorry for that. I'm a complete beginner in DSP.

Let's consider the following script:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, np.pi * 2, 256)
y = np.sin(x) + np.sin(2*x + 1)
Y = np.fft.fft(y)

freq = np.fft.fftfreq(x.shape[-1])
plt.plot(freq, Y.real)
plt.savefig('Y_real.png')
plt.close()

plt.plot(freq, Y.imag)
plt.savefig('Y_imag.png')
plt.close()

y2 = np.fft.ifft(Y)
plt.plot(x, y2.real)
plt.savefig('y2_real.png')
plt.close()

The signal in time domain (y2_real.png):

The signal in frequency domain, real part (Y_real.png):

Imaginary part (Y_imag.png):

I was expecting that the real part will represent amplitudes of two sinusoidal signals and the imaginary part will represent a phase shift of them. However the amplitude 109 doesn't make much sense to me, nor does the phase shift of +/-128.

Where did these numbers came from? Also I'm not sure whether everything is right with OX axises in the frequency domain. How exactly can I figure out which frequency represents one peak or another?

These questions will probably seem silly to you, sorry for that.

No question is "silly". My perception is that the questions are trivial and can be answered by looking things up in a textbook but if this is to get you out of a "difficult corner", here it goes:

  ...I was expecting that the real part will represent amplitudes of two sinusoidal signals and the imaginary part will represent a phase shift of them.

The real and imaginary part of the complex Discrete Fourier Transform are one number which carries many pieces of information. Specifically, the magnitude $|X[k]|$ of some $k$ complex sinusoid returns information about its amplitude.

$$|X[k]| = \sqrt{\mathcal{R}(X[k])^2 + \mathcal{I}(X[k])^2}$$

Or in more practical terms Y_amp = np.abs(np.fft.fft(y)) following the signals established in your example.

The "phase shift" is the angle $\angle X[k]$ of a given $k$ complex sinusoid.

$$\angle X[k] = \arctan\left(\frac{\mathcal{I}(X[k])}{\mathcal{R}(X[k])}\right)$$

Or in more practical terms Y_phs = np.angle(np.fft.fft(y)).

Where $\mathcal{I},\mathcal{R}$ denote the imaginary and real part of a complex sinusoid at discrete frequency $k$.

  However the amplitude 109 doesn't make much sense to me, nor does the phase shift of +/-128.

That's because of conventions used in the evaluation of the DFT integral. So, if you look at the documentation of the FFT, you will notice that there is a $\frac{1}{N}$ term (where $N$ is the length of the signal) that can be applied either on $x[n]$ before the transform, or on the $\sum$ at the inverse transform. In this case, the sum is evaluated without scaling in the forward transform and with scaling when performing the inverse transform. If you divide your amplitude spectrum values by 256, your strengths across the spectrum will be adjusted to the values you probably expect. Again, in practical terms that would be: Y_amp_normed = np.abs(np.fft.fft(y))/256.

And since we are on the topic of phase, sooner or later you are going to need unwrap to deal with phase wraping.

  How exactly can I figure out which frequency represents one peak or another?

$k$ denotes discrete frequency. To convert between discrete frequency and natural frequency use:

$$f = \frac{k}{N_{FFT}} \cdot Fs$$

For more information about this you can have a quick look at this link

Hope this helps.


Linear system stability criteria Suppose we have a closed-loop system $H(s)=\frac{A(s)}{1+A(s)f(s)}=\frac{A(s)}{1+T(s)}$. I've seen the stability of the system stated a couple ways:

If $H(s)$ has any poles in the RHP, then it is unstable since the system response contains a term of the form $e^{s_{RHP}t}$, where $S_{RHP}$ has a positive real part and so will grow without bound over time.
If $1+T(s)=0$ for any $s$ (or equivalently $T(s)=-1$ for some $s$) then the system is unstable  since there exists some $s$ such that the transfer function blows up.

What confuses me is how these definitions are supposed to be equivalent. In particular, the second definition seems to say that if $H(s)$ has any poles, then the system is unstable. Why is this the case? Take for example a simple low-pass filter, $H(s)=\frac{1}{1+RCs}$, which has a pole at $s=-1/RC$. This pole is in the LHP and so by (1) is stable, however (2) seems to consider it unstable.
What is the relationship between (1) and (2)? In what situations, if any, are they equivalent/does one imply the other?

Only the first criterion is correct when referring to input-output stability. The second criterion is just a way to compute the poles, as mentioned in a comment. All realizable continuous-time systems using lumped elements have poles, which doesn't contradict the stability requirement.
Poles in the right half-plane (RHP) make a causal system unstable. Note that if we let go of the causality requirement, RHP poles don't cause instability but they make the system non-causal.
Note that poles on the imaginary axis cause transients that don't die out, which is undesirable in almost all practical applications. So in practice you want all poles to lie on the open left half-plane, i.e., $\textrm{Re}\{s_{\infty,k}\}<0$, where $s_{\infty,k}$ is a pole of the system.


Is my solution to the DTFT of a delta function correct? I am newer to signal processing math and just figured out something cool (hopefully). I was trying to see how the DTFT of the delta function is 1, because thats what my book says. I could only find youtube tutorials that didnt really explain it mathmatically, so i tried to do it on my own. Is my process correct? See pic.
Thanks!
-Dom

You are absolutely correct.  For more cool, calculate the DTFT of $\delta \left(n - k \right)$.


There is an implant in me, is there a way to identify the frequency it is receiving signals on? There is an implant in me, is there a way to identify the frequency it is receiving signals on?
I don't want to share the evidence and my reasoning for believing there is an implant in me, you can either take my word for it or not.
I want to find the frequency it is receiving signals on, so far I have used dvb-t dongle to use it as an SDR, and I have manually scanned the entire frequency range provided by it, which is 24MHz to 1700MHz, which was very tiring, but I couldn't identify anything suspicious, it could be possible that as I was tuning, the frequency it was receiving signals on was not the tuned one.
I want to know if there a way to identify the frequency of the signal? How? And what equipment would I need?

It is very hard to detect a device that is passively receiving a signal, and even more so when it is designed to operate undetected. Some random comments and ideas:

In the UK one has to pay taxes on TV ownership. I read somewhere that they employed a technique to detect unlicensed TVs, which are otherwise passive receivers. I don't know how the technique works, but you may want to research it.

Scanning the spectrum around you is useless. There are hundreds of different signals around us all the time; it will be impossible to detect if any one of them is intended for your implant.

Any receiver must use metal components, so you could use a metal detector on yourself. Get a very small receiver (such as one of the tiny WiFi dongles) and verify that your metal detector can detect it on the other side of your head; it should then be able to detect any device inside your body.


What do we mean by Total Factor Space in Audio Processing? I am referring to the article for understanding Ivectors which says:

  Over recent years, i-vector-based framework has been proven to provide state-of-the-art performance in speaker verification. Each utterance is projected onto a total factor space and is represented by a low-dimensional feature vector.

I tried understanding on my own, but could not get anywhere.
Can someone please help me understand what does 'Total Factor Space' mean?

I read the paper, let me try to summarize what I understand.

1) For speaker identification tasks, systems need to be robust to  "session variability". This means that from one spoken utterance to another, the way a speaker pronounces the same passphrase & the environmental/background noise impacts adds a bit of variance to the captured audio signal. I.e. it is not exact.

2) Historically, speaker-verification techniques based on GMM-UBM->MAP are not robust to channel distortion (i.e distortion during the recording) or phonetic variability (i.e. differences in pronunciation) or other "nuisance factors"

3) A single space named Total Factor Space is constructed to model speaker and channel variability jointly. Joint factor analysis (JFA) is one technique which decomposes a speaker's "mean supervector"  (which is a linear combination of the set of eigenvoices) into speaker component S and channel component C

So instead of doing separate modeling of speaker and channel spaces, all factors (i.e. both speaker and channel) are combined into one space. The JFA model from the slideshare below looks like:

Have a deeper look at this slideshare. Good luck!


Erosion and dilation in images I don't know if I am posting in the correct place. I wanted to ask about dilation and erosion in bitmaps. I have as a homework to implement in a programming language these morphological operations but I fail in understanding how they work. I know that erosion will shrink an object and dilation will enlarge it.

For example I have
$$A=\matrix{0 &1 &0 &1 & 1\\
0 &1 &0 &1 & 0\\
0 &1 &1 &1 & 0\\
0 &1 &1 &0 & 0\\
0 &1 &1 &1 & 0},\;
se=\matrix{1 &0 &0\\
1 &1 &0\\
0 &0 &1}$$

What would be the eroded and dilated output and how to compute it step by step?

Here's a step-by-step procedure for erosion/dilation by hand:

Print out A and se on two sheets of paper
Place the se paper on every pixel of the A sheet in turn
At each position: Take the pixel values of A at the respective positions where se is 1. For the first top-left position, this would be 0,0,1,1 as I have tried to illustrate here:
For an erosion, the result for the current pixel is the logical AND of the values you just wrote down. For dilation it is logical OR. (Or min/max for grayvalue images)

Note: You will either get a result image that is smaller than A or you have to add "padding" pixels to A (typically 1 for erosion and 0 for dilation)


Adaptive volume normalization on a stream of PCM data I am writing a program that receives a continual stream of PCM data from various sources. Some sources have different volumes than others. The PCM data is always a live continuous stream of 16 not signed data single channel.

How can I dynamically make sure that all samples are not too low in volume or high in volume to effectively normalize its output realtime?

What you want is called automatic gain control or dynamic range compression. It can be implemented in a dynamics processor that takes the input audio, uses an envelope follower to analyze the volume of the input as it changes over time and modifies the audio so that the envelope of its output corresponds to the input envelope put through a non-linear function such as:

Adobe Audition Dynamics Processing

I don't know what is the best practice in how the envelope should be modified to match the target envelope, but usually there are separate parameters for attack and decay that control how fast the applied gain should be increased or decreased, respectively, to make the resultant envelope approach the target envelope. You could try the potentially unstable:

$$\text{gain}(t) = \text{target_envelope}(t)/\text{input_envelope}(t)$$

where $t$ is time. The resultant $\text{gain}(t)$ should be smoothed (lowpass filtered) before applying it, to get the wanted attack and decay behavior. The envelope detector can have additional controls for changing the (usually root mean square or peak) analysis window, and a lookahead time setting that enables the dynamics processor to respond to volume changes before they take place, by delaying the input audio after envelope detection. The root mean square calculation can be done by calculating the square of the input audio, lowpass filtering that and taking the square root of the result.

Adobe Audition Dynamics Processing attack, decay and lookahead controls

These parameters and the non-linear function should be tuned (by listening) for minimal audible artifacts such as pumping of the volume and low-frequency distortion.


Is this how the lowpass butterworth analog filter can be implemented? Is this how the lowpass butterworth analog filter can be implemented?

Yes but there are a lot of details; like merging R4 with a termination impedance and using pre-distortion to compensate for finite Q in the inductors.  I would suggest reading some wikipedia articles:
https://en.wikipedia.org/wiki/Analogue_filter
https://en.wikipedia.org/wiki/Network_synthesis_filters
https://en.wikipedia.org/wiki/Butterworth_filter
https://en.wikipedia.org/wiki/Electronic_filter_topology#Ladder_topologies
To get a feel for you want.
Then you can decide how far you want to go along this path.  Whether to learn how to do it yourself or find some prepacked program.  But, unless you  feel lucky today, proceed with caution and some knowledge.  There are a lot of in's and out's; pitfalls.
I have always found the field fascinating; but that's partially because it is complicated and fundamentally mathematical :)  Your mileage will vary.
BTW: this is a DSP forum not an analogue design forum.


Nyquist Theorem - Why unique frequencies upto Fs/2 and not Fs? f+Fs is start of Aliasing If any frequency, f, displays an alias at f + Fs, This shows that unique frequencies have a range of Fs. Why does Nyquist theorem say that actually there is only half of this with unique frequencies below Fs/2?

Has it got something to do with negative frequencies?

No clear answer yet only confused :( please help

Since you stated you are still confused, the following graphics may help further understand the other two good answers:

First be sure that you understand what the sampling process (which is multiplying an analog signal by a stream of impulses in the time domain) looks like in both time and frequency; a stream of impulses in the time domain is a stream of impulses in the frequency domain. Multiplication in the time domain is convolution in the frequency domain, so if we sample by multiplying by an impulse every T seconds, we convolve our waveform with a frequency impulse spaced every 1/T Hz:

An example of this is sampling a 3 Hz cosine (where a 3 Hz cosine is represented in the frequency domain by an impulse at + and - 3 Hz) with a sampling clock of 20 Hz as shown in the graphic below. The top spectrum represents the input (analog) spectrum of the cosine wave. The middle spectrum is the spectrum of our sampling process (an impulse every 20 Hz including DC), and the bottom spectrum is the digital spectrum resulting from convolving the previous two.

Note that the analog spectrum was replicated in this process; so everything from $-F_s/2$ to $+F_s/2$, what we call the "first Nyquist zone", is duplicated centered on every multiple of $F_s$. (Where $F_s$ is the sampling rate). Thus the unique spectrum is completely defined by what is between $-F_s/2$ to $+F_s/2$, but we could easily show instead the spectrum for 0 to $F_s$, so it is somewhat arbitrary what slice of the spectrum we choose.

Note for real signals, the negative frequency components are the complex conjugate of the positive frequency components, so in that case not unique: they are completely defined by the positive frequency components so in that case that is all we need to present to provide all relevant information (the spectrum from DC to $F_s/2$).

That is not the case for a complex signal, as explained by hotpaw2 and shown in the graphic below.


Eye Diagram of QPSK Modulated Signal I need to plot eye diagram for RRC filtered QPSK modulated signal. The specifications are given as follows

Sampling frequency : $1.024$ MHz

Symbol Rate : $256$ Kbps (Hence, Upsampling rate = $4$)

Roll off factor = $0.2$

So I wrote a matlab code, but the resolution of eye diagram plotted is very low. I used object commscope.eyediagram.

Increasing the sampling frequency gives a better resolution, but it also increases the number of taps in the RRC filter which would not be a viable option for hardware implementation.

Kindly suggest what changes I need to make to get a better resolution of eye diagram without increasing the sampling frequency?

Kindly suggest what changes I need to make to get a better resolution of eye diagram without increasing the sampling frequency?

There is none, since the resolution of your eye diagram is directly linked to the sampling rate of your observer. In your case, that observation rate is identical to the sampling rate of your signal-producing system.

However, in reality, you'd use an oscilloscope that has a sampling rate much higher than your transmitter.

You can simulate the same – just use a proper interpolator between your modulator's output and your eye diagram's input. Pay attention to the fact that an insufficiently designed interpolation filter will distort your signal – and that is exactly the reality of digital transmitters: the reconstruction filters after the DAC distort the signal, too, and so does the RF channel between a transmitter and the oscilloscope (and of course, the filters and active components inside an oscilloscope do that, too).


How to measure the distance between two circle-like-lines over 360 degree I am looking for a slim and with that quickly programmable image processing strategy to measure the distance between two circles over 360°. It is important to notice that the circles are not perfect, see Fig. 1.

Fig. 1: Sketch of the circles

The application of this, so you can imagine the question a bit better, is: I have a cylindrical specimen with an radius r1(alpha) over the 360°, lets call the angle alpha. In the sketch it is the black circle.
Now I would like to measure the thikness "delta r" over alpha of a slim layer on this specimen. In the sketch the thin layer is painted in blue.

I draw the circles in kind of paint style on purpose as in the real image I would draw the line where I think both rings would be and would then start with calculating/image processing. I know that it would be possible to detect the edge with some kind of edge detection algo. (eg. Canny), but I am looking for a slim and fast way of getting to my result.

The perfect result could look like this: A function which is an approximation of the black line over alpha and then another function for the blue line.

Here a slim idea I had: Drawing two perfect circles over both cirles and substract them -> delta r, but the made error might be big.

So my question (1) is:
How would you solve this to get delta r over alpha? Thinking outside the box is welcome of course!

Question (2) is:
I know that the circles I drew are - from a mathematical point of view - not described at each alpha -> so I guess an approximation of the circle is somehow needed. Any ideas how to approximate?

Edit 1: Here you can see a real pictures. For now I don't have the optical system set up, so I just took a photo with my smartphone. Assume the length of the picture in x direction as 40 mm so you have a reference.

Fig. 2: Picture of the real specimen.

As @A_A said, it's a relatively straightforward problem once you've got good images.

Assuming you do have that:

Use the Hough transform to find the two main circles in the image
From the center of the inner circle, for each alpha, get a segment (by rotating your image and extracting a triangle) which should contain both your blue and black line
within that triangle, use appropriate signal processing (a simple "average per triangle column" or something, followed by a low pass filter, followed by a threshold, for example) to detect the exact position of the blue and black line


Modulation of data bits into complex symbols in MATLAB In communication theory, a modulator takes the input sequence and and maps it to some complex numbers. In MATLAB the following commands can do that.

data = randi([0 1],96,1);
hModulator = comm.PSKModulator(16,'BitInput',true);
modData = step(hModulator, data);

My question, is regarding the last step. Why is a step function able to do this mapping. In particular this step function represents the step response.

Thanks for any advice.

The method step(H,X) used here has nothing to do with the step function. It modulates the input data X with the modulator object H, as described here. See also here.


Is there a way to separate 2 signals recorded on the same channel? I have 2 low frequency signals on the same line that look like this:

I would like to separate them so that I obtain the two signals that can be seen here:

The sampling is not perfectly intercalated, so it's not as simple as taking one point each in alternation.

I looked at ICA but my understanding is that I would need 2 recordings of the 2 mixed signals in order to separate them.

Any ideas?

You should look into a generalization of ICA called Independent Subspace Analysis (ISA). Practically it is Principal Component Analysis (PCA) followed by ICA.

There are threads on dsp.stackexchange that discuss this as well as github repos (and here) that provide various implementations.


When to use EKF and when Kalman Filter? I'm learning Kalman Filter for a week now. I just discovered that EKF (extended Kalman Filter) might be more appropriate for my case.

Le't suppose I'm applying KF/EKF for variometer (the device that tells planes and parachuters what's their vertical position and velocity). In my case I've generated some sample data: the first few seconds he (the parachuter for e.g.) is falling (the velocity is positive) then he is rising (velocity is negative).

As far as I can tell this system is linear. So should I use KF or EKF?

The answer is simple: if your system is linear, then a (regular) Kalman filter will do just fine. A very brief summary of the differences between the two:

The extended Kalman filter (EKF) is an extension that can be applied to nonlinear systems. The requirement of linear equations for the measurement and state-transition models is relaxed; instead, the models can be nonlinear and need only be differentiable.

The EKF works by transforming the nonlinear models at each time step into linearized systems of equations. In a single-variable model, you would do this using the current model value and its derivative; the generalization for multiple variables and equations is the Jacobian matrix. The linearized equations are then used in a similar manner to the standard Kalman filter.

As in many cases where you approximate a nonlinear system with a linear model, there are cases where the EKF will not perform well. If you have a bad initial guess of the underlying system's state, then you could get garbage out. In contrast to the standard Kalman filter for linear systems, the EKF is not proven to be optimal in any sense; it's merely an extension of the linear-system technique to a wider class of problems.


Convolution to do FIR with initial and final conditions! I need to apply an FIR filter on chunks of very long signals. I can use MATLAB's filter function as follows to preserve the state of the filter when a new chunk of data is available:

[y, state] = filter(b, 1, x, state);

where b includes the coefficients of my FIR filter. However, b is fairly long and even each chunk of data is long and it appears to me that conv performs much faster than filter for this case. The following is a simple benchmark to compare the execution time for these two functions:

xin = rand(1, 1e6);
b = rand(1, 20e3);
tic; y = filter(b, 1, xin); toc
Elapsed time is 9.533244 seconds.

tic; y = conv(b, xin); toc
Elapsed time is 0.634058 seconds.

Unfortunately I can't preserve the state of my filter using conv function. I already tried to overlap the chunks of data and mimic the "state" for the convolution but haven't been successful yet. All suggestions and advice are greatly appreciated!

Your benchmark is actually quite surprising! I would've expected filter to choose an appropriate implementation of convolution to filter the input data, but thinking about this, since filter needs to pass the state in a format that isn't dependent on the convolution implementation actually used, the whole process of generating the state might take longer than the actual filtering; note that this really is just speculation. (Assumption: filter with long input data uses FFT FIR implementation internally, but the state that filter returns are taps of a direct-form FIR implementation – you basically need to reverse quite a bit of the filtering you've done to get that)

  I already tried to overlap the chunks of data and mimic the "state" for the convolution but haven't been successful yet.

What you probably would like to have is a well-implemented FFT FIR! Since b doesn't seem to change, this would be a very efficient way of implementing the filtering operation.

The basic idea would be described as overlap-add filter (wikipedia), and the rough idea is this:

Pad your tap vector $b$ with zeros to length $N$. Get the $N$-point discrete Fourier transform of that using the FFT. Call the result $H$.
Take segments of length $L$ of the input signal, and pad them to $N$, transform as above.
point-wise multiply, since multiplication in discrete freq domain is equivalent to circular convolution in time domain.
inverse DFT (via IFFT)
By zero-padding enough, we now that a $L$ long part of the output signal can be generated by adding up the $N-L$ last samples of the last iteration to the start of the samples of this iteration.
save the $N-L$ last output samples for next iteration

personal addition: If you already know what filter you want, and you're actually after speed and don't care as much that you might need to leave Matlab to do the filtering (I personally cherish every moment that I don't have to use Matlab), look at this explanatory blog post that I've written, which actually does a lot of of band pass filters in parallel – you'd only need one of these, and you could choose a FFT FIR block to do the work; I've been doing 160 megasamples per second with these, on a "normal" PC. This is, as I'll frankly admit, based on GNU Radio and the VOLK acceleration library, because that's what I prefer for streaming samples DSP (which you seem to be doing in Matlab).


Is there an existing feedback elimination algorithm that can be implemented in iOS? The project I'm working on allows a user to speak through a PA system with their iphone as the microphone. Unfortunately, the existing audio units in iOS don't eliminate the feedback between the PA speakers and the iphone microphone and there is instant and debilitating feedback. Is there a way to overcome this issue using existing algorithms? If software feedback elimination is not an option, would it be possible to tune the sensitivity of the iphone microphone to be more directional?

I've never implemented this myself, but here are my (possibly incorrect) thoughts:

The general method of feedback elimination is to

Detect the frequencies that are feeding back
Tune a filter to eliminate those frequencies

Step 2. is relatively easy, especially if you are using iOS and have an audio unit that can do it for you. There are some things to watch out for -- because you need a fairly steep filter you need to make sure it is stable even at extreme settings -- but hopefully the folks at Apple have taken care of that stuff for you.

Step 1. is trickier -- how do you know something is feedback and not a genuine part of the audio? There are some advanced techniques related to echo cancellation and adaptive filtering for doing this, but I'm not aware of anything widely available as code or libraries. There might be some shortcuts for your specific use-case, e.g. If you are dealing with speech only, you might be able to simply track pitch, and if a single pitch is strong enough assume that it's a result of feedback, but even doing that quickly enough may not be trivial.

Note that this entire process usually needs to be done continuously/adaptively, at least to some extent to make sure movements of the mike and speaker don't cause new feedback.

There are several patents associated with this, mostly because it's an important aspect of designing hearing-aids. Perhaps some of them are expired and you use them:

https://patents.google.com/patent/US4783818A/en?q=method&q=of&q=and&q=means&q=for&q=adaptively&q=filtering&q=screeching+noise&q=caused&q=by&q=acoustic+feedback
https://patents.google.com/patent/US4677677A/en?q=acoustic+feedback&q=cancelation&page=1

(Note, wikipedia mentions another method based on pitch shifting. I'm not familiar with this, and it sounds pretty silly, but it might be worth investigating.)