'A priori' and 'a posteriori' SNR I am reading this paper about speech enhancement using a minimum mean square error short time spectral amplitude estimator for our speech enhancement project. Basically in the paper they mix a clean speech with noise and try to recover the clean speech.

The paper keeps mentioning a priori SNR and a posteriori SNR. I am confused by what does it mean. A priori SNR is the SNR from before the mixing? But that can't be right. If we knew the clean speech signal from the beginning this exercise would be pointless.

This is the paper if you are wondering:
Y. Ephraim and D. Malah, “Speech Enhancement Using a- Minimum Mean Square Error Short-Time
Spectral Amplitude Estimator”, IEEE Transactions on Acoustics, Speech, and Signal Processing, Vol. ASSP-32,
No. 6, Dec. 1984.

thanks

The a priori SNR is the ratio of the power of the clean signal and of the noise power. The a posteriori SNR is the ratio of the squared magnitude of the observed noisy signal and the noise power. Both SNRs are computed for each frequency bin.

Of course, the only signal we have is the observed noisy signal. The noise power as well as the power of the clean signal must be estimated. The noise power can be estimated during absence of the desired signal, but in practice it can be difficult to detect the presence or absence of the desired signal, at least when the SNR is bad. There are so-called noise tracking algorithms that estimate the noise power without relying on detecting the presence or absence of the desired signal. They are based on the assumption that the desired signal is "more non-stationary" than the noise. The power of the desired clean signal can be estimated using a decision-directed approach, as described in Ephraim and Malah's paper. In that approach, the estimated clean amplitude, i.e. the output or 'decision' of the algorithm, and the a posteriori SNR are used to estimate the a priori SNR.


Deriving the Fourier transform of cosine and sine In this answer, Jim Clay writes:

  ... use the fact that $\mathcal F\{\cos(x)\} = \frac{\delta(w - 1) + \delta(w + 1)}{2}$ ...

The expression above is not too different from $\mathcal F\{{\cos(2\pi f_0t)\}=\frac{1}{2}(\delta(f-f_0)+\delta(f+f_0))}$.

I have been trying to obtain the later expression by using the standard definition of the Fourier transform $X(f)=\int_{-\infty}^{+\infty}x(t)e^{-j2\pi ft}dt$ but all I end up with is an expression so different from what's apparently the answer.

Here's my work:

\begin{align}
x(t)&=\cos(2\pi f_0t)\\
\Longrightarrow \mathcal F\left\{x(t)\right\}&=\int_{-\infty}^{+\infty}\cos(2\pi f_0t)e^{-j2\pi ft}dt\\
&=\int_{-\infty}^{+\infty}\frac 12 \left(e^{-j2\pi f_0t}+e^{j2\pi f_0t}\right)e^{-j2\pi ft}dt\\
&=\frac{1}{2}\int_{-\infty}^{+\infty}\left(e^{-j2\pi f_0t}e^{-j2\pi ft}+e^{j2\pi f_0t}e^{-j2\pi ft}\right)dt\\
&=\frac{1}{2}\int_{-\infty}^{+\infty}\left(e^{-j2\pi t\left(f_0+f\right)}+e^{-j2\pi t\left(f-f_0\right)}\right)dt\\
&=\frac{1}{2}\left(\int_{-\infty}^{+\infty}\left(e^{-j2\pi t(f_0+f)}\right)dt+\int_{-\infty}^{+\infty}\left(e^{-j2\pi t(f-f_0)}\right)\right) dt
\end{align}

This is where I'm stuck.

Your work is OK except for the problem that the Fourier transform of
$\cos(2\pi f_0 t)$ does not exist in the usual sense of a function of $f$,
and we have to extend the notion to include what are called distributions,
or impulses, or Dirac deltas, or (as we engineers are wont to do, much
to the disgust of mathematicians) delta functions.  Read
about the conditions that must be satisfied in order for the Fourier
transform $X(f)$ of the signal $x(t)$ to exist (in the usual sense)
and you will see that $\cos(2\pi f_0 t)$ does not have a Fourier transform
in the usual sense.

Turning to your specific question, once you understand that impulses
are defined only in terms of how they behave as integrands in an
integral, that is, for $a < x_0 < b$,
$$\int_{a}^{b} \delta(x-x_0)g(x)\,\mathrm dx = g(x_0)$$
provided that $g(x)$ is continuous at $x_0$, then
it is easier to deduce the Fourier transform of
$$\cos(2\pi f_0 t)
= \left.\left.\frac{1}{2}\right[e^{j2\pi f_0 t} + e^{-j2\pi f_0 t}\right]$$
by musing on the fact that
$$\int_{-\infty}^\infty \delta(f-f_0)e^{j2\pi ft}\,\mathrm df = e^{j2\pi f_0t}$$
and so it must be that $\cos(2\pi f_0 t)$ is the inverse
Fourier transform of $\displaystyle \left.\left.\frac{1}{2}\right[\delta(f-f_0) + \delta(f+f_0)\right]$.


Calculate transfer function of two parallel transfer functions in a feedback loop I have rational transfer functions:

$$H_1(z) = \frac{ b_0 + b_1z^{-1} + b_2z^{-2} }{a_0 + a_1z^{-1} + a_2z^{-2}}$$

$$H_2(z) = \frac{ q_0 + q_1z^{-1} + q_2z^{-2} }{p_0 + p_1z^{-1} + p_2z^{-2}}$$

And they are being used in the following feedback loop:

Sometimes this system blows up, and other times it is stable depending $H_1$ and $H_2$. I'm guessing that when it blows up, the poles of the combined system are outside the unit circle.

How do I calculate the transfer function of the system so I can start re-arranging it?

Your two transfer functions are in parallel, i.e they simply add up. So your feedback transfer function is simply $G(z) = H_1(z)+H_2(z)$. You want to makes sure that the magnitude of $G(z)$ is smaller than one.

and the overall closed loop transfer function is

$$H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1+H_1(z)+H_2(z)}$$


Autocorrelation of sinc function I'm having trouble on computing the autocorrelation of the sinc function
I want to compute
$$R_{hh}(\tau)=\int_{-\infty}^{\infty}\operatorname{sinc}(t) \ \operatorname{sinc}(t-\tau) \ \mathrm{d}t$$
where
$$\operatorname{sinc}(t) \triangleq \begin{cases}
\dfrac{\sin(\pi t)}{\pi t} \qquad & t \ne 0\\
\\
\quad \ \ 1 & t = 0
\end{cases}$$

Is the result another sinc function?

Because of the symmetry of the $\operatorname{sinc}(t)$ we have
$$R_{hh}(\tau)=\int_{-\infty}^{\infty}\operatorname{sinc}(t) \ \operatorname{sinc}(t-\tau) \ \mathrm{d}t=\int_{-\infty}^{\infty}\operatorname{sinc}(t) \ \operatorname{sinc}(\tau-t) \ \mathrm{d}t = \operatorname{sinc}(\tau) \ast \operatorname{sinc}(\tau) $$

Analyzing in the frequency domain we have
$$\mathcal{F}(R_{hh}(\tau)) = \mathcal{F}(\operatorname{sinc}(\tau))^2=\operatorname{rect}(f)^2=\operatorname{rect}(f)$$

$$\Rightarrow R_{hh}(\tau)=\mathcal{F}^{-1}(\mathcal{F}(R_{hh}(\tau))) = \mathcal{F}^{-1}(\operatorname{rect}(f))=\operatorname{sinc}(\tau)$$

where $\operatorname{rect}(f)$ denotes the rectangular function in the frequency domain


QPSK using Root Raised Cosine Filter The implementation of QPSK requires 2 orthogonal basis.

So when using root raise cosine filter to implement QPSK, what is the 2 orthogonal signals?

The way I prefer to look at it is to think of the QPSK signal as a complex envelope; that is, as a baseband signal. In this representation, you have a single basis, but the amplitudes are allowed to be complex.

So, let $p(t)$ be a root-raised cosine pulse such that $$\int_{-\infty}^\infty p(t-lT_p)p(t-mT_p) = \begin{cases}1\textrm{, if $l=m$}\\ 0\textrm{, if $l\neq m$,}\end{cases}$$ where $R_p=1/T_p$ is the pulse rate. Now you can think of the QPSK signal $s(t)$ as $$s(t)=\sum_k a_kp(t-kT_p)$$ where $a_k$ are complex symbols taken from a QPSK constellation. As you can see, you only need one set of real orthonormal pulses. The complex amplitudes $a_k$ can be recovered by matched filtering.

Another way to look at it is that you have two basis, one real and the other imaginary (the set $jp(t)$).


Triaxial accelerometer to single signal I have some smartphone recorded accelerometer data, so I have acceleration along 3 axes: $$a_{x}, a_{y}, a_{z}$$

The accelerometer data was collected by securing a phone to a walking subject. I'm primarily interested in being able to identify the time between each step they took. Normally, I could just search for peaks in the vertical acceleration vector, but in this case I don't know the orientation of the phone during the walk so I cant be sure which axis the vertical movement is coming in on

Individually, the 3 signals dont tell me anything useful because of the orientation problem, so I'm trying to condense them into a single vector. Presumably walking steps will cause the largest change in acceleration, so I can just time the peaks in this combined signal. Ultimately I'm looking to create a vector that will tell me the amount of total acceleration at each time point of the walk.

I've seen some people apply the following:

$$A = \sqrt{a_{x}^2 + a_{y}^2 + a_{z}^2}$$

A few questions about this:

What is the correct name for the resulting vector? I've heard it referred to as magnitude of acceleration, but when I google that term I see that its calculated as the following and unsure if they're referring to the same thing: $$\Delta v/\Delta t$$
If I wanted to calculate that value in Python, would it be as simple as taking each acceleration vector, element-wise square each one, element-wise sum the 3 vectors together, then take the root of each individual element? I know this is just the Pythagorean theorem, but unsure how to apply it in the case of vectors. Should the result be a vector?
How would you interpret the resulting value? Would it be the total amount of acceleration at each time point?
Given my goal of identifying time between steps, is this an appropriate approach? Or is there a more preferred way of combining the 3 acceleration signals?

Your expression $A=\sqrt{a_x^2+a_y^2+a_z^2}$ calculates the length of the 3-dimensional vector with coordinates $(a_x, a_y, a_z)$. So, it is the magnitude of the acceleration. Your expression $a=\Delta_v/\Delta_t$ also gives magnitude of acceleration, but in a slightly different context. This comes from measuring speed at different times, and then you average the acceleration between the time steps. You can ignore this expression for your case.
given ax, ay, az are $N$-dimensional vectors that contain the measurements of $a_x, a_y, a_z$ over time, you would do A=np.sqrt(ax**2+ay**2+az**2) to calculate the magnitude over time.
Yes, it's the overall acceleration at every time point. Though, it is only a relative value, its unit will not be exactly $m/s^2$, but there is some scaling involved to transform your measurement $A$ to the common $m/s^2$ units. Though, you don't need to do it, since you are only interested in relative changes of the acceleration over time.
Looks reasonable. Just plot $A$ vs time for some steps and decide on your own, if your measurement would contain sufficient information.


DSSS direct sequence spread spectrum - encoding each raw data bit (question) For implementing DSSS, it is understood that each raw data bit needs to be encoded as a new recognizeable/identifiable bit pattern.

In many literature involving DSSS, there is mention of using PN sequences involving LFSR (linear feedback shift registers).

It appears that maximal length linear shift registers have a maximal length that is always an odd-valued number, such as 15. So, while LFSR generators can be relatively fast and hard-ware implementable, there could be a hardware implementation problem elsewhere.

A pseudorandom sequence from a maximal length LFSR has bit pattern that has an odd number (such as 15 bits). So if (for example) the raw data bit has a bit-rate of 1 bit per second, then to encode each raw data bit with this 15 bit pattern would require the PN generator to operate at 15 bits per second (ie. 15 times higher clock rate than the raw bit rate).

Will this mean that we have to find hardware that generates a clock rate (for the PN generator) having 15 times the frequency of the raw bit rate?

Initially, I was thinking that a high frequency clock could be used as a reference. However, conventional flip flop frequency dividing methods divide frequencies by even amounts, like divide by 2, divide by 4 etc. There is no option for 'divide by 15' (or divide by an odd number).

So if I want to use LFSR PN sequence patterns for encoding raw data bits, then what kind of method should be used? Is this the usual problem faced when using maximal LFSR generators for encoding raw data bits?

My assumption is that a maximal LFSR produces a repetitive cyclic pattern.... such as 0100101 (eg. 7 bit sequence for a 3 bit LFSR), which keeps repeating as 0100101. So at the beginning of every repetition cycle, we can use this pattern to encode our raw data bit.

So a raw data bit such as '1' needs to be translated to this long pattern 0100101. And raw data bit '0' needs to be translated to the negative version of the pattern, such as 1011010. The problem is .... 7 bits is an odd number. So if the raw data bits need to be processed at 1000 bit per second, then we would need a PN clock rate that is at 7000 bits per second.

Due to the odd-valued number of bits, is there generally going to be an issue with coming up with a system for generating the clock rate that is 7 times higher than the raw data rate? I can understand that generating 8 times data rate would be relatively easy to do, because 8 is a 2^N value.

I'm thinking it would have been extremely nice if nature had allowed for maximal length LFSRs to be produced with even-values for maximal lengths.

Would fractional DDS type frequency generation methods be used instead?

Thanks for any comments in advance!!

In your title and your question overall, you confuse bits and symbols. Not a big deal – we can just work with the idea that we use a constellation where one symbol is equivalent to one bit (e.g. BPSK). But remember that there's no restriction that says you can't DSSS e.g. a 64-QAM modulated symbol stream. (In fact, that's what IEEE802.11b does, iirc).

I'll hence use "symbol" in place of "bit" where I cite you below, whenever it doesn't make a difference. When it does, I'll try to explain.

  Will this mean that we have to find hardware that generates a clock rate (for the PN generator) having 15 times the frequency of the raw symbol rate?

Yes! We call that higher rate chip rate, by the way.

  Initially, I was thinking that a high frequency clock could be used as a reference. However, conventional flip flop frequency dividing methods divide frequencies by even amounts, like divide by 2, divide by 4 etc. There is no option for 'divide by 15' (or divide by an odd number).

Well, you can certainly build a digital counter that counts to 15 and resets – your 1/15 clock divider. So, I think you might be imagining problems where there aren't problems. It's not as trivial as a $2^N$ counter, but honestly, a $2^N-1$ counter is just either a $2^N$ counter where you reset all bits to 0, but the lowest bit to 1 (giving you a single count of headstart), or really just a counter with a $N$-input NOR gate. Or really, just a counter that can compare to an arbitrary value. All these aren't hard to build. This isn't the 1960s!

  So if I want to use LFSR PN sequence patterns for encoding raw symbols, then what kind of method should be used? Is this the usual problem faced when using maximal LFSR generators?

The thing is even easier: You've got a device that repeats itself every 15 states, and changes at the high clock rate. Why not simply use that exact device to build your frequency divider? You could just compare all the state of SR to the one state on which you output a pulse. Done!

But again, integer frequency division really isn't that hard to implement at the rates you're probably worried about.

As another matter of fact: you'd probably not even use the higher frequency as reference, but a lower frequency, from which you synthesize the symbol and the chip rate, and quite possibly also the RF local oscillator. Since the relationship between chip and symbol rate is a rational one, a relatively simple integer-N synthesizer might be used, but frankly, aside from spur performance, there's little reason not to go to fractional-N synthesizers – which, by the way, are especially made for such systems.

  Due to the odd-valued number of bits, is there generally going to be an issue with coming up with a system for generating the clock rate that is 7 times higher than the raw data rate? I can understand that generating 8 times data rate would be relatively easy to do, because 8 is a 2^N value.

There's no hard tasks here. Count to 7, emit a pulse, reset your counter. Counters are relatively easy to build capable of high clock rates.

  Would fractional DDS type frequency generation methods be used instead?

Not quite sure what fractional DDS is, but a DDS method involves jumping through a table at high rates, to produce lower-frequency tones. I don't even think that's easier than a odd-number-counter to build a frequency divider. DDS is generally used to produce tones, not digital clocks.


Time delay estimation of two sensor signals I have two sensors which are placed around 30  to 40 cm apart , both of them measure the same   signal say 'lambda signal'. so there is a delay in the measurement of the signal from the second sensor. During the travel there is some noise added to the signal. What are the procedure that I can use to calculate the time delay between them.
According to my knowledge if I cross correlate both signals, the peak gives me the delay value of the second, but this is not working. Is there any other methods to find the time delay between eachother. Please guide me as I am new to signal processing.![Part of the signals][1]![Another part signal][2]

link: https://docs.google.com/open?id=0BzT9Q3HfpmDKZERPSGJoTEFyTzA
link : https://docs.google.com/open?id=0BzT9Q3HfpmDKeVhMSnZLa3BfMjQ
Thanks in advance.
Naveen

Comment: As a new user I am not allowed to add a comment to your question - hence this approach.

You mention that you do get a peak from the cross-correlation but that the lag estimation is incorrect.  This implies that you know the correct result - is from a calibration signal or ..? And what is the lag roughly?

Could you explain what the signals look like (sinusoidals for example)?  Also what is the SNR, sequence length, sampling frequency? Have you filtered the signals before sampling - if so how?

If the offset is incorrect it might be related to insufficient time resolution. Using a sampling frequency $f_s$ means that your time resolution is $T_\mathrm{s}=\frac{1}{f_\mathrm{s}}$.  If this time resolution is not sufficiently small you may get an estimate that does not agree with your expectations.

There may be other techniques but this depends on the signals you have - and also if you for example know that the lag is always a fraction of a period time of a periodic signal. Some more details could help with this.


What is sampling frequency in this case? I am using Matlab fdatool to implement highpass IIR butterworth filter
But 'Fs' value cell is gray and i am not able to change sampling frequency, what is the actual value in this case of attached photo?

Cutoff frequency for a discrete time filter can be expressed in different ways.  It can be expressed in radians, which is mathematically useful when analyzing the transfer function in the z-domain. It can be expressed normalized to the sample rate, as it is here, which is handy because it is a simple number to work with and is not dependent on the sample rate. Or you can express it as a true frequency, which requires a sample rate and is handy because it’s the actual frequency in the implementation.
It looks like you just need to change the units in the drop down box to set the sample rate. Probably something like Hz, but I’ve never used that tool.


Can the magnitude of DFT of a triangular wave be max for at 0? I'm taking an online class on Signal Processing, and one of the slides of from the class shows the magnitude of DFT for a simple triangular wave to have the maximum at 0.

How can this be possible? This means it only has presences of the component form a sine wave with frequency 0 right?

What you see in the red graph is not an impulse ($\delta$) at zero frequency (in such case we would have a constant zero frequency which is equivalent to a constant dc signal $x[n]=c$).

What you see is actually a squared $\mathrm{sinc}$ function. You can easily verify that the frequency representation of a triangular function (both in continuous and discrete time) is in the form of $\mathrm{sinc}^2$. Look at here, for instance for the continuous triangular signal to see that

$$\mathcal{F}\{\mathrm{tri}(t)\}=\mathrm{sinc}^2(f)$$

and $\mathrm{sinc}$ has a maximum at zero.


How to compute Laplace Transform in Python? I am trying to do practicals for signal processing where I need to Laplace Transform a function. Used 'fft' of numpy before. Nothing of Laplace is found in the documentation. Do we have any other alternative?

Please go through the notebook to understand the problem (would love to get suggestions/contributions]

https://github.com/sachinmotwani20/Raw-Signal-Processing-Python/blob/master/LaplaceToFourier.ipynb

Update : Given the approach started in the OP's Github code I have this suggestion :
Observe that the unilateral Laplace Transform given as:
$$X(s) = \int_0^\infty x(t)e^{-st}dt$$
Is just the Fourier Transform of a causal function with a weighting exponential:
$$X(s) = \int_0^\infty x(t)e^{-(\sigma+j\omega)t}dt$$
$$X(s) = \int_0^\infty [e^{-\sigma t}x(t)]e^{-j\omega t}dt$$
Which is the Fourier Transform of $e^{-\sigma t}x(t)$.
So to proceed, learn how to produce surface plots in python, and then index through $\sigma$ within the Region of Convergence (see below) and compute the FFT of $e^{-\sigma t}x(t)$  to create your complex surface values given $\sigma$ and $\omega$.

The challenge with what you are trying to do is that the Laplace Transform is a function of the complex variable "s", so for each possible value of "s" (which is simply the set of all complex numbers) the Laplace Transform would have a complex result with a magnitude and phase.
So when you try to plot this, you could plot the magnitude and phase as two separate plots but these would be 3D surface plots showing either the result of the magnitude or the phase over a 2D surface, similar otherwise to how we plot the magnitude and phase of the Fourier Transform as a 2D function over the line representing all frequencies.
For the graphical representation of the Laplace Transform, we typically just show the locations where that function goes to infinity (poles) or is zero (zeros). In fact every other sample on the surface is uniquely defined by the pole and zero locations, so that is all we need to show to define it.
This is best shown with an example. Since I already have the graphic for this particular case, consider the time domain function of a decaying sinusoid given by the formula below and the plot below that where we see in the dashed red line the envelope for the decaying function $2e^{-0.2t}$. $u(t)$ is the step function which is $0$ for time $t<0$ and $1$ for time $t \ge 0$.
$$ x(t) = u(t)2e^{-0.2t}sin(0.5t)$$

To get the Laplace Transform (easily), we decompose the function above into exponential form and then use the fundamental transform for an exponential given as :
$$\mathscr{L}\{u(t) e^{-\alpha t}\} = \frac{1}{s+\alpha}$$
This is the unilateral Laplace Transform (defined for $t = 0$ to $\infty$), and this relationship goes a long way since we can describe the response of any causal linear system using such exponential forms.
So the equation above, assuming $t>0$, and using Euler's identity becomes :
$$ x(t) = u(t)2e^{-0.2t}sin(0.5t)$$
$$ = 2e^{-0.2t}\frac{(e^{+j0.5t}-e^{-j0.5t})}{2j}$$
$$ = je^{-0.2t}(e^{+j0.5t}-e^{-j0.5t})$$
$$ = je^{-(0.2-j0.5)t} -je^{-(0.2+j0.5)t}$$
Which we can then easily take the Laplace Transform to get a function of the complex variable $s$ :
$$ X(s) = \frac{j}{s+0.2+j0.5}-\frac{j}{s+0.2-j0.5} = \frac{1}{(s+0.2+j0.5)(s+0.2-j0.5)}$$
The graph of the magnitude of this is shown below, which is a surface plot for all values of $s$, where $s$ is a complex variable with real and imaginary components traditionally described as $s = \sigma +j\omega$.  So we simply choose a particular complex value $s$, and plot the magnitude of the result of $X(s)$.

Caveat : This is indeed the magnitude of the Laplace Transform if the Laplace Transform exists. For causal functions, the Laplace Transform does not exist for all values of s whose real part is to the left of the leftmost pole. Why? Because for any values of $s$ in that region (left of the leftmost pole), the Laplace Transform (given by the integral) will not converge (grows to infinity). While for all values of s at and to the right of the leftmost pole, the transform will converge to be the magnitude on the surface plotted above. So the proper graphic would only be valid to the right of the right most pole. This can be confusing as $X(s)$ certainly converges as was done to make this plot, but it is in the transformation itself that the result cannot be obtained.
This is very clear if we consider the envelope in our example function $x(t)$ which was given by $2e^{-0.2t}$ for all $t>0$, and the Laplace Transform $X(s)$ for $s = -1$ :
$$X(s) = \int_0^\infty x(t)e^{-st}dt$$
$$X(s= -1) = \int_0^\infty 2e^{-0.2t}e^{t}dt = \int_0^\infty 2e^{0.8t}dt$$
which will not converge since the function $e^{0.8t}$ continuously grows larger for larger $t$.
Hopefully after reading this you no longer feel the need to plot the Laplace Transform, and in practical application a plot of it is never used beyond showing the pole and zero locations. However what is very useful is knowing that the Fourier Transform is the Laplace Transform when $s = j\omega$. You can see this if you compare the two equations, and the small breakout in upper right-hand corner of the plot above is also showing this, which is the Frequency Response specifically. The one precaution is that the Fourier Transform is often given as a bilateral function (t extending from $-\infty$ to $\infty$) so to be truly equivalent unless the function is declared to be causal, we must be using the bilateral Laplace Transform for the two to be exactly identical (which is also seldom used). This explains why the Fourier Transform of a sine wave would appear as two impulses at $\pm \omega_c$ while as you might deduce from the plot above, the Laplace Transform of the sine wave would appear more as tent poles on the $j\omega$ axis. You would see this as well with the Fourier Transform if you did the transform of $u(t)cos(\omega t)$ specifically.


Design Stable IIR Filter in MATLAB I need to design a stable IIR filter in MATLAB.  Specifically, I am trying to use a butterworth filter to correct the phase and gain of a 2 pole RC filter:

d = fdesign.arbmagnphase('n,f,h', n, f, h);
myFilter = design(d,'iir','butter');

The butterworth part is not really important, there are several options to choose from.  Is it possible to constrain MATLAB to produce an IIR filter that is stable?

After working through some possible solutions, we figured out it would be easier to design argitrary magnitude and group delay separately in MATLAB:

fdesign.arbmag
fdesign.arbgrpdelay


How Can an Object Move by Less Than a Pixel Between Two Consecutive Frames? I think video is a series of picture frames. Consider a black & white scenario (no RGB). A picture frame is a 2D array of pixel. Each pixel has a single color.
A pixel can't be half black or half white. It will be filled completely with one shade of color.

Consider 2 consecutive frames of video - let's say frame 1 and frame 2 .
Suppose frame 1 has a square-shaped object with coordinates $(x_1,y_1)$, $(x_2,y_2)$, $(x_3,y_3)$ and $(x_4,y_4)$ (top-left, top-right, bottom-right and bottom-left, respectively).

Now, if the object moves to the right then in frame 2 its coordinates can be $(x_1+1,y_1)$, $(x_2+1,y_2)$, $(x_3+1,y_3)$ and $(x_4+1,y_4)$ (top-left, top-right, bottom-right and bottom-left, respectively).

How is less than that possible?

The frame is 2D Grid which samples the continuous world in Discrete way.

Leave alone video for a second.
Let's talk about taking a photo of a white paper laid on a black surface.
The photo is taken from above and let's assume there is no lens issues.
Do you expect the end of the paper to end at edge of a pixel of the camera perfectly?

It won't.
The end of the paper will fall in the middle of a pixel.
Let's say that now this pixel (I leave out the blurring due to limited Frequency Sampling) has 70% White and 30% Black.
Now the camera moves the smallest movement and now the pixel is 30% White and 70% Black.

This is a less than 1 pixel movement.


OFDM Simulation process I'm trying to understand OFDM by making a simulation. Are these steps correct?

generate M random complex QAM symbols. example: (1+j,1-j,-1-j,1-j....)
Split my M samples into a N 2048 sized arrays
Take IFFT of each array individually
Add cyclic prefix to each array
Make N square shaped signals by oversampling my IFFT arrays
Pulse shape the individual signals with raised cosine
Mix up the signals with frequency separation = 1/symbol_period by multiplying them by $e^{j2\pi ft}$

generate M random complex QAM symbols. example: (1+j,1-j,-1-j,1-j....)
  Split my M samples into a N 2048 sized arrays
  Take IFFT of each array individually
  Add cyclic prefix to each array

Exactly! In real-world OFDM systems, not all carriers are used; typically, you leave out one or two DC carriers and leave a few guard carriers at the band edges, so you'd have something like

  split my M samples into 1950-sized chunks, and map them to the elements of a 2048-vector, leaving the center ones and edge ones free

The edge guard carriers make the fifth step unnecessary:

  Make N square shaped signals by oversampling my IFFT arrays

so scratch that, usually.

  Pulse shape the individual signals with raised cosine

Don't do that! That's wrong. This is OFDM, and it's only used when you encounter multipath channels; matched filtering does nothing good here; it only makes things worse by making the channel more complicated and reducing average power.

Remember why you'd do raised cosine with single carrier systems: You want things to be

nicely pulse-shaped, which simply doesn't make sense in the context of OFDM, and
you want matched filtering because you want to maximize SNR in an AWGN channel, and OFDM is never used in channels that are purely AWGN, and
because after convolution with the conjugate pulse shaping filter in the receiver, you might have an easier time recovering timing, but that doesn't apply to OFDM either, because timing recovery must, in the wideband channel scenario, be done using other means (keywords: Schmidl-Cox, autocorrelation timing recovery methods).

  Mix up the signals with frequency separation = 1/symbol_period by multiplying them by $e^{j2\pi ft}$

No; your subcarrier separation is what the IFFT does; you just mix up the baseband signal generated by step 4 to bandpass with your $e^{j2\pi ft}$.

Since you're only simulating things, and we know that bandpass-region channels have an equivalent representation in baseband (in fact, only because we know that works is that we can use OFDM at all): don't do this at all, but work with the complex baseband signal directly.


Is multiplying by Ts faster than dividing by fs? Say fs = 1000 and Ts = 0.001. Would it be faster to compute Ts at the beginning and subsequently multiply by 0.001 instead of dividing by 1000 when computing frequency-dependent quantities?

Generally, it makes sense to ensure that your code is logically correct, that it is numerically well-behaved, intuitive to read and tested. That is hard enough. Only when you observe that some innerloop or library call is a real hotspot, affecting the functionality of your software does it make sense to rewrite code for speed, and then you should always profile before and after.
If a constant is known compile time, the compiler may apply the inversion to substitute division for multiplication, if this is within precision constraints and runs faster for a given target. If ppssible, I would rather outsource that complexity to the compiler.
Edit:
It is not "technically the same operation". To see why, have a look at this MATLAB snippet:
a = single(10.0)
b = 1/a
c = 42/a
d = 42*b

c-d

ans =

  single

 -4.7684e-07

Since floating-point is operating with finite precision and intermediate rounding, the order of operations does matter. Depending on compiler flags, the compiler may be allowed to re-order floatingpoint arithmetic even though the result will differ to some degree.
If we look at the binary representation we see that they differ in the lsb:
dec2bin(typecast(c,'uint32'),32)
dec2bin(typecast(d,'uint32'),32)

ans = '01000000100001100110011001100110'
ans = '01000000100001100110011001100111'

-k


Continuous time double exponential filtering in state space form? I'm trying to determine the continuous time formulation of the double exponential filter so that I can adapt it more flexibly for my particular problem.

Typically, this model is expressed as a pair of discrete, recursive equations as such in @denis' post here)

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

Most references on the subject of double exponential filtering and state space revolve around discrete time kalman filtering approaches. Question: How should I put this model into a noiseless, state-space ODE form?

Using slightly different notation, here is my current approach. The base model, with continuous time input $u$ and continuous time states $l$ and $s$, is

$\dot{l} = r_1(u - l)  + r_1 s$

$\dot{s} = r_2(\dot{l} - s)$

$\ \ = r_2(r_1(u - l)  + r_1 s - s)$

$\ \ = r_2(r_1(u - l)  + (r_1 - 1)s)$

$\ \ = r_2 r_1(u - l)  + r_2(r_1 - 1)s$

with $\dot{l}$ and $\dot{s}$ representing time-derivative state dynamics with system constants $r_1, r_2 \in [0, \infty)$.

Assuming I made no mistakes there, I attempt to put this into state space form by

$$ \dot{\bf{x}} = A{\bf x} + Bu$$
$$ A = \begin{bmatrix}
-r_1 & r_1\\
-r_1r_2 & r_2(r_1 - 1)\\
\end{bmatrix} \ \ \ \ B = \begin{bmatrix}
r_1\\
r_1r_2\\
\end{bmatrix} \ \ \ \ \bf{x} = \begin{bmatrix}
l\\
s\\
\end{bmatrix}
$$

The output of this system is just the state of the level variable.

$y(t^{'}) = l(t^{'}) = C{\bf x}$.

where $C = \begin{bmatrix}1 & 0\end{bmatrix}$

We can make future predictions of $y$ based off of some assumption of the form of $u$, we can project $y({t'})$ using the system response formula for continuous time state space systems

$y({t^{'}}) = C\bigg(e^{A(t^{'} - t)}x(t) + \int_t^{t^{'}}e^{A(t - \tau)}Bu(\tau)d\tau\bigg)$

which can then be discretized as necessary later.

However: I am finding that the system matrix $A$ is not stable.

By my assumptions system matrices for all positive values of $r_1, r_2$ should be stable but this is not the case. Empirically I just tested a couple settings and the real part of the eigenvalues of $A$ must be negative, which is satisfied are for $r_1 = 1, r_2=1$ but not for $r_1 = 3, r_2 = 3$

Where dd I go wrong?.

I think I figured it out, dumb mistake in writing the system in continuous time.

$\dot{l} = r_1(u - l)$

$\dot{s} = r_2(\dot{l} - s)$

$\ \ = r_2(r_1(u - l)  + r_1 s - s)$

$\ \ = r_2(r_1(u - l)  + (r_1 - 1)s)$

$\ \ = r_2 r_1(u - l)  + r_2(r_1 - 1)s$

So $$
A = \begin{bmatrix}
-r_1 & 0\\
-r_1r_2 & -r_1r_2\\
\end{bmatrix}
$$

which is stable for all positive $r_1, r_2$!


Use polyphase representation to construct a perfect reconstruction filter bank I'm posting this question in response to Question about zero-phase filter bank

This is a homework question for which I found an answer, but I didn't do it the way I was supposed to. I'm wondering what the correct approach is. I've been told to use polyphase representation to solve the problem but I couldn't come up with a way to do it.

Let $H_0(z) = 1 + 2z^{-1} + 3z^{-2} + 2z^{-3} + z^{-4}$ and $H_1(z) = H_0(-z)$. Find causal FIR filters $F_0(z)$ and $F_1(z)$ such that $\hat{x}(n)$ agrees with $x(n)$ except for a possible delay and (nonzero) scale.

I have found a solution with no delay and unity scale (see the linked question) but I didn't use polyphase representation. I tried a few times to work it using polyphase but I kept ending up with IIR solutions. What approach should be used that takes advantage of polyphase?

This is a problem from Multirate Systems and Filter Banks by P. P. Vaidyanathan. In Section 4.6.2.E, there is a discussion on what he calls Euclidean Complementary Functions. Specifically, he discusses one of the results of Euclid's theorem:

  ... if $H_0(z)$ and $H_1(z)$ are relatively prime, there exists polynomials $F_0(z)$ and $F_1(z)$ such that
  $$H_0(z)F_0(z) + H_1(z)F_1(z) = c, c \neq 0.$$

In the text, he uses $H_k(z)$ and $F_k(z)$ to denote the analysis and synthesis filters, respectively. He uses $E_k(z)$ and $R_k(z)$ to denote the Type I and Type II polyphase components of the analysis and synthesis filter banks, so I'll do the same here.

If you take the polyphase decomposition of $H_0(z)$, we find that $E_0(z) = 1 + 3z^{-1} + z^{-2}$ and $E_1(z) = 2 + 2z^{-1}$. Because $E_0(z)$ and $E_1(z)$ are relatively prime, we can find an $R_0(z)$ and an $R_1(z)$ to fulfill

$$E_0(z)R_0(z) + E_1(z)R_1(z) = c, c \neq 0.$$
We can recognize that $c$ is the inevitable delay for $\hat{x}(n)$, so we should probably make $c = 1$. To cancel out the $z^{-1}$ terms so that $c = 1$, the order of $R_0(z)$ needs to be one less than that of $R_1(z)$, by inspection. So, we can represent $R_0(z) = \alpha$ and $R_1(z) = \beta + \gamma z^{-1}$, then solve for the coefficients. If we multiply everything out, we'll see that we want

$$\alpha + 2\beta = 1,\\
3\alpha + 2\beta + 2\gamma = 0, \text{and}\\
\alpha + 2\gamma = 0,$$
which is just solving

$$\begin{bmatrix}
1 & 2 & 0 \\
3 & 2 & 2 \\
1 & 0 & 2
\end{bmatrix}
\begin{bmatrix}
\alpha \\
\beta \\
\gamma
\end{bmatrix}
=
\begin{bmatrix}
1 \\
0 \\
0
\end{bmatrix}
$$
to get $\alpha = -1$, $\beta = 1$, and $\gamma = 0.5$. Now, we can solve the problem by composing $F_k(z)$ from the Type II polyphase components: $F_0(z) = z^{-1} R_0(z^2) + R_1(z^2)$ and $F_1(z) = z^{-1} R_0(z^2) - R_1(z^2)$.

There is actually also a nonzero scaling factor of 2 that comes from the DFT filter bank representation, because $\mathbf{W}^*\mathbf{W} = M\mathbf{I}$, where $M = 2$ because we have two channels.


Which channel coding technique is suitable for wireless SPIHT compressed image: Reed-Solomon or convolutional code? I want to transmit a SPIHT-based (Set Partitioning In Hierarchical Trees) image compressed image over a wireless fading channel, which will be suitable as channel coding technique: Reed-Solomon or convolutional code?

The channel coding used has little to do with your source coding. So, this depends on how you model your channel (and noise, and interferers), and what kind of errors you're willing to tolerate. Saying "SPIHT-compressed image" answers neither!


Good Language for Implementing Computer Vision System / Algorithm / Application? I am trying to implement a content based image retrieval system but before I do so I would like to get an overview of some programming languages suitable for this task (having good libs and such).

Does anyone know some good languages and libs for that kind of task? What about Python or Java?

Best

Maybe you can be more specific about the scope and scale of your work (academic project? Desktop or Mobile commercial product? Web-based commercial project?).

Some recommendations and comments:

Matlab is common in the academic world, and quite good for sketching/validating ideas. You will have access to a large body of code from other researchers (in CV and machine learning); prototyping and debugging will be very fast and easy, but whatever you will have developed in this environment will be hard to put in production. Depending on what your code is doing, you might have memory/performance problems (there are situations where you can't describe what you want to do in terms of Matlab's primitives and have to start looping on pixels and Matlab's being an interpreted language is not helping in this context). Interaction with databases, web servers etc is not easy, sometimes impossible (you won't get a Matlab program to become a Thrift server called by a web front-end). Costs $$$.
C++ is what is used for many production-grade CV systems (think of something at the scale of Google's image search or Streetview, or many commercial robotics applications). Good libraries like OpenCV, excellent performance, easy to put into a production environment. If you need to do machine learning, there are many libraries out there (LibSVM / SVMlight, Torch). If you have to resort to "loop on all pixels" code it will perform well. Easy to use for coding the systems/storage layers needed in a large scale retrieval system (eg: a very large on-disk hash map for storing an inverted index mapping feature hashes to images). Things like Thrift / Message Pack can turn your retrieval program into a RPC server which can be called by a web front-end. However: not very agile for prototyping, quite terrible for trying out new ideas, slower development time; and put in the hands of inexperienced coders might have hard to track performances and/or instability problems.
Python is somehow a middle ground between both. You can use it for Matlab style numerical computing (with numpy and scipy) + have bindings to libraries like OpenCV. You can do systems / data structure stuff with it and get acceptable performances. There are quite a few machine learning packages out there though less than in Matlab or C++. Unless you have to resort to "loop on all pixels" code, you will be able to code pretty much everything you could have done with C++ with a 1:1.5 to 1:3 ratio of performance and 2:1 to 10:1 ratio of source code size (debatable). But depending on the success of your project there will be a point where performance will be an issue and when rewriting to C++ won't be an option.


A question about sampling rate of cosine signal Given $$c(t) = \cos(2\pi\cdot 30 \cdot t) $$

If we sample this signal at the Nyquist rate 60 Hz and at a higher rate of 80 Hz, we get the following:

There is no aliasing as $f$ = 30 Hz is less than or equal to the folding frequency, 30 Hz and 40 Hz, respectively.

Why does sampling at a rate of 80 Hz result in this "distorted" cosine signal, although it is higher than the Nyquist rate? (I know if we go way higher than 80 Hz the signal will look better.)

It is actually not distorted, it is sampled at high enough rate. What fools you is the straight lines drawn between sample points, it gives you a false impression of the waveform. It shows you a linear interpolation of the signal. It does not represent how the signal would actually look like. A sampled signal exists only at the sample points, and to convert that into continuous analog representation that also exists between the sample points, you should use sinc interpolation which limits the bandwidth to half of the sampling rate.


Time delay estimation of oscilloscope signals using cross correlation I have recorded 2 signals from an oscope. They look like this:

I want to measure the time delay between them in Matlab.
Each signal has 2000 samples with a sampling frequency of 2001000.5.

The data is in a csv file. This is what I have so far.
I erased the time data out of the csv file so that only the voltage levels are in the csv file.

x1 = csvread('C://scope1.csv');
x2 = csvread('C://scope2.csv');
cc = xcorr(x1,x2);
plot(cc);

This gives this result:

From what I've read I need to take the cross correlation of these signals and this should give me a peak relating to the time delay. However when I take the cross correlation of these signals I get a peak at 2000 which I know is not correct. What should I do to these signals before I cross correlate them? Just looking for some direction.

EDIT: after removing the DC offset this is the result I am now getting:

Is there a way to clean this up to get a more defined time delay?

EDIT 2: Here are the files:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv

@NickS

Since it is far from certain that the second signal in the plots is in fact a solely delayed version of the first, other methods besides the classical cross-correlation have to be attempted. This is because the cross-correlation (CC) is merely a maximum likelihood estimator if your signal(s) are delayed versions of each other. In this case, they are clearly not, to say nothing about the non-stationarity of them either.

In this case, I believe what may work is a time-estimation of the significant energy of the signals. Granted, 'significant' can or cant be somewhat subjective, but I believe that by looking at your signals from a statistical point of view, we will be able to quantify 'significant' and go from there.

To this end, I did the following:

STEP 1: Compute the signal envelopes:

This step is simple, as the absolute value of output of the Hilbert-Transform of each of your signals is computed. There are other methods to compute envelopes, but this is pretty straight forward. This method essentially computes the analytical form of your signal, in other words, the phasor representation. When you take the absolute value, you are destroying phase and only after energy.

Furthermore since we are pursuing a time-delay-estimate of the energy of your signals, this approach is warranted.

STEP 2: De-noise with edge-preserving non-linear Medial Filters:

This is an important step. The objective here is to smooth out your energy envelopes, but without destruction or smoothing out your edges and fast rise times. There is actually an entire field devoted to this, but for our purposes here, we can simply use an easy to implement non-linear Medial filter. (Median Filtering). This is a powerful technique because unlike mean filtering, medial filtering will not null out your edges, but at the same time 'smooth' out your signal without significant degradation of the important edges, since at no time is any arithmetic being performed on your signal (provided the window length is odd). For our case here, I selected a medial filter of window size 25 samples:

STEP 3: Remove Time: Construct Gaussian Kernel Density Estimation Functions:

What would occur if you looked at the above plot sideways instead of the normal way? Mathematically speaking, that means, what would you get if you projected every sample of our denoised signals onto the y-amplitude-axis? In doing this we will manage to remove time so to speak, and be able to study the signal statistics solely.

Intuitively what pops out of the figure above? While noise energy is low, it has the advantage that it is more 'popular'. In contrast, while the signal envelope that has energy is more energetic than the noise, it is fragmented across thresholds. What if we considered the 'popularity' as a measure of energy? This is what we will do with (my crude) implementation of a Kernel Density Function, (KDE), with a Gaussian Kernel.

To do this, every sample is taken and a gaussian function constructed using its value as a mean, and a pre-set bandwidth (variance) selected a-priori. Setting the variance of your gaussian is an important parameter, but you can set it based on noise statistics based on your application and typical signals. (I only have your 2 files to go off on). If we then construct the KDE Estimation, we get the following plot:

You can think of the KDE as a continuous form of a histogram so to speak, and the variance as your bin-width. However it has the advantage of guaranteeing a smooth PDF that we can then perform first and second derivitave calculus on. Now that we have the Gaussian KDEs, we can see where the noise samples peak in popularity. Remember that the x-axis here represents  the projections of our data onto the amplitude space. Thus, we can see which thresholds the noise is the most 'energetic' in, and those tell us which thresholds to avoid.

In the second plot, the first-derivative of the Gaussian KDEs is taken, and we pick the abscissa of the first sample after first-derivative after the peak of the mixture of Gaussians to attain a certain value close to zero. (Or first zero-crossing). We can use this method and be 'safe' because our KDE was constructed of smooth Gaussians of reasonable bandwidth, and the first derivative of this smooth and noise-less function was taken. (Typically first-derivatives can be problematic in anything but high SNR signals since they magnify noise).

The black lines show then at what thresholds we would be wise to 'segment' the image at, such that we avoid the entire noise floor. If we then apply to our original signals, we attain the following plots, with the black lines indicating the start of energy of our signals:

This thus yields a $\delta{t} = 241$ samples.

I hope this helped.


LTI, causal, discrete time system output Consider a discrete time LTI causal system $S: y = S(u)$, with its impulse response $h:{Z} \rightarrow R:h(n)=3^{n+1}{H(n)}$ with $H$ the Heaviside function. We know the values of input:

$$u(0) = 2$$
$$u(1) = 1$$
$$u(2) = 1$$
$$u(3) = 1$$
$$u(4) = 2$$

Also the $y(2) = 66$.
I'm trying to find the value of $y(3)$.

I may say, for a LTI, causal system the output is the convolution of the input with the impulse response:

$$y(n) = \sum^\infty_{k=0}3^{k+1}u(n-k)$$

Then

$$y(3) = 3u(3)+9u(2)+27u(1)+81u(0) = 201$$

Which is the correct answer! But I don't understand, why should I stop at $u(0)$ or why the negative values of the input $u(n), n< 0$ are zero here ?

Alternative solution:

recursion gives us this form:

$$y(n) = 3y(n-1)+3u(n)$$

Unless explicitly stated, you shall not assume that the input is zero for $n<0$, hence the convolution sum should be

$$y(n) = \sum^\infty_{k=-\infty}3^{k+1}u(n-k)=\sum^\infty_{k=0}3^{k+1}u(n-k).$$

In this situation, without making any assumption on $u[n]$ for $n<0$ you can still compute $y[n]$ for $n \geq 0$ based on the recursive relation as you have outlined :

$$y(n) = 3y(n-1)+3u(n) \implies y[3] = 3y[2] + 3u[3] = 201$$

What can be said about $u[n]$ for $n<0$ ? Actually nothing but, given that $y[2] = 66$ we may try the folowing :

$$
\begin{align}
y[2] &= \sum^2_{k=0}3^{k+1}u(2-k) + \sum^\infty_{k=3}3^{k+1}u(2-k) \\
&= 66 + \sum^\infty_{k=3}3^{k+1}u(2-k) \\
66 &= 66 +   \sum^\infty_{k=3}3^{k+1}u(2-k) \\
0 &= \sum^\infty_{k=3}3^{k+1}u(2-k)
\end{align}
$$

That the sum in the parenthesis equals zero, has a trivial solution of $x[n] = 0$ for $ n <0$.

However, we may also ask whether another, nonzero, solution is also possible that the sum be zero? Consider, for example, the input  $$u[n] = (-1)^{2-n} ~~3^{n-3}$$ for $n<0$ which yields for the sum

$$\sum^\infty_{k=3}3^{k+1}u(2-k) = 1-1+1-1+1-...$$

This is an alternating sequence and it does not converge. However, if it's possible to restrict the length of the input to be even for $n<0$, then the sum will be zero too...

Nevertheless, in order to uniquely specify $u[n]$ for $n<0$ a single sum is not enough. We must know values of output $y[n]$ for all $n<0$ for that. Hence you cannot conlcude that $x[n]=0$ for $n<0$ under the given situation.


Skew map to match aerial image I have the following aerial image (apparently taken from an airplane):

I also have a map of the area shown in the image:

On the image and the map I can mark the corresponding streets, street intersections and landmarks (I have shown only three here, but I can provide more points for precision).

I would like to skew the map to match the image (camera position, angle, etc).  How can I do this?

(I can program OpenCV and AForge.Net, and can use Gimp and Blender, all on Windows.)

Well you could try to estimate a transformation matrix that maps between the coordinates in the two images. The idea is that

\begin{equation}
\begin{bmatrix}
x_a\\
y_a
\end{bmatrix}
=
\begin{bmatrix}
A_{xx} & A_{yx}\\
A_{xy} & A_{yy}
\end{bmatrix}
\begin{bmatrix}
x_m\\
y_m
\end{bmatrix}
\end{equation}

where $x_a, y_a$ are the coordinates (in pixels for example) of the aerial image and $x_m, y_m$ are the coordinates of the photograph. Writing this out leads to the system of equations:

\begin{equation}
x_a = A_{xx}x_m + A_{yx}y_m
\end{equation}
\begin{equation}
y_a = A_{xy}x_m + A_{yy}y_m
\end{equation}

That way you can represent scaling, rotation and shearing operations between the images. You need at least four points per image to solve that system but I recommend more and then calculate the least squares solution using linear regression. This will also give you an indication on how well everything works. If there are large residual errors after the fit, then the linear transformation might not be sufficient.

If you then want to overlay the images, just apply the transformation matrix to your map and the result should be an image that has the reference coordinates on top of each other and the areas in between are linearly transformed accordingly.

There are many things that can go wrong here. For once, it only works for roughly flat terrain  (even though that effect becomes weaker if the plane is further away). Also, if there is lens distortions in the photograph, that will affect the accuracy of the fit (usually close to the image borders). And finally, map inaccuracies might throw you off.

In the end, it's just a matter of trying. If the 2D-transformation model is not capable of delivering accurate results, you need to model more effects. But estimating those will probably also require more referene points.

One other idea would be to put the image on a plane (like a texture) and place the plane in a 3D-space and match it with the plane of the photograph. That would end you up with a 3-dimensional transformation matrix for this plane. Might work as well.


Kalman Filter EM Estimation of Covariances The question might be very simple, but I get a strange result from Kalman Filter. Let us consider the simplest state-space model, the random walk plus noise:
$$
y_{t} = x_{t} + \varepsilon_{t}\\
x_{t} = x_{t-1} + \mu_{t-1}
$$
Let us generate $N_{data} = 3000$ data points and with parameters: $Var(\mu_{t}) = 1$ (which is transition variance) and $Var(\varepsilon_{t}) = 1$ (which is observation variance).

Using Kalman EM in Python, I try to estimate the transition and observation variance.

First, I would like to estimate "observation_covariance" and I assume that "transition_covariance" is known:

kf = KalmanFilter(
    transition_matrices=[1.], observation_matrices=Data['x'],
    transition_covariance=[1.], initial_state_mean=3,
    initial_state_covariance=[1.],
    em_vars=['observation_covariance']
)

kf = kf.em(Data['y'])

the result is $0.79$.

Second, I do the opposite, i.e. I estimate "transition_covariance" and I assume that "observation_covariance" is known:

kf = KalmanFilter(
    transition_matrices=[1.], observation_matrices=Data['x'],
    observation_covariance=[1.], initial_state_mean=3,
    initial_state_covariance=[1.],
    em_vars=['transition_covariance']
)

kf = kf.em(Data['y'])

the result is $0.001$, which is far from $1$!

Next, I assume that both "transition_covariance" and "observation_covariance" are unknown:

kf = KalmanFilter(
    transition_matrices=[1.], observation_matrices=Data['x'],
    initial_state_mean=3, initial_state_covariance=[1.],
    em_vars=['transition_covariance', 'observation_covariance']
)

kf = kf.em(Data['y'])

the result is "observation_covariance" = 0.77 and "transition_covariance" = 0.0002.

What do I do wrong? Shouldn't those estimated variances be close to $1$?

The all code I use is below:

#%%
import pandas as pd
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt
import pylab as pl
from pykalman import KalmanFilter

N_data = 3000
Data = pd.DataFrame(columns=['NoiseTrs', 'NoiseObs', 'x', 'y'], index=range(N_data))
Data['NoiseTrs'] = np.random.normal(loc=0.0, scale=1.0, size=N_data)
Data['NoiseObs'] = np.random.normal(loc=0.0, scale=1.0, size=N_data)

# random walk

x_f_v = 3 # first value for the simulation
for i in range(N_data):
    if i == 0:
        Data.loc[i, 'x'] = x_f_v
    else:
        Data.loc[i, 'x'] = Data.loc[i - 1, 'x'] + Data.loc[i, 'NoiseTrs']

# plt.plot(Data['x'])
# plt.show()

Data['y'] = Data['x'] + Data['NoiseObs']

# plt.plot(Data[['x', 'y']])
# plt.show()

F = [1.]                                               # transition matrix
H = Data['x'].values.reshape(N_data,1,1).astype(float) # observation matrix
Q = [1.]                                               # transition noise
R = [1.]                                               # observation noise

#%%
init_state_mean = [3.]
init_state_cov = [1.]

kf = KalmanFilter(
    transition_matrices=F,
    observation_matrices=H,
    # transition_covariance=Q,
    # observation_covariance=R,
    initial_state_mean=init_state_mean,
    initial_state_covariance=init_state_cov,
    em_vars=['transition_covariance', 'observation_covariance']
    # em_vars=['observation_covariance']
    # em_vars=['transition_covariance']
)

kf = kf.em(Data['y'].values.astype(float))

print("observation_covariance R :", kf.observation_covariance)
print("transition_covariance Q :", kf.transition_covariance)

```

Actually, I have found something. Indeed, it is a complicated problem to estimate transition and observation covariances.

In Wikedia it says:

https://en.wikipedia.org/wiki/Kalman_filter#Estimation_of_the_noise_covariances_Qk_and_Rk

  Practical implementation of the Kalman Filter is often difficult due to the difficulty of getting a good estimate of the noise covariance matrices $Q_{k}$ and $R_{k}$. Extensive research has been done in this field to estimate these covariances from data. One practical approach to do this is the autocovariance least-squares (ALS) technique that uses the time-lagged autocovariances of routine operating data to estimate the covariances. The GNU Octave and Matlab code used to calculate the noise covariance matrices using the ALS technique is available online under the GNU General Public License.

Also, similar things are discussed here:

Kalman Filter Covariance

and here

https://quant.stackexchange.com/questions/8501/how-to-tune-kalman-filters-parameter


Creating a Phantom Sound through 5.1.2 Speaker Configuration I am trying to use a 5.1.2 speaker configuration (left front, center, right front, left rear, right rear, left front height, right front height and Low-Frequency) to create a phantom sound at a certain location. I am looking for an algorithm which takes the 3D coordinate and current velocity of the emitter to calculate sound level at each speaker.

I have been doing research on this and what I found was use Head Related Transfer Function to calculate the sound pressure at right / left pinna, and use cross-talk cancellation to eliminate cross-talk effect from multiple speakers. However, I suppose HRTF is meant for creating virtual sound from headphones which does not help from my perspective. Can someone give some ideas on this? Thanks!

I believe that Crosstalk Cancellation is not the method you are looking for (based on the fact that you are using more than 2 loudspeakers for sound reproduction).  You have not described exactly what the problem is that you are currently experiencing, and I therefore provide you with the following as general advice and guidance:

A binaural recording is generated utilising HRTFs that describe a particular listeners head's filtering (due to head, pinna, torso etc) of the original sound source. When played back over a pair of loudspeakers (or in your case many) the crosstalk paths from loudspeakers to the actual listeners ears degrade the virtual image, which is why Crosstalk Cancellation filters are typically employed (to overcome this problem of "crosstalk" from each loudspeaker to each ear).  A great, simple, explanation of this can be found here http://pcfarina.eng.unipr.it/Aurora/crostalk.htm, however if you wish to delve deaper into Crosstalk Cancellation theory then I would recommend this paper for you: http://asa.scitation.org/doi/abs/10.1121/1.3206660

If you are extending your playback system to a larger number of loudspeakers, then depending on your desired listening position/area you may wish to consider alternative soundfield reproduction methods such as Ambisonics or Wave Field Synthesis.  You therefore don't have to rely on an individuals HRTF, nor Crosstalk Cancellation etc...  Since you have an irregular loudspeaker placement you may wish to consult research such as https://www.researchgate.net/profile/Yukio_Iwaya2/publication/268328126_High_order_Ambisonic_decoding_method_for_irregular_loudspeaker_arrays/links/54f5044e0cf2ba6150645b70.pdf for more information on how to design your playback algorithms.

I hope this is the help you were after?!


Estimating velocity from known position and acceleration I am stuck at modeling a system model, i.e. getting my state vector and input vector. My guess is that position and velocity are state vector and acceleration is input vector. My 2nd guess is that all three quantities are in state vector and none in input vector.

So... what is state vector and what is input vector in my case?

--

Additional info:

I get measurements from position sensor and acceleration sensor. Everything is happening in 1D, for example on a straight line. I want to merge these readings (and remove the noise) to get an estimation of velocity for each timestep.

These equations describe the system; I am not sure if they're modeled right though. If I understand correctly it's safe to predict that acceleration is constant (even though in reality it changes) - because process covariance matrix fixes this assumption (right?).

I have also some sample data to work with (input values aren't noised here for simplicity):

 time    pos     acc      what I should get as output (velocity)
[0.0s]  0.000, -0.000  | 18.850
[0.1s]  1.885, -0.113  | 18.850
[0.2s]  3.768, -0.227  | 18.839
[0.3s]  5.650, -0.340  | 18.816
[0.4s]  7.528, -0.452  | 18.782
[0.5s]  9.401, -0.565  | 18.737

ADDITION 2:

For better communication I'm creating a new answer but should be treated as a comment to the first answer. Jason you've already helped me tremendously and I really am grateful for your time. I still have problems with this though - the results from Kalman Filter are not as expected. May you find the time please review the following, thanks. I already owe you a beer or two (or coffies if you like) - if you have paypal contact me on primoz[a t]codehunter.eu :)

I've implemented the model Jason had proposed in first answer. I added the jerk as 4th state variable. After hours of reviewing I decided to come back here for help. The values I get out of KF aren't as expected. Table below represents the data from first 10 iterations of algorithm. Notice how jerk is increasing each time step thus making other estimates wrong. After one second the difference between real acceleration and estimated is more than 1m/s² (see table, last row)!

           real           measured                   estimated                    real
time   pos     acc       pos    acc        pos      acc    jerk    vel[!]       velocity
0.0   0.000  -0.000    -0.040  0.030  |  -0.300   -0.060   0.000   18.850  <-->  18.850
0.1   1.885  -0.113     1.965 -0.153  |   1.585   -0.061  -0.006   18.844  <-->  18.844
0.2   3.768  -0.227     3.778 -0.247  |   3.469   -0.066  -0.035   18.835  <-->  18.827
0.3   5.650  -0.340     5.750 -0.370  |   5.351   -0.090  -0.122   18.815  <-->  18.799
0.4   7.528  -0.452     7.358 -0.452  |   7.228   -0.152  -0.291   18.769  <-->  18.759
0.5   9.401  -0.565     9.251 -0.555  |   9.094   -0.282  -0.574   18.673  <-->  18.708
0.6   11.269 -0.677    11.309 -0.717  |   10.938  -0.518  -1.006   18.494  <-->  18.646
0.7   13.130 -0.788    13.260 -0.758  |   12.752  -0.840  -1.490   18.233  <-->  18.573
0.8   14.983 -0.899    15.043 -0.949  |   14.520  -1.286  -2.096   17.854  <-->  18.488
0.9   16.827 -1.009    16.977 -1.089  |   16.235  -1.838  -2.770   17.362  <-->  18.393
1.0   18.661 -1.118    18.831 -1.168  |   17.890  -2.477  -3.476   16.762  <-->  18.287

My matrices are here:

What is causing this addition in each timestep for jerk? Is any of my matrices wrong?

The same goes with the first solution (only 3 state model) - the acceleration isn't changing as it should.

LAST EDIT:

I've finally managed to make it work. I am not sure if there was an implementation error or wrong P&Q matrices.

You would not make a constant-acceleration assumption in this case. You would typically do that if you had no means of measuring the system's acceleration, but you say that it is observable for your case. The most obvious way to model this system would be using the state vector

$$
\mathbf{x_k} = \left [
\begin{array}{c}x_k \\
\dot{x}_k \\
\ddot{x}_k
\end{array} \right ] = \left [
\begin{array}{c}x_k \\
v_k \\
a_k
\end{array} \right ]
$$

where $x_k$ is the position, $v_k$ is the velocity, and $a_k$ is the acceleration of the system at time instant $k$. Since you say that you can measure the position and acceleration of the system, the measurement vector $\mathbf{z_k}$ would be:

$$
\mathbf{z_k} = \mathbf{Hx_k} + \mathbf{v_k}
$$

with

$$
\mathbf{H} = \left[
\begin{array}{cc}
1 && 0 && 0 \\
0 && 0 && 1
\end{array}
\right]
$$

resulting in the measurement model:

$$
\mathbf{z_k} = \left[
\begin{array}{c}
x_k \\
a_k
\end{array}
\right] + \mathbf{v_k}
$$

where $\mathbf{v_k}$ is the (Gaussian) noise inherent in the measurement. Now, unless you are adding some known force to the system that would affect its dynamics, then the input vector $\mathbf{u_k} = \mathbf{0}$. See this recent question for more detailed discussion of the various terms in the Kalman model, which you need to know, and which you don't.

Edit: Regarding the below comment on how to handle the state transition matrix for the acceleration term: there are a couple different ways to handle that. You could assume a zero-jerk (jerk is the time derivative of acceleration) model; this is equivalent to assuming that $a_{k+1}=a_k$. Or, you could add a fourth element to your state vector:

$$
\mathbf{x_k} = \left [
\begin{array}{c}x_k \\
\dot{x}_k \\
\ddot{x}_k \\
\dddot{x_k}
\end{array} \right ] = \left [
\begin{array}{c}x_k \\
v_k \\
a_k \\
j_k
\end{array} \right ]
$$

where $j_k$ is the system's jerk, or time rate of change of acceleration, at time instant $k$. Since you aren't measuring jerk directly, you would have a measurement matrix of:

$$
\mathbf{H} = \left[
\begin{array}{cc}
1 && 0 && 0 && 0\\
0 && 0 && 1 && 0
\end{array}
\right]
$$

And the state transition matrix would become:

$$
\mathbf{F} = \left[
\begin{array}{cccc}
1 && dt && \frac{1}{2}dt^2 && \frac{1}{6} dt^3 \\
0 && 1 && dt && \frac{1}{2}dt^2 \\
0 && 0 && 1 && dt \\
0 && 0 && 0 && 1
\end{array}
\right]
$$

This is a "constant-jerk" model (similar to constant-velocity or constant-acceleration models you might see in simple examples); your state transition matrix makes an implicit assumption that the jerk is constant for all values of k. This is not likely to be true. To handle that aspect of the problem, you would introduce a process noise term in the jerk component of the state transition equation.

Qualitatively, this allows you to express that you aren't sure how the jerk term is going to change from time step to time step, but that you expect the changes to be random with some Gaussian distribution. This is a tool often used to "tune" adaptive filter models like this one until you find a set of parameters that works well for your application.


what is the best technique to detect a given signal : energy detection, autocorrelation, and matched filter based sensing, CRC? In order to detect the presence of a signal impaired by noise, what is the best technique (best metrics) one can use ?

energy detection,
autocorrelation,
matched filter based sensing,
CRC
data aided or not ..
an explanation of the underlying algorithms would also be appreciated

There is no best rule, it depends on what we know about the signal.

Energy Detector: The energy Detector makes a decision on presence or absence of a signal based on sum of squared samples. This rule comes from the fact that the signal to be detected is inherently assumed to be random following a wide sense stationary process with a known PDF, the likelihood ratio test then simply determines this rule.

Matched filtering and autocorrelation: The matched filtering is nothing but the correlation of the. Received samples with the known signal itslef, so they are inherently the same, the difference compared to energy Detector is that here the signal is assumed to be a known deterministic quantity. The matched filter is the optimal detector in low SINR scenarios.

CRC data: A check on CRC is an agreement between the receiver and tarsmitter in whether an error has occurred or not while transmitting data on a wired/optical/Wireless channels. By itself it is not a signal detection scheme. It only detects whether an error has occurred or not. Once an error is detected the receiver can then request the trasnmitter to send the data again by sending a NACK (not acknowledge) signal to the trasnmitter (there are other ways to handle it as well, like soft combining etc.)


When are two signals orthogonal? The classical definition of orthogonality in linear algebra is that two vectors are orthogonal, if their inner product is zero.

I thought this definition might be applied to signals as well, but then I thought about the following example:

Consider a signal in the form of a sine-wave, and another signal in the form of a cosine-wave. If I sample both of them, I obtain two vectors. While sine and cosine are orthogonal functions, the product of the sampled vectors is almost never zero, nor does their cross-correlation function at t=0 vanish.

So how, then, is orthogonality defined in this case? Or is my example off?

As you may know, orthogonality depends on the inner product of your vector space. In your question you state that:

  While sine and cosine are orthogonal functions...

This means that you have probably heard of the "standard" inner product for function spaces:

$$\langle f,g\rangle=\int\limits_{x_1}^{x_2} f(x)g(x) \ \mathrm{d}x$$

If you solve this integral for $f(x)=\cos (x)$ and $g(x)=\sin(x)$ for a single period, the result will be $0$: they are orthogonal.

Sampling these signals, however, is not related to orthogonality or anything. The "vectors" you obtain when you sample a signal are just values put together that make sense to you: they are not strictly vectors, they are just arrays (in programming slang). The fact that we call them vectors in MATLAB or any other programming language can be confusing.

It's a bit tricky, actually, since one could define a vector space of dimension $N$ if you have $N$ samples for each signal, where those arrays would indeed be actual vectors. But those would define different things.

For simplicity, let's suppose we are in vector space $\mathbb{R}^3$ and you have $3$ samples for each signal, and all of them are real-valued. In the first case, a vector (i.e. three numbers put together) would refer to a position in space. In the second one, they refer to three values a signal reaches at three different times. In this example it is easy to spot the difference. If you had $n$ samples, then the notion of "space" would be less intuitive, but the idea still holds.

In a nutshell, two signals are orthogonal if the inner product between them (namely, the integral I wrote above) is $0$, and the vectors/arrays obtained by sampling them tell us nothing about their being orthogonal.


Time-varying waveform My situation is as follows: i am trying to generate a waveform the hard way, by constructing the samples one by one and then saving the result to a .wav file using Python.

When the frequency is constant, everything is fine: i use $y(t) = \sin(2 \pi \cdot f \cdot t)$.
However, if i change the frequency to be a function of time, things go wrong. If the function is a linear one of the form $f(t) = a + bt$, it still works. But if i choose, for example, $f(t) = 40 + 10 \sin(t)$, the max frequency increases over time, reaching a maximum higher than the expect 50Hz.

I have read something about instantaneous frequency, namely, this: Why does a wave continuously decreasing in frequency start increasing its frequency past the half of its length?. But doing the integral evaluation makes the sound even weirder. And the method i currently have works for linear function of time, so i figured there must be something else wrong.

I also tried to generate chunks of sound in the frequency i need, at each time, and then glue them together. I calculated the period of the oscillation, so that a chunk has as many samples as necessary to make a whole period on that frequency, so that there are no "jumps" between different frequencies. But generates a cracking sound in the sample.

Here is an example of the kind of a JavaScript implementation of the kind of sound i need: http://jsfiddle.net/m7US6/4/.

Any help would be appreciated. Thanks.

Consider the function $\sin(2\pi f t)$. When $2\pi f t$ goes from $0$ to $2\pi$ you get oscillation of the sine wave, $2\pi$ to $4\pi$, another, and so on. So every time the argument changes by $2\pi$ you get one oscillation.

Now lets plot $2\pi f(t) t$ where $f(t) = 10 + 10 \sin(2 \pi t)$. (I've modified it a bit to show the effect more):

As $t$ increases, this value changes faster and faster, meaning $\sin(2\pi f(t) t)$ will have higher and higher frequency.

Another way of thinking about it - the frequency is the derivative of $f(t) t$. For the example this is $20 \pi (\sin(2 \pi t)+2 \pi \mathbf{t} \cos(2 \pi t)+1)$. Note the $t$ multiplier in there - $t$ increases thus the frequency increases.

So what you really want is
$\tfrac{d}{dt} f(t) t = 40 + 10 \sin(t)$ which gives $f(t) t = 40t - 10\cos(t) + c$ and for your example try $ \sin(2\pi(40t - 10 \cos(t)))$


Do waterfalls generate pink (1/f) or brownian (1/f^2) noise? Does anybody know if natural waterfalls generate pink noise or brownian noise? I have tried to search about this but the matches are non-scientific and contradictory.

For example

Livescience.com mention Brownian noise to resemble big waterfall.
Soundonsound.com magazine says Pink noise to occur in nature, e.g. waterfall.

Even though it feels weird, it could be that the type of noise is dependent on the size of waterfall, for example by affecting to alpha of the power spectral density of the broad definition of pink noise: $S(f) = 1 / f^{\alpha}.$

Thank you in advance.

Interesting question. I juts happen to live right outside a waterfall so I can relate some hands on experience.

The sound depends A LOT and water flow and weather conditions. Water flow is the main engine but wind can generate significant modulations
Depending on conditions you can get strong harmonics in the signal. Sometimes you can actually see standing waves running up and down the fall. These tend to be very low in frequency: maybe around 1 Hz or so.
The fall is pretty wide and has strong directivity. So the spectrum heavily depends where you are.
Overall it sounds more brown than pink
Fun fact: if the fall hits the right resonance, it noticeably shakes the house and the dishes in the kitchen cabinet start to rattle.

This is just for one specific fall and I'm sure others behave quite differently. So all in all any simple model would be a gross over-simplification.


Blind source separation on images - with a known source A camera outputs an image $A$ that is a linear combination of two uncorrelated images $B$ and $C$:

$$A = B + kC$$

Images $A$ and $C$ are known; the unknowns here are $B$ and a scalar coefficient $k$ that varies the intensity of $C$.

Is there a way to retrieve $B$? I could use $A$ and $C$ as separate sources and apply the ICA algorithm, but I think that it's an overkill in this case, and there should be something simpler. Any suggestion? (I'm using OpenCV to process the images, but I can implement non-OpenCV algorithms if I need to.)

You already know one of the sources here, so you don't need ICA.

There are infinitely many couples $k$ and $B$ that verify:
$$
A = B + kC,
$$
 but in your case if you assume that $B$ and $C$ are uncorrelated, we can choose the solution that minimizes the correlation between $B$ and $C$.

Using the previous equation and the linearity of the correlation operator, we have:
$$
\text{corr}(A,C) = \text{corr}(B,C) + k \text{corr}(C,C)
$$

So you need to find $k$ that minimizes $(\text{corr}(A,C) - k \text{corr}(C,C))^2$, which is simply:

$$
k = \frac{\text{corr}(A,C)}{\text{corr}(C,C)}.
$$
And finally:
$$
B = A- \frac{\text{corr}(A,C)}{\text{corr}(C,C)}C
$$


What is the relationship between poles and system stability? I see two notions that describe the relationship between poles and system stability. But they are not the same from my understanding

The system is BIBO stable if and only if all the poles are in the
left half of the complex plane
A LTI system with a rational system function H(z) is stable if and only if all of the poles of H(z) lie inside the unit circle.

Why these two notions are different? Is that in different conditions?

I agree with Peter K.'s answer, but I would like to add one important point: the two statements in the question are only true for causal systems. The most general statement about stability for LTI systems described by rational transfer function is:

  An LTI system with a rational transfer function is stable if the region of convergence (ROC) of its transfer function includes the $j\omega$-axis (for continuous-time systems), or the unit circle $|z|=1$ (for discrete-time systems).

Since for causal systems the ROC is to the right of the pole with the largest magnitude (or outside the circle with radius equal to the largest pole radius in the $z$-plane), the above statement coincides with the statements in the question. However, non-causal systems can have poles in the right half-plane (or outside the unit circle of the $z$-plane) and still be stable.

Take as an example the discrete-time transfer function

$$H(z)=\frac{z^2}{(z-a)(z-b)},\quad a<1<b\tag{1}$$

which has one pole inside the unit circle and one pole outside the unit circle. You can't say that the corresponding system is unstable, because it depends on which system you mean. The transfer function (1) has 3 possible ROCs, corresponding to 3 different systems:

$|z|>b$
$|z|<a$
$a<|z|<b$

The first one corresponds to a causal system which is unstable, because one of its poles is outside the unit circle (the ROC does not include the unit circle). The second one corresponds to an anti-causal system which is also unstable, because it has one pole inside the unit circle (again, the ROC does not include the unit circle). Finally, the third ROC includes the unit circle, i.e. it corresponds to a stable system. It is a non-causal system with a two-sided impulse response which decays exponentially as the time index $n$ goes to $\pm \infty$.


Estimating an Image Filter from Examples I am sure this question must be fairly easy and could have been answered on this site but perhaps I am not finding the right keywords to search for it.

I would like to know what are the common techniques to estimate an image filter from a set of input (unfiltered) and post-filtered images. From a machine learning point of view, this seems easy to do under some assumptions and constraints. Assuming the filters are linear, we might have to estimate the parameters of a convolution kernel using an optimization algorithm such as gradient descent. A similar approach may be done for non-linear kernels assuming a suitable non-linear model. However, I have still not found lots of literature discussing this.

As such, my questions are the following: How I can perform this operation in practice, how it is typically referenced in the image processing literature, and where can I find more information about it?

I will talk about the Linear Spatially Invariant (LSI) Model.

Since the Filter (System) is LSI it can be described in the domain of Fourier Transform:

$$ Y \left( \omega \right) = H \left( \omega \right) X \left( \omega \right) $$

Where $ Y \left( \omega \right) $ is the Blurred (Output) Image, $ X \left( \omega \right) $ is the Input Image and $ H \left( \omega \right) $ is the filter model.

Now, assuming the filter doesn't vanish in the Fourier Domain, then:

$$ H \left( \omega \right) = \frac{Y \left( \omega \right)}{X \left( \omega \right)} $$

Now, to practical issues with this model:

Noise.
Vanishing Frequencies.

There are many methods to deal with those issues (Pseudo Inverse, Wiener Filter, Lucy Richardson) all are having different assumptions (White Noise, Gaussian Noise, Linear Model, Poisson Noise, etc...).

Till now I described the classic "Signal Processing" Model.

What you can do is build this model in the Spatial Domain.
The Filtering Operator will be a very large Toeplitz Matrix which you can estimate using LS Model or any solver for any cost function you'd want.

Practically in the world of Image Processing, the quality is in the regularization yo enforce.

P.S.
There is a new community proposal in the SE network:

http://area51.stackexchange.com/proposals/86832

It would be great if you could follow, up vote question with less than 10 points and propose your own questions.

Thank You.


Remove Frequency Range I got a signal and I would like to remove a specific frequency range.

I am a frequency noob, so I will explain it by images:
I have a frequency filter (butter bandpass):

I have a signal and it's filtered signal:

The frequency spectrum of the signal looks something like this:

And the extracted frequencies like this:

So far everything is fine. But if I try to subtract the extracted signal from the original signal in time domain, then I get this frequency spectrum:

Instead, I would like to get the red marked spectrum.

What do I need to calculate in order to remove the specified frequencies from a signal?

I got the code from here and added FFT computations.

Kind regards,
Thomy800

The problem is most likely that you aren't accounting for the group delay of the filter. As an LTI system, your bandpass filter does two things to each frequency component present in its input signal:

It multiplies the frequency component's magnitude by some value $A(\omega)$.
It shifts the frequency component's phase by some phase angle $\theta(\omega)$.

Together, these two components make up the filter's frequency response. The second of these two effects is what yields the filter's group delay. That is, it takes some time for each component of the input signal to propagate through the filter. If you are looking to do some operation on the filter input and output that relies upon them being time-aligned, you have to account for this effect first.

With that said, doing such a thing with an IIR filter such as yours can be tricky. In general, IIR filters do not have linear phase. This means that your filter will delay some frequency components more than others. Again, if you're relying on time alignment of the filter input and output to do some operation, then this is going to work against you.

An alternative would be to use a symmetric FIR filter. These filters do have linear phase, so it will delay all frequency components by an equal number of samples (the delay $D$ is equal to half of the filter order, so if you have $N$ total coefficients, the delay is $\frac{N-1}{2}$  samples). If you did this, then you could accomplish what you want by doing the following:

Pass your input signal $x[n]$ through the filter to get an output $y[n]$.
Calculate the difference signal that you're looking for by adjusting your indexing according to the filter's group delay:

$$d[n] = x[n] - y\left[n + \frac{N-1}{2}\right]$$

Addendum: All of this begs the question of precisely what it is you're trying to do anyway. It seems like maybe you're trying to implement a bandstop filter by using a bandpass filter to isolate one frequency band and then subtract it from the original. It might make more sense to use a bandstop filter directly. For example, techniques exist to take a lowpass prototype filter and transform it to other types, like highpass, bandpass, or bandstop. This process might have been used to design your original bandpass filter anyway.


Why Does 2D FFT of Gaussian Looks More Sharper than Gaussian Itself? I am trying to understand why 2D FFT is done on a Gaussian process in a particular code. From my understanding from these posts:

https://www.researchgate.net/post/Why_would_we_want_to_do_Fourier_transform_of_a_Gaussian_signal
Why use a gaussian low pass filter, when we can convolve with gaussian kernel
Zero padding effect on a FFT of gaussian noise

It appears that it is used to

reduce noise
reduce the sharpness on the edges

I am working with MATLAB, and the code snippet is as follows. Please note the peak has been moved to the top-left in the gaussian_shaped_label function, with wrap-around:

%window size, taking padding into account
window_sz = floor(target_sz * (1 + padding));

%   we could choose a size that is a power of two, for better FFT
%   performance. in practice it is slower, due to the larger window size.
%   window_sz = 2 .^ nextpow2(window_sz);

%create regression labels, gaussian shaped, with a bandwidth proportional to target size
output_sigma = sqrt(prod(target_sz)) * output_sigma_factor / cell_size;
yf = fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

I was trying to plot the two results.

The Gaussian (left) zoomed on the top-left edge(right):

Plotting the code:

  gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

The absolute value of the 2D FFT of the Gaussian mentioned in point 1 (left) zoomed on the top-left(right):
Plotting the code-

fft2(gaussian_shaped_labels(output_sigma, floor(window_sz / cell_size)));

The plot of 2D FFT seems sharper than the Gaussian. Shouldn't 2D FFT actually blur the edges and make it smoother or something? I think I am going wrong somewhere and any help will be appreciated.

P. S.
If you can help me understand the mathematically too, it would extra helpful but not a hard requirement here.

Why Does 2D FFT of Gaussian Looks More Sharper than Gaussian Itself?

Have a look at the Fourier Transfrom of a Gaussian Signal.

$$ \mathcal{F}_{x} \left\{ {e}^{-a {x}^{2} } \right\} \left( \omega \right) = \sqrt{\frac{\pi}{a}} {e}^{- {\pi}^{2} \frac{ {\omega}^{2} }{a} } $$

First, Gaussian Signal stays Gaussian under Fourier Transform.

As you can see, the parameter which multiplies the variable is inverted.
Let's say $ a = 5 $, then it means that in time we will have very sharp and thin Gaussian while in frequency we will have very smooth and wide Gaussian.

This is related to other property of Forier Transform.
In simple words, what's thin on Time / Spatial Domain is Wide on Frequency Doman and vise versa.
This is basically the Uncertainty Principle of Fourier Transform.


What is the driving Voltage? While I was looking into ultrasonic transducer in order to transmit and receive a beam, I came across specification of a transducer stating driving voltage as 20Vrms. I am not sure what actually does it mean and where is the application of driving voltage in transducer(I mean is it used at transmitting or receiving end? or what is its function?). I wonder what parameters of the transducer will be affected due to application of driving voltage?
Below is the specification of the transducer.

20Vrms is the maximum voltage you can apply across the transmitting transducer without the risk of immediately damaging it.
The amplitude of sound it produces is determined by the driving voltage. The transmitter is characterized at 10Vrms, so about 28Vp-p assuming a sine wave, probably where you would prefer to use it for reliability and long life.


Extra energy when altering frequency amplitudes How can I convert a time signal (audio) into the frequency domain, multiply frequency amplitudes by a vector of scalars, and convert it back into the time domain, without messing up timing information (ringing?)?

Some context:

I am calibrating ultrasonic speakers. I have a frequency roll-off curve of the attenuation dB for each frequency, from a reference point. I now want to compensate my output using this curve.  Here is an example that was gathered experimentally:

I can then multiply my output signal in the frequency domain by this attenuation vector, and convert back into the time domain. For the most part this works, however, there is some extra energy at the wrong times, and I am not sure why. Here is what it looks like:

Here is some python code to demonstrate what I have. To make this a compilable example I try to approximate the attenuation vector with a sawtooth for variation.

    import numpy as np
    from scipy.signal import chirp, sawtooth

    def smooth(x,window_len=11):
        """smooth the data using a window with requested size."""

        s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]

        w = np.hanning(window_len)
        y=np.convolve(w/w.sum(),s,mode='valid')
        return y[window_len/2:len(y)-window_len/2]

    fs = 5e5
    duration = 0.2
    npts = duration*fs
    t = np.arange(npts).astype(float)/fs

    # my output stimulus
    x = chirp(t, f0=5000, f1=1e5, t1=duration)

    X = np.fft.rfft(x)

    # stand-in attenuation vector
    variation = sawtooth(np.arange(npts/2 + 1)/500)
    attendB = np.linspace(0, 30, npts/2 + 1) + variation

    # attempt to minimize ringing?
    attendB = smooth(attendB)

    H = 10**((attendB).astype(float)/20)

    # signal adjusted to calibration
    A = X * H
    a = np.fft.irfft(A)

    # ==== plot ====
    import matplotlib.pyplot as plt

    plt.subplot(211)
    plt.specgram(x, Fs=fs)
    plt.title('original')
    plt.subplot(212)
    plt.specgram(a, Fs=fs)
    plt.title('adjusted')
    plt.show()

I can eliminate the extra energy by increasing the smoothing. However, this also reduces the accuracy of the calibration, as assesed by comparing the desired and recieved frequency response. Is this a trade-off I have to deal with, or am I going about this the wrong way? Alternatively, is it a feasible solution to implement this as a filter in the time domain?

Try zero-padding your data by a large amount (approximating the length of the impulse response associated with your frequency domain compensation), and then use a longer FFT so that circular convolution artifacts don't wrap around in time.

A smoother frequency compensation curve will usually have a shorter impulse response, thus requiring less zero-padding for a given noise floor.


What is the best estimate of a signal given multiple noisy realisations I assume that the noise corrupting the sought-for signal is [Edit: Gaussian] white noise of zero mean and unknown variance. Is the optimal solution simply the mean over the realisations (the ensemble mean), or, we can have a better estimate, say, in terms of a smaller root-mean-square-error?
Edit: I found that my question have been asked already. Although the practical context is different, as in my case there is no concern that the different observations on different sensors could be correlated or that the corrupting process is slightly different. The answer to that question suggests that superposition, i.e., ensemble mean, is the best (“efficient"?) estimate. However, it occurs to me that it might not be the case, as a moving-window temporal average of the ensemble mean would yield an even smaller RMSE if the signal is sufficiently smooth. Unfortunately, there was no published reference given there. It would be great to find the answer in a textbook or paper, as it would surely clarify the conditions under which the result/theorem applies.
In any case, even if the moving average of the ensemble mean yields a better RMSE than the ensemble mean alone, I wonder if there is an even better estimate. Intuitively -- and it might be a common failure of intuition -- there might be value (information) in having the relaisations available "separately" as opposed to only their superposition, the ensemble mean. Any further help would be much appreciated.

The term you're looking for is efficient estimator. That's the least-variance (== least RMS) estimator for a given amount of observation.
Whether or not that's the average depends on your noise distribution.
For example, consider the two noises:

continuous-valued uniformly distributed noise
Noise that only takes values $\in \{-\sigma, 0, \sigma\}$ with equal probability.

Both can have the same variance, both are zero mean. Both can be white. But while the average sounds like a good idea for the first one, for the second one you can get a 0 RMS estimator as soon as you've observed three different values by picking the middle one.


Filter long audio data using cascading biquad filter I would like to lowpass filter a long pcm wav file. I designed a 20 order butterworth lowpass filter in Octave and transform the zero pole gain data to a series of second order sections. Then I am applying these SOS filters to each processing block on the long audio. The script is like this:

clear all; close all;
[y, fs, bit] = wavread('sweep.wav');

# design butter
[z, p, k] = butter(20, 1200 / (fs / 2), 'low');
[sos, k] = zp2sos(z, p, k);

# frame index
framesPerBuffer = 512;
currentFrame = 1;

# final condition buffer
zbuf(2) = 0;

while currentFrame + framesPerBuffer < size(y, 1)
    # read one frame
    data = y(currentFrame:currentFrame + framesPerBuffer - 1, 1);
    for i = 1:size(sos, 1)
        data = filter(sos(i, 1:3), sos(i, 4:6), data);
    end
    y(currentFrame:currentFrame + framesPerBuffer - 1, 1) = data;
    currentFrame += framesPerBuffer;
end
# apply gain to filtered signal
y = y * k;

The "sweep.wav" input is a mono pcm file with 48kHz sample rate and 16 bit depth with duration of 40 seconds. It is a step sweep signal with each sine tone duration about 1 second, frequency range from 200-12000 Hz. The overall wave shape is like this:

After the filtering process above the output signal has the wavform of this:

Since butterworth filter has the flat passband/stopband, the waveform is clearly not correct. When I listen to the filtered wave it is like a filtered wav plus some low frequency modulation.

I have also tried to add final condition z buffer for the filter function but it is still perform bad. Is there anything wrong for me filter process? How can I get the correct filtered data out of this long PCM? Thank you very much!

You are not storing the filter states between buffers.

In this case, each filter is a 2nd order filter, which has 2 state variables (the values stored in the 2 delays implementing each filter).

So for each new frame, you are starting with state 0, introducing a glitch between frames.

You must store the resulting state for each filter, and apply it for each next frame.

Built-in Function: [Y, SF] = filter (B, A, X, SI)

 If the fourth argument SI is provided, it is taken as the initial
 state of the system and the final state is returned as SF.  The
 state vector is a column vector whose length is equal to the
 length of the longest coefficient vector minus one.  If SI is not
 supplied, the initial state vector is set to all zeros.

So for each filter, you must store the corresponding SF, and provide as SI on the next frame.


Modeling Load Modulation I'm trying to create a simple NFC/RFID software simulator and am having a hard time understanding the math fully. I know that NFC tags respond by load modulating the 13.56 MHz wave that the reader is sending it. I'm pretty sure NFC tags use a subcarrier wave with a frequency of ~847 kHz. Where I'm getting confused is how to model the result of the load modulation.

Can I simply take the 13.56 MHz wave and do a point-wise multiplication with the 847 kHz wave? Using the Manchester Encoder, this nets me a signal that looks like this:

After measuring a NFC tag with a USRP and making a couple of simulations in octave and python, it seems like I was incorrect regarding the point-wise multiplication. As it stands, I believe the best simulation is to do an addition of the carrier wave (13.56 MHz) and the subcarrier wave (847 kHz), with the subcarrier wave's amplitude set at 0.1x the amplitude of the carrier wave.


Count the number of times the phone rings before its answered I have a dataset of phone recordings that I'd like to do some characteristic-level analysis on.  The audio is sampled at 44khz stereo, with each channel containing the signal from one party in a conversation.

I would like to start by counting the number of time the call rings before the phone is picked up.  I believe this to be simple, but am a little lost, so any pointers in the right direction would be much appreciated.  I believe I need some sort of a matched filter, or maybe a bandwidth filter.  The signal is generally clean -- I have attached a couple of graphs of the signal.

General solutions are appreciated, but I will be doing this in python so anything specific to that would be extra appreciated.

Additionally, if you have other ideas for extracting features using signal characteristics, I'm all ears (eyes).

If you can identify one phone ring, then you can extract it and use it to construct a matched filter. The beginning of the wiki article is math heavy, but check out the example section. Additionally, Matlab has a page describing the matched filter from one of their toolboxes. You should really give it a try; it works like magic.

In the time domain, the matched filter is a cross-correlation of the signal against a known pattern (which is in this case a recording of the ringing). As a cross correlation it can efficiently be implemented using a DFT. The output of the matched filter is a time series.

The time series output will have a series of peaks whose locations correspond to the beginning of each ring. All you need to do is count the number of peaks which exceed some threshold.

If the rings could be different but you are guaranteed that they're selected from a finite set of rings (imagine different phone carriers or something) then you could use a bank of matched filters.