Basic PLL theory: voltage-controlled oscillator confusion I am trying to write my first PLL in Python, and I've been making use of this excellent tutorial:

http://www.pythonpowerelectronics.com/contents/tutorials/tutorial1/pll_tutorial.pdf

Reading through the math on pages 7 and 8 (very minimal- it's just a few lines), I have one point of confusion:

Let's assume, like they do in the tutorial, that the reference (power grid) voltage signal to be tracked is $V_{r}\sin(\omega_{r}t)$ and that our VCO voltage can be expressed with the two components $\sin(\omega_{v}t)$ and $\cos(\omega_{v}t)$. We integrate the incoming grid voltage to get a cosine term of $V_{r}\cos(\omega_{r}t)$, and then combine our four sinusoidal terms together in a way that reduces to $V_{r}\sin(\omega_{r}t-\omega_{v}t)$. According to the tutorial, we feed that function into the VCO, which then changes $\omega_{v}$ until the frequencies are equal and the term in the sine is zero for all $t$.

That all seems well and good when the initial phase of the incoming reference wave is 0. In that case, $V_{r}\sin(\omega_{r}t-\omega_{v}t)$ is positive when $\omega_{v}$ is less than $\omega_{r}$ and it is negative when $\omega_{r}$ is less than $\omega_{v}$. But what if the grid voltage is actually at some initial phase offset relative to the VCO? Let's say the true equation of our input to the VCO is, for example, $V_{r}\sin(\omega_{r}t-\omega_{v}t+\pi)$. Then $\omega_{v}$ could be less than $\omega_{r}$ and $V_{r}\sin(\omega_{r}t-\omega_{v}t+\pi)$ would still be negative. The VCO would falsely interpret this as a "you're ahead of the reference frequency, slow down" command, $\omega_{r}$ and $\omega_{v}$ would separate even further over the course of the wave's negative half-cycle, and I don't think a PI frequency controller would ever recover from that.

I'm not seeing a way around this problem. We don't have access to the actual frequency difference, only to the value of the combination sine function we've created. We also can't make any claims about the initial phase angle of the reference voltage, since we could have started our PLL at any point in its cycle and it could have any initial phase relationship to our VCO. If anyone could help me understand how we get past this issue, it would be a huge help. Thanks in advance!

EDIT: PLL block diagrams such as the one in the link below have a low-pass filter on the input to the VCO, but I don't believe that has anything to do with my question. Unless I'm mistaken, that LPF is just there to remove the high frequency component $\omega_{r}+\omega_{v}$ that results from the phase detector, i.e. the "combine VCO and reference signal" block. The LPF then leaves us with a pure (in the ideal case) sine wave $V_{r}\sin(\omega_{r}t-\omega_{v}t)$. If I'm wrong about that, though, please let me know!

https://www.google.com/search?q=pll+diagram&tbm=isch&source=iu&ictx=1&fir=BxBgiyvYxsuLiM%253A%252CgtPdFL-0ZhSHZM%252C_&usg=AI4_-kQbgUQu5HbOzaiSAkInzcAP8kbYQw&sa=X&ved=2ahUKEwju9dr_ztbdAhWOFTQIHRobD64Q9QEwAHoECAAQBA#imgrc=ZTkMCllE8hHDZM

EDIT 2: Had a helpful exchange with Dan, who gave me some pointers on more efficiently modeling a PLL. He also pointed out that it would be helpful clarify that my confusion is regarding how to implement a PLL. I am assuming that I don't know the phase or the frequency of the incoming reference signal ahead of time. In fact, that's why I'm building a PLL in the first place: I have a sinusoid of unknown phase and frequency, and I am trying to determine both of those quantities by locking a VCO to it and then just querying the VCO.

Here are some helpful tips to proceed with an all digital PLL implementation specific to your desire to determine the phase of your reference by observing the phase of your (cleaner) VCO.

First I would recommend using an NCO instead of a VCO. Here is a previous post where I provided hopefully helpful information on implementing an NCO.

Numerically Controlled Oscillator (NCO) for phasor implementation?

The advantage with this is that you can directly read the instantaneous phase from the phase accumulator (it also implies that you have some master reference of time since the phase accumulator will be continuously counting consistent with the output frequency as $d\phi/dt$.

Next assuming your "capture range" is narrow enough, you can use a simple multiplier as the phase detector, followed by a wide low pass filter for the purposes you describe; a real multiplier will create the sum and difference of the two signals; you would be multiplying the NCO digital output with a digitized representation of your reference signal resulting in a component at the multiplier output that is at twice the reference frequency (sum of reference and NCO specifically) that you want to filter out and the difference of the two which results in the phase error between the two. For more details on such phase detectors see the "Phase Detector" section of this post:

Phase synchronization in BPSK

Also at the bottom of the above link is an all digital carrier recovery loop that would be similar in structure to such a PLL implementation for your purpose. Specifically after the error signal shows the implementation diagram for the Proportional-Integral Loop filter that can be used after the wider filter mentioned above.  The output of the Loop Filter is the control word as the input to the NCO (The NCO is also shown in this diagram as a complex implementation providing both sine and cosine outputs, to simplify your first approach you can simply use a real NCO and real phase detector which is a single multiplier).

Note depending on phase detector used, the loop will likely lock in quadrature to the input signal (if you use a real multiplier specifically this will happen as the phase detector output is 0 when the two inputs are in quadrature).

As far as "capture range" when using the simple multiplier as a phase detector; typically this is referred to the frequency range in which a PLL locks within a single beat-note between the frequencies compared, and the tuning range of the NCO should be limited to stay within this range.   Here is a reference that develops the capture range for a simple XOR gate phase detector (note that an XOR gate is a hard-limited multiplier so has a similar operational range).

http://www.cypress.com/file/60101/download


Why is it necessary to average signals with filters? As someone who is new in working with filters, why do signals need  to be averaged? What are the advantages, by performing filters on signals? As I know, there are different types of filters for signals, is there some sort of schema, which filter for which signal do I have to use?

I am sorry, if my question sounds somewhat ignorant, but I am having a hard time to understand it, especially since it is in another language. I am asking here on this forum, since I know there many experts here, who know more than me. If someone is able to explain it in a compact and easy way to explain, then you have my eternal gratitude.

An average is a summary of different values: if you have several signal processing exams, you know that both your answers and the grading person may incur some variability. On top of that, it is complicated to compare student, or the same student at different years, based on a series of grades.

So an average of grades both:

compacts different grades in one single number,
balances exam variability.

The latter is of first importance in signal processing: a signal sample may be disturbed by noise, sensor fluctuations, quantization, etc. While the signal is supposed to have some regularity, the disturbances are often chaotic, incoherent, from one sample to the other. If you average the distances alone, their amplitudes tend to diminish (not in all cases, but let us say that for now). So if close sample values are averaged, the ones from the signal don't change a lot (like a regular working student), while the disturbances reduce. Thus, the signal looks smoother, more pertinent. This is the basis for the linear filters: further, you can weight the samples differently.

The summarizing property is also important. You can average three samples, and keep the average only, and switch to the other sample triplet.  This reduces the quantity of information, and is the basis for subsampling, data compression, etc. JPEG uses a lot of filtering and subsampling.

Both combined are the root of filter banks, putting several filters in parallel with subsampling. Then, you have nonlinear filters, but this is a whole new story.


Electret microphone RC filter design contradiction among information sources I'm building a circuit for an electret microphone and I want to build a bandpass filter around the op-amp.
I was using this source until I have found that most sources ( e.g. here and here) indicate that RC filter needs to be grounded.
I can't think of a way check this (FFT doesn't seem to indicate whether the first source is correct or not). Any ideas on how to check or information which source is correct?

First, the youtube source shows an 1st order active filter, maxim source shows 2nd order active filter and last source shows first order passive filter. That's three different designs.

A 1st order active filter can be implemented either in the negative feedback branch of the op amp (as in youtube source), this does not need to be grounded.

Or it can be done as a classical RC circuit followed by impedance converter or amplifier. Here, either the R goes to ground (highpass) or the C goes to ground (lowpass)

If you combine these two approaches you get a 2nd order filter. So, wether "the filter needs to be grounded", as you ask, depends strongly on where its implemented.


Impulse Response of a Moving Average System According to Discrete Time Signal Processing by Al Oppenheim 3rd Edition,  the impulse response of a moving average filter is as below

My question is that shouldn't  h[n] be having the dirac delta function instead of the unit function u[n] (See Eq 2.9).
This is mentioned in example 2.3 of the same book

Taking Equation $2.74$ from Example $3$ and setting $M_1=0$ gives:
$$ h[n] = \frac{1}{M_2+1} \sum_{k=0}^{M_2} \delta[n-k] $$
Now lets take a closer look at this: $h[n]$ is non-zero only for certain values of $n$. Its a good a starting point as any, so lets look at $n=0$. We have:
$$ h[0] = \frac{1}{M_2+1} \bigg( \delta[0] + \delta[-1] + ... + \delta[-M_2] \bigg) = \frac{1}{M_2+1}$$
For $n=1$, we have:
$$ h[1] = \frac{1}{M_2+1} \bigg( \delta[1] + \delta[0] + ... + \delta[1-M_2] \bigg) = \frac{1}{M_2+1}$$
For $n=2$, we have:
$$ h[2] = \frac{1}{M_2+1} \bigg( \delta[-2] + \delta[-1] + \delta[0]+ ... + \delta[2-M_2] \bigg) = \frac{1}{M_2+1}$$
And this pattern continues, the only non-zero value, $\delta[0]$, percolates through until $n=M_2$, then we have:
$$h[M_2] = \frac{1}{M_2+1} \bigg( \delta[-M_2] + \delta[-M_2+1] + ... + \delta[0] \bigg) = \frac{1}{M_2+1}$$
And for $n = M_2 + 1$:
$$ h[M_2+1] = \frac{1}{M_2+1} \bigg( \delta[-M_2-1] + \delta[-M_2] + ... + \delta[-1] \bigg) = 0$$
And you can see, for all $n > M_2$, $h[n]=0$. So, altogether we can write $h[n]$ as equal to $\frac{1}{M_2+1}$ for $n \in [0, M_2]$ and zero otherwise.
All Equation $2.90$ does is treats the summation of impulses (from Equation $2.74$) as unit step functions, a very useful trick! Since $h[n]$ starts at $n=0$ and ends at $n=M_2$ we do: $u[n]-u[n-M_2-1]$ to get a length $M_2$ box car, then scale it by $\frac{1}{M_2+1}$ to get the right amplitude. Hope this helps!


What does energy of the difference between two DFT's $X$ and $Y$ being 0 tell us about their corresponding signals $x$ and $y$? I think that would mean that the energy of difference between the signals $x$ and $y$ should also be 0, but how can I prove this?
Does this also mean $x$ = $y$?

If the energy of the difference is zero, than the difference is zero as well
If the difference is zero, than the signals are identical
Two signals with the same power spectral density can be completely different. A simple example: both white noise and a delta impulse have the same power spectral density


Can I use IIR or FIR to get a sliding mean? My output signal shall be computed by (meta code):

o = 0  # or some arbitrary initial value
for i in input:
    o = (o * 99 + i) / 100
    print o

I call this a "sliding mean", but maybe another term is established for this.  (If so, mentioning this could help me researching this better ;-)

Is there a way to achieve this using the a / b coefficient arrays of IIRs or FIRs?  What would they look like?

I'm aiming at a solution in Python / scipy.signal.

Or is there another established way to achieve this which I can find in such libraries like numpy, scipy, pandas, etc.?

It corresponds to an exponentially weighted moving average (EWMA) filter with parameter $\alpha = 0.01$.

The difference equation can be expressed as:

$$y(n)-(1-\alpha)y(n-1)=\alpha x(n)$$
$$y(n)-0.99y(n-1)=0.01 x(n)$$

Thus it's an IIR filter. You can easily get the coefficients from the difference equation.

$$Y(z)-0.99z^{-1}Y(z)=0.01X(z)$$
$$\frac{Y(z)}{X(z)}=\frac{0.01}{1-0.99z^{-1}}$$

Therefore, a = [1 ; -0.99] and b = 0.01.


Why use windowing function instead of truncating the signal to integer amount of periods? According to this source of information if the signal is consisted of integer amount of periods then FFT would work ideally as if the signal was infinite in time. The same source mentions that windowing is not completely eliminating the spectrum leakage but just attenuating it. Therefore, my question is why perform windowing instead of truncating the signal to integer amount of periods, shouldn't the truncation to integer amount of periods give better results?

When an integer number of periods is known, truncation when followed by FFT (assuming periodicity) is equivalent to windowing it with a uniform window. Using another window may induce some leakage.

But, if your signal was sampled this way, this either means that you were very lucky, or that you already have a lot of knowledge about your signal. And thus that you have dmaller need for Fourier analysis.

The point of using specific windows is indeed to help analysis when periods are unknown.


What Does an Adaptive Filter Do? I studied a bit about adaptive filter on internet and found that its a special filter which keep on updating its filter value as soon as it proceeds. It finds out the difference between input and output and using the error function and previous coefficients finds out the new filter coefficients.

But this doesn't make any sense. It always tries to minimize the difference between input and output. So, how is it of any use, if it tries to pass all the signals as it is.

Can anyone explain me how it is being used in real day applications.

It will also be good if you can help me through some links which could help me in implementation of adaptive digital filter.

please comment if I am unclear in expressing my doubt !

The key concept that you are missing is that you are not just minimising the difference between input and output signals. The error is often calculated from a 2nd input. Just look at the Wikipedia example related to the ECG.

The filter coefficients in this example are recalculated to change the notch frequency of a notch filter according to the frequency extracted from the mains signal. One could use a static notch filter, but you would have to reject a wider range of frequencies to accommodate the variability in the mains frequency. The adaptive filter follows the mains frequency and so the stop band can be much more narrow, thus retaining more of the useful ECG information.

EDIT:

I have looked at this again and I think I understand your question a little better. The LMS algorithm needs an error term in order to update the filter coefficients. In the ECG example that I paraphrase above, I give the error term as a second input from a mains voltage. Now I'm guessing that you are thinking, "Why not just subtract the noise from the signal-plus-noise to leave the signal?" This would work fine in a simple linear system. Even worse, most examples given online tell you (correctly but confusingly) that the error term is calculated from the difference between the desired signal and the output of the adaptive filter. This leaves any reasonable person thinking "If you already have the desired signal, why bother doing any of this!?". This can leave the reader lacking motivation to read and comprehend the mathematical descriptions of adaptive filters. However, the key is in section 18.4 of Digital Signal Processing Handbook, Ed. Vijay K. Madisetti  and  Douglas B. William.

where:

x=input signal,
y=output from filter,
W=the filter coefficients,
d=desired output,
e=error

  In practice, the quantity of interest is not always d. Our desire may
  be to represent in y a certain component of d that is contained in x,
  or it may be to isolate a component of d within the error e that is
  not contained in x. Alternatively, we may be solely interested in the
  values of the parameters in W and have no concern about x, y, or d
  themselves. Practical examples of each of these scenarios are provided
  later in this chapter.

  There are situations in which d is not available at all times. In such
  situations, adaptation typically occurs only when d is available. When
  d is unavailable, we typically use our most-recent parameter estimates
  to compute y in an attempt to estimate the desired response signal d.

  There are real-world situations in which d is never available. In
  such cases, one can use additional information about the
  characteristics of a “hypothetical” d, such as its predicted
  statistical behavior or amplitude characteristics, to form suitable
  estimates of d from the signals available to the adaptive ﬁlter. Such
  methods are collectively called blind adaptation algorithms. The fact
  that such schemes even work is a tribute both to the ingenuity of the
  developers of the algorithms and to the technological maturity of the
  adaptive ﬁltering ﬁeld

I will keep building on this answer when I get the time, in an attempt to improve the ECG example.

I found this set of lecture notes to be particularly good also: Advanced Signal Processing Adaptive Estimation and Adaptive Filters - Danilo Mandic


Mininum phase question I know in a minimum-phase system, any poles or zeros are reflected.  How do I show that a minimum phase system is unique, or that only one system with that magnitude response can be minimum phased?

If you reflect any of the poles, the system becomes unstable (but even lower in phase before it blows up).  If you reflect any of the zeros, the phase increases, so no other placement of zeros can be minimum.  If the system response is a rational polynomial, then the factorization is unique, so no other (finite?) set of poles and zeros works the same (magnitude response).


Turning a MPG video into a set of coordinates I have a video in MPG format showing a white mouse with a strong black dot on the head moving around on an off-white background, so effectively a black dot moving around on a 2-D plain surface.  From this video, I want to generate a set of coordinates that describe the position of the black dot as it varies with time, that is $(X_1, Y_1),...,(X_t,Y_t)$ relative to some arbitrary point on the surface, such as the centre of the image.  I also want to create a set of arrays representing a grid across the plain, with $0$s across the grid, except for a $1$ for where the mouse is at time $t$.  That is, the motion of the mouse down in a small surface in 3 time steps/frames would be represented by
$$
        \begin{matrix}
        0 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0
        \end{matrix}
$$
$$
        \begin{matrix}
        0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0 \\
        0 & 0 & 0 & 0
        \end{matrix}
$$
$$
        \begin{matrix}
        0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 \\
        0 & 1 & 0 & 0
        \end{matrix}
$$
or $(2,3), (2,2), (2,1)$.

I think the first stage is to somehow convert the video into a set of still images, so effectively split the video into the frames that make it up.  I have a mac, and any suggestions on a piece of software to do this would be much appreciated.

The second stage requires some kind of image processing package, but again, I'm not sure what to use, so advice would be much appreciated.

The project I am working on requires the final data (i.e, matrix and coordinate) forms, and isn't that concerned with how I get the data in the first place, so I am looking for simple software packages that lets me get this data.  However, if anyone has some code in matlab or python that lets me analyse the still images from scratch, that would be appreciated too.

Thanks!

You can use opencv library with python for image processing.

Use Videocapture object to capture video from video MPG file.
Do cap.read() to read frame by frame.
Convert your frame to grayimage by cvtColor().
Use threshold() function with threshold type-4(threshold with invert option) and with max-binary value as 1. Now you got Mat array with 1 for your mouse with zeroes for background.
Now you can search for 1 in the array and get the coordinates.


Obtain motion between image features by means of the homography matrix I have been reading several really good answers about how to obtain camera pose (rotation, translation, etc), such as this one and this one, but I am not sure how this applies to my problem.

My problem: I have two images of a planar scene (terrain image from above), and I want find the motion that the camera has experienced between the two images $(dx, dy)$.
What I have done: I have image1 and image2 of a terrain with their features matched, and obtained its homography via OpenCV's C++ function:

cv::Mat H = cv::findHomography(image1points, image2points, CV_RANSAC);

Which looks like it is indeed an affine transformation:
$$\pmatrix{a&b&c\\d&e&f\\0&0&1}$$
What I need help with:

In this matrix, I don't know what is the reference frame of the $dx$, $dy$ parameters. $(dx, dy) = (c, f)$
How can I extract the rotation angle and scaling factor between the 2D axes of the two images. (Right now I am using this suggestion, but results don't seem reliable).
Once I have all the four parameters: $dx$, $dy$, rotation angle, scaling. How can I actually compute the motion of the camera (between the images) using the first image's axes as my reference frame.

First, calibrate your intrinsics: The focal lengths and the principal point. Homography is not really a rigid transformation and rather a mapping of a plane onto another one. What you really need is a 6DOF parametrization, which is the camera pose. If you have the intrinsics, you could transfer your coordinates to the normalized coordinates and estimate the essential matrix using 8 point algorithm - or a better one would be five point algorithm (check this, this and this). The operation is the same as solving for fundamental matrix in the uncalibrated setting. But, just because you have your $K$, you could recover essential matrix. Note that, such estimation doesn't assume planar scenes and can in fact operate on 3D structures.

You could then decompose your essential matrix as

$$E=R[t]_x$$
where $[t]_x$ is the skew symmetric form. To obtain the components, proceed as described in here. Note that, this problem has 4 possible solutions and identifying the best one is usually found via checking all solutions.

If your cameras are uncalibrated, you could maybe approximate the calibration using some known dimensions or a grid, within the scene.

If you are using simulated images, then you know the FOV and image dimensions. If this is the case, you could get the camera matrix. Or even better, if you have your $K$ matrix (you could artificially make up one), then you could compute the OpenGL equivalent. Check here and here for that.

To address your questions one by one:

They are in pixel coordinates and denote the principal point (computed center of the image).
Unfortunately, $E$ is defined up to a scale and you could only recover the scale, if you have a known object with size within the scene. You just need to know one distance to do this. It could be a marker, the length of some existing furniture, or the object you inspect and etc.
Use the explanation I already provided (above).

And additionally, it is not an affine transformation, it's a projective mapping (perspective transformation). In fact homography is nothing more than the transformation of one plane to another.


Associative Property of Convolution: Consider the following sequences:

$x_{1}(n) = A$ (a constant), $x_{2}(n) = u(n)$, $x_{3}(n) =\delta(n)-\delta(n-1)$.

($\circledast$ stands for linear convolution)

If I perform the operation $x_{2}\circledast (x_{3}\circledast x_{1})$, the value I am getting is $0$, where as if I perform $(x_{2}\circledast x_{3})\circledast x_{1}$ the value I am getting is A.

Since Convolution is associative, why are the answers different?

My approach: In general, assuming $x_{1}(n)$, $x_{2}(n)$, $x_{3}(n)$ are of infinite lengths,

$$x_{2}\circledast (x_{3}\circledast x_{1})=\sum_{k=-\infty}^{\infty}x_2(k).x_{3}\circledast x_{1}(n-k)=\sum_{k=-\infty}^{\infty}x_2(k)\sum_{l=-\infty}^{\infty}x_3(l)x_1(n-k-l)$$

Let $m=n-k-l$,

$$=\sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}x_2(k)x_3(l)x_1(n-k-l)=\sum_{k=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}x_2(k)x_3(n-m-k)x_1(m)$$

$$=\sum_{m=-\infty}^{\infty}x_{2}\circledast x_{3}(n-m)x_1(m)=(x_{2}\circledast x_{3})\circledast x_{1}$$.

Hence I feel Convolution is associative even if the sequences are of infinite lengths. Where did I go wrong?

The proof of associativity of discrete convolution relies on the assumption that multiple infinite sums can be evaluated in any order. This is not true if some of the involved sequences do not converge absolutely, which is the case for the given sequences $x_1[n]$ and $x_2[n]$. Note that the convolution sum $x_1\star x_2$ does not converge, i.e., $x_3\star (x_1\star x_2)$ gives yet another (infinite) result.

In continuous time you have the same problem. Associativity of continuous convolution relies on Fubini's theorem for switching the order of integration. If the assumptions of Fubini's theorem are not satisfied for the given functions, convolution is not associative.


How do derivative masks work for finding edges in image? As I know, masks like follow are derivative mask,

but I don't get it why they can detect edges and which one of them can find vertical edges and which one can find harizontal edges?

Basically, the triplet 1,0,-1 (and similarly 2,0,-2) is a discrete, 3-point derivative operator. It can be applied column- (left) or row-wise (right). Let us try first with a row: a flat one, f = 0,0,0,0,0,0,0,0, a step-like s = 0,0,0,1,1,1,1,1 and one hat-like 0,0,0,1,1,0,0,0. 0 here denotes black, and 1 white.

The results will be, on the $6$ central points for each:

f: 0,0,0,0,0,0
s: 0,1,1,0,0,0
h: 0,1,1,-1,-1,0

The same reasoning applies to columns. The Sobel operators are composed of a second filter in the opposite direction: 1,2,1 that smoothes the data. So the global effect: it differentiates in one direction, smoothes in the other. THe resulting value will be high or low, positive or negative, depending on the type of contour in the image.

The (close-to-)zero values are features not detected by the derivative (hopefully non-edges). High amplitude values and their signs are hopefully the strength of the edge and the grading (positive: black-to-white, negative: white-to-black).


Phase transitions in GPS C/A signal I'm reading "Fundamentals of Global Positioning System Receivers" section 7.3 "Maximum data length for acquisition" and it's talking about phase transition:

  Theoretically, if there is a navigation data transition, the
  transition will spread the spectrum and the output will no longer be a
  cw signal. The spectrum spread will degrade the acquisition result.
  Since navigation data is 20 ms or 20 C/ A code long, the maximum data
  record that can be used is 10 ms. The reasoning is as follows. In 20
  ms of data at most there can be only one data transition. If one takes
  the first 10 ms of data and there is a data transition, the next 10 ms
  will not have one.

  In actual acquisition, even if there is a phase transition caused by a
  navigation data in the input data, the spectrum spreading is not very
  wide. For example, if 10 ms of data are used for acquisition and there
  is a phase transition at 5 ms, the width of the peak spectrum is about
  400 Hz (2/ (5 × 10−3)).

My understanding was that the C/A signal is a pseudo-random Gold code that is XORed with the data (each data bit being 20*1023 chips long) and the resulting binary signal is BPSK-modulated.

So I would think phase transitions in the C/A signal (caused by the BPSK modulation) would occur much more often, not just between data bits but between chips, every few chips.

What is wrong in my understanding of the modulation of the text?

That text is confusing on its own as the spectrum spreading as transmitted is NOT affected as described. This is suggesting that if you acquire on any arbitrary 20 ms sequence of received data, your acquisition results will vary since a data transition can occur equally likely anywhere in that 20 ms given you select a random starting point for your acquisition. You avoid this whole issue by effectively performing acquisition using a 20 ms sliding correlation on 40 ms of received signal (assuming you want to maximize processing gain over the full 20 ms available).

The BPSK modulation includes the C/A Code specifically making the carrier change back and forth 180°. The C/A code is 1023 chips long; when a C/A code chip is = 1, the carrier prior to data modulation (such as L1 at 1575.42 MHz) is transmitted with 0° phase reference. When a C/A code chip = 0, the carrier is transmitted with a 180° phase reference. The C/A code is at a 1.023 MSps rate, so there are 1023 phase transitions in each 1ms C/A code sequence.

Then there are 20 C/A code sequences in one data bit:  Each data bit additionally toggles the carrier between 0 and 180°: So the C/A code is sent 20 times to send a "1" data bit, and the negative of the C/A code (0110... instead of 1001) is sent 20 times to send a "0" data bit. Thus the data rate is 50 bps.

Below is an example of a correlation to SV24 for a single 1023 length sequence associated with that satellite from actual data captured from a GPS antenna with a digital sampling scope. The correlation was done by sliding the sample-matched 1023 length sequence through the received signal that in every other way appeared as white noise. With frequency offsets removed, the result appears as in the plot below showing the correlation peaks every 1 ms at every point in time that the local C/A code aligned with the transmitted C/A code. What we see in the plot is three data symbols (either 1 1 0 or 0 0 1 since phase ambiguity hasn't been removed).

This is the same correlation result plotted on a complex plane (the above plot is the real part of this versus time). This shows how the carrier offset was removed and suggests carrier tracking approaches.


Making a low pass filter for irregular samples We are interested in reconstructing the red channel of a scene.

The red channel is slowly varying.

Unfortunately the camera does not have a Bayer sensor, and the red samples are irregularly spaced over the grid:

How do you design a 2D low-pass filter (convolution kernel) that:

only depends on the samples of red that you have
has as close to constant frequency response as possible?

I was thinking of something along the lines of solving an optimization problem relating to the NDFT. Would this idea work? :

  Take the NDFT of the red components we have. Then in continuous space, form the sum of all the sinusoids from the frequency information in the NDFT. Re-sample the function along the whole grid.

I realize that this is not a filter, but I do not believe any filter will give you a reasonable result given the strange sample locations. This seems to be borne out by several references on the web. This one for example and this PDF one.

I suggest something something along these lines, using projections onto convex sets (POCS).

Algorithm Descroption

Assume that the known red pixel values $R_i$ are at $(x_i,y_i)$ so that your start red image is:
$$
I_{\tt red}(x,y) = \sum_{i=1}^N R_i \delta(x-x_i, y-y_i)
$$

Set $k=0$, and initialize your estimate of the fully sampled red image with $I_{\tt red}(x,y)$:
$$
\hat{I}^k_{\tt red}(x,y) = I_{\tt red}(x,y).
$$

Take the 2D FFT of $\hat{I}^k_{\tt red}$:
$$
F_{\tt red}(n,m) = {\tt FFT}\{\hat{I}^k_{\tt red}(x,y)\}
$$
Because the red channel is slowly varying, remove all frequencies above a cutoff:
$$
F^{lp}_{\tt red}(n,m) = F_{\tt red}(n,m) L(n,m)
$$
where $L$ is a low pass mask.

Form the inverse FFT to get $\hat{I}^{k+1}_{\tt red}(x,y)$
$$
\hat{I}^{k+1}_{\tt red}(x,y) = {\tt FFT}^{-1}\{F^{lp}_{\tt red}(n,m)\}
$$

Set the known pixels of $\hat{I}^{k+1}_{\tt red}(x,y)$ from $I_{\tt red}(x,y)$.

If the difference between $\hat{I}^{k+1}_{\tt red}(x,y)$ and $\hat{I}^{k}_{\tt red}(x,y)$ is small enough, stop, otherwise reiterate with $k=k+1$.

Results

An example run with just a $6 \times 6$ matrix is below. The left image shows the original pixels. The right image shows the POCS version.  The red dots show the original pixel locations in both images.

Different Initial Conditions

I tried setting the initial image to

Ired <- array(100000,c(ImageSize, ImageSize))

and

Ired <- array(128,c(ImageSize, ImageSize))

and the results are below, so I do not believe the initial values for the unknown pixels matter.

The plot scales the greyscale values, which is why the 100,000 example looks wrong in the left-hand plot.

R Code Below

# 32018
# Non-Bayer 2D interpolation

# Model pixels as:
# P(x,y) = (a_0 + a_1 * x + a_2 * x^2) * (b_0 + b_1 * y + b_2 * y^2)
#        = [a_0 a_1 a_2 b_0 b_1 b_2] * [1 x x^2 1 y y^2]^T
# And so we need to find the a_i and b_i from specific pixel and coordinate values.
#
# install.packages('ptw')

require(ptw)

pad <- function(x,N)
{
  p1 <- ptw::padzeros(x,N,side="right")
  p2 <- t(ptw::padzeros(t(p1),N,side="right"))

  return(p2)
}

Npoints <- 8
xy <- array(c(c(1,1), c(1,5), c(2,3), c(3,6), c(4,2), c(4,4), c(5,6), c(6,3)), c(2,Npoints))
xy_norm <- (xy - min(xy))/(max(xy) - min(xy))
Pxy <- array(runif(N, 0, 255), c(1,Npoints))

ImageSize <- max(xy)

Ired <- array(0,c(ImageSize, ImageSize))

for (idx in 1:Npoints)
{
  Ired[xy[1,idx], xy[2,idx]] = Pxy[idx]
}

Ihatred <-array(0,c(ImageSize,ImageSize))
Ihatred2 <- Ired
count <- 0

while (sum(abs(Ihatred - Ihatred2)) > 1)
{
  Ihatred <- Ihatred2

  Fred <- fft(pad(Ihatred, 2*ImageSize))

  Fredsize <- dim(Fred)[1]
  cutoff1 <- Fredsize / 8
  cutoff2 <- 7 * Fredsize / 8

  Fred[4:16,] <- 0
  Fred[,4:16] <- 0

  IFred <- fft(Fred, inverse=TRUE) / Fredsize / Fredsize

  Ihatred2 <- Re(IFred[1:ImageSize, 1:ImageSize])

  for (idx in 1:Npoints)
  {
    Ihatred2[xy[1,idx], xy[2,idx]] = Pxy[idx]
  }

  count <- count + 1
  if (count > 100)
  {
    break
  }
}

par(mfrow=c(1,2))
image(Ired, col= grey(seq(0, 1, length = 256)))
for (idx in 1:Npoints)
{
  points(xy_norm[1,idx] , xy_norm[2,idx] , col="red", lwd=5)
}
title('Original with only known red pixels')
image(Ihatred2, col= grey(seq(0, 1, length = 256)))
title('POCS solution')
for (idx in 1:Npoints)
{
  points(xy_norm[1,idx] , xy_norm[2,idx] , col="red", lwd=5)
}


Filter Design in Matlab I want to design a filter to remove all frequency components from a .wav file except those within $\pm 25\,\text{Hz}$ of $523\,\text{Hz}$ as well as its harmonics (up to the Nyquist frequency).

Attempt

I tried to suppress frequencies below 523-25 Hz using a high-pass filter, and suppress the ones higher than 523+25 Hz using a low pass filter. Here's my code:

[s, Fs] = wavread('x.wav');

for N=2:1:(Fs/(2*(523+25)))
    % High-pass filter
    FsNorm =  (523-25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'high');
    sHigh = filtfilt(b, a, s);

    % Low pass
    FsNorm = (523+25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'low');
    sLow = filtfilt(b, a, s);
end

wavwrite(sLow, Fs, 'x_filtered');

The code doesn't work, the resulting file is mute. My intention is to apply each filter to the output of the previous filter, but I don't know how to implement this in code.

I used $(523+25) N_{max} < \frac{F_s}{2}$ to include all the harmonics up to the Nyquist frequency.

So what do I need to do to get this filter to work?

Edit:

Here's my new code:

[s, Fs] = wavread('x.wav');

sNew = zeros(size(s));

for N=2:1:(Fs/(2*(523+25)))
    % High-pass filter
    FsNorm =  (523-25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'high');
    sHigh = filter(b, a, s);

    % Low pass
    FsNorm = (523+25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'low');
    sLow = filter(b, a, s);

    Now sum the results
    sNew = sHigh + sLow;
end

sNew

wavplay(sNew,Fs);
wavwrite(sNew, Fs, 'x_filtered');

I still get sNew=0 and the sound is still mute. What's wrong?

You just need to make a simple modification:

Send the original signal through each filter and sum the output of each filter together to get your final result.

At the moment you are passing the filtered data on to each subsequent filter. If you think carefully about that, the output will be silent i.e. you filter out everything but 523 Hz  +/- 25 which leaves you just frequencies around 523 Hz, then you take that and filter out everything but 1046 Hz  +/- 25 which will strip out the frequencies around 523 Hz and then you have silence after just two iterations of your loop.

s = s(:).';
sNew = zeros(1,length(s)); % The eventual output

for N=2:1:(Fs/(2*(523+25)))
    % High-pass filter
    FsNorm =  (523-25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'high');
    sHigh = filter(b, a, s);

    % Low pass
    FsNorm = (523+25).*N / (Fs/2);
    [b,a] = butter(10, FsNorm, 'low');
    sLow = filter(b, a, sHigh); % Note we are filtering sHigh, not 's'

    % Now sum the results
    sNew = sNew + sLow;
end


Laplace transform of a time domain sampled data MATLAB I have two sets of one second voltage data sampled with 4000Hz and I can plot all the voltage points vs time points in MATLAB. So it means I have a data matrix with with length of 4000 one column for the time in seconds the other for the voltage.

Now I have simultaneously sampled two data matrix in time domain with this way. One is input to a filter (Vin,t) the other is output (Vout,t). I want to find the transfer function both for amplitude and phase shift.

How can I take the Laplace transform of this data? Can I do it without converting it to a expressible function like poly-fit?

Since you are sampling a real signal, you can just use the a Fourier transform, the fft function in MATLAB. You can divide the FFT of the output by the input, and then fit a curve to the result to approximate the transfer function


Getting audio samples from power spectrum Let's say I have a signal sig in a form of a .wav file, I've read that file, and now I have an array of samples of that audio, after that I've performed FFT on that array and got an array of complex numbers, that'll be fft, which consists of positive as well as negative values. After that I've computed power spectrum, ps using this formula:

ps[a] = sqrt(fft[a].real^2 + fft[a].imag^2).

So, basically process looks like this: sig$\rightarrow$ fft$\rightarrow$ ps. Now, the question is: can I go backwards, meaning: ps$\rightarrow$ fft$\rightarrow$ sig?

No, the power spectrum only holds half of the information that the DFT'ed signal held, which held all the information that the original audio signal held.

Sadly, that half is important, and not reconstruction is possible. You might come up with an estimator that "guesses" the phase information you're losing by doing the magnitude square, but it will not recover the original signal in the general case. I'm not aware of any such estimators that do something useful.

In short: you can reverse the DFT, it's a reversible operation. You cannot revert calculating the magnitude square of complex numbers – that is a irreversible operation, and hence, you can't recover the original signal from a PSD estimate.

To illustrate this:

Imagine you record one hour (let's say from $t=-18000$ to $t=1800) of perfect silence, and in the middle of that is one second of a pure 100 Hz sine tone.

The DFT of that is pretty simple: In time domain, it's a sine (which is infinite) multiplied with a unit rectangular window. We know the spectrum of that; it's the convolution of the Fourier transform of the sine (which is dirac with "area" $\frac 12j$ at the positive frequency, and $-\frac12j$ at the negative frequency) convolved with the Fourier transform of a rectangular window (which is always some scaled $\text{sinc}$ function). Hence, we know that the DFT of that recording must be the sum of two shifted sincs, one of those being "mirrored" on the frequency axis.

Now, we estimate the PSD by just calculating the magnitude square point-wise, and we end up with something that looks like two positive lumps at the positive and negative frequencies. For obvious reasons, it's symmetric to the magnitude axis. We call the DFT $X(f)$, and that leaves us with the PSD estimate $|X(F)|^2$

We now decide that the 1s of tone shouldn't be at the $t=0$, but let's say at $t=900$; now, instead of recalculating the whole spectrum, we just apply the time-shift property of Fourier transform: a $\Delta t$ shift in time domain is a multiplication with $e^{j2\pi f \Delta t}$ in frequency domain. That way, we just need to multiply with this function in frequency domain.

And lo! Let $X(f)\cdot e^{j2\pi f \Delta t}$ be our spectrum, and $|X(f)\cdot e^{j2\pi f \Delta t}|^2$ be our PSD estimate. Observe that the modulation with the complex $e$xponential doesn't change magnitude square – hence, you can't tell time shifts from your PSD alone.


Forward Linear Predictor I just need some advice on a project. I have an unwanted sinusoid that has time-varying frequency (its not sinusoidally varying) which is relatively constant over a small section. I've been advised to use a Linear Predictor in order to deal with it and remove it.

I just need some advice on how to go about implementing it. I understand the theory, my main issue is with this equation for a forward-linear predictor

$$
\hat{x}(n) = \sum_{i=1}^p a_i\,x(n-i)
$$

Now, what's the point in estimating the next value of the signal if it's also going to contain the same unwanted sinusoid component? I understand how to use this (finding coefficients using Yule-Walker/Levinson-Durbin)

EDIT: More info, Sinusoid is added on top. I can almost hear the telephone conversation clearly when the sinusoid dips down. I understand that for a small enough section, I could just use a notch filter but of course it's not feasible to do this for each small section manually. I need a way of predicitng it's frequency hence the linear predictor The sinusoid has a significantly higher amplitude than the conversation and it's frequency varies from 200Hz to 2kHz

So, the algorithm I'd propose is the following:

Filter your input signal with a 200 Hz – 2 kHz bandpass filter
Since in that band your sine is the dominant oscillation, a PLL would lock onto that. So, use a second-order PLL to track that wave; there's many PLL implementations that'll output the reference as tone with the correct phase.
Typically, PLLs lose amplitude information, so you can't simply remove the sinusoid by subtracting it from the audio signal. Instead, you'd use the reference tone to multiply the signal with, which mixes down the original tone to 0 Hz. A low-pass filter (in this case, even a single-pole IIR) might be trivially used to extract the instantaneous amplitude.
Now that you have the right tone with the right frequency and phase, you can subtract it from the input
To predict how the tone is going to be developing its frequency, you could use a linear predictor, but honestly, a second-order PLL is a linear predictor for phase (a linearly changing phase is a sinusoid; the coefficient of the linear term is the frequency). If you need to work with a sinusoid which might be undergoing a linear change in frequency: third order PLL.

Even more elegant: Use the output of 2. Mix your signal down with this. Apply 0 Hz Notch filter (again, single-pole IIR might do, or other kind of narrow high-pass filter). Shift up in frequency again.


Tracking a Sine Wave with Kalman Filter - How to Account for Offset (DC Signal)? I am attempting to create a Kalman filter to track a sine wave (I am using a linear Kalman filter example assuming I already know the frequency of the sine wave) - the example I am using is derived on pages 194-196 of "Fundamentals of Kalman Filtering: A Practical Approach" 2nd edition by Paul Zarchan and Howard Musoff.

It is working to track the AC part of the signal, however the offset of the sine wave from the $x$-axis is not correct, it seems to be tied to my value of $R$, I'm guessing it is because in the derivation of the model no offset is included but I cannot see where an offset term should be inserted, how should I go about this?

The matrices describing my system are like so:

$$\mathbf \Phi_k = \begin{bmatrix} \cos(\omega T_s) & \dfrac{\sin(\omega T_s)}{\omega} \\ -\omega \sin(\omega T_s) & \cos(\omega T_s) \\\end{bmatrix}$$

Where the state matrix is:

$$\mathbf X = \begin{bmatrix} x \\ 0 \end{bmatrix}$$

And I have set $\mathbf Q$ and $\mathbf P$ like so:

$$\mathbf Q = \begin{bmatrix} \dfrac{T_S^3}{3} & \dfrac{T_s^2}{2} \\ \dfrac{T_s^2}{2} & T_s \end{bmatrix}, \quad\mathbf P =  \begin{bmatrix} 9999999999999 & 0 \\ 0 & 9999999999999 \end{bmatrix} $$

I am now attemping to use this to track a generated signal of a noisy sine wave of amplitude and frequncy $1$ with an offset of $300$.

If $R$ is set to $0.1$ I get the following output of my filter:

Which is offset by $300/20$, as can be seen from this plot where I add an offset to the Kalman filter output:

Changing $R$ changes this offset (as one might expect as $R$ denotes the uncertainty in your data).

How can I add this offset into my model such that my Kalman filter can tracking the sine wave correctly regardless of offset?

EDIT: I've been thinking more about my question and realised that in a real signal with noise this offset could be considered a separate constant signal with frequency $0$ which could be extracted by using a Kalman filter fitting a constant value (or simply a moving average filter) and that the offset is not inherently part of the sine wave.

In which case how might I remove the offset in the Kalman filter output entirely? (e.g. if I used a bandpass filter no offset would be present).
Since my model does not include an offset I might expect no offset at all but I do see an offset of $300/20$, as discussed. Why does this occur if there is no such offset parameter in my model?

Well, in continuous time, a sinusoid with a bias can be seen as the output of the linear system
\begin{align*}
 \begin{bmatrix}\dot x_1\\\dot x_2\\\dot x_3\end{bmatrix} &= \begin{bmatrix}0 & 1 & 0\\-\omega^2&0&0\\0 &0 &0\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}\\
 y &= \begin{bmatrix}1&0&1\end{bmatrix}
\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}
\end{align*}
As the initial condition $(x_1(0),x_2(0),x_3(0))$ vary in $\mathbb R^3$, the system generates all the possible $y(t)$ of the form
$$
y(t) = a\cos(\omega t+\phi) + b
$$
The system is observable for all $\omega\ne 0$,therefore a Kalman Filter is possible. With $x:=col(x_1,x_2,x_3)$ and
$$
A=\begin{bmatrix}0 & 1 & 0\\-\omega^2&0&0\\0 &0 &0\end{bmatrix},\qquad C=\begin{bmatrix}1&0&1\end{bmatrix}
$$
The filter has the form
\begin{align*}
\dot{\hat x} = A\hat x -L(C\hat x-y)
\end{align*}
Where $L$ is chosen according to the Kalman theory, on the basis of the covariance matrices of the noises and the solution of the Riccati Equation. For discrete time systems the procedure is the same.
EDIT:-----
If you only use a model for the oscillator you are actually restricting the $A$ map to the $x_1$-$x_2$ subspace. That's equivalent to take $C=[1\;0\;0]$ in the above system. In this latter case you loose observability (actually you do not have neither detectability) and therefore you cannot observe the bias.


Filters applied separately to right and left channel of stereo microphone yielding inaccurate results I must estimate the direction of arrival of an acoustic source, and I hope to distinguish between a voice input and musical input using a frequency cutoff. In order to do this, I apply the same bandpass elliptic filter on the left and right channels of the stereo microphone. I have applied the filter to white noise and verified using a spectrogram that the filter does indeed retain the frequencies that it is designed to retain. I calculate the power after filtering, and use a power threshold to decide whether the frame of data is noise or not. I then use cross-correlation between the left and right channel, to estimate the time delay between left and right channels, and use that to calculate the angle of arrival.

I had first simulated the situation with two signals in numpy, with random noise added to it, and delayed by a certain number of samples. I then applied the same filter to the 'left' and 'right' channels, where the 'right' channel was simply delayed by a certain number of samples from the 'left', and random noise added. I was then able to compute the correct delay between the signals by giving the filtered 'left' and 'right' channels as input to the cross-correlator.

In real time, however, after applying the filter to the left and right channels, the time  delay is neither accurate not consistent. Could it be possible that the process of computing and applying the filter (in python) is interfering with the acquisition of signals, thereby adding some random time delays and making my time delay measurements entirely inaccurate? If so, how would I go about resolving the issue?

I hope to distinguish between a voice input and musical input using a frequency cutoff.

The use case scenario is not exactly clear here. Are you trying to get the vacuum cleaner to be called to a location by responding to a command, or are you trying to get the vacuum cleaner to follow a person that is talking in general. (A specific person perhaps?).

The second use case scenario ("follow me as I am talking to you") is probably a "voice" recognition application and localisation could be added in two stages. In stage 1, a classifier that is trained to the voice of a specific person is fed frames from the sound card. Once the classifier detects that a frame contains a "voicing", then stage 2 is activated. In stage 2, you run a frequency analysis and try to locate the fundamental tone of the person's voice. Once you locate that, you then base your direction finding on the phase difference of that particular harmonic (or band). Obviously here, because the voice is processed in frames, the vacuum cleaner's reaction time depends on the length of the frame and their overlap. Very briefly, if you are processing 1 second long frames (with zero overlaps), then the robot makes 1 direction decision every 1 second.

The first use case scenario ("come towards me when I call you") is probably a pattern recognition problem. If the vacuum cleaner is to respond to a voice command, like "Robot, clean up", then, at the very least, you could use a matched filter approach. You can then rectify and integrate the output of the matched filter and pass it through a threshold. The combined output of this will be a pulse whose duration will be, approximately, equivalent to the duration of the voice command. Running this process on two microphones with some separation between them will give you two pulses (one on each channel) with a slight time difference from which you can then estimate direction.

A variation of this technique would be to call the robot via a whistle pattern (sounding a bit like R2D2 but obviously, at human speeds). The whistle pattern is attractive because:

You can filter your input to the expected range of the whistle (and it's very hard to whistle at 20kHz)
A whistle is very near to a sinusoid, therefore you can use very simple means of detection
Being a whistle and very close to a sinusoid, allows you to use the established techniques of phase difference to detect direction. In fact, for a given whistle pattern, you could choose to focus on the high frequency parts of it to get more accuracy out of the direction estimation.

So, given an audio stream which can even be processed on a sample-by-sample basis, establish a filter bank of filters that are atuned to a given frequency. Rectify and integrate their output (so that you essentially measure power at that band) and run that through a threshold. Once the user hits the right notes, you get a sequence of pulses and the times at which the outputs crossed the threhsold. Analysing the timings within one channel you can infer if the user just whistled the "secret pattern" and then analysing the timings across channels (or simply using a phase detector) you can infer direction.

Hope this helps.


How to add pitch-shifted version of signal to itself? I want to write an algorithm to add pitch-shifted version of audio signal to itself (like Electro-Harmonix POG does). And I need to do it in realtime. How to implement it?

My guess is: break signal into pieces, get FFT of each piece, copy each of the frequencies in the original signal (for example, for 1 octave up shift copy 20hz to 40hz, 40hz to 80hz and so on), get the iFFT, play it back. But I don't know if I got the right idea. Will it work?

okay, try googling "pitch shifting AMDF" or "pitch shifting TDHS" and see what you hit there.  also check out DSP Dimension and maybe music-dsp and the mailing list.

you need three concurrent processes:

a good pitch detector that gives you a good estimate of the period (to a precision of a fraction of a sample).  this period estimate should be updated as least as often as every 15 or 20 ms.
a real-time resampling and interpolation alg with a pointer that has precision at the fractional sample.  this pointer will have both integer and fractional parts, the integer part tells you where in your delay buffer you are interpolating, the fractional part tells you how you're interpolating (like the interpolation coefficients).  try linear interpolation at first, it's easier and you can improve it later.  now this pointer will advance faster than the write pointer for upshifting and slower for downshifting.  either way, that pointer will run into the edge of the delay buffer (and glitch or click or pop when it wraps around) unless you do something about it.
a splicing and crossfading alg that allows you to splice or seamlessly "jump" from one instantaneous pointer position in the delay buffer to another pointer position.  the difference between those two pointer positions should be one period as determined by your pitch detector.

if i hint more than that, i might be spilling beans (internal trade secrets) from one manufacturer or another.  but this is the main idea.  think of splicing as necessary to time-scale the input (say time-expand if you're upshifting), and then step 2 is playing that time-expanded audio back faster which both returns it to the original tempo, but bumps it up in pitch.


Hilbert Transformer processing block by block I am processing audio and trying to obtain a 90º phase shift when processing block by block.

I have implemented a Hilbert Transformer (from Lyon's 3rd. Ed. & MATLAB). MATLAB code below (I believe this is similar to what MATLAB hilbert does).

function result = my_hilbert(s)
    fprintf("Using My Hilbert\n");
    n = length(s)

    %fftLength = 2 ^ nextpow2(2 * n - 1);
    fftLength = n

    H = zeros(1, fftLength);
    H(1) = 1;
    H(2:n / 2) = 2;
    H(n / 2 + 1) = 1;

    x = zeros(1, fftLength);
    x(1:n) = s
    X = fft(x)

    Y = X .* H
    y = ifft(Y)

    result = y(1:n);
end

This implementation seems to assume one is processing the complete audio stream in a single call.

I need to be able to achieve a phase shift but will be processing one block (block size between 1024 and 2048 samples) at a time.

When I simply call this implementation on each block I receive audible pops around the block boundaries.

I have considered overlap+add but when using an FFT length other than length of the input, the results are no longer the same as MATLAB's hilbert() function.

How can one obtain a 90º phase shift when processing block by block without these discontinuities?

The key to block (or stream) filtering in Matlab is using the filter's command ability to return and reuse its internal state.

To illustrate, first create a filter that implements the Hilbert transform.

hilb = designfilt('hilbertfir','FilterOrder',60, ...
       'TransitionWidth',400,'SampleRate',1e4);

This example is taken directly from the documentation (type doc hilbert in Matlab).

Now, we need this form of the filter command (see doc filter):

[y, zo] = filter(b,a,x,zi)

which uses initial conditions zi and returns final conditions in zf.

To put everything together, say you have incoming data and store it in variable block. Then, you can do something like:

block = acquire_block_of_data();
[output, zo] = filter(hilb, 1, block);
while 1
    block = acquire_block_of_data();
    zi = zo;
    [output, zo] = filter(hilb, 1, block, zi)
end

Note that we need an initial call to filter without state; alternatively, you can set let zi have an initial value equal to a vector of zeros with the same length as the filter order.


What can be done with I/Q data which would otherwise be impossible? As I understand it, I/Q data is the complex representation of any waveform. It enables us to understand two things at any arbitrary point in time:

The projected amplitude of the waveform.
The direction of the waveform.

However, it's still unclear to me what makes I/Q data so special.
Logically, I can discern the direction of a waveform implicitly: If I'm talking (or transmitting), I know a signal is travelling away from me. On the other hand, if I'm listening, I know any signal I hear must have been travelling toward me.
And if I want to know the amplitude of a waveform, I can just wait until I see a peak.
This in mind, what sorts of modulations/demodulations would be impossible with real data instead of I/Q data? Respectfully, I'd appreciate an example based on reality as I understand that I/Q data makes the result of a FFT transform deterministic, but I don't really understand why that's relevant.

The main motivation for using quadrature (I/Q) signals is spectral efficiency.
Most signals are naturally baseband: they have no carrier, their energy is concentrated around DC, and they are limited to bandwidth $B$. An example would be an audio signal, which has bandwidth of roughly 20 kilohertz.
To transmit such a signal, often it must be modulated (or "upconverted") using a carrier. For example, in AM broadcasting, the audio signal is first limited to a bandwidth $B=10\,\text{kHz}$, and then upconverted to a frequency of roughly one megahertz.
The issue is that, when upconverted, a signal occupies a bandwidth $2B$, which is double the original bandwidth $B$. So, an AM signal occupies $20\,\text{kHz}$ "on the air". The extra 10 kilohertz do not have any additional information; they are simply a loss of half of the available bandwidth.
Since bandwidth is a very limited and expensive resource, such a waste is unacceptable.
One way around this issue is to use single side-band modulation (SSB), which eliminates the waste. However, SSB signals are more difficult to generate and have other disadvantages.
The most popular solution is to use quadrature signals. This allows transmitting two signals of bandwidth $B$ over a chanel of width $2B$, restoring the average spectral occupancy to $B$ hertz per signal.
To see more of the math behind quadrature transmitters and receivers, see this answer.


Choosing Drive Signals for System Identification? (I'm just learning a little about system identification so apologies in advance if this question is badly worded)

How do you go about choosing drive signals for system identification? I've seen PRBS signals used but it seems like that's going to work well for frequencies around the chip rate but not really low frequencies; I've also seen frequency sweeps.

If I have a SISO system that I know is close to a 2nd order linear system with poles in a certain range, and I can drive it with an arbitrary signal up to some amplitude A for up to some time length T, how do I pick a signal that would give me the best responses for determining the accuracy of the transfer function?

I tried googling for "system identification drive signals" but I don't see anything that pertains to my question.

edit: one particular type of SISO system I've dealt with is an (input=power dissipation, output=temperature) system for power semiconductor thermal behavior, and it seems very hard to model because there's usually a dominant pole at very low frequencies (<1Hz) and the next one might be 100 times higher, so any high-frequency drive signals just get very heavily attenuated.

For linear systems, you can completely characterize the transfer function using its frequency response, so a frequency sweep would be one possible choice. However, you would need to ensure that at each test frequency, you allow time for the system's transient response to die out before measuring its steady-state amplitude/phase response.


Removing marker pen artifacts from scanned microscope slide images I am developing an image processing algorithm / plugin for Fiji (ImageJ) that processes scanned microscope slide images (20x magnification) and extracts features from the slides.

The slides contain histological tumor sections that are stained immunohistochemically with an antibody that identifies nuclei positive with the Ki-67 protein, along with a counter-stain that identifies all nuclei. The so called proliferation index, which measures how fast the cells are dividing, is then calculated as the ratio of positive-to-negative nuclei. For the feature extraction, I basically need to extract two color components from the RGB image; the one corresponding to the positive nuclei (diaminobenzidine; brown color) and the other corresponding the negative cells (methyl green; light green color).

I am using Fiji's Color Deconvolution plugin (http://www.mecourse.com/landinig/software/cdeconv/cdeconv.html) for separating the two stains. This works remarkably well, except when the slide has been drawn onto with a marker pen. Pathologists sometimes mark regions of interest on the slide with a marker pen so that it would be easier to inspect the slide later on. My data set contains slides that have been drawn onto with a red/blue/green or a black marker. The problem is that the presence of colors other than the two stains fails the color deconvolution algorithm totally.

 Above is an example of a slide at 20x magnification showing the edge of a line drawn with a black marker pen on the top-half of the image. The lower part of the image shows the brown positive nuclei and the light green negative nuclei.

Is there any chance for filtering out these marker pen artifacts? The effect is similar to looking through a welder's mask, and I'd wan't to somehow subtract the mask from the image. If I can't filter out the marker, I'll probably just need to manually crop out the portions having these artifacts.

Edit:

I've given high-pass filtering a try as suggested by Aaron. I used ImageJ's High-pass filter plugin, which creates a Gaussian blurred image from the source and then subtracts this blurred image from the original. The image is here (the radius / blur length or standard deviation of the Gaussian used was 20 pixels):

The artifact is almost gone and the nuclei are still recognizable, which is nice. The yellowish background color in the original (due to a calibration problem with the scanner) is also gone.

However, I'd also like to understand how this works. :) What does the Gaussian do to the RGB values of each pixel? How should I select the pixel radius for the Gaussian? Wikipedia's Gaussian Blur article wasn't all that clear to me. Are there other ways for implementing a high-pass filter for image data?

You could also try local contrast enhancement or similar. Fiji has a function called Enhance Local Contrast (CLAHE) which does this. Looks like high pass filtering works reasonably well for you though.

As for your second question. The Gaussian Blur wikipedia you mention refers to Low pass filter. For the low pass filter each pixel is average over the surrounding pixels, weighted by a Gaussian centred on the pixel. In your case you are doing a high pass filter. This can can be given by the original minus the low pass filter, which is probably the easiest way of thinking about it (IMO).


Possible to Derive Receive Signal Power from RSSI? My Set-Up

I am sending packets at a specific frequency and with a specific gain from one N-210 USRP to another of the same model.

NOTE: I have the transmitting USRP connected to a series of attenuators (so that I don't overwhelm the spectrum analyzer) and then to an RF splitter - one output is connected to a spectrum analyzer and the other output is connected to the receiving USRP. This is how I am reading the receive signal's magnitude.

My Problem

I am attempting to calibrate the receive signal power measurements (in dBm) that I can visually record from the spectrum analyzer with the RSSI measurements (in dB) that are output into a file by a script (I read them from the text output of the ofdmflexrframe_rx function). My plan is to come up with a linear correction factor to relate RSSI measurements to the receive signal power measurements from the spectrum analyzer.

I read on Wikipedia:

  RSSI is the relative received signal strength in a wireless environment, in arbitrary units.

I also read this from a response to a USRP user's question about RSSI measurements:

  [The] Received Signal Strength [Indicator is] always relative to some
  signal model, incorporating considered bandwidth, assumptions on the
  modulation scheme, duration of transmission, generally: It's a
  estimation of received signal strength based on some property of
  the received signal.

Now, I think that comparing the receive signal power levels on the spectrum analyzer with the RSSI measurements that function provides is not the right way to come up with a correction factor.

What is the difference between RSSI measurements and the dBm values that I record from the spectrum analyzer?
Is it possible to relate these RSSI measurements to the receive signal power levels that I record from the spectrum analyzer?
If I cannot relate these two measurements by a correction factor, is there a way to obtain the receive signal power measurements (in dBm) programmatically (GNUradio perhaps)? (NOTE: I don't believe there is any way to obtain the receive signal power measurements through the library that contains the function I linked above)

I also read this from a response to a USRP user's question about RSSI measurements:

    [The] Received Signal Strength [Indicator is] always relative to some signal model, incorporating considered bandwidth, assumptions on the modulation scheme, duration of transmission, generally: It's a estimation of received signal strength based on some property of the received signal.

Yes, I wrote that, and I still think it applies, so:

  What is the difference between RSSI measurements and the dBm values that I record from the spectrum analyzer?

Well, the RSSI is calculated based on some signal-model based estimator, whereas the dBm value on your spectrum analyzer is calculated with another metric – namely, the power measured going through a filter with the filter bandwidth you configured at your spectrum analyzer.

These two might, or might not be related, taking into account that even a signal with very low power spectral density values might be sufficiently good for something that has e.g. matched filters to extract the signal from noise, or any other form of processing gain.

So, your only option is to understand both estimators, with the Spectrum Analyzer estimator being relatively simple, and the flex estimator being free and open source (so you can just read the code ), and relate the two; so basically, you'll have to find, given a (class of) signal(s) $s(t)$, the functions $f_{SA}$, $f_{flex}$ and finally $g$ from

$$\begin{align}
P_\text{spectrum analyzer}\left(s(t)\right) &= f_{SA}\left(s(t)\right)\\
P_\text{flex frame sync RSSI}\left(s(t)\right) &= f_{flex}\left(s(t)\right)\\
\leadsto
P_\text{flex}'\left(P_\text{spectrum analyzer}\right) &= g(f_{SA}\left(s(t)\right))\text{ .}
\end{align}$$

I've just skimmed the liquiddsp code for the RSSI estimation, and the interesting part seems to happen in qdetector_cccf_execute_align.

However, this can never work without making a lot of assumptions. For example, if your signal happens to have a very narrow, but very strong single tone interferer, it will probably not break flex reception, and if that RSSI estimator is any good, not change RSSI very much, but it will drastically change the dBm value displayed by your spectrum analyzer. So if the question was

  Can I find an estimator that always gives me the same reading as the liquidDSP flex RSSI based on the dBm value displayed by a spectrum analyzer?

then the answer is a

  No. See above proof by counterexample (single tone interferer).

However:

  If I cannot relate these two measurements by a correction factor, is there a way to obtain the receive signal power measurements (in dBm) programmatically (GNUradio perhaps)?

Sure. Why not? You got the raw samples, don't you? So do a magnitude squared of them, and find a linear correction factor. That will work (have done so multiple times), because now you're actually measuring to things that relate to each other (signal amplitude as seen by your SDR frontend's ADC and signal power as seen by the spectrum analyzer simply have a quadratic relation).

Note that you will have to re-calibrate for different frequencies, sampling rates, if settable, analog frontend widths (only applies to B2xx, E3xx, and a few very old daughterboards for your N2xx currently; by the way, the TVRX* boards have always-on AGC, so they can't be used for this).

The relationship between input voltage and digital number amplitude is linear for quite a large range for the Ettus frontends – however, as soon as your input amplifier and mixers start reaching saturation, this will no longer apply, so you did good to use a few dB of attenuation.

If you happen to have a WBX, SBX, CBX or UBX daughterboard, there's IP2/IP3 measurements (along with effective vs set gain) available that will help you estimate the region of linearity on http://files.ettus.com/performance_data/.


Given a plot of both the magnitude $|H(\omega)|$ and its angle, How can you find the $H(\omega)$? I'm specifically trying to use an inverse Fourier Transform to find $h(t)$, but I'm finding it difficult to get $H(\omega)$ in the first place.

I'm under the impression from my textbook that $H(\omega) = |H(\omega)|e^{j\arg\{H(\omega)\}}$, and that $h(t)$ would just be the inverse Fourier Transform of that.

Yes you are right. Let's practically implement it with Matlab / Octave.

N = 8;
h = [1,2,3,4,4,3,2,1];    % just an impulse response (FIR)

H = fft(h,N);             % DFT H[k] of h[n] as samples of H(w) at w = 2*pi*k/N

Hm = abs(H);              % extract the magnitude of H[k]
Hp = angle(H);            % extract the phase angle of H[k]

Hr = Hm.*exp(j*Hp);       % Reconstruct H[k] from its magnitude and phase

hr = ifft( Hr, N);        % reconstruct h[n] from the reconstructed Hr[k]


Cascaded One-Zero Filter IR given the one-zero digital filter,

$$y[n] = x[n]\cdot \frac{1}{2} + x[n-1] \cdot \frac{1}{2}$$

Which has the transfer function:

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

and taking $M$ of these filters in series as a cascade, giving a filter $G$ with

$$G(z) = H(z)^{M}$$
Is there a known closed form of G's impulse response?

I have calculated one and am not sure if it's documented somewhere.
I know we can take the inverse fourier transform of $G(z)$, but we arrive at an integral/sum. I wonder about a closed form.

Thank you very much.

This is called a binomial filter. Its transfer function can be written as

$$H(z)=\frac{1}{2^M}\sum_{n=0}^M{M\choose n}z^{-n}\tag{1}$$

which gives

$$h[n]=\frac{1}{2^M}{M\choose n},\qquad 0\le n\le M\tag{2}$$

for its impulse response.

As mentioned in a comment, such a filter can be used to approximate a Gaussian filter. Also take a look at this related question and its answers.


Detecting Arabic verses spanning multiple lines I need to build a Quran app and I want to read out the verses when a user touches it. The problem I'm facing is that some verses may expand to one and half lines (highlighted red verse) or just fit in a quarter of a line (highlighted green verse). So adding each verse to textview or some other view wont work it seems.

I want to detect verses like the red ones in the second image. I have audio files for the verses so no need of text to speech conversion

This can be solved fairly straightforwardly with simple template matching. I don't know exactly how you have it set up, so I'll just describe the algorithm generally and use illustrations.

Observe that the verse numbers have a distinctive border that can easily be used to detect the start and end of a verse. So create a binarized template for that pattern and store it. Something like this:

Since the number of lines in a screen are known in advance (you're formatting the page) and each verse has a constant height, you can easily infer (algorithmically) where the Y coordinates for the centerlines of the verses should be on the screen. This demonstrates the idea:

When the user touches a verse, get the X-Y coordinates and snap the Y coordinate to the nearest verse center.
Then starting with the X coordinate, perform a simple template matching (cross-correlation) across that row. The first match (peak in the cross-correlation) in the forward direction (to the left), will be the end point for the verse. If there are no matches in the reverse direction (to the right), then move up one verse (which you can do, because you know the Y coordinate of the centerline) and repeat. The first match from the left end will be the start point of the verse. Similarly, if there is no forward match on the line, move down one line and repeat.

Here's a short illustration of the idea. The yellow box is where the user touches the verse. You then do the cross-correlation with your template and the blue circles will be the match.

I also use template matching in this answer, if you're interested in seeing it in action.
Once you've determined the start point for the verse, then use an Arabic text recognizer to infer the verse number inside that border and play the corresponding audio file.

Simpler solution:

A simpler solution, if you don't want to go through this is to store the X-Y coordinates of the verse starting points (keep it simple and use the center points) and once you get the coordinates of the user input, you can again snap it to the centerline and then walk backwards to see where the verse starts. This might have the advantage of being faster.

I didn't put this forward as the first solution because you seemed to reject a similar idea in the comments. In the end, it depends on your constraints — would you rather do computational work (template matching — which, by the way, also requires you to store the template) or using memory (storing coordinates).

If I were you, I'd probably go with this one, but the image processing solution can be fun to try.


Eigenvalues and eigenvectors of a 3D image Laplacian I need to compute the eigenvalues and eigenvectors of a 3D image Laplacian. I'm trying to evaluate the heat kernel on the 3D uniform grid (the uniform structure generated by the voxelized image) at different time values, to implement a Volumetric Heat Kernel Signature (please see the "Numerical computation" section). The domain I'm working on is not rectangular, so I have 1s on some grid points and 0s on other points, producing a binary 3D shape. I've been reading about this for 2D: in 2D, the laplacian of every grid point is calculated and the eigenvalues and eigenvectors of the resulting matrix are then calculated. How can I do this for 3D?All the information and examples I have read are for 2D images. The specific questions are:

How can I obtain a matrix from the 3D image I obtain after applying it a (3x3x3) discrete Laplace operator?
If there is no way to obtain a matrix, then I should compute the eigenvalues and eigenvectors of a 3D structured grid. What is the simplest way to do this? Could you recommend me any soft for doing this?

I would really appreciate any hint.

Thanks in advance.

This answer comes a little late, but I think that it's necessary to clear up some of the confusion about what the eigenstructure of a Laplacian is and how it is calculated.

First of all, it's important to stress that this is not about properties of the local kernel used for calculating discrete derivatives. Instead you have to understand the Laplacian as linear operator on a vector space whose elements are in your case the volumetric data sets.

That means applying the Laplace operator is a linear map that maps one vector (data set) to another vector (data set). In this context, the eigenvectors of the Laplacian are again vectors (data sets) in that very same vector space. So the answer to your question would be a set of volumetric data sets, in fact as many as the dimension of your vector space (i.e the number of independent voxels).

Let's consider a very simple example. Take a 1-dimensional image, i.e. a single row of pixels, and let's also only use very few pixels, namely 4. Then the two center pixels have to direct neighbours, while the first and last pixel only have one neighbour each.

$$1 - 2 - 3 - 4$$

With this pixel geometry, we can give the discrete Laplace operator's result for the two center pixels 2 and 3 as linear functions of the pixels' values:

$$l[2] = p[1] - 2*p[2] + p[3]$$
$$l[3] = p[2] - 2*p[3] + p[4]$$

The other two pixels 1 and 4 don't have enough direct neighbours to calculate the discrete second derivative. We could fix that by assuming that pixels 1 and 4 are direct neighbours, closing the topology to a circle and imposing what is called a circular boundary condition. Or we just take the second derivatives at the boundaries to vanish. Let's do both, but start with the cyclic boundary condition. So we have:

$$l[1] = p[4] - 2*p[1] + p[2]$$
$$l[4] = p[3] - 2*p[4] + p[1]$$

This map is linear and we can write it as a matrix equation, mapping the column vector $p := (p[1],p[2],p[3],p[4])$ to $l = (L[1],L[2],L[3],L[4])$ by multiplication with the matrix $M$.

$$ L = M \cdot p $$

We call this matrix the discrete representation of the Laplace operator, and for our case it is
$$M=
\left(
\begin{array}{cccc}
 -2 & 1 & 0 & 1 \\
 1 & -2 & 1 & 0 \\
 0 & 1 & -2 & 1 \\
 1 & 0 & 1 & -2 \\
\end{array}
\right)
$$

The eigenvectors of this matrix are
$$ v_1 = (-1,1,-1,1) $$
$$ v_2 = (0,-1,0,1) $$
$$ v_3 = (-1,0,1,0) $$
$$ v_4 = (1,1,1,1) $$
with the associated eigenvalues
$$ \lambda_1 = 4
, \lambda_1 = 2
, \lambda_1 = 2
, \lambda_1 = 0 $$

You may recognise these vectors as the basis vector of the discrete Fourier transform on this vector, and the eigenvalues as their discrete frequencies. This is true in general, and in fact the decomposition of a vector (or more generally, a function) into the eigenspectrum of a Laplace operator generalises the idea of the Fourier transform.

Now let's investigate what happens if we use the alternative boundary conditions where $l[1] = 0$ and $l[4]=0$. The Matrix $M$ is then

$$M=\left(
\begin{array}{cccc}
 0 & 0 & 0 & 0 \\
 1 & -2 & 1 & 0 \\
 0 & 1 & -2 & 1 \\
 0 & 0 & 0 & 0 \\
\end{array}
\right)$$

Cearly, this matrix has a different set of eigenvectors and eigenvalues. They're not as intuitive and insteresting, so I won't list them explicitly. It's worth noting however that we now get the eigenvalue $0$ twice, that means the eigensubspace where the laplacian vanishes is two dimensional.

So how do things change if we have a proper image instead of a single row of pixels? Not much. We only need to write down the Laplacian for each single pixel while considering the direct neighbour relationship, or topology, of the image. To make things a little trickier, let's go with an irregularly shaped two dimensional image.
$$
\begin{array}{ccccccccc}
    && 1 &-& 2 &-& 3 &&  \\
    && | & & | & & | \\
 4 &-& 5 &-& 6 &-& 7 &-& 8 \\
 | & & | & & | & & | & & | \\
 9 &-& 10 &-& 11 &-& 12 &-& 13 \\
    && |  & & |  &&  | & &  \\
    && 14 &-& 15 &-& 16 && \\
\end{array}
$$

Obviously now we have to take a 2-dimensional Laplacian by summing the second partial derivatives in horizontal and vertical direction. For that we require the two direct neighbours in each direction. The inner points $5,6,7,10,11,12$ therefore have a full 2-d laplace expansion. For point $5$ it looks like this for example:

$$l[5] = p[4] - 2 p[5] + p[6] + p[1] - 2 p[5] + p[10] \\= p[4] + p[6] + p[1] + p[10] - 4 p[5]$$

For the corner points $1,3,4,8,9,13,14,16$ we cannot construct a discrete second derivative so we use a boundary condition, e.g. $l[1]=0$

There are two points left, $2$ and $15$. Both have two direct neighbours in horizontal direction, but not in vertical direction. We can therefore apply a boundary condition that only affects the vertical direction, by setting the vertical second derivative to zero, while we evaluate the discrete second derivative horizontally and get $l[2]=p[1]-2p[2]+p[3]$ and likewise for $l[15]$.

Following this construction we get a linear equation for every point that relates it to the pixel values. Again, we write it as a matrix equation $L=Mp$, where the matrix in this case has $16\times 16$ entries. To be specific, it's
$$M=\left(
\begin{array}{cccccccccccccccc}
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 1 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 1 & -4 & 1 & 0 & 0 & 1 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\
 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
\end{array}
\right)$$

And again we can solve for the eigensystem of this matrix. A nice physical interpretation of the images that you get as eigenvectors is that they represent the vibrational modes of a membrane shaped like your image, with frequencies given by the eigenvalues.

You can easily step this game up to any number of dimensions, as long as you know the neighbourhood relations between your voxels. Simply formulate the individual linear equation like above, construct a matrix, find the eigensystem.

With the information gained from the eigenvectors of the Laplacian solving difference equations involving the discrete Laplacian can be greatly simplified. Once the eigenstructure is found, depending only on the geometry of the region, all data sets can be easily decomposed into the eigenbasis and the difference equations become trivial.


Necessary Conditions for stability in z domain? What are the necessary(must) conditions for stability in z domain?
I am sure about one(ROC must include unit circle)
Is there any other such condition which states that there shouldn't be any poles in ROC?

The one and only condition for BIBO stability of a 1D discrete-time system, in the z-domain, is that its transfer functions's ROC (region of convergence) should include the unit circle : $|z| =1$. Therefore, it's a necessary and sufficient condition for BIBO stability of a 1D SISO system.

There are no other conditions (to my knowledge).

EDIT [based on comments] :
By definition of ROC (region of convergence) there cannot be any poles inside the ROC region. This's related not with stability but with convergence. And also, causality is not related with stability too.

Properties of ROC is a different question to be answered.


Does MSAA make sense in audio / 1d data streams? In graphics, MSAA is a technique where  you take multiple samples of something per pixel and then average the values to get the final pixel value. The image quality improves as a result, and its considered a form of anti aliasing.  Does this technique have an equivalent in audio samples or other 1d data streams?

No. The reason is that MSAA works in the visual domain. Optically, we have hardwired edge detectors. MSAA eliminates unintended edges in generated images.

Acoustically, we do not have such edge detectors, and therefore there's no reason to remove them. Also, loudspeakers already filter out such edges mechanically.


Why Fourier transform is not sufficient and we have to use Laplace transform? Is there an easy way to explain the motivation behind the use of Laplace transform instead of Fourier transform?
Isn't that any periodic function can be represented by sines and cosines? - Why to introduce exponential idea?
Why not using differential equations with Fourier transform? An example would help.

*Asked the same question a while ago at math.stackexchange but no answers given.

The Laplace Transform is more representative of real systems that have a starting point, which is why the integral starts at 0, and also why the unit step function is generally talked about alongside the Laplace Transform. With the Laplace Transform, we can examine the transient and steady-state behavior of a system.

Using $e^{st}$ instead of $e^{iwt}$ allows us to examine different aspects of a physical system. The variable $s$ is complex, and if the real part was set to 0, it would reduce to a truncated Fourier Transform. The real part of $s$ is related to the amount of damping in the system. Also, with the Laplace Transform, a system's stability can be considered.

In short, Laplace is used to consider damping, stability, transient and steady-state behavior of a physical system (represented by a differential equation).


Why is Cheby 2 filter uncommon? I've been reading the Chebyshev filters chapter in The Scientist and Engineer's Guide to Digital Signal Processing and it says

"The Chebyshev filters discussed in this chapter are called type 1 filters,
meaning that the ripple is only allowed in the passband. In comparison, type 2 Chebyshev filters have ripple only in the stopband. Type 2 filters are
seldom used, and we won't discuss them. "

Why is it that Type 2 filters are uncommon? It seems to me that it's much more useful to have the ripple in the stopband (the bit you don't want) than in the passband (the bit you do want).

Have a look at this plot from this page at the Mathworks site:

The main issue is, as that page says:

  The frequency input to the Chebyshev Type II design function sets the beginning of the stopband rather than the end of the passband.

As you can see, the Chebyshev Type I filter has controllable passband gain out to the cutoff frequency.  That is not true for the Type II: Only the stopband gain is controllable there.  The passband roll-off for Type II makes it not very interesting for most purposes.


Negative dB signal to noise ratio Is it common to find systems having signal to noise ratio (SNR) less than 1 i.e. noise is more than signal? I am currently working with SNR of -5db to -20db. And how useful is it to work with sensor readings having negative SNR? I am asking it because the algorithm I am using to solve a problem works well with negative SNR but at positive and higher SNR, it is numerically unstable. This instability arises because of computation of hypergeometric function $_pF_q\left(;;\frac{z}{σ^2}\right)$. At high SNR, i.e. low $σ^2$, the argument $\frac{z}{σ^2}$ becomes very high, making the convergence of hypergeometric function very slow and the algorithm becomes numerically unstable.

Yes, GPS signals are below the noise floor. Link


Gibbs phenomenon in Hamming's digital filters In 'Digital Filters' by Hamming there is a cryptic section where he describes how the Gibbs phenomenon can be viewed as the displacement between the centers of two functions as they are convolved together. This is on pages 112 - 113 of the 3rd edition.

In the process of this he shows that truncating the Fourier Series is the same as multiplying the Fourier coefficients with a rectangle function. He then goes on to show that the function that has 2N+1 coefficients which are 1 is:

$h(\theta)$ = $\frac{sin(2N+1/2)\theta}{sin(\theta /2)}$

I'm confused: I thought the frequency response of a rectangle was a sinc function, and another page he shows this (when he derives the Lanczos smoothing factors).

Could anyone please clear this up for me?

I think the main issue is you are jumping ahead of yourself.  You probably remember or read somewhere that the Fourier transform of a $rect$ function is a $sinc$ function.  This is true; however, no where in this section does he mention Fourier transform!  In fact, what he is doing is not Fourier transform.

What he does in this section is to represent any periodic function as a Fourier series:

(1)    $f(\theta)=\frac{a_0}{2}+\sum_{k=0}^\infty{(a_k \cos k\theta + b_k \sin k\theta)}$

The key here is that this function doesn't contain every frequency.  It contains frequency 0 (the DC) and $\frac{k}{2\pi}$ where $k=1,2,3,...\infty$.  This is actually very important.  This function is always periodic with a period of $2\pi$.  The coefficients $a_k$ and $b_k$ are almost sampled values of the Fourier transform (but not quite, why? homework exercise! :p).  You have countably many of these coefficients.  Later in the book, you will learn that when you sample data in one space (frequency or time) it necessarily makes the counter part in the other space (time or frequency) periodic.

In the first case, he does the Lanczos smoothing derivation where he averages the function by running a rectangular window through it (convolving with rect).  What he shows is, not surprisingly, that the coefficients get multiplied by this term:

(2)    $\frac{\sin k\pi/N}{k\pi N}$

which should look very familiar to you, of course, because it is the $sinc$ function.  However, what you are missing is that $k$ is discrete!  It is actually a sampled version of the $sinc$ function.

Effectively, he convolves a function with a $rect$ and shows that the coefficients of the resulting Fourier series (read loosely as: sampled Fourier transform) is a sampled $sinc$ function.  No surprise there.  Convolution thm says convolution in time turns into multiplication in frequency.  Fourier transform, which you know, of $rect$ is $sinc$, so convolution by $rect$ is multiplication by $sinc$ in frequency space.

In the next section, he does something different.  He takes the Fourier series (read: sampled Fourier transform) and removes all the higher frequency coefficients.  In effect, he is taking the Fourier transform, multiplies by a $rect$, and then samples it.  For simplicity, he sets all the Fourier coefficients that did not get discarded to $1$.

What he's left with is this:

(3)    $h(\theta)=\frac{sin(N+1/2)\theta}{sin(\theta/2)}$

And you ask, why isn't this a $sinc$ function?  Can you answer it now?

The quick answer is because what is applied in frequency domain is not just a truncation, it's a truncation and a sampling operator.  What you know is that when you truncate (i.e. multiply by rect) in frequency domain, the time domain gets convolved with a $sinc$ (by Convolution thm and Fourier transform of $rect$), but this is without sampling.

As for why the formula looks the way it does, there are two ways to look at it.  The first, which he shows, is that you can just sum the Fourier series from $-N$ to $N$, and that's what you get.

The second, which is more profound and may come up later, is that (*) when you sample a function in one space (say frequency), the corresponding function in the other space (say time) becomes the sum of shifted versions of itself.  In fact, it probably won't come up exactly like that.  The typical scenario is that sampling time domain creates periodic replication in the frequency domain (btw: this is the reason for what people call aliasing).  However, you can apply the duality property of Fourier transform to get (*).

Does equation (3) make sense now?  It is periodic.  The closer you get to 0, the more it looks like a $sinc$ function.

So, an exercise for you is to derive equation (3) by sampling and truncating in frequency space and applying inverse transform.


Eight Sub Band Filter Design I have to process a RF signal band with bandwith between 75MHz and 100MHz to divide it in up to 8 contiguous subbands, each one with a bandwith programmable between 50kHz and 25MHz, in 25kHz steps.

The objective is to have the capability to attenuate each sub-band differently.

Transition skirts should be aprox. 6dB at 25kHz, on-band ripple less than 0.5dB and out-of-band attenuation better than -85dB at > 500kHz of the transition frequencies.
I would prefer a linear phase filter to ensure a constant group delay (constant, but as low as possible).

My first thoughts were to downconvert the RF signal band to a lower band (or to 0Hz to create a complex I+Q signal), then use 8 separate FIR filters, and finally upconvert the summing result.
However, I guess that the contiguity between adjacent filters may lead to some ways to reduce the computation cost.
May you give me some suggestions?

This sounds like the classical job of a polyphase filterbank, and a cascade of polyphase synthesizers. Yes, do this in complex baseband.

First, you use your filterbank to divide the 25 MHz in 25 kHz wide channels. Then, you use polyphase synthesizers to combine subsets these channels back to the channel widths you need, and then you use another one to combine these back to 25 MHz of bandwidth.

Another option would be to simply have a prepared table of filter coefficients for your potential subband widths, and just shift these in frequency (which is rather trivial) accordingly to get out your subbands of choice, then combine these back together with a bog-normal adder.

I'm assuming none of this requires that group delay in all these subbands is identical, for example. Otherwise, the design approach would have to amount for different filter delays.


What is the significance of a spike like autocorrelation function I have some questions and doubts to which the answers are difficult to extract from text books. Hence I have posted this question here. Help would be extremely beneficial. Thank you.

(1) The normalized ACF of white noise is similar to delta impulse function. What is the purpose and benefit of this spike nature in signal processing and statistics?

(2) If the length of the spike for white noise is higher in comaprison to the length of the spike for a random process then is there a significance and what soes this convey? Does the length of the spike have a special name and a role in applications to signal processing?

(3) A delta function is not a random process. Then why is the ACF of white noise, which is closer to the delta function preferred?

Consider a linear time-invariant system with unknown impulse response $h(t)$. If the input is a zero-mean wide-sense-stationary random process $\{X(t)\}$, then, with $\{Y(t)\}$ denoting the output process, we have that
$$R_{Y,X}(\tau) = E[Y(t)X(t+\tau)] = h\star R_x\big|_\tau.$$
So. if $\{X(t)\}$ is a white noise process wth autocorrelation function
$K\delta(t)$, then $R_{Y,X}(\tau) = Kh(t)$. Why is this at all important? Well, most processes are assumed to be ergodic, and so $R_{Y,X}(\tau)$ can be estimated by cross-correlating a sample path of the output with a sample path of the input with a delay of $\tau$. Thus, the spikiness of the autocorrelation function has a role in system identification of unknown systems.