Gradient of Total Variation of Magnitude of Complex Function for Denoising Say I have a complex function $f^*$ (e.g. a MRI image) that has a near piece-wise constant magnitude, but a non constant phase.

If I have an optimization problem to find $f^*$ and set up an objective function with a total variation term (e.g. for denoising or compressed sensing) it usually has the following form:

$$
obj_1(f) = \ldots + \text{TV}(f)
$$

However, since I assume that $f$ has a piece-wise constant magnitude, I think it might be better to use:

$$
  obj_2(f) = \ldots + \text{TV}(|f|)
$$

However, for a gradient based solver, one would have to know the gradient of obj2. The gradient for $obj_1(f)$ is: $\text{TV}'\left(TV(f)\right)$. What is the gradient of $obj_2(f)$?

Update:

Intuitively I would assume something like the following (since the phase has no influence on $obj_2$, leave the phase "untouched"):

$$
\text{TV}'\left(TV(|f|)\right)* e^{i \arg(f)}
$$

However, my knowledge in complex analysis is very limited and I am not sure if this makes sense.

The problem with $\left|f\right|$ is that since is not analytic the standard definition of complex derivative does not apply. A solution is to use Wirtinger derivatives:

http://en.wikipedia.org/wiki/Wirtinger_derivatives

A detailed account of Wirtinger calculus for signal processing problems is

http://arxiv.org/abs/0906.4835

Another (probably simpler) option is to treat the complex image as a two-channel (real,imag) image and use the definition of derivative for vector fields. This paper has a very clear explanation on how to do this:

Lee, H.-C.; Cok, D.R.; "Detecting boundaries in a vector field" (IEEE Transactions on Signal Processing, vol.39, no.5, pp.1181-1194, May 1991)


Linear Prediction of AR Process A discrete signal x is generated by the recursive process
$$
x_n = x_{n-1} - 0.2 x_{n-2} + w_n
$$

where $w_n$ is white noise with zero mean and unit variance. What is the optimum order of a linear predictor for this signal? What are the values of the prediction coefficients? What is the average power of the residual?

I would really appreciate if someone could help me with this question.

It's a past paper question not homework.

Thanks

From the definition of the process you know that

$$x_{n+1}=x_n-0.2x_{n-1}+w_{n+1}\tag{1}$$

Since $w_n$ is white you can't predict it, so the best linear predictor for the given process is the filter

$$P(z)=1-0.2z^{-1}\tag{2}$$

which is a first order filter. It estimates the future sample $x_{n+1}$ by computing

$$\hat{x}_{n+1}=x_n-0.2x_{n-1}\tag{3}$$

The residual error is equal to $w_n$, which has an average power of $1$.


How to detect Premature ventricular contractions (PVCs) in a ECG-signal I want to know how to detect a  Premature ventricular contraction (PVC) in a ECG-signal. A PVC is best seen from the V1-lead in an Electrocardiogram.

Figure 1: A) PVC coming from the left ventricle, B) PVC coming from the right ventricle
The definition of a PVC can be found here

I'm an audio guy and not a biomed guy, but if the ECG is decently periodic except for an occasional artifact or anomaly like a PVC, you can do "pitch detection" (estimating the period or fundamental frequency) of the waveform and tune a comb filter to null out the periodic component.  Then the anomalies will be the only thing left after subtracting or nulling out the periodic component.

Then you can perhaps have a dictionary of waveshapes for possible anolamies and cross-correlate (sliding correlation) each of the waveshapes against this remaining waveform.  That operation is, essentially, what we call a "matched filter", which is used to recognize "grains" of signals.


Isolate g(t) starting from the definition of power spectrum Starting from the following definition of power spectrum of a signal $g(t)$ in time:
$$S(f)=\int_{-\infty}^{+\infty}\rho(\tau)e^{-i2\pi f \tau} d\tau$$
where
$$\rho(\tau)=\int_{-\infty}^{+\infty}g(t)^{.}g(t+\tau)dt$$
Is it possible to isolate $g(t)$?

Is it possible to isolate g(t)?

No. At least not if you mean by "isolate", "can I calculate $g(t)$ from $S(f)$"
You can recover $g(t)$ from it's Fourier Transform $G(f)$. This is related to the power spectrum through
$$S(f) = |G(f)|^2$$
So you can get the magnitude of $G(f)$ but not the phase. Information is lost and you cannot recover $g(t)$ uniquely. There are many different time domain signals that have the same power spectrum.


Spectral estimation in Spectrum analyzer I understand that there are various techniques for spectral estimation of the signals.
Methods such as parameteric(AR,ARMA etc.), non-parameteric( periodogram , welch etc..).
But I fail to understand which one is used by spectrum analyzer.
I am unable to relate the theoretical methods with spectrum analyzer.
Can anyone explain please?

The methods you're mentioning are methods for estimating spectra from a (typically discrete) time-domain signal – what the "classical" spectrum analyzer (SA) does is not that!
The classical sweeping SA consists of four parts:

A tuneable oscillator
a mixer
a good filter
a power detector

Here, 1. and 2. are used to mix down different parts of the spectrum - either by sweeping (typical for completely analog SA) or stepping through a set of frequencies.
The output of the sweeping mixer is then filtered by 3., which effectively means 1.-3. are a method of getting a varying piece of RF passband down to a frequency where 4. can be used to measure how much power is in that piece of RF spectrum.
You'd use the instantaneous frequency of the oscillator as x-coordinate, and the power "seen" by the power detector as y-coordinate. That's your point!
Note how you're not looking at the whole spectrum at once. You'll find "probability of intercept" considerations in literature, i.e. how likely is it that a bursty signal is actually active when the SA "tunes by" and observes the frequency they're on.
There's nothing here that really deals with the time-domain signal. For all that matters, the signal simply "has" some energy at some frequency. No need to estimate the spectrum from some time signal – they're actually measuring the energy at a single frequency, then move on to the next.
If you want to know more, there's really many introductions to using a SA, and you'll find the same components (oscillator, mixer, filter, power detector) again, and how you can configure them (sweep speed/dwell time, amplification, filter bandwidths, sensitivities).
Now, I didn't use "classical" with quotation marks without reason: For the last ~20 years, SA typically do have the possibility to sample the mixed-down signal and analyze it to improve the resolution, to implement better filters, or to do signal analysis beyond estimating a spectrum. Typically, these measurement methods are pretty limited, even in upper-class SA; they are literally labeled "FFT" without specifying exactly what they do, but you can expect things to be Bartlett's method for as long as the SA dwells on one frequency. How the SA pieces together different spectral estimates from different tuning frequencies is something that I haven't even found in manuals. Armed with an SDR and time you can figure that one out, but in all honesty, that's rarely what you really need an SA for: SDRs and even oscilloscopes can give you a better spectral estimates (though I can tell you, as DSP guy, the nomenclature of some oscilloscope vendors describing their spectral estimation methods drives me nuts).


Are noncausal filters ever used in practice? Basically what it says in the title; I have just started reading about these things and find noncausal filters pretty interesting in concept, but also they do not seem like they would have any advantage worth sacrificing real-time processing. Since "I have just started reading about these things" I feel as if I should make sure, and would also be interested in hearing: does anyone know if there are noncausal filters commonly used in practice? Why are they preferred? Thanks.

Jeff Boucher.

Yes.

The problem with a system that operates in (near) real time is that you can't look into the future. One way you can deal with this is if you only need a finite amount of look-ahead is to put some delays and then delay the output so you're still causal.

However, many filtering problems have non-causality allowed, e.g. filtering a file on a disk (which occurs a lot. Audio or image or video files you download, time series such as finance data, or histories of systems). For example, you can collect a time series and smooth it with the rauch-tung-striebel filter rather than running a fixed lag smoother.


Where does homomorphic filtering stand in regards to DSP applications? I am studying my Oppenheim and Shaffer book, (New Edition), and the last chapter deals with something called homomorphic filtering. I have read the wiki and some other websites about it, but they do not really delve into details on applications. Moreover, I have almost never heard of it be brought up in discussions or recent papers, etc.

So what are some applications of this processing, and why/how is it so quiet out there in regards to it? Is it because it has been overtaken by more advanced techniques, or is it just too clunky or some other reason? Have you used it before, if so, for what purpose?

The cepstrum, or cepstral analysis, which is one form of homomorphic processing, has been used in music and speech analysis to estimate the deconvolution of source (exciter) pitch and vocal tract (and/or other) resonances; sometimes just used as a music/speech pitch detection/estimation method, but also useful for formant tracking as well.


Decoding Matrix Encoded Surround Sound Through Convolution I wish to create a finite impulse response filter to decode a dolby digital matrix encoded 2 channel signal into 5 channels. These filters would then be used in a realtime pipeline on a Linux machine using jack and jconvolver to do the convolution.

The equation is as follows:

The reference for the above equation is: http://en.wikipedia.org/wiki/Matrix_decoder

It would seem this should be able to be done but I haven't the faintest how. Any help, and in particular working code that could generate the filters (Octave or other open source preferred) would be much appreciated!

Suppose that we had two channels encoded in some form in two signals according to a matrix.  Recovering the original signals would be straightforward using the matrix inverse (or impossible, if the inverse did not exist).  The situation is no different here, except that the matrix is not square, so you cannot use the conventional matrix inverse.  You can, however, use a pseudoinverse or generalized inverse to reconstruct the original signals.

Other matrices which are not equal to an inverse may work also succeed, and may in fact sound better.  I don't know the psychoacoustic situation, but it may more sense to do some transform based on the position of speakers for example. Certainly, the inverse seems like a good place to start.

I rewrote your encoder table as a 2x5 input matrix $ E $, and used Wolfram Alpha to determine the inverse, $ D = E^{-1} $.  Here is the 5x2 result it provided, using $ i $ as the imaginary unit, with $ j = i $ and $ k = -i $.

$$
D = \frac{1}{18(46+6\sqrt{2})} \begin{bmatrix}
405 & -3(27 - 36\sqrt{2}) \\
-3(27 - 36\sqrt{2}) & 405 \\
6(18+27\sqrt{2}) & 6(18+27\sqrt{2}) \\
-i(27\sqrt{3} + 99\sqrt{6}) & i(63\sqrt{3} + 27\sqrt{6}) \\
-i(63\sqrt{3} + 27\sqrt{6}) & i(27\sqrt{3} + 99\sqrt{6})
\end{bmatrix}
$$

As one would expect, the decoding matrix shows left/right symmetry.

Now we can express the encoding/decoding operations as:
$$
\begin{bmatrix}
L_{encoded} \\
R_{encoded}
\end{bmatrix} = E \space{} \begin{bmatrix}
L \\
R \\
C \\
L_{rear} \\
R_{rear} \\
\end{bmatrix}
$$
$$
\begin{bmatrix}
L' \\
R' \\
C' \\
L_{rear}' \\
R_{rear}' \\
\end{bmatrix} = D \space{} \begin{bmatrix}
L_{encoded} \\
R_{encoded}
\end{bmatrix}
$$

Now you can easily obtain formulas to express the three front components as functions of the encoded left and right signals.  For the two rear components, you have an equation in the time domain with a complex component.

Looking at the left rear component, and converting constants from the matrix $ D $ to numeric form:
$$
L_{rear}'[n] = -i0.2949 L_{encoded}[n] + i0.1787 R_{encoded}[n]
$$

The Hilbert transform can be used to remove the imaginary unit, because this has the same effect on the frequency domain of performing a ±90º phase shift.  The Hilbert transform only needs to be applied once, since convolution ($ * $) is distributive and commutative and all other good things.
$$
L_{rear}'[n] = (-0.2949 L_{encoded}[n] + 0.1787 R_{encoded}[n]) * H[n]
$$

The operating principle, then, is that the encoding technique functions in part by encoding the left/right rear components in quadrature.

The Hilbert transform is infinitely long, so it must be truncated and expressed as either a FIR filter, or depending on your circumstances, a FFT-IFFT operation.  A causal Hilbert transform will delay the rear channels, so your other filters will need to be delayed by some amount to compensate.

For example, if $ k $ is the group delay of a causal approximation of a Hilbert filter $ H_k[n] $ then you delay must delay the three front channels by $ k $ as well:
$$
C'[n] = (0.3437L_{encoded}[n] + 0.3437R_{encoded}[n]) * \delta[n-k]
$$

In practice, all of the output signals depend on the two input signals.  That is, you take two input channels and produce one output channel five times.  Your JACK pipeline will need to accommodate this.


Why does the Fourier series not include phase information? From my understanding, the Fourier series is a way to describe an arbitrary continuous signal in terms of sinusoids of varying frequencies and amplitudes, as shown here in Equation 1.

Why are the phases of the sinusoids not considered?

(The following is the answer to the original version of the question asking about the Fourier Transform. But exactly the same argument also applies to Fourier series).

They are. Assume that $F(\omega)$ is the Fourier transform of the signal $f(t)$:

$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt$$

In general, $F(\omega)$ is a complex-valued function and can be written in polar form

$$F(\omega)=|F(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase of $F(\omega)$.
The signal $f(t)$ can be written in terms of $F(\omega)$ by virtue of the inverse Fourier transform:

$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{j\omega t}d\omega=
\frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|e^{j(\omega t+\phi(\omega))}d\omega\tag{1}$$

From (1) it can be seen that $f(t)$ is represented by complex sinusoids with varying amplitudes $|F(\omega)|$, phases $\phi(\omega)$, and frequencies $\omega$.


Ways to differentiate three respiratory signals I have three 1-min respiratory signals I derived from ECG signals. I am wondering if you guys have any suggestion on differentiating those three? The red line represents the experimental signals and the blue ones are the computational results. Clearly, the second one is more representable than the other ones but I am not sure how I can differentiate it programatically or mathematically. I am not really asking for a specific way to process it. More like specific features that can help differentiate different waves below.

Any suggestions??

You're talking about ECG-derived respiratory (EDR) signals, but you didn't detail how you obtained them. Please detail your methods. Maybe then it will be easier to understand what each of your signals mean.

By 'experimental signal', I take it that you mean the reference/gold standard annotation of when respiration occurs, am I correct? But you did not detail if this was the onset or end or peak of the breathing cycle, as well as what kind of measure it was (tidal volume, airflow, other measure?)

Edit #1:
As for features, frequency is a good starting point. EDR signals may not evolve the same way as the reference signals but they should evolve at about the same period. If the breathing is controlled during a significant amount of samples, you can compare the frequency of the whole signals. If not controlled, then take segments of more than 2 seconds. You can apply FFT to get the frequency.

Phase will also be of interest, but it will difficult to comment further without knowing the algorithm for calculation of the EDR as well as the reference signal.


Application of UKF on quaternions I'm trying to perform a state estimation on quaternions to predict the future orientation of a human head. The only sensor data I can obtain (from the AR headset) is the current orientation of the head, sampled at 200 Hz i.e. I don't have access to any gyroscope or acceloremeter data. After getting the state estimations, I would like to reuse the process/motion model (constant angular velocity) to make predictions further into the future, e.g. 20 ms to 100 ms.

Since the process model is nonlinear (due to quaternions), one option is to use an Unscented Kalman Filter (UKF). However, as discussed in this paper, quaternions cannot be directly used in the UKF and some conversions need to be made to obtain "quaternion sigma points" (Sec. 3.2 of the paper).

My question is, does it make sense at all to use this kind of method if I only have attitude measurements (quaternions) and no gyro or acceloremeter data? In this case, my state vector would be 7D: four quaternion and three angular velocity components. However, the paper (and most other works I encountered) always have some kind of gyro/acceloremeter measurements which makes me wonder if it's feasible to just have the attitude information and have the filter estimate the angular velocity.

Yes, it's perfectly possible.  All that you'll need is to model how you think the angular velocity components of the state will evolve.  Usually simple brownian (random) motion is enough, at least to start with.  If you know more about how the angular velocities are constrained, then you can include that in the model.

All that it means is there are no outputs from your signal model directly from those elements of the state.

For example, in this answer to another question the state is made up of the $x$ position and the $x$ velocity, but there is no velocity measurement.

$$
\begin{align*}
\mathbf{x}_{k}  & =\left(\begin{array}[c]{c}x_{k}\\
\dot{x}_{k}\end{array}
\right)
\end{align*}
$$

All that means is the output matrix $H$ zeros out that part of the state when producing the measurement.


In DSP - is it true that the only type of interfering wave occurs when frequencies match? Therefore, any other frequency, no matter the amplitude, won't interfere or degrade original wave in the DSP frequency domain systems?

This is a question where theory and practice will diverge somewhat.

Theory

I assume that by "frequencies" you mean sinusoids of varying frequency. In that case, yes, you're correct, for any amplitudes $a_1$ and $a_2$ and frequencies $f_1$ and $f_2$ (where $f_1 \ne f_2$), the signals

$$
x_1(t) = a_1e^{j2\pi f_1 t} \\
x_2(t) = a_2e^{j2\pi f_2 t}
$$

are orthogonal to one another. The exact meaning of orthogonal varies from context to context, but in this sense, you can take it to mean that you could isolate the two signals from one another. If you had infinite time and infinite resources (and perhaps the ability to bend physics to your will if you want to realize it in hardware), you could build a filter that would pass $x_1(t)$ and eliminate $x_2(t)$, or vice versa. In theory.

In our theoretical utopia, the relative amplitudes of $x_1(t)$ and $x_2(t)$ don't matter either. Your theoretical processor can have unlimited dynamic range and/or numeric precision, so having one really large and one really small sinusoid isn't a problem.

Practice

Practical systems are hamstrung by several limitations that can invalidate the above claims. Among them:

If $f_1$ is really close to $f_2$, it can be difficult to resolve them as two separate sinusoidal components. If you were trying to do so digitally using short-time Fourier analysis, for example, you need a longer observation time in order to properly separate the two frequencies in your spectral estimate.

Another similar way of looking at it: if $f_1$ and $f_2$ are close to one another, and you wanted to build a filter that would pass one and eliminate the other, you would need a filter with a very sharp transition band. It can be shown that, qualitatively, the width of a filter's transition band is inversely proportional to its overall delay. So, as your transition band gets smaller and smaller, the delay required to implement the filter gets larger and larger. Furthermore, such a sharp transition band would require higher- and higher-order filters, possibly beyond the point of being practically realizable.
Practical systems do not have unlimited dynamic range. There will be some finite ratio between the maximum and minimum amplitudes that a system can reliably measure. This can limit your ability to resolve sinusoids that have widely-varying amplitudes, even if their frequencies are far apart. A very strong signal can desensitize your observer to the point where you can't see smaller signals.

Try shining a bright flashlight and a small LED at your face at the same time. Even if they are comfortably spread apart in space, you'll have a hard time seeing the LED due to the bright flashlight; your visual system has limited dynamic range.


Can a linear reconstruction in compressive sensing perform well? I am trying to implement compressive sensing for grayscale 2D images, then reconstructing them using a multi-layer perceptron(MLP). It seems to perform well no matter how many layers I add or remove, even without activation functions which is basically a linear combination of the inputs. there must me something wrong with my implementation. this is what I'm doing:
Encoding

Load image (PIL)
normalize image by dividing pixel values by 255
partition image into $B \times B$ blocks(patches) to get a matrix of size $P \times B \times B$ where $P$ is the number of blocks.
Flatten the blocks to form a matrix of size $P \times B^2$ let's call it $X$
Take M measurements from each block by getting the dot product of $\Phi$, an $M \times B^2$ Gaussian matrix, and $X^T$, the matrix containing the flattened, normalized blocks. The result is the measurement matrix $Y$ of size $P \times M$

Decoding
From what i understand, this is a sparse approximation problem. meaning my images(or blocks) are supposed to be sparse in some basis $\Psi$ in order to get an acceptable construction. I assumed the basis to be 2D-DCT. Therefore for each block, i try to reconstruct $S_i$ which is a $k$-sparse vector representing the flattened block in frequency domain. So my optimization problem becomes as follows:
$$\min|S_i|_1   s.t.   |\Phi\Psi^{-1} S_i-Y_i|_2 < \mu$$
$\mu$ being a small quantity. Now this applies to per-image optimization methods. What i am trying to do is an MLP which learns a universal function $f(Y_i)$ for reconstruction. So after training, my decoding process is as follows:  $$\hat{S_i} = f(Y_i)$$ following the acquisition of $\hat{S_i}$ i apply $\Psi^{-1}$ to it in order to get $\hat{X_i}$. finally, $\hat{X_i} * 255$ gives me the flattened block pixel values. reshaping each $\hat{X_i}$ into a $B \times B$ matrix then arranging them gives me the full image $X$.
Training
My dataset consists of $50,000$, $256 \times 256$ images. I just use the first batch (1000 images) since adding more data does not affect my model at all somehow. I do the aforementioned encoding process to every image then concatenate and shuffle all the blocks (i.e instead of images, i now have a matrix of flattened shuffled blocks from multiple images). My model inputs are the measurements taken from these images, of size $1000*P \times M$. To get my training labels, i apply 2D-DCT to each unflattened block before taking the measurements (right after encoding step 3 above), then flatten the sparse labels $1000*P \times B^2$, let's call the matrix $S$. Loss is calculated by taking $MSE(\hat{S}, S)$. Number of layers/nodes-per-layer do not seem to affect the model accuracy at all.
My questions

A Linear(no activations) MLP with a single hidden layer performs the same or sometimes better than a Multi-layer model with ReLU activations. What might be the cause of this?
In the literature, to guarantee construction $X$ has to be sparse in some domain. But i tried removing all transformations from my model and just worked in spatial domain and the model performs exactly the same. shouldn't it struggle to reconstruct the signal?

Output example
A model with 1 hidden layer having $B^2$ nodes takes the input of shape $M$. then the intermediate output is given to the output layer having $B^2$ nodes. in short, $M \rightarrow B^2 \rightarrow B^2$, no activations. training $MSE$ is $0.002$. $10$ of the blocks are reserved for validation, validation loss is also $0.002$. $B = 8$, M = $ 0.25 * B^2$. Now trying the trained model on an unseen images like "lena.bmp" $(512 \times 512)$ gives the following results compared to OMP:

N.B. The same $\Phi$ used in training is used in testing.
Any help is appreciated.

This is a good idea, I was pondering upon awhile ago. I drop some thoughts.

A Linear(no activations) MLP with a single hidden layer performs the
same or sometimes better than a Multi-layer model with ReLU
activations. What might be the cause of this?

The dimension reduction transform, i.e. taking CS measurements is a linear operation,  a measurement is weighted linear combination of input "sparse" vector. Are you reconstructing the same image that is used for training the network? The results could be different if the training image and test image are drastically different.

In the literature, to guarantee construction X has to be sparse in
some domain. But i tried removing all transformations from my model
and just worked in spatial domain and the model performs exactly the
same. shouldn't it struggle to reconstruct the signal?

There is actually no need to add that. Transformation is like any other function that network can "learn". It does not matter if you map measurement to "sampling" domain or the transform domain.
Maybe you might be interested in seeing Baraniuk's work on similar thing as well:
https://patentimages.storage.googleapis.com/7d/cc/57/9f8ebb657d6fb7/US10985777.pdf


Convolution that outputs a unit impulse Im thinking whether any convolutional operation can output a unit impulse, an example to further explain:
where a convolution between system $h[n]$ and unknown system $g[n]$ would output $\delta[n]$.

conv(h[n], g[n]) = δ[n]

where

$$h[n] = \delta[n] + \delta[n-1]/2 + \delta[n-2]/4 + \delta[n-3]/8 + \delta[n-4]/16$$

In regards to Z-transform - it can be concluded that

$$H[z] = 1 + 1/2z^{-1} + 1/4z^{-2} + 1/8z^{-3} + 1/16z^{-4}$$

Moreover, it is know that in z domain conv becomes multiplication therefore:

$$H[z] G[z] = 1$$

And in conclusion $G[z] = 1/H[z]$

$$G[z] = \frac{1}{1 + 1/2z^{-1} + 1/4z^{-2} + 1/8z^{-3} + 1/16z^{-4}}$$

If the above is somehow correct, then when trying to achieve the inverse z-transform we will get:

$$g[n] = \delta[n] + 2/\delta[n-1] + 4/\delta[n-2] + 8/\delta[n-3] + 16/\delta[n-4]$$

We get $1/\delta[n]$, which does not bode well for anybody.

Im I in the wrong? Can it be possible?

This problem is related to deconvolution and equalization. You are basically undoing the effect of a filter by another filter, such that the total system has a flat response, i.e., has a unit impulse as its impulse response.

From $$(h\star g)[n]=\delta[n]\tag{1}$$ it follows that $$H(z)G(z)=1\tag{2}$$ must be satisfied. So the solution to the problem is
$$G(z)=\frac{1}{H(z)}\tag{3}$$

as you've already found out.

However, the problem with $(3)$ is that the corresponding filter might not be causal and stable. Note that the poles of $G(z)$ are the zeros of $H(z)$, and the pole locations determine whether $G(z)$ is causal and stable. The inverse filter $G(z)$ is stable if all its poles are inside the unit circle. I.e., all zeros of $H(z)$ must be inside the unit circle for the inverse filter $G(z)$ to be causal and stable. Such a filter $H(z)$ with all its zeros inside the unit circle is called minimum-phase.

As a simple example take $h[n]=\delta[n]+\frac12\delta[n-1]$. We have $H(z)=1+\frac12z^{-1}$, and the zero of $H(z)$ is $z_0=-\frac12$. We can conclude that the inverse filter $(3)$ is causal and stable, because its pole is inside the unit circle:

$$G(z)=\frac{1}{1+\frac12z^{-1}}\tag{4}$$

The inverse $\mathcal{Z}$-transform of $G(z)$ is

$$g[n]=\left(-\frac12\right)^nu[n]\tag{5}$$

It's straightforward to show that the sequences $h[n]$ and $g[n]$ indeed satisfy $(1)$.

In practice, we usually assume that an exact inversion $(3)$ is not possible, so we try to approximately invert a given system. And we usually allow a certain delay, so $(1)$ is changed to $$(h\star g)[n]\approx\delta[n-N]\tag{6}$$ with some appropriately chosen $N>0$.

As a final remark, you have to work on your $\mathcal{Z}$-transform skills. Your attempt at inverting $G(z)$ was not successful.


Stabilizing an unstable system and removing oscillations from step response I have a system with the transfer function given in this MATLAB code
num = [-2.341e07 -2.478e10 -1.393e12 -2.58e13 - 1.027e14 1.689e15 1.474e16 2.227e16];
den = [1 2474 1.451e06 2.622e08 1.544e10 6.874e11 2.656e13 5.294e14 4.09e15 8.13e15 3.771e15 1.122e13];
G2 = tf(num,den);

The system has 2 poles and 1 zero in the RHP. The poles are a complex conjugated pair.
ans =

   1.0e+03 *

  -1.7134 + 0.0000i
  -0.4939 + 0.0000i
  -0.2000 + 0.0000i
   0.0032 + 0.0424i
   0.0032 - 0.0424i
  -0.0375 + 0.0000i
  -0.0211 + 0.0000i
  -0.0117 + 0.0000i
  -0.0020 + 0.0000i
  -0.0007 + 0.0000i
  -0.0000 + 0.0000i

This obviously makes the system unstable, but how do I stabilize the system? I have read that you can cancel out poles in RHP by adding a zero in the exact same spot, however it was not recommended and I am not sure how you would do that when the pole is complex.
Looking at the root locus plot, however, shows me that if the unstable poles will move to the LHP with appropriate gain (the cyan and purple are the unstable poles):-

But the problem is, that if I add a gain high enough such that the unstable poles move to the left, then some of the poles in the LHP will move to the RHP in the process:-

Edit
I made it stable! I added this term $G_d = \frac{s+0.002}{0.002s+1}$ and also a gain of $Kp = 1.6 $ to the forward branch. The closed loop step response with unity feedback looks like this (the 1 at the top indicates a stable system):-

The next step is just to remove the oscillations and improve settling time.

The OP has a mistake in the numerator vector as there is an additional space between the minus sign and the 1.027e14 coefficient which will cause that to operate as an expression (-2.58e13 - 1.027e14) instead of negate the sign.
Had the system been correct, the resulting closed loop system with the compensator and gain chosen is underdamped. This can be resolved in this case by lowering the gain such that the closed loop poles are equidistant from the real and imaginary axis. Further details below.
The  resulting ringing is because the closed loop poles are too close to the $j\omega$ axis relative the the distance to the real axis, an underdamped system. The root locus shown by the OP shows all the possible closed loop pole locations depending on what loop gain is used, and thus as is with the current controller chosen the OP can lower the gain to select a condition with closed loop poles equidistant from both axis (which is an ideal trade of settling time and ringing). The general considerations for closed loop pole locations are summarized in the chart below.


How does mic combined with loudspeaker work? There is variety of loudspeakers on the market that have built-in microphone so that they can produce sound and perceive it at the same time (and do it surprisingly well).

How they managed to overcome audio feedback problems?

By using Acoustic Echo Cancellation. Since the signal that goes to the speaker is known you can subtract it out at the microphone, provided you know the exact transfer function from speaker to microphone. The determination of the transfer function is typically done through an adaptive filter. Similar techniques can be used to reduce reflections, reverberation and echoes from the room.

In addition, there are often multiple microphones using beam forming. The beams can be steered so that the pickup from the speaker is minimized (basically steering a null in the polar pattern to the speaker) and it can also be used to maximize pickup of the desired sound.

All of these things are really difficult, complicated and hard to get right. Cancellation and adaption often create artifacts that significantly reduce speech intelligibility. The quality of commercial implementation varies all over the place. In my name brand laptop, it's so bad that it's worse than useless and I simply turned it off. I have tried multiple USB speakerphones and also found a huge variability in quality. Finally I found one that actually works :-) but product recommendation are out of scope of this site. Let's just say that the one pictured is fairly decent.


how to simulate multipath fading? I want to simulate a multipath fading channel, different amplitude and delay is given. I want to use delay in padding style. I read that the different values for delay like [0 .6 3.1 ...] represent the delay between the first and last signal arrival? Giving the number of paths, How to simulate it?
regards

@MohammedFatehy If you put up -exactly- what you currently know and have, we can help you more.

Generally speaking, let us say you have a signal x[n]. And lets say your sampling rate is 1 Hz. So you take one sample every second.

Now you want to construct a channel for multipath. Right off the bat, your multipath channel is going to be an FIR filter btw. Let us say someone tells you, "your first echo is going to be 5 seconds away, and attenuates the amplitude by 50%. The second echo is going to be 8 seconds away, and attenuates the amplitude by 70%".

Now, your channel is going to be simply a vector, (call it channel), such that: channel =  [1 0 0 0 0 0.5 0 0 0.3]. (You have the 1 in the beginning, because that is the 'line of sight' co-efficient. That is, the signal without any echos).

Now, you simply filter your signal x[n] with channel, and you have a multipath response. Does that make sense?


What is the unit-power information signal? I am learning signal processing for my subject and have some problems with the unit-power information signal. I refer to this paper: https://arxiv.org/pdf/1906.03949.pdf
In this paper, they formulate the received signal as: $y=h_{sd} \sqrt{p}s+n$ with $h_{sd}$ is the channel, $p$ is the transmit power, $s$ is the unit-power information signal, and $n$ is the noise.
What is the unit-power information signal? As I understand, the unit-power signal is the signal that has the power of 1. However, how can I formulate it (e.g., using Matlab)? For example, with BPSK signal, the signal is formulated as:
$s_1(t) = A_c \cos(2\pi f_ct)$ and $s_0(t)=A_c \cos(2\pi f_ct + \pi)$ for bit 1 and 0, respectively.
With a given transmit power $p$, how can I turn $s_1(t)$ and $s_0(t)$ into unit-power signal?
Is there any suggestion or reference?

The power of a deterministic signal $s(t)$ is given by
$$P_s=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^Ts^2(t)dt\tag{1}$$
A unit power signal is obtained by normalizing $s(t)$ by the root of its power:
$$\hat{s}(t)=\frac{s(t)}{\sqrt{P_s}}\tag{2}$$
In the case of $s(t)=A\cos(2\pi f_0t)$ you obtain from $(1)$
$$P_s=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^TA^2\cos^2(2\pi f_0t)dt=\frac{A^2}{2}\tag{3}$$
So the corresponding unit power signal is
$$\hat{s}(t)=\frac{\sqrt{2}}{A}s(t)=\sqrt{2}\cos(2\pi f_0t)\tag{4}$$


What is the zero in this transfer function? What is the zero for the following transfer function, $-1/2$ or $-2$?

$$H(s) = \frac{2s +1}{(s + 3)(s + 2)}$$

This appears to give a zero of $-1/2$.

I can transform this into the standard form that includes the gain $K$:

$$H(s) = \frac 12 \left(\frac{s + 2}{(s + 3)(s + 2)}\right)$$

Where the leading half is the gain. The zero now appears to be $-2$. What is the conceptual error I am making? Should I take the gain out or leave it in the numerator?

UPDATE: I think this question is worthy of deletion given the nature of the error but I am unable to. Instead, see my answer below.

$$2s+1=2\left(s+\frac12\right)$$
That's all I can say.


Unexpected peaks in power density following downsampling and filtering I first applied a 100 Hz lowpass filter to my data, which was recorded at 30000 Hz:

import numpy as np
from scipy import signal as ss
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt

def butter_lowpass(low_cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_low_cutoff = low_cutoff / nyq
    b_low, a_low = butter(order, normal_low_cutoff, btype='lowpass', analog=False)
    return b_low, a_low

def butter_lowpass_filter(data, low_cutoff, fs, order=5):
    b_low, a_low = butter_lowpass(low_cutoff, fs, order=order)
    y_low = ss.lfilter(b_low, a_low, data)
    return y_low

# Filter requirements.
order = 5
fs = 30000      # sample rate, Hz
low_cutoff = 100 # desired cutoff frequency of the filter, Hz

print('Filtering data')
filtered_array = butter_lowpass_filter(array, low_cutoff, fs, order)

Then I down sampled my data from 30000 Hz to 250 Hz

updated_array=ss.decimate(filtered_array, 12, ftype = 'fir')
newarray = ss.decimate(updated_array,10, ftype = 'fir')

And then apply the decimated result to a spectrogram:

frequencies, time, Sxx = ss.spectrogram(newarray,sampling_rate, ss.get_window('hamming', bin_size), noverlap=0, nfft=sampling_rate*4)

However, the resulted plot showed some signal higher than 100Hz even though I already used a lowpass 100Hz filter earlier.

For reference, this is the plot if the filter was not sed before decimation:

Is there any explanation for this?

Thanks a lot!

2 comments:

First, scipy signal decimate already low-pass filters the signal. It first filters the signal and then subsamples. So there is no need to pre-filter your signal before decimation. In your case, the signal is filtered at 125 Hz, which is half the resulting sampling frequency. (decimate doc here)

Second, a Butterworth filter is not an ideal filter. This filter will have a 3dB attenuation at the cutoff frequency, no matter the order of the filter. So you are expected to have components remaining between 100 Hz and 125 Hz.

If you want nothing detectable above 100 Hz, then design a filter with stop-band at 100 Hz and sufficient attenuation. This will result in a cutoff frequency below 100 Hz.

Depending on your requirements, an elliptic or Chebyshev filter will result in a narrower transition than a Butterworth of the same order.


Pitch correction detection? What type of signal analysis (if any) might be performed on a recorded music vocal track to determine if pitch correction might have been applied to a performance?

someone recently pointed me to a paper that purports to detect and count duplicated fragments.  i presume that if the density of duplicated fragments exceeds some threshold, one might suspect that up-shifting may have been done to that piece of audio.  i don't see how it could detect down-shifting (how does it detect a removed fragment?).

if, by "pitch correction", you mean pitch quantizing or autotuning™ on a monophonic sound, then i suppose you might run a good tracking pitch detector (maybe not such an easy problem) and maybe a histogram on the detected pitches.  we would expected pitch corrected vocals to have, in the histogram of pitch, a higher frequency of occurrence for pitches more closely located adjacent to pitch semitones (if autotune was set to correct to equal-temperament) or against some note identified as the tonic and then other notes that are "just intoned" (frequency ratios of small integers) with that tonic.

pitch detection/estimation is a signal processing issue, but determining what the tonic is, might be sorta a musicological question.  also something could be pitch "corrected" to a pitch standard that is not A-440.  perhaps all of the semitone locations are scooted down a quartertone from A-440.

if one is not going for the "Cher effect", pitch correction can be very subtle when it is done well.  one need not clamp the pitch hard against the semitone targets.  one can cause the amount of pitch correction to "slew" a little, which makes it slower to correct.  one can separate the vibrato out of the pitch detector output (with a LPF applied to the pitch contour signal), apply pitch correction to the remaining pitch contour, then add the extracted vibrato back.  that would be very hard to detect.  maybe the singer just sang it well.

not an easy problem.  and not a well-solved problem, at least with the current state-of-the-art.


DC to phase relationship with Low pass filtered sine and square wave multiplication and it's pitfalls. IQ Demodulation From what I have read about signal (de)modulation I understand that multiplying two sine waves, both with the same frequency, gives a two frequency component sine wave. One at DC and one at $2\omega$. Using a low pass filter will let us remove the $2\omega$ component. Now DC amplitude is proportional to the phase difference between the two sine waves.

So:
$$
A \sin(\omega t) \sin(\omega t + \theta) = \frac{A}{2}\cos(\theta) + \frac{A}{2}\cos(2\omega t + \theta) \rightarrow {\rm LPF} \rightarrow \frac{A}{2}\cos(\theta)
$$
where $\omega$ is the frequency, LPF is Low Pass Filter, $A$ is the amplitude and $\theta$ is the phase difference between the sine waves.

I made a Matlab program to plot this, then I tried with a square wave multiplied with a sine wave. I can see that there is still a connection between DC and phase by doing a FFT as well as in the time domain. I would like to see the math for the relationship between DC and phase.

For sine $\times$ sine, I have read that sampling with a frequency four times the modulated signal frequency will make it possible to get the I and Q values by subtracting samples in the following order:

$$
I(t) = s(t) - s(t-2)
$$
$$
Q(t) = s(t-3) - s(t-1)
$$
where
$$
s(t) = A\sin(\omega t + \theta)
$$

This I do not quite understand yet and therefore leads me to be unable to determine if this also is valid for a square $\times$ sine. Can anyone try to explain this to me?

I will appreciate any feedback on how to effectively do IQ demodulation on a micro-controller using a square reference waveform.

Additional information:

I wish to use a micro-controller to sample the analog sine and square (reference) signal and do IQ demodulation. I almost only find information about sine $\times$ sine and therefore want to know if there are any pitfalls using a reference square wave that I am not seeing yet.

So if your sampling rate is four times the modulated signal frequency (and sampling synchronously) then $\omega= \pi/2$. So then
$$
\begin{array}
 II(t) &=& s(t) - s(t-2)\\
&=& A\sin(\pi/2 t + \theta) - A\sin(\pi/2 (t - 2) + \theta)\\
&=& 2A\sin(\pi/2 t + \theta)
\end{array}
$$
and
$$
\begin{array}
 QQ(t) &=& s(t-3) - s(t-1)\\
&=& A\sin(\pi/2 (t-3) + \theta) - A\sin(\pi/2 (t - 1) + \theta)\\
&=& A\sin(\pi/2 t + \theta - 3\pi/2) - A\sin(\pi/2 t + \theta - \pi/2)\\
&=& 2A \cos((\pi/2 t + \theta - 3\pi/2 + \pi/2 t + \theta - \pi/2)/2) \\
&&\sin((\pi/2 t + \theta - 3\pi/2 - (\pi/2 t + \theta - \pi/2))/2)\\
&=& 2A\cos(\pi/2 t + \theta)
\end{array}
$$

The trick will be to think about the square wave case, and see what this derivation tells you about that.


New capture card YUV 4:2:2 leads to different images I have just bought a new video capture card whose output format is YUV 4:2:2 (as my old one) and surprisingly when I capture images with both cards I get different results. Differences are very small (have to zoom on edges if you want to see something) but I would like to know where do they come from? I was wondering if it was the transformation:

YUV 422 -> YUV 444

That could change from one driver to another?

Here is the original image:

Here is the capture with my old card:

Here is the capture with my new card:

4:2:2 means the color information has half the horizontal sample rate compared to the luminance.  Depending on the input, this lower sample rate may require having a low pass filter in the color channels.  Different cards may use different filters in phase and frequency response, resulting in slightly different results near any sharp transitions.


Instantaneous phase calculation problems I am currently struggling with a problem when calculating the instantaneous phase of a wavefield using a 1D Hilbert transformation.

The scheme is as follows.

1D Hilbert transformation of wave field $U(x,z)$
$$
q_z(x,z) = H_z[U(x,z)] = \int_R U(x-\xi,z)\frac{d\xi}{\xi}
$$

Calculation of the instantaneous phase $\phi(x,z)$

$$
\phi(x,z) = arctan\left(\frac{q_z(x,z)}{U(x,z}\right),
$$
where $x$ and $z$ are the coordinates and $H_z$ is the 1D Hilbert transformation over the vertical direction.

My problem is, that $\phi(x,z)$ which is implemented using the intrinsic function $atan2()$ contains phase shifts [$\pi \rightarrow -\pi$].

Using the instantaneous phase for further calculations in which its gradient is derived hence leads to imaging artifacts as shown in the attached pictures:

TOP LEFT: input wavefield $U(x,z)$

TOP RIGHT: Hilbert transformed wavefield $q_z(x,z) = H_z[U(x,z)]$

BOTTOM LEFT: instantaneous phase $\phi(x,z)$

BOTTOM RIGHT: derived propagation angle with artifacts

.

My question: I am able to overcome this problem using Gauss filtering or Median filtering. Yet this does not wotrk properly for more complicated wavefields and higher frequencies. Is there another way to overcome those sharp transitions in the area of my wavefield? (not the noise around)

Thanks alot!

One way to overcome this problem is phase unwrapping, which is also implemented by the Matlab function unwrap. But I believe that the best option for obtaining the derivative of the phase is to avoid a direct computation, but instead start out from the complex function itself. To simplify notation, let's use $U$ (without arguments):

$$U=re^{j\phi}=a+ib$$
where $r$ is the magnitude and $\phi$ is the phase. The derivative of $U$ is given by

$$U'=r'e^{j\phi}+rj\phi'e^{j\phi}\\
\frac{U'}{U}=\frac{r'}{r}+j\phi'$$

Since $r'/r$ is real-valued, we get

$$\phi'=\text{Im}\left\{\frac{U'}{U}\right\}=\frac{ab'-a'b}{a^2+b^2}\tag{1}$$

where $a$ and $b$ are the real and imaginary parts of $U$, respectively. Using (1) instead of a direct derivative of the phase should solve the problem.


What is the significance of Butterworth filter in image processing? Filters like mean,median,gaussian,bilateral and adaptive forms of these filters are used commonly in image processing but not Butterworth filter.
The reasons  I could think of are

a) Butterworth filtering can be done only in frequency domain(AFAIK there is no spatial filtering equivalent for Butterworth filter, correct me if I am wrong)

b) One of the main advantages of Butterworth filter is since there is no sharp cutoff frequency, there won't be any ripples, but since Gaussian also has this property there is no need for Butterworth.

Are there any specific applications Butterworth filter offers advantage over other filters? Do Butterworth have any significance in image processing?

Some important differences:

In Image Processing (IP), there is no causality like in Signal Processing (SP), hence there is not a tradeoff between filter quality and sampling sequence.
In IP, the FIR versions of SP are preferred instead of the IIR version (which are rare as you pointed). A possible relevant cause for this is FIR are designed as linear phase, unlike IIR which cannot be linear phase. Linear phase in FIR means the delay is constant (some pixels or fractions of pixels) while nonlinear phase in IIR means distorsion (like a blur) over the image toward the axis in which the filter is applied.
In IP there is not frequency spectrum, you do not seek to have a cutoff at x Hz and a flat passband. In a Butterworth this is the main design focus.
A butterworth or any other SP filter IIR version can be readily be extended to IP. The implementation is not a limitation.
A FIR like the Gaussian and others, dont have ripples in their impulse response. Butterworth, Chebychev, eliptical and ALL other IIR filters have an infinite number of ripples! This is another way to see the distortion which are produced by them.
The butterworth is ripple-less in their frequency response (the spectrum of the inpulse response) which we already agreed it do not have a direct sense in IP.
IIR means infinite impulse response: hence if you apply a Butterworth,  one convolution affects all other axis pixels.
Nevertheless, there are frequency figures on IP. You can define and extend 2D noises distribution based on an spatial frequency, and thus define frequency cutoffs such as the SP case. But this normally solved with a FIR filter, without paying the cost of the distortion an IIR implies.

Which is the significance? Most of their good properties in SP are not applicable in IP, with a major demerit added.


Low pass equivalent of bandpass white noise Let
$$n(t) = Re\{z(t)e^{j2\pi f_ct}\},\ $$ where $z(t)$ is the low pass equivalent of white noise $n(t)$. We know the autocorrelation of $n(t)$ is
$$R_n(\tau) = \frac{N_0}{2} \delta(\tau),\ $$
$\delta(\tau)$ being the Dirac delta function and the PSD is $S_n(f) = N_0/2$.
From Proakis' Digital Communication book, I learned that the autocorrelation of $z(t)$, $R_z(\tau) = N_0\delta(\tau)$, and PSD is naturally, $N_0$.
Now consider
$$z(t) = x(t) + jy(t)$$
I also know that $$R_z(\tau) = R_x(\tau) + jR_y(\tau)$$ So what are $R_x(\tau)$ and $R_y(\tau)$? Proakis' book mentions $R_z(\tau) = R_x(\tau) = R_y(\tau)$. But that seems a bit off to me.
Since, $R_z(\tau)$ is real then shouldn't $R_y(\tau)$ be zero being as it is the imaginary term in a real quantity?
Reference: Proakis, Digital Communication, Third edition, Chapter 4: Characterization of communication signals and systems, section 4-1-4, Representation of Bandpass Stationary Stochastic Processes.

Note that for the complex noise envelope $z(t)=x(t)+jy(t)$, the autocorrelation $R_z(\tau)$ is defined by (cf. Eq. $(4.1.47)$ in Proakis)
$$R_z(\tau)=\frac12E\big\{z^*(t)z(t+\tau))\big\}=\frac12\big[R_x(\tau)+R_y(\tau)\big]+j\frac12\big[R_{xy}(\tau)-R_{yx}(\tau)\big]\tag{1}$$
As shown in the chapter you refer to, for the real-valued bandpass noise $n(t)$ to be stationary, the following must be true:
$$\begin{align}R_x(\tau)&=R_y(\tau)\tag{2}\\R_{xy}(\tau)&=-R_{yx}(\tau)\tag{3}\end{align}$$
Plugging $(2)$ and $(3)$ into $(1)$ we obtain
$$R_z(\tau)=R_x(\tau)+jR_{xy}(\tau)\tag{4}$$
Note that this expression is different from the one you suggested in your question.
If the power spectrum of $z(t)$ is even, as is the case for bandpass white noise, the autocorrelation function $R_z(\tau)$ must be real-valued, and, consequently, $R_{xy}(\tau)=0$. But note that we don't require $R_y(\tau)=0$.
As a final note, defining a complex envelope is only useful for bandpass processes, such as bandpass white noise, but not for ordinary white noise, because the latter is no bandpass process.


Deconvolution of Synthetic 1D Signals - How To? I convolved a square wave with a Gaussian wave using linear convolution. Can I get the original square wave back by deconvolving my output with the Gaussian function?

I took the FFT of both signals, divided and then took the IFFT to get back the square wave. The output I got looked like random noise. I thought this might be due to the fact that I am doing a division and some denominator values might be very low, hence causing this error. I tried to set up a threshold and obtained the result depicted below.

How can I improve the result?

Edit 1 : Adding the code used to generate Gaussian signal. Was using matlab.

fs      =  200;           % Sampling frequency

t1 = 1/fs : 1/fs : n1/fs; % Where n1 is the length of the square wave signal

s2   =  4*gaussmf(t1, [ 0.4 4.5 ]);   % Generating Gaussian signal

Edit 2   :  Adding code used to generate square wave. I am afraid, it's a clumsy code!

fs      =  200;                           % Sampling Frequency
t       =  1/fs : 1/fs : 3;
n       =  length(t);

sq(1 : round(n/3) )  = eps;
sq( (round(n/3) + 1) : 2*round(n/3) ) = 3;
sq( (2*round(n/3) + 1) : n ) = eps;

sq( (3*round(n/3) + 1) : 4*round(n/3) ) = 3;
sq( (4*round(n/3) + 1) : 5*round(n/3) ) = eps;
sq( (5*round(n/3) + 1) : 6*round(n/3) ) = 3;

sq( (6*round(n/3) + 1) : 7*round(n/3) ) = eps;
sq( (7*round(n/3) + 1) : 8*round(n/3) ) = 3;
sq( (8*round(n/3) + 1) : 9*round(n/3) ) = eps;

You cant't recover the original signal through deconvolution.

A Gaussian kernel is in essence a lowpass filter, i.e. it will remove information at higher frequencies from the signal. Once it's gone, it's gone and you can't recover it.

This problem shows up as "divide by zero" or "divide by a very small number", which then amplifies numerically noise of the original convolution to a massive amount.

In order to recover the original signal the system must have "sufficient signal to noise ratio in the bandwidth of interest". What exactly that is, depends heavily on your specific application.

You can de-convolve from a low-shelf or high-shelf filter, but not from a low-pass or anything with a stop band or lots of attenuation at specific frequencies.

You'd also lose some information at the beginning and end of the signal.


Why do butterworth filter sometimes make signals diverge? When using Butterworth filter in MATLAB, if I chose an order too high, sometimes the filtered signal would grow without bound to infinity, and if I decrease the order, the problem disappears, has anyone encountered the same problem, and why?

Basically, this is bound to happen with higher order filters and finite precision. The filters' poles are getting closer and closer to the unit circle for sharper transition bands. But poles very close to the unit circle lead to unstable filters for some level of precision.

Things you could do:

Increase the precision of your numeric type (e.g. double instead of single precision)
Factor the filter into a series of 2nd order sections

Matlab and Octave allow you to directly compute the poles (p), zeros (z) and the overall gain (g) instead of the filter coefficients b and a:

[z,p,g] = butter(...);

and the result can be converted to a SOS matrix (second order sections):

[sos,g] = zp2sos(z,p,g);

Each row in the sos matrix represents a biquad, a 2nd order IIR filter. Using this second order section factorization -- meaning: a couple of 2nd order IIR filter in series -- is less subspectible to rounding errors as far as I know. Actual filtering can be done with sosfilt instead of filter:

y = sosfilt(sos,x)*g;

And if you still have trouble getting it stable, you should probably replace butter with ellip and allow a bit of passband ripple and/or do multiple passes of the same lower order filter.


Nyquist criterion for zero ISI in the frequency domain Can someone explain in simple words /graphs the necessary and sufficient Nyquist criterion for zero ISI in frequency domain, namely:

The constant folded Fourier Transform.

The short answer is you will see that when you follow that symmetry, the resulting transform (the impulse response in time) will all have nulls spaced at T where T is the symbol duration which by definition is zero-ISI. Work through those transforms as well as ones without such symmetry and I believe it will become very clear to you what is going on. The simplest case to start is obvious: a rectangular function in one domain is a Sinc function in the other. If it is a rectangular function with width 1/T in frequency, it is a SInc in time with nulls at integer multiples of T.

Another very simple example is to then convolve the two rectangular fucntions in frequency to get a triangular fucntion which will have this "folded frequency" relationship. This is simply multiplying in time, so although we changed the Sinc response, we have not changed the location of the nulls.


Deep learning model for phone recognition - issues with dimension the model I seem to have some problems understandind how the model described in this paper has been designed

This is what is written about the model dimension..

  ...In these experiments we used one convolution ply, one poolingply
  and two fully connected hidden layers on the top. The fullyconnected
  layers had 1000 units in each. The convolution andpooling parameters
  were: pooling size of 6, shift size of 2,filtersize of 8, 150 feature
  maps for FWS..

So according to ^ does the model consist of

Input

Convolution

Pooling

Input being the 150 feature maps (each with shape (8,3)

Covolution being 1d as kernel size is 8

and pooling is with size 6 and stride 2.

What was expected of output would be a shape of (1,"number of filters), but what i get is (14,"number of filters)

Which I understand why i get, but I don't understand how the paper suggest this can give an output shape of (1,"number of filters")

when using 100 filters I get these outputs from each layer

convolution1d give me (33,100)

pooling (14,100)..

Why i expect the output to be 1 instead of 14

The model is supposed to recognise phones, it takes in a 50 frames (150 including deltas) as input, these being a context frame, meaning that these are used as support to detect one single frame... That usually why context windows are used.

Since next layer is fully connected it does not really matter what shape your pooling output would be. You have 14x100, you can rearrange them as 1x1400 as input for next layer, 1000 elements as output. The paper says that they select the size of the fully-connected layer of 1000 to be close to the output size of the pooling layer.


How to detect time of arousal with excitation signal in speech signal? I have several seconds of speech signal and I have to detect moments in time of arousal with excitation signal, which is creating by larynx. I'm looking for useful reference since several days. I know that a speech signal is a convolution of excitation signal with a filter response and I've read about cepstral analysis and linear prediction, but I still do not know how to take advantage of these methods for my puspose. I think that I have to separate these two signals in some way.

I found the fundamental frequency by create autocorelation function and reading maxima, which gave me the period of fundamental tone. But actualy I don't know how to use it.
I need to do it in matlab. Please help me.

Here is good example what I need.
http://www.iba.zut.edu.pl/dydaktyka/pliki/referat_js/image002.jpg
In the top plot these impulses are marked as short lines and this is what I need to do.

You are looking for the instants of glottal closure, or GCIs.  There are many freely Matlab algorithms to do such.  Check out, for example, the DYPSA algorithm.  A paper explaining how these algorithms generally work can be found here.


Analogous Library to OpenCV for Audio Processing / Analysis I understand OpenCV is the de facto library for programming image processing in C/C++; I'm wondering if there is a C or C++ library like that for audio processing. I basically want to filter raw waves from a microphone, and analyze them with some machine learning algorithms. But I may eventually also need:

Multiplatform audio capture and audio playback
DSP - Audio filters
Tone detection
Tonal property analysis
Tone synthesis
Recognition given some recognition corpus and model
Speech / music synthesis

Any advice would be appreciated.

Consider the following:

clam-project.org:

  CLAM (C++ Library for Audio and Music) is a full-fledged software framework for research and application development in the Audio and Music Domain. It offers a conceptual model as well as tools for the analysis, synthesis and processing of audio signals.

MARF:

  MARF is an open-source research platform and a collection of voice/sound/speech/text and natural language processing (NLP) algorithms written in Java and arranged into a modular and extensible framework facilitating addition of new algorithms. MARF can run distributedly over the network and may act as a library in applications or be used as a source for learning and extension.

aubio:

  aubio is a tool designed for the extraction of annotations from audio signals. Its features include segmenting a sound file before each of its attacks, performing pitch detection, tapping the beat and producing midi streams from live audio.


Audio silence on 1.0 I have audio files with silence on 1.0 (max amplitude). I attach an example of audio waveform. I don't know why it is on 1.0? How can I convert it to have silence on 0 (as usual)? I noted that at the beginning of the waveform it is starting from 0 and goes to 1.0 - any ideas why? The files are telephone recordings, they are sampled with 8kHz and coded with GSM. I don't have more information about recording procedure.

For transmission of audio by amplitude modulation of for example light, the audio can be made unipolar in order to make silence correspond to no transmission which for some receivers means no noise during silence. Just a guess but your audio may have gone through such processing. From the look of the waveform it may have worked like this:

Start with bias = 0. For each time step:

bias = bias + constant, where constant determines the slope seen during silent segments in your graph
If bias + input > 1, then set bias = 1 - input to avoid clipping.

It is not possible to exactly recover the original audio from the output of such processing, but you could use a similar method that keeps shifting the bias towards the direction that centers the audio, and changes the bias appropriately any time the audio would be clipped.


How is quadrature mixing done? I've a Gaussian modulated signal whose quadrature and in phase components are separated. I'm multiplying these components with quadrature and in phase components of a signal of 10MHz frequency, respectively. Then i added the two products.This sum is quadrature mixed data.Right?

If your signal's in-phase and quadrature components are $x_I(t)$ and $x_Q(t)$, respectively, then your quadrature mixed signal is

$$s(t)=x_I(t)\cos(\omega_ct)-x_Q(t)\sin(\omega_ct)\tag{1}$$

where $\omega_c$ is the carrier frequency in radians per second. You can also use a '$+$' sign in $(1)$, that's just a matter of convention.


How do I go about detecting whistles, pops and other sounds in live audio input? I've read many questions on SO, and frankly, each one of them is not describing any particular way to go about it. Some say "do FFT" and some say "zero crossing" etc. But I've only gone as far as understanding that the digital audio input consists of an array of amplitudes for a particular frequency and well.. Don't really know much beyond it.

Now I do know Nyquist theorem, frequency, amplitude, Fourier series etc, and that's from 2-3 years ago when I did that in my university program in some semester. But back then we really weren't taught the real world usage of Fourier much, and I didn't bother to dig up more on the subject beyond learning just enough to pass the subject. But now I'm gonna have to use all that stuff.

Here's a snapshot of the sounds I'm trying to detect:

Clearly the sounds have unique graphs. I just want to understand how to extract their particular unique characteristics for their unique spiky lines in the graph. Like what amplitude, frequency, etc. And for how much time - although that is trivial I guess.

I want a simple step by step, un-vague instruction list - I can google the terminology I won't understand.

Like perhaps this? --

Get input audio data
Plot spectogram
Get spectogram graph for the sound you want to detect in a noiseless environment
Study that graph - draw unique characteristics of that sound
Make some kind of function that can detect those particular characteristics in live audio feed, using the characteristics of the sound found in (4)
If a match is found, well - the job's done.
Polish the algorithm to weed out false negatives.

Where I am confused is - how do I go about 3,4,5. Particularly 3 and 4.
Because I'm making a software where in the user would be able to record any sound as a gesture input to be used later, I would want my system to be able to record any sound and match it against live audio feed to detect that sound gesture. Like say, coughing, or finger snapping, or whistling, or clapping. Or saying consonants - $A$ or $B$ or $X$ or $Z$. Just about any sound.

I was thinking of making the user record the sound they want to store as a gesture in a quite environment. And the user would make the sound only in between the quiet time padding; 3 secs in the beginning and in the end of the recording.

So say, during the first 3 secs, my system would establish that the current input is the normal quiet background sound. And then a sudden change in the graph would be the sound input starting. And when that stops, the recording would go on for another 3 secs, the trailing quiet time padding. This would be done manually by the user. It would then automatically store the characteristics of only that part during which the sudden change in the graph lasted - somewhere between the padding times.

And thus the characteristics of that part would be saved as that sound's gesture data, which would be used to detect that particular sound in live audio feed later on.

Problem is, that I'm thinking all this in plain English. I need to think in maths and physics, to be able to implement it efficiently in my code. I'm just so damn clueless about what to write and where to write it in my code - even with so many libraries and questions on SO at my disposal.

And sorry if this was long.

While I agree with one of the comments that using speech recognition techniques might be a good start, these sounds are different and I'm not aware of anyone having done any research into categorizing them (the paper cited by Nathan appears only to distinguish between speech and noise), so unless someone presents something to contradict that, you are going to have to invent your own technique, and that's going to take a lot of learning and a lot of work. The best I can do is get you started.

First off, don't expect anyone to be able to produce some magic formula. There's no equation to get from sound to what the sound is. Humans and computers must learn by incorporating data to make educated guesses about what sounds are. The reason people say in their answers "use FFT" or "use zero crossing" is because these are some of the basic DSP building blocks used in speech recognition and related algorithms. But the FFT and zero crossing rate are usually just first steps in creating a set of parameters describing the sound. These parameters are then analyzed statistically (not through some magic function) to determine which category they most likely belong to. Note that I said "most likely": even the best speech detection (and the human brain!) has a high error rate at this stage and further processing using knowledge of the language and grammar is used to correct these mistakes.

So, some parameters you could look for include:

zero-crossing rate
frequency centroid
envelope (this is actually a set of parameters, including, eg attack time)
spectral envelope
even/odd harmonic balance
voiciness
fundamental pitch

Once you have a set of parameters that you think will allow you to distinguish your sounds, you will need to use some statistical method to categorize them. Hidden Markov Model is often used in speech. You can also look into logistic regression, K-means, and I'm sure there are other choices, but I think HMM is tried and true.


How have common PRBS polynomials been chosen, and what are they? Why are most common or standard (is there actually a standard?) PRBS uneven in their number of state bits, e.g. PRBS7, PRBS15, PRBS23, if there are plenty of generator polynomials that have even, byte-multiple numbers of bits, like 8, 16 and 24?

More generally, how have the common sequences, denoted with the aforementioned names, actually been chosen? Do they have any favorable properties compared to other LFSR maximum length sequences, or was it arbitrary?

Finally, where can I actually find what common polynomials are used? There are precomputed lists of maximum length sequences, but no mention of what the common/standard ones are. If I instead search for "PRBS63" in Google, I can find mentions of that sequence, but not the actually used generator polynomials.

The number of maximal-length sequence generators increases rapidly with the degree and so tabulating all of them is impracticable for large degree. A standard table of all irreducible binary polynomials of degree 16 or less (and a few selected polynomials of degrees up to 34) takes up 17 pages in in an Appendix of Peterson and Weldon's Error-Correcting Codes, (2nd ed.), MIT Press 1972 and they are much more densely packed than the prolix notation used in the OP's reference.  As to which ones are in common use or which have desirable properties, it depends on what is meant by desirable; mathematically, they are all equivalent, practically, there might be good reason to choose one over the other. One preferable feature might be small number of nonzero coefficients (which reduces the number of XOR gates (or instructions in software implementations).  Another might be smaller overall complexity of implementation: the primitive polynomial of degree 8 used in the Reed-Solomon code implementation in the NASA system was carefully chosen to minimize the overall complexity of the decoder (and no, it is not the first one in the Peterson&Weldon table).


Is imadjust a Matlab-only library? If so, can it be replaced? I want to apply a homomorphic process to an image, and one of the steps involves the imadjust function. The source of this code is http://www.csse.uwa.edu.au/~pk/research/matlabfns/

One of the steps in the homomorphic process is the histogram truncation. This is where toward the end, imadjust is applied like this

newim = imadjust(im,gv,[0 1]);

Unfortunately Octave doesn't seem to have this function, so I'm guessing it's a Matlab Image Processing Toolkit function. Is there a way to replace it? I haven't managed to find one that uses just two arguments, in addition to the image.

I found this imadjust.m but it has 4 arguments (image, in, out, and gamma)
Not sure how it's related to the one I'm using which is 3 arguments (image, gamma?/interpolated value, and two values array).

The imadjust in Matlab includes 4 arguments as well. You always need to specify the shape of the curve (gamma), and map the values in I to new values in J such that values between low_in and high_in map to values between low_out and high_out.

You can try the link you found yourself, or try this version.


Get FIR Hilbert coefficients from Matlabs envelope() function How can I get the used Hilbert filter coefficients from the Matlab function envelope()?

filter_length = 10;
[up1,lo1] = envelope(input_data, filter_length,'analytic');

From the mathworks documentation of the function envelope():

  The filter is created by windowing an ideal brick-wall filter with a Kaiser window of length fl and shape parameter $\beta = 8$.

So without hacking the function you can't directly get the filter coefficients, but you can easily find them yourself by just doing what they do, i.e., windowing the ideal impulse response by a Kaiser window. The ideal impulse response of a Hilbert transformer is given by

$$h[n]=\begin{cases}\frac{2}{n\pi},& n\textrm{ odd}\\0,&n\textrm{ even}\end{cases}\tag{1}$$

You just have multiply these coefficients with the window of the desired length.


"River" detection in text Over on the TeX stackexchange, we have been discussing how to detect "rivers" in paragraphs in this question.

In this context, rivers are bands of white space that result from accidental alignment of interword spaces in the text. Since this can be quite distracting to a reader bad rivers are considered to be a symptom of poor typography. An example of text with rivers is this one, where there are two rivers flowing diagonally.

There is interest in detecting these rivers automatically, so that they can be avoided (probably by manual editing of the text). Raphink is making some progress at the TeX level (which only knows of glyph positions and bounding boxes), but I feel confident that the best way to detect rivers is with some image processing (since glyph shapes are very important and not available to TeX). I have tried various ways to extract the rivers from the above image, but my simple idea of applying a small amount of ellipsoidal blurring doesn't seem to be good enough. I also tried some Radon Hough transform based filtering, but I didn't get anywhere with those either. The rivers are very visible to the feature-detection circuits of the human eye/retina/brain and somehow I would think this could be translated to some kind of filtering operation, but I am not able to make it work. Any ideas?

To be specific, I'm looking for some operation that will detect the 2 rivers in the above image, but not have too many other false positive detections.

EDIT: endolith asked why I am pursuing a image-processing-based approach given that in TeX we have access to the glyph positions, spacings, etc, and it might be much faster and more reliable to use an algorithm that examines the actual text. My reason for doing things the other way is that the shape of the glyphs can affect how noticeable a river is, and at the text level it is very difficult to consider this shape (which depends on the font, on ligaturing, etc). For an example of how the shape of the glyphs can be important, consider the following two examples, where the difference between them is that I have replaced a few glyphs with others of almost the same width, so that a text-based analysis would consider them equally good/bad. Note, however, that the rivers in the first example are much worse than in the second.

I have thought about this some more, and think that the following should be fairly stable. Note that I have limited myself to morphological operations, because these should be available in any standard image processing library.

(1) Open image with a nPix-by-1 mask, where nPix is about the vertical distance between letters

#% read image
img = rgb2gray('http://i.stack.imgur.com/4ShOW.png');

%# threshold and open with a rectangle
%# that is roughly letter sized
bwImg = img > 200; %# threshold of 200 is better than 128

opImg = imopen(bwImg,ones(13,1));

(2) Open image with a 1-by-mPix mask to eliminate whatever is too narrow to be a river.

opImg = imopen(opImg,ones(1,5));

(3) Remove horizontal "rivers and lakes" that are due to space between paragraphs, or indentation. For this, we remove all rows that are all true, and open with the nPix-by-1 mask that we know will not affect the rivers we have found previously.

To remove lakes, we can use an opening mask that is slightly larger than nPix-by-nPix.

At this step, we can also throw out everything that is too small to be a real river, i.e. everything that covers less area than (nPix+2)*(mPix+2)*4 (that will give us ~3 lines). The +2 is there because we know that all objects are at least nPix in height, and mPix in width, and we want to go a little above that.

%# horizontal river: just look for rows that are all true
opImg(all(opImg,2),:) = false;
%# open with line spacing (nPix)
opImg = imopen(opImg,ones(13,1));

%# remove lakes with nPix+2
opImg = opImg & ~imopen(opImg,ones(15,15));

%# remove small fry
opImg = bwareaopen(opImg,7*15*4);

(4) If we're interested in not only the length, but also the width of the river, we can combine distance transform with skeleton.

   dt = bwdist(~opImg);
   sk = bwmorph(opImg,'skel',inf);
   %# prune the skeleton a bit to remove branches
   sk = bwmorph(sk,'spur',7);

   riversWithWidth = dt.*sk;

(colors correspond to width of river (though color bar is off by a factor of 2)

Now you can get the approximate length of the rivers by counting the number of pixels in each connected component, and the average width by averaging their pixel values.

Here's the exact same analysis applied to the second, "no-river" image:


Conceptual definition of entropy What is the conceptual difference between Kolmogorov-Sinai entropy, Shannon entropy and Boltzmann entropy? Are they interchangeable and mean the same? Where can I find a good lucid explanation about Kolmogorov entropy and why it is needed and its significance. Thank you.

Let me give you a quick answer, focused on ways you can improve and research for a deeper knowledge by yourself.

1 - Kolmogorov and Shannon  Entropies are related to measures of information and strongly related to one another. (See the wikipedia page for Kolmogorov complexity for further info.)

2 - Boltzmann entropy, which is an older concept, is a completely differenty and unrelated concept from the theory of Thermodynamics. (Again, wikiepdia...)

3 - The formal definition of kolomogorov entropy uses the Shannon entropy in the formula.

If you are interested in the information entropy concept, i suggest you read the Book by Cover &Thomas, which has an entire chapter on Kolmogorov Complexity.


stationary vs. undecimated wavelet transform I have a little bit confused on the difference between stationary wavelet transform and un-decimated wavelet transform.

So, can anyone tell me, if there is a difference between them?

From what I understand un-decimated wavelet transform (UDWT) involves the elimination of the decimation step to enable DWT to have translational invariance property. The stationary wavelet transform is one such algorithm which is also known as algorithme à trous. There are a number of other algorithms of UDWT a quick google scholar search can you show that but they all use similar ideas to the stationary wavelet transform.


What is the general form of a transfer function I repeatedly see two representations of the general transfer function in the literature. The first is the following which is factorization of the numerator and denominator polynomials:

$$\textbf{G}(s) = k \frac{ (s - a_{1})(s - a_{2})(s - a_{3})\ldots
                              }{ (s - b_{1})(s - b_{2})(s - b_{3})(s - b_{4})\ldots
                                }$$

The alternative, below, eludes me and I'll checked many sources but none seem to justify it. The difference is a lowly $s$ term in the denominator raised to a power.  Is this general form related the open-loop transfer function?

$$\textbf{G}(s) = k \frac{ (s - a_{1})(s - a_{2})(s - a_{3})\ldots
                              }{s^m (s - b_{1})(s - b_{2})(s - b_{3})(s - b_{4})\ldots
                                }$$

My question is what is the difference, if any in these two forms?

You may write the second equation on the first form by using $b_i = 0$ corresponding to those poles appearing at $s = 0$. Hence, the only difference is that in the second form you know that there are $m$ poles at $s=0$ while in the first form they may still be there, but one will have to check the values of $b_i$ to determine if they are there.

Note that a system with poles at $s=0$ is a marginally stable system and that the poles are typically introduced by integrators, which have a transfer function $1/s$.


Transform RGB image to *look like* Infrared Context: I'm trying to improve a pose estimation model so that it works better when my camera is in Infrared mode. Unfortunately I only have RGB images to train on.

I realize that you can't convert RGB to IR directly, but my hypothesis is that converting the RGB images to look more like IR, and then training on a dataset of combined RGB and IR images, will lead to better performance.

Are there any libraries that have tried to implement a function like this? I'm essentially looking for a function that something like this ("IR effect") - http://funny.pho.to/infrared/

Without very much additional info, I'd presume the red channel of your camera has the highest correlation with the infrared spectrum.

Since you only have RGB and no knowledge of how that was calculated from the color sensor pixels: Take the R channel, it's as good as it gets.

The "IR effect" is for artistic purposes only. I'd doubt it has anything to do with what you see in an IR picture.

Notice that I also think you haven't tried to properly understand the physics of your problem: While I'd presume that a visible light photograph is highly correlated with a near-IR photo, that's very likely not the case for mid- and long-wavelength IR and certainly not for far infrared.

So, define your requirements better. And: ML really depends on having data. A lot of it. you might start with wrong data, but as DNN really don't lend themselves to understanding of the weights, it's very questionable that anything that you do with "emulated" data transfers to "real" IR data well.


Minimum phase systems with pole at infinity If a system is given by a transfer function in the $z$ domain that has all poles and zeros inside the unit circle except for a factor of $z^{-1}$ in the denominator (pole at infinity), can it still be considered minimum phase?
If not, how would I create an all pass system to neutralize this pole and obtain a minimum phase representation of the system without also adding a zero at infinity (which i assume would also cause the system to not be min-phase?

The transfer function:

$$H(z) = \frac{1 - \frac 12 z^{-1}}{z^{-1}}$$

If I see it correctly, your transfer function is

$$H(z)=\frac{1-\frac{1}{2}z^{-1}}{z^{-1}}=z-\frac{1}{2}$$

This system is not a minimum-phase system because its pole is at infinity (as you've already pointed out). You could use a trivial all-pass filter with impulse response $h[n]=\delta[n-1]$, i.e. a simple delay by one sample to change your original (non-causal) transfer function to a stable and causal transfer function

$$H_m(z)=z^{-1}H(z)=1-\frac{1}{2}z^{-1}$$

which is also minimum-phase, because it has a pole at $z=0$ and a zero at $z=1/2$. This means that all the system's poles and zeros are inside the unit circle and the minimum-phase requirement is met.


Why eigen values and poles of a system are equivalent? In control systems engineering, the stability of a system (modeled in the form of Transfer Function) is determined by the poles of the system in the right or left hand sides.
When the model is represented using State Space approach, the eigen values of the (A) state matrix are equivalent to the poles in the Transfer Function approach.

why are both eigen values and poles equivalent? why do they give the same information?

Let's consider a discrete-time state space model (the derivation for a coninuous-time system is completely analogous):

$$\begin{align}\mathbf{q}[n+1]&=\mathbf{Aq}[n]+\mathbf{b}x[n]\\
y[n]&=\mathbf{c}^T\mathbf{q}[n]+dx[n]\tag{1}
\end{align}$$

where $x[n]$ is the input, $y[n]$ is the output, and $\mathbf{q}[n]$ is the state vector. Taking the $\mathcal{Z}$-transform of (1) gives

$$\begin{align}z\mathbf{Q}(z)&=\mathbf{AQ}(z)+\mathbf{b}X(z)\\
Y(z)&=\mathbf{c}^T\mathbf{Q}(z)+dX(z)\tag{2}
\end{align}$$

From the top equation of (2) we can express $\mathbf{Q}(z)$ in terms of $\mathbf{A}$, $\mathbf{b}$, and $X(z)$:

$$\mathbf{Q}(z)=(z\mathbf{I}-\mathbf{A})^{-1}\mathbf{b}X(z)\tag{3}$$

where $\mathbf{I}$ is the identity matrix of appropriate dimension. Plugging (3) into the second equation in (2) gives

$$Y(z)=\mathbf{c}^T(z\mathbf{I}-\mathbf{A})^{-1}\mathbf{b}X(z)+dX(z)\tag{4}$$

from which follows the expression for the transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\mathbf{c}^T(z\mathbf{I}-\mathbf{A})^{-1}\mathbf{b}+d\tag{5}$$

The matrix $(z\mathbf{I}-\mathbf{A})^{-1}$ can be written as

$$(z\mathbf{I}-\mathbf{A})^{-1}=\frac{1}{\text{det}(z\mathbf{I}-\mathbf{A})}\mathbf{B}\tag{6}$$

where $\mathbf{B}$ is the adjugate matrix of $(z\mathbf{I}-\mathbf{A})$ whose elements are polynomials in $z$.

Plugging (6) into the expression for the transfer function (5) shows that the denominator of $H(z)$ is given by $\text{det}(z\mathbf{I}-\mathbf{A})$, which is the characteristic polynomial of $\mathbf{A}$ the roots of which are the eigenvalues of $\mathbf{A}$. Since the roots of the denominator of $H(z)$ are the system's poles, the eigenvalues of $\mathbf{A}$ and the poles of $H(z)$ are equivalent.


Equations for particle filter The particle filter is based on the state and observation model equations

$x_{t+1}=f_t(x_t, v_t)$

$y_t=h_t(x_t, u_t)$

The idea is to randomly generate some particles then propagate them through the equations, resample, and normalize to get an estimate of the state $x_t$.

My question is: we need to know the state and observation model equations so why do we need the particle filter? I hope I am missing something because it seems like if this knowledge is required then why bother with particle filter, can't we just do ML estimation?

why do we need the particle filter?

When your model is non-linear there is no analytical solution to obtain the PDF of the state. The Kalman filter/smoother can't be applied in the sense that it is not optimal. The particle filter is one method to obtain an estimate of state PDF in a sequential setup.

  can't we just do ML estimation

Dependent on the model you may find that even an ML estimator can't be found without going through some approximations. But the particle filter is also used in a Bayesian contexts where other estimators are relevant.


Why do the OFDM training (and payload) symbols have silent sub-carriers (aka virtual carriers)? I am examining an OFDM signal.

Training symbol:

The first two OFDM symbols are training symbols, they are 2PSK modulated (within the subcarriers).
In the first OFDM symbol the $(0+18N)^{\rm th}$ subcarriers are silent ($N$ is integer).
In the second OFDM symbol the $(9+18N)^{\rm th}$ subcarriers are silent.

Payload:

In the first payload symbol, i.e. in the $3^{\rm rd}$ OFDM symbol, the $(6 + 18N)^{\rm th}$ subcarriers are silent.
In the second payload symbol, i.e. in the $4^{\rm th}$ OFDM symbol, the $(6 +9 + 18N)^{\rm th}$ subcarriers are silent.
In the third payload symbol, i.e. in the $5^{\rm th}$ OFDM symbol, the $(3 + 18N)^{\rm th}$ subcarriers are silent.
In the $4^{\rm th}$ payload symbol, i.e. in the $6^{\rm th}$ OFDM symbol, the $(3 +9 + 18N)^{\rm th}$ subcarriers are silent.
The payload symbols are QPSK modulated.

Then this pattern seems to repeat.

Pattern:

So to summarise, the offsets are : $0, 9, 6, 6 + 9, 3, 3 + 9, 0, 9, \ldots$

Questions:

How does this make sense ?
What purpose does such a unique silencing of subcarriers serve ?
I suspect it has to do with equalizations , but why is silencing useful ?
And why this odd way of silencing the subcarriers ?
Is this (way of silencing subcarriers) somehow connected to the Schmidl and Cox method ?

EDIT: Google gives this on "OFDM silent subcarriers"

EDIT2: The signal has been emitted by a MIMO system.

EDIT3: This is the Google search for "OFDM virtual carriers".

EDIT4: This is a 1024 carrier OFDM signal.

From the first to the sixth OFDM-Symbol you have spread virtual carriers equally across [0,..,51] with a distance of 3. With this way you have spread your 'energy' loss over a wide spectrum. In difference to having a single carrier serving as virtual carrier all the time you can transmit information over every 'data' carrier and maybe have a constant mean energy level across all subcarriers.

You could use this virtual carries to detect carrier frequency offset (CFO) since you know the pattern. In the worst case you have to store 6 OFDM-Symbols to find the beginning of your pattern.

If you have an OFDM-Symbol with 64 subcarriers you can set 6 carriers at each end as subcarrier and use the other 52 with this virtual carrier pattern to transmit information with equalization.

I think this virtual carrier pattern is not connected to the Schmidl & Cox synchronsiation pattern.


Is it possible exctract sinusoids from non periodic signal? Digital signal UT1-UTC is not periodic but is including many sinusoids (periodic elements in IERS nomenclature) that are not multiples of some fundamental. For example tidal sinusoids are not multiples of yearly seasonal sinusoid because lunar month is not submultiple of the year. Then, is it possible these sinusoids be extracted from UT1-UTC by Discrete Fourier Transform?

Fourier's theorem say almost any (non-pathological) waveform, periodic or not, can be decomposed into sinusoids (or complex exponentials).  Whether, or how well, those sinusoids correspond to any underlying pseudo-periodic phenomena or not is another issue.

Note that with a DFT, you may need to interpolate a periodic signal between the DFT result bin sinusoids.


Designing a feature vector for discriminating between different sonic waveforms Consider the 4 following waveform signals:

signal1 = [4.1880   11.5270   55.8612  110.6730  146.2967  145.4113  104.1815   60.1679   14.3949  -53.7558  -72.6384  -88.0250  -98.4607]

signal2 = [ -39.6966   44.8127   95.0896  145.4097  144.5878   95.5007   61.0545   47.2886   28.1277  -40.9720  -53.6246  -63.4821  -72.3029  -74.8313  -77.8124]

signal3 = [-225.5691 -192.8458 -145.6628  151.0867  172.0412  172.5784  164.2109  160.3817  164.5383  171.8134  178.3905  180.8994  172.1375  149.2719  -51.9629 -148.1348 -150.4799 -149.6639]

signal4 = [ -218.5187 -211.5729 -181.9739 -144.8084  127.3846  162.9755  162.6934  150.8078  145.8774  156.9846  175.2362  188.0448  189.4951  175.9540  147.4631  -89.9513 -154.1579 -151.0851]

We notice that signal 1 and 2 look similar and that signal 3 and 4 look similar.

I am looking for an algorithm that take as input n signals and divide them into m groups, where the signals within each group are similar.

First step in such an algorithm would usually be to calculate a feature vector for each signal: $ \mathbf{F}_i $.

As an example we might define the feature vector to be: [width, max, max-min].
In which case we would get the following feature vectors:

$ \mathbf{F}_1 = [13,146,245] $

$ \mathbf{F}_2 = [15,145,223] $

$ \mathbf{F}_3 = [18,181,406] $

$ \mathbf{F}_4 = [18,189,408] $

The important thing when deciding on a feature vector is that similar signals get feature vectors that are close to each other and dissimilar signals gets feature vectors that are far apart.

In the example above we get:

$| \mathbf{F}_2 - \mathbf{F}_1 | = 22.1, |\mathbf{F}_3 - \mathbf{F}_1 | = 164.8 $

We could therefore conclude that signal 2 is much more similar to signal 1 than signal 3 is.

As feature vector I might also use the terms from the discrete cosine transform of the signal. The figure below shows the signals along with approximation of the signals by the first 5 terms from the discrete cosine transform:

The discrete cosine coefficients in this case are:

F1 = [94.2496  192.7706 -211.4520  -82.8782   11.2105]

F2 = [61.7481  230.3206 -114.1549 -129.2138  -65.9035]

F3 = [182.2051   18.6785 -595.3893  -46.9929 -236.3459]

F4 = [148.6924 -171.0035 -593.7428   16.8965 -223.8754]

In this case we get:

$| \mathbf{F}_2 - \mathbf{F}_1 | = 141.5, |\mathbf{F}_3 - \mathbf{F}_1 | = 498.0 $

The ratio is not quite as large as for the simpler feature vector above.
Does this mean that the simpler feature vector is better?

So far I have only shown 2 waveforms. The plot below show some other waveforms that would be the input to such an algorithm. One signal would be extracted from each peak in this plot, starting at the nearest min to the left of the peak and stopping at the nearest min to the right of the peak:

For instance signal3 was extracted from this plot between sample 217 and 234. Signal4 was extracted from another plot.

In case you are curious; each such plot corresponds to sound measurements by microphones at different positions in space. Each microphone receive the same signals but the signals are slightly shifted in time and distorted from microphone to microphone.

The feature vectors could be sent to a clustering algorithm such as k-means that would group together the signals with feature vectors close to each other.

Do any of you have any experience/advice on designing a feature vector that would be good at discriminating waveform signals?

Also which clustering algorithm would you use?

Thank in advance for any answers!

You want only objective criteria to seperate the signals or is it important that they have some sort of similarity when listened to by someone? That of course would have to restrict you to signals a bit longer (more than 1000 samples).


Is it necessary to scale a signal for quantization noise analysis? I am doing a series of tests to analyse the round-off (quantization) noise performance of single and two stage biquad cascades (direct form 1), at different sampling rates. The system I am simulating for this analysis uses 8.24 fixed point arithmetic.

At the moment the expected data for the filter will only range between 1 and -1, and no scaling will be implemented for the time being, to make use of the upper bits. The analysis I am doing uses a white noise signal.

Is it necessary to scale the stimulus up to the full bit depth of the filter, in order to get an accurate measure of the round-off noise of the system?

As I do not expect the filter to overflow in practice due to the massive headroom available, I am not sure if there is a more appropriate amplitude to scale the signal to, that exercises the most bits without causing an overflow.

Many Thanks

No, you don't need to scale the stimulus up to the full bit depth. That would actually end up underestimating the effect that roundoff noise would have on your actual system, as you would be artificially increasing the signal power (thus lowering the predicted signal to quantization noise ratio).

It makes sense to simulate the filter with a range of inputs that covers the space that you're designing the system to operate with. If the inputs will always be in the range $[-1, 1]$, then simulate using that range. On the other hand, if it's very likely that the data will be in that range, but possible that it might be larger, then it would be appropriate to test at higher levels. This requirement would likely be driven by higher-level considerations in your design.


How do I add more resolution / taps to an analog -> digital filter? (I am trying to create an IIR audio filter that adds reverb to an initial sample)

Say I designed an analog filter to model acoustic attenuation based on the following mathematical model:

$$
I = I_0 e^{pt},
$$

Where $p$ is some constant such that $-1 \leq p < 0$.

In the laplace domain, this is simply

$$
\mathfrak{L}\{I\} = \dfrac{1}{s-p}
$$

Using the impulse invariant transform, I get the digital equivalent

$$
Z\{I\} = \dfrac{z}{z-e^{pT_s}},
$$

Which inherently has one pole.

How do I add more resolution to this filter? Implementing this as a digital filter makes my audio sample sound like garbage. Would adding arbitrary taps keeps it from modeling the original analog decay?

Using experimental attenuation data and frequency sampling for an FIR equivalent, I can easily obtain 50000 poles and create a very clear reverb-adding filter.

Both IIR and FIR methods were tested in MATLAB, with the FIR using MATLAB's built-in church impulse response.

(Sorry if the question is unclear, I'm a bit new at this.)

Doing a decent sounding audio reverb with IIR filters is difficult. You need way more poles that you can generate with a normal IIR structure and you need a way to do it efficiently.

There is a significant trade off between sound quality versus MIPS, memory, latency & controllabilty
A good starting point is the original Schroeder reverb, which uses comb filters and warped allpass filters. See for example: https://ccrma.stanford.edu/~jos/pasp/Example_Schroeder_Reverberators.html
You want to get lots of poles cheaply. For example, wrapping  a feedback loop around a long delay will create lots of poles with just a single multiply + add.
Ideally the poles are randomly distributed and have no regularity to them. The simple Schroeder algortihm doesn't quite do this, so it tends to sound "ringy" or "metallic".
A better alternative are feedback delay networks. See for example https://ccrma.stanford.edu/~jos/pasp/FDN_Reverberation.html. These work by having multiple delay lines of different lengths with a feedback matrix: i.e. each delay is fed back into all delays. This creates a very large number of relatively low Q poles.
Finally you can dial in all the frequency dependencies and also add some early reflections as sparse FIR filters (with IIR frequency shaping on the spares tabs)