Design causal FIR filter to approximate hilbert transform I'm attempting to design an FIR filter that approximates the Hilbert transform in order to get a $90^{\circ}$ phase shift and unity gain. However, I'm having trouble adjusting the filter to make it causal. I've already read this answer, which was very helpful. Here's what I've got so far.
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hamming

N = 50
H = hamming(2 * N + 1)

def coeff(n):
    if n == 0:
        return 0
    return H[int(n + N - 1)] * (2 / np.pi) * (np.sin(np.pi * n / 2) ** 2 / n)

def freq(w):
    return np.sum(
        [coeff(n) * np.exp(-1 * 1j * w * n) for n in range(-N, N + 1)]
    )

if __name__ == "__main__":
    freq_vals = np.linspace(0, np.pi, 1000)
    freq_resp = [np.abs(freq(w)) for w in freq_vals]
    ang_resp = [np.angle(freq(w)) for w in freq_vals]
    # plot frequency response
    plt.plot(freq_vals / np.pi, freq_resp)
    # plot phase shift
    plt.plot(freq_vals / np.pi, ang_resp)
    plt.show()

coeff computes an ideal hilbert transform coefficient using the equation
$$
h[n] = \begin{cases}
\frac{2}{\pi}\frac{\sin^2(\pi n/2)}{n} && n\neq 0,\\
0 && n = 0,\\
\end{cases}
$$
which I've taken from Discrete-Time Signal Processing (3e) p.959. I've restricted $n\in[-50,50]$. Additionally, it applies a Hamming window to reduce the Gibbs phenomenon (see the linked answer above).
freq then computes the frequency response for a provided frequency value, using the equation
$$
H(\omega) = \sum_{n=-\infty}^{\infty} h[n]e^{-in\omega}.
$$
Now, when I plot this (plt.plot lines), I get the magnitude and phase response I expect:
mag (x range is 0 to $\omega$, which I've normalized to 1):

phase:

However, I computed the impulse response values for $h[-50],h[-49],\ldots,h[50]$. So, this filter is not causal. I believe this filter should be time-invariant, so I attempted to simply shift it so that $n$ would now range from 0 to $2N+1$. Here's the adjusted freq function:
def freq(w):
    return np.sum(
        [coeff(n - N) * np.exp(-1 * 1j * w * n) for n in range(0, 2 * N + 1)]
    )

This produces the same frequency magnitude response, but with the following "incorrect" phase response:

I expect I've misunderstood something fundamental about this process but I'm not sure what. How can I make this filter causal with the desired phase response? Why does my attempt to make the filter causal keep the correct magnitude response but not phase response?

Basically, Hilbert transformers are, by definition, non-causal, always. Any zero-phase filter is non-causal. We can shift the impulse in time but the phase then becomes linear and not flat (a function of frequency).
To be implemented, the input signal must be delayed by the amount that is half the length of the filter minus one (N in this case). Or one can convolve the input signal with the impulse of the Hilbert (windowed, as you did, to minimize Gibbs). Never done the latter so check yourself.
BUT! There is a configuration that uses a pair of 4th order allpass filters that turn out to be pretty close to 90 degrees out of phase with one another. Unfortunately, the phase of each output differs from the original but depending on what you need, it can be hacked...
Reference: https://www.dspguide.com/ch19/4.htm, Design Hilbert filter at matlab (pure 90 phase and magnitude 0db)


detrending in real time I have a certain signal that has a linear trend that should be removed:

The problem is that the signal needs to be processed in real time(sampling rate is 20Hz).
If I only apply the detrending per window of say, 20 samples and then concatenate detrending results, the signal is disrupted.

I wonder if there is some efficient technique to detrend the data per window in real time and then concatenate it so the final signal will make sense? Or should I just detrend every new incoming data with the older data which is far from being efficient.

I would like to note that taking the entire signal and processing it offline using MATLAB's detrend works correctly. The problem occurs when breaking the signal into certain fixed size windows and then applying detrend on each.

Thanks,
Mike

What is the nature of the trend?

If the trend is due to DC drift (perhaps due to the signal acquisition apparatus), you may be able to remove it with a simple high-pass filter around 0.1 - 0.5 Hz. This would eliminate the need for any real-time post-processing.

DC drift gives a slowly rising/falling DC component, such as in the first figure here:


Filterpy Kalman Filter batch processing with multiple measurement sources In pythons module for kalman-filtering, filterpy, there is a function batch_filter() to batch filter a list of measurements that then can be used for RTS-smoothing. See the documentation here.

I want to smooth my measurements with a RTS-smoother like:

(mu, cov, _, _) = kf.batch_filter(list(np.array(centroids)))
(x, P, K, Pp) = kf.rts_smoother(mu, cov)

The problem is now that I have two measurements from two different sources with different measurement noises. The function batch_filter() can only process one source of measurements. When I calculate and save x and P for every time step with the sequence "predict, update with R1, update with R2", is this the same thing batch_filter() would do?

I am not use how you can success multiple sensors using batch filter, but, in sequential filtering you need to perform time update before each measurement. You can do it using the timetag information of each sensor measurement. You can modify the batch filter code however, I strongly recommend to write your own function depending on what you need.


x(t) and x(-t) are flipped (left to right) signals of each other $x(t)$ and $x(-t)$ are flipped (left to right) signals of each other.  but $x(2t-3)$ and $x(-2t+3)$ do not share that kind of relation as I thought incorrectly.  Can some one comment on this?  Is there a name for the relationship of this kind of signals?  Is it important enough to study this relation and give a name?

$x(t)$ and $x(-t)$ are flipped left-to-right about the point $t=0$.  that point of reflection at $t=0$ is determined by asking: "When is $t$ and $-t$ equal to the same value?"

$$ t = -t \qquad \qquad \implies t=0 $$

so it turns out that $x(2t-3)$ and $x(-2t+3)$ are also flipped left-to-right, but the point of reflection is not $t=0$.  So then, what value of $t$ is the point of reflection and how do we determine that value?  Like above, it is when the two expression of argument to the common function are equal to each other:

$$ 2t-3 = -2t+3 $$

solve for $t$ and you'll find your point of reflection.


Obtaining the frequency and phase spectra from fourier transform? For a question on my homework, I am given f(t) = e^(-a*abs(t)) and asked to find the Fourier transform of the function. I understand how to do that part, but then the answer requires me to "Compute the amplitude and phase spectra" of the function. My question is, what are these spectra and how do I obtain them from the given function or its Fourier transform?

They simply are the magnitude and angle of the Fourier transform.


Normalization of signal against reference I am trying to figure out if there is any standard way of normalizing a sequence such as a time series against a reference sequence. I am not an expert on signal processing but I was hoping someone would be able to point out a good way of doing this. Specifically, I am working on image analysis where I measure pixel intensities along linescans. This gives me two intensity profiles along space (from 2 channels). I am trying to normalize 1 channel against the 2nd one (reference).

I assume your reference image data y_reference is acquired when there is no light on the image sensor - that's the standard thing to do. Then the usual normalized signal is
y_normalized = (y_signal - y_reference) / y_max
where y_signal is the signal image data and y_max is the maximum possible pixel output. y_max can be measured by increasing the light level on the image sensor until the output no longer increases.


Band limited signals that are sparse I am looking for examples of signals that are band limited but also sparse. For instance, the spiking of the neurons can be modeled as banded in frequency domain and sparse in time domain. What other examples exists for such type of signals?

If a signal is sparse in the time domain, it has infinite support in the frequency domain, thus is not band-limited in the mathematical sense. However it might be so close to band-limited that the out-of-band spectrum disappears under the quantization or other noise floor.


Confusion regarding lacunarity of koch curves with same fractal dimension but different orders I'm trying to use the concept of lacunarity to classify images, as per how a few research papers have done it. To verify whether my equations and implementation are correct, I'm trying to verify the output of my program - for images with same fractal dimension, lacunarity is inversely proportional to the denseness of an image - denser the image, lower is its lacunarity.

For example, if I plot two triangular Koch curves, one with order 4 and another with order 6(order meaning iterations), both would have the same fractal dimension. Then you would expect the image with order 4 to have comparatively higher lacunarity since its less dense.

Now, after computing the lacunarity, there are two possible plots you can get - non-normalized and normalized lacunarity. When I plot the non-normalized lacunarity, I get the correct output - the order 4 curve has higher lacunarity. Whereas when I plot the normalized lacunarity, the plots exchange their positions and the order 6 curve ends up having higher lacunarity which seems incorrect.

The research papers that I'm referring are only using normalized lacunarity in their applications. I'm not sure where am I making an error.

Lacunarity is defined as follows -
For a m x n binary image, you define box size r, where r ranges from, for example, from 1 to 100.
Lacunarity is defined for every r. For a given r, it is computed as the ratio of (A+B)/B
Where A = Variance of X
B = Square of the mean of X

Where X = no. of 1s within a box of size r x r. So you compute mean and variance over all possible boxes of size r x r in the m x n matrix. Then r is varied from 1 to 100 for example, and hence lacunarity value is computed for 100 values of r, after which you can plot lacunarity.

Normalized lacunarity is obtained by dividing Lacunarity[r] by Lacunarity[1] for all r.

So I need help with 2 questions:
[1] Am I making an error in the concept of lacunarity being inversely proportional to denseness, or is it that I'm making some mistake while plotting normalized lacunarity itself?

[2] Is there a way I can verify that I am implementing lacunarity calculation correctly? For example, if there are standard images whose lacunarity plots are well known. I can plot lacunarity for those standard images and verify that I am getting a similar lacunarity plot.

I think this is a conceptual issue. Your images are not fractals, they are approximations to fractals. Lacunarity is a property of the fractal curve, not of its approximation. The Koch curve has a specific lacunarity.

Your program computes a measure that is inspired by lacunarity, from a given approximation to a fractal. This measure determines what fraction of the image can be considered a hole at a given scale.

Your normalized lacunarity I hadn’t seen before. But it normalizes the curve by the value at r=1, the smallest possible scale in the image. Obviously at this scale the difference between your two images is the largest, since at a larger scale they should be identical (the difference is in the fine details). By normalizing by this value, it is expected that the ratio elsewhere, where the differences are smaller, will reverse.


What is Faster Than Nyquist signaling? Faster than Nyquist signaling is used to improve the spectral efficiency by reducing the time spacing (relaxing the orthogonality constraint) to pack more data in the same channel while tolerating a certain intersymbol interference (ISI).

What is exactly faster than Nyquist? Is it by increasing the sampling rate higher than the Nyquist rate? But this doesn't create ISI. What is the relationship between the Nyquist rate and the spacing in the time domain?

Harry Nyquist made so many contributions that it's easy to get confused.

Related to sampling, Nyquist proved that a signal $s(t)$ bandlimited to $B$ Hz can be reconstructed from samples taken at a rate larger than $2B$. This is unrelated to Faster than Nyquist (FTN) signaling, though.

Related to communications, Nyquist also showed that the maximum  ISI-free pulse rate in a channel of bandwidth $B$ is $2B$ pulses per second, using a technique known as "orthogonal signaling".

Note that Nyquist didn't say that communication at a rate faster than $2B$ was impossible; he just said that you will have to handle the ISI somehow. For a long time, people didn't know how to do that, so $2B$ became engrained as a sort of actual bound.

One exception is partial response signaling, which introduces ISI in a controlled manner. However, in this case the main purpose is spectral shaping more than an increase in data rate.

In the 1970s, Mazo published a paper showing that, in some cases, the ISI can not only be handled, but it actually doesn't detract from a system's error rate performance. Start from an orthogonal signaling system operating at pulse rate $R=1/T$; that is, find an orthonormal set of pulses $p(t-kT)$ (for integers $k,l,m$) such that
$$
p(t-lT)p(t-mT) = \begin{cases}
1,\text{when $l=m$} \\
0,\text{otherwise.}
\end{cases}
$$
Now let's say that you want to transmit the data sequence $\lbrace a_0, a_1, a_2, \ldots, a_M \rbrace$. You can do that by transmitting
$$
s(t) = \sum_{k=0}^M a_k p(t-kT)
$$
and the receiver (ignoring noise for simplicity) can recover any $a_k$ by calculating
$$
a_k = \int_{-\infty}^{\infty} s(t) p(t-kT)\, dt.
$$

Mazo's proposal was to convert this orthogonal-signaling system to an FTN system by introducing a time-compression factor $0 < \tau < 1$, so that the transmitted signal becomes
$$
s(t) = \sum_{k=0}^M a_k p(t-k \tau T)
$$
Note that the bandwidth of $s(t)$ is still $B$, but now we transmit pulses at rate $1/(\tau T) > R$. Mazo's main result was to prove that, for certain pulses $p(t)$ and certain time-compression factors $\tau$, the bit-error rate performance of the accelerated system is the same as the orthogonal system.

At the time, the problem was that handling the ISI was impractical. The ISI acquires a certain structure that can be represented as a trellis, but in any non-trivial case the trellises were huge and well beyond the computational capabilities available.

It is only recently that some new techniques, along with faster computers, have made FTN somewhat practical. I recommend this tutorial paper:

J. B. Anderson, F. Rusek, and V. wall, “Faster-Than- Nyquist Signaling,” Proceedings of the IEEE, vol. 101, no. 8, pp. 1817–1830, Aug. 2013.

Mazo's original paper is:

J. E. Mazo, “Faster-Than-Nyquist Signaling,” Bell System Technical Journal, vol. 54, no. 8, 1975.


My Power Law (Gamma) Transformation Code Doesn't Compatible With Gamma Value I'm trying to apply Power Law (Gamma) Transformation. Formula is simple: s = c*rγ.
(c=1 , r= intensity , s= outputintensity). Here is detailed information: https://theailearner.com/2019/01/26/power-law-gamma-transformations/.
I tried code (witch C#) below. But it doesn't work properly.

 public void gammatransform(Bitmap CikisResmi)
    {
        int ResimGenisligi = Img.Width;
        int ResimYuksekligi = Img.Height;
        int Greyscale;
        double y =0.4;

        for (int i = 0; i < ResimGenisligi; i++)
        {
            for (int j = 0; j < ResimYuksekligi; j++)
            {
                OkunanRenk = Img.GetPixel(i, j);
                Greyscale = Convert.ToInt16(OkunanRenk.R * 0.3 + OkunanRenk.G * 0.58 + OkunanRenk.B * 0.12);
                Greyscale = Convert.ToInt32(Math.Pow(Greyscale, y));
                if (Greyscale > 255)
                {
                    Greyscale = 255;
                }
                CikisResmi.SetPixel(i, j, Color.FromArgb(Greyscale, Greyscale, Greyscale));

            }
        }
        return;

}

Gamma correction can be applied to values ranging between 0 and 1. Based on your code, it looks like there is either no conversion from 0..255 range to 0..1 range and back. If there range really is 0..1, the intermediate result is stored in an integer so it is an issue anyway.


position of null carriers in N=256 subcarier OFDM can anyone tell me the position of null subcarrier in $N=256$ subcarrier OFDM. since edges re assigned null to prevent interference from an adjacent symbol, i assuming that the subcarrier $-127$ to $-100$ then $0$ and then $100$-$128$ are null, is that correct assumption, please respond!

Positioning of null carriers in an OFDM signal is a protocol-dependent piece of information. Specifications will often call for some null carriers around the edge of the band (to allow for rolloff in the antialiasing filters in the front end). Another common place you'll find a null carrier is at the band center. This may be found in specifications that are targeted at low-cost (i.e. zero-IF) receiver technology. Zero-IF receivers' baseband DC offsets can cause large spurious responses in the tuned band center, hence the use of a null carrier there.

So, the short answer to your question is that it's unanswerable: it's dependent on the specifics of the system you're working with.


Help in understanding the formula of Signal-to-Noise-Ratio (SNR) - Part 1 Question 1:

Consider an Autoregressive model :
\begin{align}
y[n] &= y[n-1] + x[n]\\
z[n] &= y[n] + w[n]
\end{align}
where $y$ is the output observation, $x$ is a random input and $w$ is Additive White Gaussian Noise (AWGN). In general, I have seen in many research articles and tutorials that the SNR = $\frac{E[y^2]}{\sigma^2_w}$

How to arrive at this formula and what is the meaning of Expectation? Is it calculating the mean?

$\mathsf{SNR}$ (signal-to-noise ratio) is a generic term whose value can be defined in different ways by different people, and as long as one states clearly what is meant by $\mathsf{SNR}$ in a particular document, there is no confusion. Thus, there is no "arriving" at the formula
$\mathsf{SNR} = \frac{E[y^2]}{\sigma_w^2}$ at all: it is the definition
of the term $\mathsf{SNR}$ as it used in that particular document.  There is, of course, a general feeling among communications systems designers and analysts that the error probability $P_e$ should be a decreasing function of $\mathsf{SNR}$, but exactly which decreasing function it is depends on the definition of $\mathsf{SNR}$ in use. For example, someone could define $\mathsf{SNR}$ as $\frac{E[|y|]}{\sigma_w}$ or as
$\frac{\sqrt{E[y^2]}}{\sigma_w}$ or as $\frac{E_b}{N_0}$ (a quantity that I like to call the BEND ratio ("bit energy to noise density ratio")), because those definitions make more sense in that particular application, or lead to more memorable formulas. If $\mathsf{SNR}$ is defined to be the BEND ratio, then the bit error probability for coherent demodulation of
binary orthogonal FSK would be $P_e = Q(\sqrt{\mathsf{SNR}}\,)$ (where $Q(\cdot)$ is the complementary standard Gaussian CDF) while for differentially coherent demodulation of differentially encoded PSK, the bit error probability would be $p_e =\frac{1}{2}e^{-\mathsf{SNR}}$, etc.
Skipping over the fine details to the executive summary, large values of  $\mathsf{SNR}$ are better than small values of $\mathsf{SNR}$.

Since the $x[n]$'s are random (variables), so are the $y[n]$'s random
variables since they are sums of the $x[k]$'s. Note that
\begin{align}
y[n] &= x[n]+y[n-1]\\
&= x[n] + \left(x[n-1]+y[n-2]\right)\\
&= x[n]+x[n-1]+\left(x[n-2]+y[n-3]\right)\\
&\ddots\\
&= \sum_{k\leq n} x[k]
\end{align}
Thus, $z[n] = y[n]+w[n]$ is the sum of two independent random variables
and so
\begin{align}\require{cancel}E[(z[n])^2] &= E[(y[n])^2] + E[(w[n])^2]+2E[y[n]]\cancelto{0}{E[w[n]]}\\
&= E[(y[n])^2] + \sigma_w^2\tag{1}\label{1}\end{align}
where $E$ denotes expectation (as the OP suspected) and the
expectation is indeed calculating the mean, not of the random variable $z[n]$ but of the random variable $(z[n])^2$, and $\sigma^2$ is the variance of the AWGN.  So, the definition of $\mathsf{SNR}$ that the
OP is reading about, viz. $\frac{E[y^2]}{\sigma_w^2}$ is the ratio of the two terms in $(1)$. The numerator is the instantaneous power of the
output signal $y$ and the denominator is the instantaneous power of the AWGN.  This definition of $\mathsf{SNR}$ is deemed by the authors of
the paper that the OP is reading to be the most convenient for their
purposes. Others may choose to carry out similar analyses using different definitions of $\mathsf{SNR}$ and their formulas will appear to be different. But, arguments about "My $\mathsf{SNR}$ is larger than your
$\mathsf{SNR}$" are readily resolved by refusing to compare apples with oranges and looking at the important parameters: which system achieves smaller error probability, or needs smaller (peak or average) transmitter power, uses smaller bandwidth, etc. The rest is just semantics.


Python's $\tt resample$ vs $\tt resample\_poly$ vs $\tt decimate$ I have multiple signals each measured with a sampling rate 300Hz. They are mostly periodic, although there are some parts of some of the signals, that seem to be random and do not follow the rest of the pattern(could be noise). I want to downsample this signal and I have found these three functions: decimate, resample_poly, and resample.

I do not know much about signal processing, so I was wondering whether someone could explain in simple terms what is the difference between the three functions and if possible suggest which would be the most appropriate to use for my case.

Note that I would like to apply machine learning methods to classify new signals to 2 categories. I need to downsample to decrease the time it takes to classify signals.

Downsampling and upsampling are operations that change the sampling rate of a signal. Each one of them is composed of two steps, changing the sampling rate and filtering. Usually, the amount of change is expressed as a ratio.

When downsampling, we are trying to take the signal from some $Fs$ to some $Fs_n < Fs$. The key problem with doing this is that we are effectively reducing the available bandwidth. Therefore, if we were to leave any frequency $f > \frac{Fs_n}{2}$ within the resampled signal, we would get aliasing.

Aliasing is a byproduct of sampling and it distorts the original signal. In DSP, we are trying to minimise it.

Therefore, when Downsampling, we first apply a low pass filter to the signal at its $Fs$ rate. This filter suppresses all frequencies that might cause aliasing. After applying the filter, we apply decimation. Decimation is a very simple process by which we are throwing away samples. For example, to downsample by $\frac{1}{2}$ you throw away every other sample (i.e. keep 0,2,4,6,8,...).

When upsampling, we are trying to take the signal from some $Fs$ to some $Fs_n > Fs$. The key problem now is not that we can't fit the amount of original data to the new bandwidth but the exact opposite. We now do not have enough data.

Sampling works by examining the value of a continuous signal at discrete times, specifically $Fs$ times per second. But, we don't really know what exactly does the original signal do in between those discrete instances we are examining its value.

Therefore, when upsampling, we first upsample, which is an operation that inserts silence in between the samples and then apply interpolation to make up the signal values during the silence bits. Interpolation is usually another low pass filter.

Obviously, we can use both of these operations to change the $Fs$ by any given ratio. But, no matter what this ratio is we first upsample and then downsample (with each upsample, downsample step composed of their own steps mentioned above). Otherwise, you are throwing away information.

So, having said this:

Decimate does downsampling.
Resample does upsampling / downsampling and the associated filtering, entirely in the frequency domain, using the Fourier Transform. Because it is using the Fourier Transform, a key assumption is that the signal that is fed in its input is periodic. In other words, if you supplied a recording of the word "Blue is a co", the Fourier Transform assumes that in reality, the signal sounds like "Blue is a coBlue is a coBlue is a coBlue is a coBlue is a coBlue is a coBlue is a co"...to infinity. This explains the Note towards the end of that page.
Resample_poly does upsampling / downsampling and the associated filtering, using a polyphase filter](https://en.wikipedia.org/wiki/Polyphase_quadrature_filter). (Notice the series of operations applied as mentioned in the documentation).

2 and 3 are essentially the same thing with different techniques for reasons of computation speed. In the techniques mentioned, filtering and decimation and upsampling and interpolation are performed in one step because they are inherent to the technique.

As you can see, changing the sampling rate involves throwing away information or "diluting" the information that was already captured in some signal $x$ during its acquisition at some $Fs$. This is why I asked you what are you going to do after the decimation because under certain conditions, your classification/ clustering might be affected if those tiny little characteristics of signals that differentiate them from each other are smoothed out.

By what you are mentioning, if your highest non-negligible frequency component is at 50Hz, this means that (theoretically), you can reduce your $Fs$ by $\frac{1}{3}$. In practice, try not to get too close to a lowest $Fs$ of 100Hz. Maybe leave some lee-way, say 120Hz.

Hope this helps.


How to find the output of a system with the unit-step response and the input? Convolution is used to find the output when the input and the impulse response is known. I am curious if there is a similar operator, like convolution, or a method that can determine the output of a system. I am asking this question, right of the top of my head, so I don't have any concrete example to apply it to.

If I am not wrong, one can use deconvolution to find the impulse response and go from there, but this is not what I am searching for.

This can indeed be done. Write the input function $x(t)$ as

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

where $x^{\prime}(t)$ is the derivative of $x(t)$, and $u(t)$ is the unit step function. From (1), the system's response $y(t)$ to the input $x(t)$ can be written in terms of the step response $a(t)$:

$$y(t)=a(\infty)x(-\infty)+\int_{-\infty}^{\infty} x^{\prime}(\tau)a(t-\tau)d\tau\tag{2}$$

So, apart from the first constant term, the output signal $y(t)$ is given by the convolution of the derivative of the input signal with the system's step response.

Equation (2) is valid because the system's response to the constant $x(-\infty)$ is $H(0)x(-\infty)$, where $H(0)$ is the system's frequency response at DC:

$$H(0)=\int_{-\infty}^{\infty}h(t)dt=a(\infty)$$

The integral on the right-hand side of (2) is the convolution

$$((x^{\prime}*u)*h)(t)=(x^{\prime}*(u*h))(t)=(x^{\prime}*a)(t)$$

EDIT: The discrete-time case is completely analogous. Write the input sequence as

$$x[n]=x[-\infty]+\sum_{k=-\infty}^n(x[k]-x[k-1])=x[-\infty]+\sum_{k=-\infty}^{\infty}(x[k]-x[k-1])u[n-k]$$

This results in the following expression for the output sequence in terms of the step response:

$$y[n]=x[-\infty]a[\infty]+\sum_{k=-\infty}^{\infty}(x[k]-x[k-1])a[n-k]$$

Again, apart from the constant first term, the output sequence is obtained by convolving the first order difference of the input sequence with the system's step response.


Mathematical model of a signal in Compressive Sampling Currently I am reading a paper on Compressive Sampling, and trying to understand each and every parts of it. When I came across the mathematical model of the signal in paper, I got confused.  I have mentioned the part of paper for your convenience

Consider the following mathematical model for a class of discrete multitone signals. Let $W/2$ be a positive integer that exceeds the highest frequency present in the continuous-time signal $f$. Fix a number $K$ that represents the active tones. The model contain each signal of the form
$$ f(t)=\sum_{\omega \in \Omega} a_{\omega} e^{-2\pi i \omega t},\hspace{5mm} for \hspace{2mm} t \in [0,1)$$
Here, $\Omega$ is a set of $K$ integer-valued frequencies that satisfies
$$ \Omega \subset \left \lbrace  0, \pm 1 , \pm 2, \ldots , \pm (W/2-1),W/2\right \rbrace$$
and
$$ \left\lbrace a_{\omega}:\omega \in \Omega\right\rbrace$$
is a set of complex-valued amplitudes. We focus on the case where the number $K$ of active tones is much smaller than the available bandwidth $W$.

Now my questions are:

I think there should be $\pm W/2$ in the set $\Omega$.
Does the signal model consider both the positive and negative value of frequency in a single signal?
What will be the signal $f(t)$ when $\Omega = \left\lbrace 0, \pm1, \pm2\right\rbrace$?

Thanks.

The model is the Fourier series of a signal up to the W/2th component, with sparse coefficients. That is, most of the set $\{a_k\}_{k=-W/2}^{W/2}$ is zero.

It is because the signal is bandlimited, and thus for a real signal, $a_{W/2}$ will be real. However, if this constraint is not met it is better to use $\pm W/2$
For a real valued signal, the coefficient of the negative frequency is just the conjugate of the positive frequency coefficient. It really just is a mathematical artifact. You could also write the model as

$$f(t) = \sum_{\omega\in\Omega} b_w \cos(2\pi \omega + \phi_\omega)$$

where $b_\omega$ is real and $\phi_\omega \in [0,2\pi)$ is a real phase value.

That is, the model represents the signal as the sum of sinusoids of different amplitudes and phases.

It will be the Fourier series approximation of the signal using just the first two frequencies.

You should take a look at this paper and the example code:

Towards a mathematical theory of super-resolution


FFT window phase syncronising I'm currently devising a solution for demodulating a DPSK signal in C#. So far I've got a working FFT with phase detection, BUT:

I now have an issue of aligning the window of the FFT to synchronise with each symbol transmitted on the audio stream - now I've done some manual testing and for a test transmission I rendered, and I have symbols with a window width of 19 samples.

I do understand that DPSK is designed to avoid a reference signal for demodulating the input signal but FFT window alignment is still needed.

Now the problem is that because this will be potentially used on amateur radio where transmission times are entirely random, the alignment of each symbol will have to be determined. The transmission I rendered has a clock timing header of 1024 symbols, but how can I be sure to get the alignment right?

For example, the header is like:

1 0 1 0 1 0 1 0 etc etc

and the phase of the sine wave rendered from this bit stream would be:

0, 90 0, 90, 0, 90 etc etc

How can I make sure that the receiver will determine the transmitted phase and not anything like

15, 105 ,15, 105, etc etc

as this would also introduce spurious phase calculations from the FFT detecting the higher-frequency harmonics due to the jump in phase.

Here are a few things that aren't right:

1) Regardless of the Carrier Synchronization:
 You are using 0 and 90 degrees for your signaling. You would need an I and Q channel in order to decode the signal if that were the case. For BPSK (Binary) you should be using 0 and 180 degrees. Only if you are using QPSK should you use 90 and 270 degrees, then again its just mimicking using an I and Q channels.

There is a system block diagram below that will make this more apparent.

2) DPSK does not require any phase carrier synchronization. You pay for it by halving your SNR.

Encoding DQPSK

(for BPSK, eliminate Q and sqrt(2)

I-inphase channel Bit values
Q-quadrature bit values
C = (I + J*Q)/sqrt(2); // recombined I and Q channels
for n = 0 to length( C )
    x_diff = x_diff(n-1) *C(n);  //where X(n) <0 == 0
//X_diff is the differentially encoded channels

Decoding DQPSK

I - retrieved I channel
Q - retrieved Q channel
C =I + j*Q;
for n = 0 to length( C )
    decoded(n) =C(n)/decoded(n-1);  //where X(n) <0 == 0

4) You should never see a phase value because of thresholding at the sampler block

Transmitter Block Diagram

Receiver Block Diagram

(Remember you can ignore anything with the Q channel path if you are using BPSK, ignore the Convolutional encoding block and the input string block this whole thing will be easier to understand) .

The Sampler block acts like a Analog to Digital Converter. It samples the output of the I and Q channels at the bitrate your transmitting at (100bits/s, 200bits/s etc) it will give you the AMPLITUDE of the signal recovered.

Because, hopefully, you've used NRZ encoding, (-1, +1) to signal 0 and 1, you must threshold your signal. All values below Zero are 0's and all values above zero are 1's.  If you've used RZ encoding that is (0, 1) instead of (-1, +1) for your signaling you should threshold at 0.5.

Therefore you should NEVER work with a phase value.

5) Additional Block Diagram Information:

This is purely to be thorough and give readers a better idea about what is going on.

 1) The Users inputs a string ("Hello World" for example) this gets converted into its
    Ascii binary value.
3)This then gets put through a convolutional encoder, this is for feedfoward error
  correction. This is to reduce the transmission errors due to noise in the
  communication channel. Like power-lines or someone talking.
4)The binary values are converted into Non-return to Zero Encoding (NRZ) this means
  that 0's are converted to -1's and 1's are converted to +1's.
5) differential encoding is applied to eliminate the need for carrier (discussed above)
    synchronization.
6)The Signal is split into the I and Q channels (all even bits go to Q. All odd bits
   go to the I channel)
7) They are then shaped into pulses by the wave generator.
8) An Inter-Symbol Interference filter is applied (root raised cosine filter).
9)The I and Q channels are multipled by cosine and sine (same frequency, just shifted
    90 degrees) then recombined and transmitted.
10)The Receiver process is the same thing as the transmitter in reverse except:
   Convolutional Decoder: Use the Viterbi algorithm.
   Detector Block - this is where the NRZ thresholding takes place, values greater
   than Zero  are 1's values less than Zero are 0's

6) Working with M-ary (like Bin-ary, or trin-ary)

If you want to work with more data points in QPSK, you need to generate your constellation diagram.

Remember that all points on a QPSK system should lie on a circle.

Think of each point as a combination of the I and Q channels and not as purely a Phase value used for signaling.


Resources on Solving Convex Optimization Problems in the Compress Sensing Field When I read papers of compressed sensing, sparse representation and whatever requiring optimization of a cost function, I just find the final results as an iterative equation or so which will converge after few iterations and solve the problem.

But, when I try to add some other constraints to the cost function, I cannot derive the iterative or whatever the solution is and I know that is for my poor mathematics.

Can anyone lead me to a book or so to read and understand how to solve these types of mathematics? What book chapter? What tutorial?

There are few options:

Stephen Boyd, Lieven Vandenberghe - Convex Optimization.
This is the classic in this field. Very well written book.
Also have a look on other papers of Boyd on similar subjects such as the The Alternating Direction Method of Multipliers (ADMM).
They also have a great MOOC Course Stanford Online CVX 101 - Convex Optimization.
Amir Beck - Introduction to Nonlinear Optimization - Theory, Algorithms and Applications & Amir Beck - First-Order Methods in Optimization.
Amir is a great teacher and his book are both deep, practical and easy to read.
CMU Statistics 36-725 - Convex Optimization: Fall 2016 by Ryan Tibshirani.
Videos are available on YouTube.

I think once you skim through those you'll be able to handle most of the cases in Convex Optimization.


Timing recovery in GNU Radio : Variable rate vs Fixed rate operation I have been wondering about the implementation of timing recovery blocks in gnuradio. In literature, for asynchronous sampling, all that is needed are error detector (MM, ML, Gardner etc), a loop filter (2nd order), interpolation control (modulo-1 counter etc) and an interpolator. An example of a timing recovery scheme with ML error detector is shown below.
I have two questions regarding the implementation in GNU Radio

All blocks I have seen (MM, Andy Walls Symbol Sync )seem to be approximating the number of samples per symbol and use it to update the rate of the block in the forecast function. My question is why do we need to use variable rate? Can't a scheme like the one below with fixed rate work? Isn't this supposed to be asynchronous sampling?
Out of curiosity: Can one use a third-order PLL instead? Or will this just be an overkill?

Regards,
M.

EDIT

Interpolation control:

Matlab code

Symbol timing offset and symbol frequency offset are two different distortions. Like a carrier phase offset and a carrier frequency offset. The difference is that a carrier frequency offset imposes itself on the demodulated signal on the scale of a symbol time quite harshly due to the high frequencies of operation. The symbol rates are comparatively low and hence symbol frequency offset manifests itself slowly within a few symbol times. This is in regards to our modelling. Any real oscillator will probably exhibit a random component over this first order approximation.

Now Michael Rice is considering a case of fixed symbol timing offset in this figure. So the loop is just finding the optimal clock phase after asynchronous sampling. Considering that he's using a second order loop, a simple mechanism of stuffing and deleting samples is needed to adapt it to a varying symbol timing phase offset. GNU radio code runs on real signals so it has to adapt to that symbol frequency offset plus any small random variations over a period of time.

Finally, a 3rd order loop would ideally generate a steady state zero error for a changing frequency offset (sometimes called drift) but not many systems encounter such a scenario with a few exceptions. Higher order loops are used less in practice due to more complex stability issues. See Gardner's book Phase Lock Techniques for details.


Keypoint orientation with the oFAST from Rublee's ORB I've just read the papers from Rublee et al. (ORB: an efﬁcient alternative to SIFT or SURF) and Rosten et al. (Machine learning for high-speed corner detection) (=fast-detector). Rublee tries to explain the oriented-FAST-Detector an its improvements.

But there is one thin I can't understand. He assumes that a detected corner is not the real centroid of this corner. He uses the image moments to determine the real centroid. He then computes the atan2 to get the orientation of the intensity patch.

The Question is: What will he do with this orientation? I don't see any sentence about a shift of a possible corner-point using the orientation and den centroid.

The point descriptor usually requires samples at some specific locations from around the detected point.

If you have point location refined and orientation assigned, you just shift and rotate the locations before computing point descriptor.

Some simple descriptors require samples at discrete locations (e.g. 15x15 pixel image patch) and applying shift+rotation will make the sample locations non-discrete (the patch is arbitrarily rotated and shifted) and you have to use bilinear or bicubic interpolation to take the pixel values.

Here is a sample of detected points - notice the corresponding patches are rotated:

To get e.g. 64x64 matrix of values for each patch, one has to sample points rotated accordingly.


What is the current "state of the art" in future audio sample prediction, given recent developments in "image inpainting" using deep learning? In the last few years, there have been many breakthroughs in the image processing world regarding repairing "damaged" images, images with corrupt pixels, or even reversing artefacts from formats such as JPEG, which have received media attention. I am trying to determine if these methodologies can similarly be used to predict future samples of audio, rather than pixels of an image.
Some examples:

https://www.engadget.com/2018-04-24-nvidia-ai-fixes-photos.html
https://www.dpreview.com/news/3422748291/this-neural-network-can-fix-damaged-photos-and-enhance-low-res-shots
https://heartbeat.fritz.ai/guide-to-image-inpainting-using-machine-learning-to-edit-and-correct-defects-in-photos-3c1b0e13bbd0

much of which derives from an NVIDIA thesis published in EECV 2018: Image Inpainting for Irregular Holes Using
Partial Convolutions, EECV 2018.
In theory, this basic method should have a direct analog to audio: where, rather than predicting missing pixels, we're trying to predict the next ms or so of future audio from the last second. So, it can really be thought of as a type of nonlinear prediction, which clearly performs very well.
QUESTION: What is the state of the art in using this method for audio processing rather than image processing? Is there any published literature on this? Does anyone have any references?
Below are some examples of the technique that I am talking about in the image processing realm. In all these situations, the missing pixel data is reconstructed from the remaining pixels, which is basically 2D sample prediction; I am hoping to get some insight into the 1D analog of this.

These Images are FAKE !!!
There is absolutely no way, no algorithm, no nothing that exist that can tell you whether the eyes of a girl are blue or black by just looking at her lips !!!
You can guess, something, good or bad, based on your prior training and the probabilities based on that, but that's just a guess.
The fact that the synthesized images look so natural is the key to deceive people into thinking that the algorithm correctly predicts the impossible to fix missing information.
No algorithms can tell whether she had a tattoo in her screwed part of the image !!!
What the algorithm is doing is that you feed it with correct information in many many other training images, from different angles, different illuminations etc. Then the algorithm will try at best (and may succeed just as painting artist would do) to replace the missing pixels from those of the available set... And if she had a tattoo in the mean time; forget about it !
Same is for audio...


Unexpected behaviour when using lfilter with initial filter delay multiple times I want to filter a very long signal in smaller parts, therefor I am currently using scipys lfilter with an initial filter delay doing multiple iterations of:

(signal_part_filtered, filter_delay) = lfilter(b, a, signal_part, zi = filter_delay)

The first few iterarions are going as expected (original is blue; filtered is red):

But after a few iterations the inital filter delay keeps getting bigger:

I think the issue might be the limited precision of the float datatype, is there any way to solve this?

The problem is not in the filtering process but already at the design stage. Your specifications are very difficult to realize because your desired band is at very low frequencies. With these specs, the butter routine returns an unstable filter (i.e., with two pole pairs outside the unit circle of the complex plane).

One thing you could try is reduce the filter order. Always check the maximum pole radius by computing the (magnitude of the) roots of the denominator polynomial. It could be that the design turns out to be OK for a lower filter order. Of course, your filter will be less steep.

Another, better approach would be to downsample your signal before filtering, so your pass band is not anymore in such a low frequency range (compared to the sampling frequency). This will avoid numerical problems in the filter design routine.


What is the best file format to save short signal samples? We have the need to save signal samples measured with an arbitary sampling rate. The samples are usually about one second in duration. Currently we save them in a text file format which consumes a lot of disk space and is very slow. Also this type of format doesn't suit for streaming applications.

I have been searching for a binary file format for this kind of data, but I'm not sure what factors I need to consider when choosing a format.

The sampling rate of our data varies from 1kHz to 152kHz and currently 16 bits would be enough to represent the data (varies roughly from -1000 to 1000). This is vibration data measured from different sensors if it makes a difference. Also we would need to save pre-calculated FFT's in the same format. Is this possible?

Just to point out, I have no education in signal processing. I have the task to implement the system to be used in the analysis.

Currently the best option seems to be WavPack, which is basically a compressed WAV file and supports sampling rates from 1 to 4.3GHz and the compression is lossless.

WavPack is released in 1998 so I'm thinking are there any newer or more efficient formats available today?

Let's write down a quick spec:

The samples are usually about one second in duration,  i.e. let's say the maximum length is 4s
The sampling rate of our data varies from 1kHz to 152kHz, i.e. max rate 152 kS/s
varies roughly from -1000 to 1000,  i.e. more than 2048, so let's head straight for 16 bit integers

That makes a maximum recording raw data size of

$$4\,\text s\cdot 152\,\frac{\text{kS}}{\text s}\cdot 2\,\frac{\text B}{\text S} = 1216\,\text{kB.}$$

In other words, by all modern standards for permanent storage and RAM of machines that you'd use to store recordings... nearly no data at all, and that's the maximum case. Even when you incorporate the pre-calculated FFTs (why?? Storage access of data in this length is usually much slower than computing them from the data when you need it), that's still not much data.

So, I would recommend not to use any specific compression at all, but simply write the data as raw int16 into files. Every, language, in this world has methods of reading raw data from files.

Regarding metadata like sampling rate, time of capture, and something like that:

That really depends on your usage scenario. If you know you'll keep the data with the metadata database, I'd quickly put up a separate table (be it in CSV file, be it in a SQLite file, or in a proper relational database).

If you know you'll distribute these files somewhere else, check whether you'd rather be using a custom header format (which usually boils down to a packed C struct in one way or another) is really advantageous, or whether a simple  standard scientific data file format like HDF5 with well-defined structures within wouldn't be better – that would allow you to skip writing your own header parser in any target language, and also, you could easily put your raw samples, and an arbitrary amount of secondary data into it. Drawback is the fact that you can't just start streaming at the first byte of the data stream without loading the file to a large extent into memory. However, again, with your very small files, that isn't really a problem. I've not worked on a PC with less than 4 GB of RAM in the last decade, and neither should someone who does signal analysis. If of those 4GB of RAM, 2GB are available, you could completely have around 2000 of your maximum size data files in RAM at once.


Issues in generating AR model with a constraint I am new to the topic of system identification and looking for a large Autoregressive (AR) model. Can anybody point out a large stable AR model which has more than 2 coefficients AND there should be a large number of zero coefficients?

This link presents an example of how to have nonzero coefficients in lag 1,4 for order, $p=4$ AR(4) model. But, the values are NaN !! So, how do I apply this technique so as to generate a large AR model with more number of zero coefficients and lesser number of non-zero coefficients ?

EDIT:

Looking for any stable AR($p$) model (preferably popular cited ones, if available) of higher order $p \ge 20$ where the number of non-zero coefficients are few (constraint) . If I generate a noisy sinusoidal wave given below, then

%Generate sine wave = A*sin(2*pi*f*t + phi)
t = linspace(0,1,1000);
A = 5;
f = 2;
phi = pi/8;
sinewave = A*sin(2*pi*f*t + phi);
noisy_sine=sinewave+0.5*randn(size(t));
subplot(1,2,1);
plot(t, sinewave)
hold on;
subplot(1,2,2);
plot(t,noisy_sine);

%Generate AR model
order =20;
ARCoeff = aryule(noisy_sine,order);

It returns AR model of 20 lags. But, all the coefficients are non-zero. How do I generate AR model of any order higher than 2 where the number of non-zero coefficients are lesser than the number of zero coefficients ?

I was able to produce the AR(p) model with much of zero-coefficients. I use Scilab rather than Matlab, therefore I'll just describe the technique:

Generate some appropriate data with the desired spectral properties - the AR model will mimic these spectral properties.
Fit the AR model of the desired order to your data.
Find the desired number of the coefficients with the lowest values and reset them to zero  (something like a(find(abs(a) < stdev(a) / 2)) = 0 - valid in Scilab).

The resulting model should (may?) still has a similar spectral properties & satisfy your requirements considering zero-coefficients.


Power Spectral Density of Finite Signals I have a basic confusion which I hope someone can help clear.

Say I have two signals:

$S_1=\sin(2\pi t) \forall t\in[0,+\infty]$
$S_2=\sin(2\pi t) \forall t\in[2,5]$, $S_2=0$ otherwise

Now let's say the signals are sampled from t=0 to t=10. What would be the correct power spectral density (PSD) at 1Hz, in discrete time?

How would the answer change if I modify the sampling time interval, assuming the sampled interval always encompasses at least a portion of t=2 to t=5?

This, hopefully, will help me understand how to get the correct PSD out of STFT.

Edited to remove the word window and replaced with interval to avoid confusion.

Assuming the systems are sampled and we are looking for the discrete power spectral density, meaning a relative indication of the power level of signal versus frequency, and as suggested in the comments, the sampling rate is an integer multiple of the 1 Hz period, then each case can be described as follows:

For case 1, the Discrete Fourier Transform, if scaled by N, would result in the coefficients of Euler's Identity for a sine wave (this is explained simply by realizing that the Fourier Transform is the correlation to each exponential frequency term $e^{j\omega t}$ and as Euler's identity illuminates, the sine wave consists of positive and negative exponential frequency components.

$$sin(2\pi t) = -0.5j e^{2\pi t} + 0.5j e^{2\pi t} $$

The each of the two bins that have a non-zero magnitude ($\pm 1$ Hz) would have a complex value of +/-0.5j.

The power for each bin is found by complex conjugate multiplication or 0.25 in each case.

If the sampling rate is NOT an integer multiple of the frequency present (1 Hz in this case), spectral leakage would result such that the power would be distributed across all the other bins. The nature of this distribution is dependent on the actual sampling rate used.

For case 2 We are now multiplying our original signal in time by a smaller rectangular window. If you take the DFT of that window (all zeros from 0 to <2, ones from 2 to 5 and all zeros from >5 to 10, with the samples dictated by the sampling rate used) the result will be an aliased Sinc function; the specific result will depend on the actual sampling rate used. However since multiplication in time is equivalent to convolution in frequency, this specific aliased Sinc function result would convolve in frequency with the result from case 1 (as a circular convolution) to the provide the modified power spectral density.


Overlap and windowing a signal for signal processing I’d like to digitally signal process an audio data. Therefore my idea was to:

Split the audio data into individual blocks of 16 samples each.
Each block has an overlap of 50% (8 samples are shared across each consecutive block).
Each block of data is multiplied with a window (hanning),
Do the FFT,
apply an Equalizer to the FFT data (Which I will implement later),
iFFT the data and recover the original signal back.

Due to the Equalizer, the output might differ a little from the input, but shouldn’t be distorted .

To this end, I implemented a MATLAB code,

% I am using an audio file sampled at 16kHz
[Input_sig_16k,Fs] = audioread(‘filename);

window_size = 16;        % individual block size
overlap_win = 8;         % 50% overlap
data_full = zeros(1,length(Input_sig));
w = hanning(16);         % hanning windowing
for ii = 1:(length(Input_sig)/win_size)-1

  data_range = ((ii-1)*overlap_win)+1:((ii-1)*overlap_win)+16;
%     [data_range(1) data_range(end)]
  data_block = Input_sig_16k(data_range);

  % windowing function
  data_block = data_block.*w;
%     plot(data_block,'-rx');
%     ii

 fft_sig = fft(data_block);                 %FFT of the signal

 % I will multiply the frequency data with an Equalizer

 ifft_sig = ifft(fft_sig);                  %iFFT of the signal

 data_full(((ii-1)*window_size)+1:((ii-1)*window_size)+16) = ifft_sig;

 end

subplot(2,1,1), plot(Input_sig_16k,'-og'), grid on;   % plot and compare both the signals
subplot(2,1,2), plot(data_full,'-xr'), grid on;

However the output data is quite distorted. I’m not sure what I need to do to recover the original signal. Could someone let me know the mistake behind what I have done.

So how do I obtain the correct output?

Don't overlap.  Don't use a Von Hann window.  Zero-pad by at least the length of the impulse response of your desired equalization filter.  Then use an overlap add or overlap save fast convolution method or algorithm to combine your zero-padded and filtered FFT/IFFT blocks.

Otherwise circular convolution artifacts will corrupt (severely distort) your attempts at filtering your signal.

If you want to analyze or visualize the audio, use separate overlapped and windowed FFT blocks, not the ones used for overlap add/save fast convolution filtering.


GPS Coarse Acquisition PRN Codes Are the actual 1023 bit sequences for each of the 32 GPS satellites available to consumers? I found the IS GPS 200 but it doesn't list the entire sequences. Are verified sequences available for download anywhere online?

I'm trying to validate my python CA PRN generator. Also, what does the last footnote (**** The two-tap coder utilized here is only an example implementation that generates a limited set of valid C/A codes.) mean? Which codes is it valid for?

I have a PRN generator that I have validated with live captured signals that is available on the Mathworks Exchange site at this address and equally runs in Octave (Update: I also pasted the core of this in a code block below):

https://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator

The two tap coder is as given in the diagram in that spec and copied in the figure below:

What they mean by “two tap coder” is the selection and addition of two outputs of the G2 generator in the diagram and that this is just one way to generate the sequences (a very convenient way). An alternate approach which I have done is to implement the generator with one LFSR. In this case you simply convolve the coefficients for the two generator polynomials. This results in a length 20 shift register and non-maximum length sequence (since it clearly can be factored), and will generate 1023 different sequences each 1023 long, depending on how you seed it. (So if you change the state of the 20 elements with any 20 consecutive chips from a PRN of interest, it will continue to run the code for that SV).    So when they mention it will only generate a limited selection of the Gold Codes, they mean it won’t generate all 1023 possible codes (because with the selection of two taps there are only $9!=42$ combinations -if I did that right, and if so then I highly suspect that the selection of which of the 1023 different codes to use was based on the assumption of using this 2-tap generator). The two-tap generator shown won’t generate all possible Gold Codes for these polynomials used but it will indeed generate all the GPS C/A codes.

This is because the Gold Code (which is the GPS C/A code) is formed by adding two maximum length sequences in GF(2) (basically x-or'ing the output). Doing that is identical to what could be done with one generator where the polynomials are multiplied together, and you multiply polynomials by convolving their coefficients.

The advantage of doing that with just one long LFSR would be the creation of a flexible code generator capable of creating any maximum length sequence up to $2^{20}-1$ including all the GPS C/A codes and many other compound codes.

The interesting point about the Gold Code is that if you delay one of the generators relative to the other by even 1 chip, the result will be a completely new 1023 length sequence. So since there are 1023 possible delays, there are 1023 total possible sequences that can be generated, each 1023 long.

Further an interesting point about maximum length sequences in general, is if you GF(2) add the sequence to a delayed version of itself, it will create the exact same sequence just at a completely different delay! Thus you see what is happening in the Gold Code Generator with the "two-tap coder" which is used to select the satellite; in the figure I gave we select taps 2 and 6 resulting in SV1 (which should match the first part of the table that you didn't post). The table you have shows the equivalent delay of the G2 coder for each two tap selection that is made.

Since a login is required for the MathExchange site, here is a copy of my code with the front-end error-trapping removed (for brevity):

function g=cacode(sv,fs)
% function G=CACODE(SV,FS)
% Generates 1023 length C/A Codes for GPS PRNs 1-37
%
%
% g: nx1023 matrix- with each PRN in each row with symbols 1 and 0
% sv: a row or column vector of the SV's to be generated
%    valid entries are 1 to 37
% fs: optional number of samples per chip (defaults to 1), fractional samples allowed, must be 1 or greater.
%
% For multiple samples per chip, function is a zero order hold.
%
%
% For example to generate the C/A codes for PRN 6 and PRN 12 use:
% g=cacode([6 12]),
% and to generate the C/A codes for PRN 6 and PRN 12 at 5 MHz use
% g=cacode([6 12],5/1.023)
%
%
% Dan Boschen 12-30-2007
% boschen@loglin.com

% table of C/A Code Tap Selection (sets delay for G2 generator)
tap=[2 6;
    3 7;
    4 8;
    5 9;
    1 9;
    2 10;
    1 8;
    2 9;
    3 10;
    2 3;
    3 4;
    5 6;
    6 7;
    7 8;
    8 9;
    9 10;
    1 4;
    2 5;
    3 6;
    4 7;
    5 8;
    6 9;
    1 3;
    4 6;
    5 7;
    6 8;
    7 9;
    8 10;
    1 6;
    2 7;
    3 8;
    4 9
    5 10
    4 10
    1 7
    2 8
    4 10];

% G1 LFSR: x^10+x^3+1
s=[0 0 1 0 0 0 0 0 0 1];
n=length(s);
g1=ones(1,n);   %initialization vector for G1
L=2^n-1;

% G2 LFSR: x^10+x^9+x^8+x^6+x^3+x^2+1
t=[0 1 1 0 0 1 0 1 1 1];
q=ones(1,n);    %initialization vector for G2

% generate C/A Code sequences:
tap_sel=tap(sv,:);
for inc=1:L
    g2(:,inc)=mod(sum(q(tap_sel),2),2);
    g(:,inc)=mod(g1(n)+g2(:,inc),2);
   g1=[mod(sum(g1.*s),2) g1(1:n-1)];
   q=[mod(sum(q.*t),2) q(1:n-1)];
end

%upsample to desired rate
if fs~=1
    %fractional upsampling with zero order hold
    index=0;
    for cnt = 1/fs:1/fs:L
        index=index+1;
        if ceil(cnt) > L   %traps a floating point error in index
            gfs(:,index)=g(:,L);
        else
            gfs(:,index)=g(:,ceil(cnt));
        end
    end
    g=gfs;
end


Relating phase difference in a signal with its frequency I have a 134 point complex signal. Using MATLAB, I have observed that the phase difference between any two adjacent samples of this signal is 2 radians i.e. the phase increases linearly from 0 to 268 radians over 134 samples.

How can I relate this phase difference between adjacent samples with the frequency of the same signal?

PS I am sampling at 422 Hz.

Frequency is the derivative of phase over time.

So, what's the time between two adjacent samples? That directly gives you the frequency.


Inverse $\mathcal Z$-transform problem $B(z)+B(-z) = 2c$, explain the structure of $b[n]$ and find the constraint of its length given that $c$ cannot be $0$.

This is a homework problem. "Explain the structure" means that $b[n]$ is zero for certain values of $n$, and has a certain shape.

I'm trying to take the inverse $\mathcal Z$-transform of $B(z)$ and $B(-z)$, but I'm not sure how the inverse $\mathcal Z$-transform of $B(-z)$ is related to $b[n]$, so I'm stuck...can anyone give me some advice?

Just use the definition of the $\mathcal{Z}$-transform:

$$B(z)=\sum_{n=-\infty}^{\infty}b[n]z^{-n}\tag{1}$$

from which it follows that

$$B(-z)=\sum_{n=-\infty}^{\infty}b[n](-z)^{-n}=\sum_{n=-\infty}^{\infty}b[n](-1)^nz^{-n}\tag{2}$$

From $(2)$, the sequence that corresponds to $B(-z)$ is $b[n](-1)^n$. Now taking the inverse $\mathcal{Z}$-transform of

$$B(z)+B(-z)=2c\tag{3}$$

gives

$$b[n]+b[n](-1)^n=2c\delta[n]\tag{4}$$

For odd $n$, the left side of $(4)$ equals zero. For even $n$ it equals $2b[n]$. Consequently, equation $(4)$ can only be satisfied if $b[n]=0$ for even $n$ except for $n=0$, where we require $b[0]=c$:

$$b[n]=\begin{cases}c,&n=0\\
0,&n\text{ even}\end{cases}$$


The history of FIR filters Does anyone know the name of the first DSP (sampled data) practitioner who discovered that real-valued symmetrical-coefficient FIR filters exhibit linear phase in the frequency domain?

Let me propose a majorizing date for this discovery, at least in written form: Rader, C. M. and Gold, B. (1967) Digital filter design techniques in the frequency domain, Proceedings of the IEEE. Such a discussion appears on page 164, under the name of "frequency sampling filters", involving the combination of "elemental" and "comb filter". They say "The phase versus frequency is exactly linear except for discontinuities of $\pi$ radians. These discontinuities occur where the magnitude response is zero".

A lower bound could be Lerner, R. M. (1964) Band-pass filters with linear phase, Proceedings of the IEEE, which seems to keep analog.


How to remove noise from this image without causing the crack to 'shrink'? I am currently embarking on an image processing project to detect cracks on surfaces using OpenCV and Visual Studio via C++. I have tried writing my program and at the moment, I am stuck because I have difficulties striking a balance between noise in an image and preventing the 'crack' from shrinking due to the dilate() operation. I am completely new to OpenCV and have only coded in C before, so I hope you will be patient and explain things slowly. :)

I am currently working on the photograph below, taken from Google Images, which has the parent URL: http://www.newindianexpress.com/states/kerala/article1330634.ece

I shall type my code below:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/photo/photo.hpp"
#include "highgui.h"
#include <cv.h>
#include <iostream>
#include <stdlib.h>
#include <stdio.h>

using namespace cv;
using namespace std;

int main( int argc, char** argv ){

    //Declarations
    double alpha; /**< Simple contrast control */ //unused
    int beta;  /**< Simple brightness control */ //unused
    float minVal;
    float maxVal;
    Mat bw_img;
    Mat contrast1;
    Mat gaussian_img;
    Mat filtered_img;
    Mat weighted_img;
    Mat adaptive_img;

    //Loads image from file and put image matrix into "original_photo"
    Mat original_img = imread ("crack.jpg", 1);

    //Convert color-to-gray image
    cvtColor( original_img, bw_img, CV_BGR2GRAY );
    //imshow("bw_img", bw_img);

    //Gaussian filtering
    int i = 7;
    GaussianBlur(bw_img,gaussian_img,Size(i,i),0,0);
    imshow("gaussian_img", gaussian_img);

    //HPF kernel
    Mat kernel = (Mat_<float>(3,3) <<
        0,  -1, 0,
        -1, 5, -1,
        0,  -1, 0);

    filter2D(gaussian_img, filtered_img,-1, kernel, Point(-1,-1),0);

    //Extract high-frequency features, e.g. edges  show up.
    addWeighted(filtered_img, 1.5, bw_img, -0.5, 0,weighted_img);
    imshow("weighted_img",weighted_img);

    //Perform adaptive thresholding
    adaptiveThreshold(weighted_img,adaptive_img,255,ADAPTIVE_THRESH_GAUSSIAN_C, THRESH_BINARY, 51, 20);
    imshow("adaptive_img",adaptive_img);

    //Construct structuring element
    Mat struc_elem = getStructuringElement(MORPH_RECT,Size(5,5));

    //Perform erosion first
    Mat eroded_img;
    erode(adaptive_img,eroded_img,struc_elem);
    imshow("eroded_img",eroded_img);

    //Perform dilation next
    Mat dilate_img;
    dilate(eroded_img,dilate_img,struc_elem);
    imshow("dilate_img",dilate_img);

    double otsu_thresh_val = threshold(dilate_img, dilate_img, 0, 255, CV_THRESH_BINARY | CV_THRESH_OTSU);
    double high_thresh_val  = otsu_thresh_val,lower_thresh_val = otsu_thresh_val * 0.5;

    Mat cannyOP;
    Canny( dilate_img, cannyOP, lower_thresh_val, high_thresh_val );
    imshow("cannyOP",cannyOP);

    //wait for long long time
    int c = waitKey(20000000000000);

    //if you get impatient, end program
    if (c == 27)
    {
    return 0;
    }
}

The images I got so far are rather noisy, with a lot of black specks after doing erosion followed by dilation. I shall post the photos later as I don't have enough reputation to post more links. Sorry.

I have tried to increase the blockSize and C values in adaptiveThreshold() function, from the current 51/20 to 101/50. The noises were mostly removed, but the crack feature got reduced in size. This image shall be posted later and referred to as Image 2.

I have also tried to maintain the adaptiveThreshold() figures for blockSize and C, and tried to perform a dilation followed by an erosion. The final image was less noisy, though there were still black specks, but the crack became slightly discontinuous. This image shall be posted later and referred to as Image 3.

May I ask, how should I remove noise via either adaptiveThreshold() function or the erode() and dilate() functions without causing the crack to become thinner or wider?

This is because if the crack is thinner, the crack will 'disappear' totally in the final output image.

Thank you! :)

As first preprocessing step use edge-aware smoothing methods before converting your image to binary. These methods do not modify the sharp boundaries that are larger than a parametric size. Following search reveals many methods:

https://scholar.google.com/scholar?q=edge+aware+smoothing

I am not quite sure if the method will work in your case, but it seems there is an edge aware smoothing filter in OpenCV:

http://docs.opencv.org/3.1.0/d3/d14/tutorial_ximgproc_disparity_filtering.html


Is there a problem with using the WAV format for non-audio signals? I have biological signals sampled at 1 kHz. The signals tend to behave fairly simply, sinusoidally. The software I've used to collect these data uses their proprietary format (*.adicht) and so I have to export to something I can work with in MATLAB, for example. I also have audio recordings from the same sessions for importing to praat, which accepts WAV files.

The formats available for export include WAV, tab-delineated text, and an unwieldy .mat conversion process that sometimes causes the program to crash.

My question is whether there is a disadvantage to exporting everything as WAVs, or a specific reason why I would want to use text over WAVs? Is there a drawback to using the WAV format for non-audio signals, or is it just unnecessary?

It will save me a lot of time if I export to a single format, and since I know I need WAV for audio analysis in praat, it would be easier just to export the non-sound signals that way as well (this program does not handle macros/batch scripting very well).

I think WAV will be a good format for your data as it supports:

Integer Hz sampling frequencies (like your 1000 Hz).
Integer sample formats (like 16-bit, 24-bit, 32-bit).
A floating point sample format (32-bit IEEE float).
Multiple channels.
Seeking (quickly reading the file at a specific time point in the audio).

Other nice things about WAV:

You can use FLAC to compress WAV files for storage and transfer.
WAV files can also contain non-PCM data (such as markers and metadata).


What does it mean for a Wavelet transform to commute with translations? Referencing this article here https://arxiv.org/pdf/1203.1513.pdf

It states "A wavelet transform commutes with translations, and is therefore not translation invariant". Now I understand why it is a problem that the result is not translation invariant, however, I'm confused as to why it is.

What does it mean for a transform to commute with translation and why does the Wavelet transform commute with translations (i.e. why is the Wavelet transformation shift invariant)?.

We can start from what is "shift invarient":

Transform G is shift invariant if - $$\forall x:\sigma^nG(x) = G(x)$$
$\sigma^n$ being shift by n. Examples for transforms that are invarient to shifts are histogram and the amplitude of Fourier transform.

Commuting with shift is - $$\forall x:\sigma^nG(x) = G(\sigma^nx)$$
So it can't be shift invariant (unless G(x) is constant).

Note: I'm actually not certain that wavelet transform commutes with shift, but it's most certainly is not shift inveriant.


Image Processing/ Computer Vision: How is separability of filters implemented? I am trying to convince myself that a separable 2D filter can be implemented via two 1D filters. So, I took the example of following Sobel filter:

$1/8\begin{bmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{bmatrix}$

Based on Szeliski

This can be separated as :

$1/2\begin{bmatrix}-1&0&1\end{bmatrix}$ and $1/2\begin{bmatrix}-1\\0\\1\end{bmatrix}$

Consider an image as shown below:

$\begin{bmatrix}1&2&3\\4&5&6\\7&8&9\end{bmatrix}$

Let us first perform a normal 2D convolution which involves:

1) flipping the kernel horizontally and vertically (i.e first interchange first and last rows, then interchange first and last columns).

2) applying the flipped kernel on the image (the border is assumed to be zero padded).

3) Normalizing values by dividing them by 8.

This yields the output:

$\begin{bmatrix}-1.125&-0.75&1.125\\-2.5&-1&2.5\\-2.625&-0.75&2.625\end{bmatrix}$

Now we can try applying 1D filters:

1) The 1D filter represented by the row vector is flipped (columns are interchanged).

2) The flipped filter is applied to the entire image.

3) The 1D filter represented by the column vector is flipped (rows are interchanged).

4) The flipped filter is applied to the result from first 1D convolution.

5) The values are normalized by dividing them by 4.

This yields the output:

$\begin{bmatrix}1.25&0.5&-1.25\\1.5&0&-1.5\\-1.25&-0.5&1.25\end{bmatrix}$

The two results are different. What's the mistake and how to perform these operations properly?

Szeliski isn't saying what you think they're saying.

The vertical part of the separable filter is:
$$
\mathbf{v_2} = 1/4\begin{bmatrix}1\\2\\1\end{bmatrix}
$$
not
$$
\mathbf{v_1} = 1/2\begin{bmatrix}-1\\0\\1\end{bmatrix}.
$$

You can see this by defining:
$$
\mathbf{h}= 1/2\begin{bmatrix}-1&0&1\end{bmatrix}
$$
and looking at

$$ \mathbf{K_1} = \mathbf{v_1} \mathbf{h} = \begin{bmatrix} 1/4 & 0 & -1/4\\
0 & 0 & 0\\
-1/4 & 0 & 1/4
  \end{bmatrix}$$

and

$$ \mathbf{K_2} = \mathbf{v_2} \mathbf{h}= \begin{bmatrix} -1/8 & 0 & 1/8\\
-1/4 & 0 & 1/4\\
-1/8 & 0 & 1/8
  \end{bmatrix}$$

You can see this effect in the image I included in your question: the vertical edges are highlighted, but the horizontal ones disappear. This is because the $\mathbf{v_2}$ filter is a low-pass filter, not enhancing edges.


Phase contribution of complex poles I am struggling with understanding the phase contribution of each individual pole. Let's say we have a system (minimum-phase system if it makes a difference) and it has poles located at:

and

What is the analytical relation between the phase contribution of these two complex poles?

This is easier to see graphically, on an s-plane, knowing that the phase is from the frequency response of the system and the frequency response is determined when we restrict s to be the $j\omega$ axis.

Consider the simple ratio of a single pole and a single zero expressed as

$$H(s) = \frac{s-z_1}{s-p_1}$$

Given s as any point on the s-plane, and $z_1$ and $p_1$ as specific points, we see that the above expression can be viewed as the ratio of two phasors, alternatively expressed as $$H(s) = \frac{A\angle\phi_1}{B\angle\phi_2}$$

To see this clearly, and how it will ultimately answer your question; first consider two arbitrary points on the complex plane, for example m and n, as shown in the graphic below:

The equation $m-n$ is a phasor with a magnitude and phase $K\angle \phi$. By subtracting n from m, we are shifting the origin to be at position n with relation to m. And thus we get a phasor with length K which is the cartesian distance between m and n, and the angle $\phi$ which is the angle of m-n with the origin at position n:

As stated initially, the phase contribution of each pole and zero is given by the frequency response, and the frequency response is given by restricting s to be the $j\omega$ (or imaginary) axis. So as we sweep s along the $j\omega$ axis, each zero in the numerator and each pole in the denominator contribute to the magnitude and phase of the overall response. The phase of the numerator will add, and the phase of the denominator will subtract (for example the simple case of 1/j = -j wich is $-\pi/2$ (-90°): As phasors, this expression is magnitude 1 at angle zero divided by magnitude 1 at angle $\pi/2$).

Thus the angle contribution for each pole and zero is as described in the graphic below. We see that as s approaches infinity, the angle contribution from the zero will approach +$\pi/2$, and the angle contribution from the pole (given it is in the denominator) will approach -$\pi/2$. So in this case at infinity the resulting angle will be angle 0, and at any given frequency we can solve for the angle analytically by using $s = j\omega$ for frequency $\omega$:

You can apply the above to your specific case, except $0.9e^{j\pi/4}$ is a pole in the right half plane, which would be unstable (will not converge). If you corrected that for a case where both poles are in the left half plane, then each pole will contribute phase as described above and the overall phase response will be the sum of the two. Thus for two poles in the left half plane, the phase will approach -$\pi$ (-180°) as s approaches infinity.


Different "zoom" after reconstruction from Radon transform using Fourier Slice Theorem I have implemented the procedure to reconstruct an image from its Radon transform involving the Fourier Slice Theorem (FST). After application of the FST, the data is transformed from polar coordinates to cartesian, by regridding.
It works fine, except the reconstructed image is "zoomed out" (see picture). What causes this difference, and how can I resolve it, so that the phantom in the original, an reconstruction are exactly the same size?

Code:
N = 128;
f = phantom(N);

% Radon Transform
theta = 0:1:180;
[g, xp] = radon(f, theta);

% Fourier Slice Theorem
ghatpol = fftshift(fft(ifftshift(g)));

% Transform to cartesian coordinates
phi = -theta*2*pi/360;
[pX, pY] = pol2cart(phi, xp);
x = -N/2:1:N/2-1;
[X, Y] = meshgrid(x, x);
fhatcart = griddata(pX, pY, ghatpol, X, Y);

frec = ifftshift(ifft2(fhatcart));

figure,
subplot(121), imagesc(f)
subplot(122), imagesc(abs(frec))

MathWorks' image processing toolbox describes the Radon transform as the projection of the image intensity along a radial line oriented at a specific angle. Now, let the set of parallel lines spaced 1 px apart is drawn in a fixed direction (let it be vertical lines) and you rotate the image bounded by the square 128x128 units (in which the phantom slice is inscribed) and for each 1 degree of rotation take the projections of the image (let it be horizontal projection of rotated image), integrating the image intensity (optical density, luminosity, etc) along the parallel lines of your fixed line set. Your object under test, in all its rotation orientations, is bounded not by the 128x128 square, but by a circle of 128·sqrt(2) diameter. To make sure that no region of the source image is missing, the Radon transform is performed for the circle of 128·sqrt(2) diameter.
Now, you perform the reconstruction of all the Radon transform data. The meaningful data are located in the circle of 128·sqrt(2) diameter, not in the 128x128 square, and you inscribe the reconstructed image into in the 128x128 square as is seen in the lines
x = -N/2:1:N/2-1; [X, Y] = meshgrid(x, x);
fhatcart = griddata(pX, pY, ghatpol, X, Y);
frec = ifftshift(ifft2(fhatcart));

of your code.
As the result, the structures of the source image are zoomed out by ~1.414.
You can adjust the meshgrid
x = -N/2:1.41421:N/2-1;

and the other parameters accordingly, or zoom in the reconstructed image in other ways, if you need to keep a 100% zoom.


What's the purpose of augmenting an image with a random "gradient"? I have the following Python code

import matplotlib.pyplot as plt
import numpy as np

def make_random_gradient(image_shape):
    x, y = np.meshgrid(np.linspace(0, 1, image_shape[1]),
                       np.linspace(0, 1, image_shape[0]))

    grad = x * np.random.uniform(-1, 1) + y * np.random.uniform(-1, 1)

    return grad

if __name__ == '__main__':
    grad = make_random_gradient((64, 80, 3))
    plt.figure()
    plt.imshow(grad)
    plt.show()

In particular, you should note the name of the function make_random_gradient. This function produces images similar to the following one.

In the past, I followed a computer vision course. I've also looked at the images in the Wikipedia article related to the image gradient. In a way, the function make_random_gradient seems to generate "random image gradients", even though they are not black and white, like the images in the linked Wikipedia article. I am looking for a rigorous and mathematical explanation of why the function above generates a "random image gradient" and precisely in which sense it generates a "gradient".

Furthermore, the "gradients" generated by this function are used to augment images. In particular, if x is an RGB image with shape (64, 80, 3), each channel i of x is augmented, in my case, in the following way

x[:, :, i] = (1 - alpha) * x[:, :, i] + alpha * grad

where alpha = np.random.uniform(0.05, 0.15) and grad is an object returned by make_random_gradient.

So, each channel is augmented with the grad using a convex combination. Why would one do this?

I am aware of certain image augmentation techniques, such as cropping, flipping, the addition of Gaussian noise, but I've never seen this technique. Has anyone ever seen it? If yes, why would one augment images in this way? I am looking for papers or, in general, works that I've used a similar technique and that provide a rationale behind it.

I've also quickly skimmed through the relevant pages of the paper A survey on Image Data Augmentation for Deep Learning, but the paper is a little bit too long.

It's not generating a gradient in the computer vision / image processing sense but in the graphics community sense. For example, if you use Gimp's gradient tool you will fill images with values interpolated between 2 ends.

My guess regarding the augmentation is that the author is trying to learn invariance to some shading effects (eg., when the sun makes an angle with object in the scene and you can see walls or other surfaces becoming gradually darker).


Claim: Given sampling time T, the hold operator is approximated at low frequency by a time delay of T/2 Can someone verify this statement?

The hold operator is assumed to be a zero order hold

Then the laplace transform of this hold operator has a well known form $(1-e^{-sT})/s$

Let w approach 0, we see that our resulting function is $(1-e^{\alpha T})/\alpha$ where $s = \alpha + j\omega$

What does it mean for the time hold operator to be approximated by time delay of T/2? How can delaying a signal be equivalent to holding that signal

A zero-order hold can be modeled by a rectangular impulse response of length $T$ (where $T$ is the sampling period) and height $1/T$, starting at $t=0$. The corresponding frequency response is

$$H(f)=\text{sinc}(fT)e^{-j2\pi f T/2}\tag{1}$$

with $\text{sinc(x)}=\sin(\pi x)/(\pi x)$. From (1) it can be seen that $H(f)$ has a linear phase response with a corresponding group delay of $T/2$. Now if the sampling frequency is very high compared to the bandwidth of the sampled signal, then the $\text{sinc}$ function can be approximated by a constant around $f=0$, or, equivalently, the rectangular impulse response can be approximated by an impulse at its center (which is $t=T/2$). So for low-frequency signals (compared to the sampling frequency), the zero-order hold operation can be approximated by a delay of $T/2$.

The following figure illustrates this. The blue line is a signal with a small bandwidth compared to the sampling frequency. The green curve is the result of a zero-order hold operation, where the signal is constant over 1 sampling period $T$. Finally, the red curve is the original signal delayed by $T/2$. If $T$ is very small, the delayed signal is closely approximated by the piecewise constant output of the zero-order hold.


Matlab - Segment cell nucleus from similar background I am trying to segment the nucleus of the cell in an image so I can extract its boundary. The nucleus and the background of the image are very similar so I am finding it difficult to do this using the something like erosion.

If anyone could suggest a way to remove the background and keep the nucleus that would be very helpful.

Image:

You need to work with different color space. Try YCbCr for instance

oimg = imread('http://i.stack.imgur.com/aYhrS.png');
img = rgb2ycbcr(oimg); % conver RGB to YCbCr color space
msk = img(:,:,3)>130;  % apply threshold on Cr channel to get segmentation mask

I picked up 130 as arbitrary threshold, you might need to adjust it a bit.
Showing the input image, the image region for which msk==1, and the background:

I'll leave it to you as an exercise to think why the Cr channel performed well for this image.


How to separate mixed signals, 3 identical gears meshing together I have 3 identical gears that are meshing together, on each shaft that is holding these gears I have an accelerometer recording vibration data.
Once I move to the frequency domain via MATLAB's fft I am unable to see any specific frequency that would describe the meshing frequency of my gears, it looks all the same.

Are there any methods I should be looking at to tackle such a challenge of identifying the gear meshing frequency and probably identify which gear is inducing what frequencies?

Also would it make any sense to run the motor while disconnecting the gears and record the vibration induced by the motor, and analyze that to filter out the motor's influence?

Yes it would make sense to analyse the sound of the vibrations without the gear interactions if you can, and then delete that sound profile from the overall sound using a noise remover that works with a sample of the noise only prior to the noise + recorded signal. You should say roughly what frequencies you expect from the gears based on the RPM and number of teeth, because it affects analysis methods. On this and many other programs there is "noise reduction using noise sample profiling" http://wiki.audacityteam.org/wiki/Noise_Reduction it is tremendously efficient.

There are no obvious correlations or patterns in the graph, it is true. Read about engine vibration analysis. I found a paper on engine sound graphs when i was researchign SPWVD, and indeed it is 100ds of times more precise than FFT and pseudo smoothed wigner ville distribution makes waterfall graphs which show the clearest patterns of sound that are physically measurable.

To seriously analyse your vibrations, you should not use coarse FFT, instead use pseudo smoothed wigner ville distribution which is available as an addon. this program can do it too, see the graphs: http://www.christoph-lauer.de/sonogram


Power spectral density and SNR for AWGN I have a GMSK modulator in MATLAB and wish to generate AWGN with a specific SNR, and I'm having some issues figuring out just how.

I have found this step-by-step process to generate zero-mean Gaussian noise with a specified power spectral density W0:

Form Gaussian distributed random variable: w = randn(1,N)
Map top zero mean: w = w - sum(w)/N
Compute average power: Pw = sum(w.^2)/N
Form w = w.*sqrt(W0*fs/Pw)

Note: As pointed out by Deve, steps 2 and 3 should not be necessary, since randn() generates a random variable with zero-mean and power 1. However, it is possible that the mean and power may be slightly off, since the vector is of finite length. For my application, neglecting these steps give a slightly less accurate SNR.

So, how do I determine which power spectral density I need to achieve a specified SNR? I understand that SNR is defined as $SNR=P_s/P_n$. Measuring the signal power using $P_s=\frac{1}{N}\sum_{n=0}^{N-1} s[n]^2$ and then isolating $P_n=P_s/SNR$ seems trivial, but how do I relate $P_n$ to $W_0$?

Thanks in advance!

Edit: Using Deves answer, I have written the following code:

s           = signal;
L           = length(s);

% Convert SNR from dB
SNRdB       = 10;
SNR         = 10^(SNRdB/10);

% Measure average power of signal
Ps          = sum(s.^2)/L;

% Calculate wanted noise power
Pn          = Ps/SNR;

% Generate random vector and ensure 0-mean
w           = randn(1,L);
w           = w-sum(w)/L;

% Scale to wanted power
Pw          = 1/L*sum(w.^2);
w           = sqrt(Pn/Pw).*w;

% Measure resulting SNR
Pw_meas     = 1/L*sum(w.^2);
SNRdB_meas  = 10*log10( Ps/Pw_meas ); % Gives 10.0000

And here's a plot in time domain for anyone interested:

It doesn't seem like a lot of noise - but I guess the calculation of SNRdB_meas is quite foolproof.

[PSD plot removed]

First a comment on your noise-generation process. The Matlab function randn() generates Gaussian noise with zero mean and mean power 1. So steps 2 and 3 are obsolete.

If you'd like to achieve a given SNR, then creating the noise signal with the wanted power is as simple as

w = w .* sqrt(P_n)

Where P_n is the mean noise power and can be calculated by the equation you've already found yourself. Actually you're doing the exact same thing in your step 4, from which you can also derive that

P_n = W0 * fs

assuming that fs is the sampling frequency of your system and thus equal to the bandwidth and using Pw = 1.


Need a similarity metric that describes these two curve as highly similar I have a large number of time series data and I need to be able to compare the similarity of the curves to a reference curve. The reference curve in question is shown in red in the figure, and I'd like to compute a similarity metric for the blue curve that tells me that these things are highly similar.

The reason that I want these things to be similar is that they are both characterized by a largely monotonic rise, followed by a near plateau at some upper value, and then followed by a largely monotonic drop back down to the initial level. I have computed some similarity measures such as the Frechet distance and DTW distances, using the following snippets of code

sim_score = similaritymeasures.frechet_dist(Z, Zref)
sim_score2, _ = fastdtw(Z, Zref, dist=euclidean)

The problem is is that the scores computed using these techniques do not tell me that these curves are similar. For example, the fastdtw score is shown in the title of the figure, and this score is actually quite large when compared to other comparison scores for curves which I would judge by eye to be much more dissimilar. For example, the following curve has roughly half the fastdtw distance, despite it not having as similar of a shape.

The issue seems to be that the temporal overlap is still being weighted much more heavily then more shape-related attributes. Now it's perhaps possible that a more nuanced understanding of what's possible with DTW would give me the outcome I need, perhaps using a more fully-featured DTW package like https://pypi.org/project/dtw-python/, but I'm unfamiliar with many of the options available in that package. Or maybe there's another method that would be better suited. Of course, I could brute force this by creating a piecewise continuous curve with a bunch of variable parameters describing a step-like shape, optimize them against the target curve, and take the best error match as the distance metric. If the optimization doesn't get trapped in a local minima, then this should work, but it's more involved and messier than I'd like to do if I can avoid it.

Would appreciate any other ideas you all might be able to suggest. Thanks in advance!

While I still don't yet understand all the options in the dtw-python package, it seems as though the relative ordering of my two examples is now as expected using the code snippet

alignment = dtw.dtw(y, yref, keep_internals=True)

This allows me to create the following alignment figure using that package.

Additionally, the relative ordering of my two examples is now as I would hope. The more similar example has a distance metric of 6 whereas the more dissimilar example has a distance metric of ~11. I will have to explore the internals of this package more to see if I could be doing something even better, but this looks promising. Curious what the differences are in implementation that make this ordering so much different from the other two packages used.


Zernike Moments' implementation in OpenCV The reason why Hu moments were implemented in OpenCV and why Zernike moments were not implemented is looking like their performance similar as stated in this paper.

As stated in the paper Zernike moments' advantage is their reconstruction facility.

If there is a reliable implementation of Zernike moments compatible with OpenCV I would like to test their shape matching performance with OpenCV's Hu moments shape matching performance and if possible reconstruction facility.

Thanks.

I haven't used the Zernike module from PeRL's myself, but check out their implementation here.

It's written as an OpenCV module and looks well optimized.

However! I have no idea what their license on this is, so definitely check before you integrate it into anything you're going to distribute or sell.


Demodulating wideband FM signal in Octave from recorded baseband IQ samples Assuming I have a baseband recording of 8 bit IQ samples at 2.4MSps centered on an FM station is it correct to say that after low pass-filtering the complex signal to 240KHz, taking discrete derivative of phase and computing the envelope of the result this envelope should be the baseband audio signal?

Will this approach work for FM stereo?

I'm trying to confirm that my understanding is correct in Octave based on IQ samples recorded from RTL-SDR tuned to an FM station, but I'm getting white noise.

pkg load signal;
% 2.4MSps IQ baseband signal of an FM station from RTL-SDR recorded in SDR#
samples = audioread("C:\\Users\\remote\\Downloads\\sdrsharp-x86\\SDRSharp_20191228_175834Z_101200000Hz_IQ.wav", "native");
csamples = arrayfun(@complex,samples(:,1),samples(:,2));
% low pass filter and decimate to get 480KSps (240KHz Widband FM)
csamples_filtered = decimate(csamples, 5);
% descrete derivatives of phase
filtered_derivative = conv(angle(csamples_filtered), [-1,1]);
% center amplitude around 0
filtered_derivative = filtered_derivative - (max(filtered_derivative)+min(filtered_derivative));
% rectify, low pass filter and decimate to get the envelope (48 KSps)
envelope = decimate(sqrt(filtered_derivative.*filtered_derivative), 10);
% scale to 0..255
sound_normalized = sound_normalized * (255/max(envelope));

player = audioplayer(sound_normalized, 48000, 8);
play(player);

The FFT plot of the original data in csamples (centered on Fc) corresponds to spectrum that is shown in SDR Sharp I recorded it in (and listened to the station), so I think the input data is good.

This should work if you further filter to 15KHz after your FM discriminator. FM stereo is actually sub-modulated as follows:

L+R: 0-15KHz

Pilot Tone: 19KHz

L-R: 23KHz-53KHz as DSB-SC on a 38 KHz sub-carrier

Where L is the left channel and R is the right channel.

Update: I found a nice graphic showing this from Keysight with further details on FM broadcasting at this link: http://rfmw.em.keysight.com/wireless/helpfiles/n7611b/Content/Main/FM_Broadcasting.htm

If the pilot tone is present, that indicates the station is transmitting stereo but with this approach mono receivers would still work as you would hear the combined left and right channels.

So the composite signal is formed by first AM modulating using Double-Sideband-Supressed-Carrier (DSB-SC) the L-R channels onto a 38 KHz sub-carrier (subtract R from L and then simply multiply this with a 38Khz tone).
The tone is frequency divided by two and this coherent 19Khz pilot tone is added to the DSB-SC L-R waveform. Finally R is added to L and this baseband signal is added.

This is the signal that is then FM modulated, so would be the signal you would expect to get if you are properly FM demodulating.

A stereo demodulator after FM demodulation locks onto the 19KHz pilot to create the AM sub carrier at 38Khz (2x) and with that easily demodulates the AM DSB-SC sub-modulation (by multiplying the two) to achieve L-R at basebnad, and then can extract L and R by summing and subtracting the L+R and L-R results.

If you want to do a simple test of audibly hearing the demodulation of the FM signal, filter out the 19 KHz pilot and higher frequencies after you have FM demodulated this composite signal and you will hear the combined L+R channels.


Histogram Equalization Ignoring zero in Python I wrote a small python code for histogram equalization (ignoring zero values) for an image but it is taking too long to run. Code is written with reference to wikipedia article on histogram equalization.

def hist_eq_updated(img):
    img_flatten = img.flatten()
    img_flat_nnzero = img_flatten[np.where(img_flatten!=0)]
    image_histogram, bins = np.histogram(img_flat_nnzero,256)
    #print(image_histogram,bins)
    cdf = image_histogram.cumsum() # cumulative distribution function
    #print(cdf)

    eq_img_flatten = np.zeros_like(img_flatten)
    cdf_min = min(cdf)
    cdf_max = max(cdf)
    for i in range(len(img_flatten)):
         val = img_flatten[i]
         ind = 0
         flag=False
         for j in range(len(bins)):
             if val==0:
                 flag=True
             if (val>=bins[j]) and (val<=bins[j+1]) or flag==True:
                 break;
             ind = ind + 1
         if flag==True:
             eq_img_flatten[i] = 0
         else:
             eq_img_flatten[i] = (((cdf[ind] - cdf_min)/(cdf_max-cdf_min))*254) + 1
    return eq_img_flatten.reshape(img.shape)

I am new to image processing in python and not familiar with any other "pythonic" way to do same. Thanks in advance !!!

EDIT - Sorry, I forgot to mention that image size is ~4000 * ~4000, image type is float32.

Actually, this is not a signal processing quesiton, and would have been better at stackoverflow. Anyway, you can use a look-up-table for this problem. The output value of a pixel does only depend on the pixel value, not its position in the image. So, it makes no sense to go through each pixel. Instead, just go through each possible pixel value. Then, later on transform your original pixels to the new pixels using this look-up table

import skimage
def hist_eq_updated(img):
    img_flatten = img.flatten()
    img_flat_nnzero = img_flatten[np.where(img_flatten!=0)]
    image_histogram, bins = np.histogram(img_flat_nnzero,256)
    #print(image_histogram,bins)
    cdf = image_histogram.cumsum() # cumulative distribution function
    #print(cdf)

    eq_img_flatten = np.zeros_like(img_flatten)
    cdf_min = min(cdf)
    cdf_max = max(cdf)
    LUT = np.zeros(256, dtype=img.dtype)
    for val in range(256):
         ind = 0
         flag=False
         for j in range(len(bins)):
             if val==0:
                 flag=True
             if (val>=bins[j]) and (val<=bins[j+1]) or flag==True:
                 break;
             ind = ind + 1
         if flag==True:
             LUT[val] = 0
         else:
             LUT[val] = (((cdf[ind] - cdf_min)/(cdf_max-cdf_min))*254) + 1
    return LUT[img]

Img = skimage.io.imread("http://upload.wikimedia.org/wikipedia/commons/thumb/4/4a/Portrait_elisha_gray.jpg/220px-Portrait_elisha_gray.jpg")

plt.figure()
plt.subplot(121)
plt.imshow(Img, interpolation='none', cmap='gray')
plt.subplot(122)
plt.imshow(hist_eq_updated(Img), interpolation='none', cmap='gray')