generating log mel spectrogram using librosa

The mel spectrogram additionally includes a step of projecting (power of) STFT bins onto Mel-frequency bins via a Mel filterbank; I don't have access to path so I made demo on exponential chirp:

You can visualize the kind of projection taking place by plotting the mel basis:

Note in general the two won't look alike unless filters.mel are carefully selected (nor do they have to).

Code
import numpy as np
import librosa
import librosa.display
import matplotlib.pyplot as plt

def echirp(N, fmin=.1, fmax=None, tmin=0, tmax=1):
    """https://overlordgolddragon.github.io/test-signals/ (bottom)"""
    fmax = fmax or N // 2
    t = np.linspace(tmin, tmax, N)

    a = (fmin**tmax / fmax**tmin) ** (1/(tmax - tmin))
    b = fmax**(1/tmax) * (1/a)**(1/tmax)

    phi = 2*np.pi * (a/np.log(b)) * (b**t - b**tmin)
    return np.cos(phi)

#%%
N = 2**16
n_fft, hop_length, n_mels, sr = 2048, 512, 128, 16384
x = echirp(N)

#%%
mel_spectrogram = librosa.feature.melspectrogram(
    x, sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mels)
log_mel_spectrogram = librosa.power_to_db(mel_spectrogram)
librosa.display.specshow(log_mel_spectrogram, x_axis="time", y_axis="mel", sr=sr)
plt.show()

#%%
X = librosa.stft(x, hop_length=hop_length)
mel_basis = librosa.filters.mel(sr, n_fft, n_mels=n_mels)
Xmel = np.dot(mel_basis, abs(X)**2)
Xdb = librosa.power_to_db(Xmel)
librosa.display.specshow(Xdb, x_axis='time', y_axis='hz', sr=sr)
plt.show()

#%%
mel_basis = librosa.filters.mel(sr, n_fft, n_mels=10)
plt.plot(mel_basis.T, color='tab:blue')
kw = dict(weight='bold', fontsize=16, loc='left')
plt.title("n_mels=10", **kw)
plt.show()

mel_basis = librosa.filters.mel(sr, n_fft, n_mels=n_mels)
plt.plot(mel_basis.T, color='tab:blue')
kw = dict(weight='bold', fontsize=16, loc='left')
plt.title("n_mels=%s" % n_mels, **kw)
plt.show()


Unambiguously detecting (not correcting) between $t+1$ and $2t$ errors in a Reed-Solomon message

The Berlekamp-Massey algorithm and the extended Euclidean algorithm (both further extended with the Chien search and the Forney calculation of error values) for decoding BCH and RS codes (hereinafter referred to as the decoder) have the following characteristic that is not well understood by many people:

  If there is a codeword $\hat{\mathbf C}$ that differs from the received word $\mathbf R$ in $t$ or fewer locations, then the decoder finds $\hat{\mathbf C}$.

Note the complete absence of any mention of $\hat{\mathbf C}$ being the correct (i.e. transmitted) codeword $\mathbf C$.  Now, if $\mathbf R$ does happen to differ from $\mathbf C$ in $t$ or fewer locations, then the   decoder does indeed find $\mathbf C$. Note that since the minimum distance is $2t+1$, there is no competing codeword $\mathbf C^\prime$ that is also within distance $t$ of $\mathbf R$ to lead the decoder astray: there is at most one codeword within distance $t$ from any given $\mathbf R$ and when such a codeword exists, the decoder finds it. But just because the decoder finds a codeword, there is no guarantee that it is the correct codeword $\mathbf C$; it might be some $\hat{\mathbf C} \neq \mathbf C$ that happens to be at distance $\leq t$ from $\mathbf R$. (The decoder is following the sailor's philosophy of "When I'm not near the girl I love, I love the girl I'm near" with the addendum that on the high seas, no girls are near).  What it is correct to say is that (with the usual assumption that fewer errors are more likely to have occurred than more errors)  with high probability the $\hat{\mathbf C}$ that the decoder finds is indeed the transmitted codeword $\mathbf C$.

But the question is what happens when there have been more than $t$ errors. Can such an occurrence be unambiguously detected? The answer is NO for the reasons described above: the decoder may well find a $\hat{\mathbf C}$ that is different from the transmitted codeword $\mathbf C$ exactly as the decoder is designed to do because $\hat{\mathbf C}$ is at distance $\leq t$ from $\mathbf R$ and the poor schmuck of a system designer will have no idea that s/he has just bought a pig in a poke.  Once again, probability comes to the rescue. In the (fairly unlikely) event that more than $t$ errors have occurred, it is very much more often the case that the decoder fails to decode in several possible ways (which may or may not be detected by the system designer depending on how many safeguards have been built into the system) than that the decoder produces an incorrect codeword.

The Berlekamp-Massey algorithm or the extended Euclidean algorithm find an error-locator polynomial $\Lambda(z)$ of degree $e$ and an error-evaluator polynomial $\Omega(z)$ where the inverse roots of $\Lambda(z)$ (as found by the Chien search) are the $e$ error locations and the Forney error-evaluation formula gives the error values at these $e$ error locations. When more than $t$ errors have occurred, it is usually the case that $\Lambda(z)$ has fewer than $e$ roots and/or that the Forney error-evaluation formula returns a zero error value at a purported error location that the Chien search has found. Both of these conditions indicate that something is awry but testing for these conditions is expensive in hardware implementations and time-consuming in software implementations. The simplest test for detecting a decoder failure is to let the Chien search and Forney formula do their jobs, and the recompute the syndrome of whatever the decoder spits out.  If the recomputed syndrome is zero, then the decoder output is a valid codeword - whether it is $\mathbf C$ or an impostor $\hat{\mathbf C}$ is unknown - but if the recomputed syndrome is nonzero, then the system designer knows that the decoder output is not a valid codeword and that the decoder has actually failed to decode.  But again, such a test is expensive to implement and very often skipped by system designers who know about it, and, of course, not even thought about by those who don't know of the possibility of decoder failure that can be detected.

For a well-designed RS decoder system, it is the case that
$$P_{\text{decoder output correct}} \gg P_{\text{decoder failure}} \gg P_{\text{decoder output incorrect}}.$$
The decoding spheres (the set of vectors or possible received words $\mathbf R$ that are at distance at most $t$ from some codeword and thus lead to a successful decoding, whether correct or incorrect) have very small total volume and so it is far more likely than not that when more than $t$ errors have occurred, $\mathbf R$ is not in any of the decoding spheres, and the decoder will fail to decode and spit out garbage. Whether it is worthwhile detecting such failures is up to the system designer.
It is sometimes the case that $1-P_{\text{decoder output correct}}$ is a perfectly acceptable word error rate and so ignoring decoder failure is an option, but in high-performance systems where very low error rates are required, detecting decoder failure and asking for a
re-transmission of a codeword that was badly munged up in transmission allows the achievement of the desired low error rates at the expense of increased delay.


Maybe including some code containing an minimum working example would help. I didn't see any problem.

I tried this because I am interested in related questions and recently learned about the Hilbert transform. I also added some noise because noise is always there, but often neglected in answers to questions like this (a pet peeve of mine), even though it sometimes completely defeats certain methods. Here's what I get:

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

# Parameters
N = 1300
t = np.arange(0,N)
f = 1.0/20

# Discontinuous amplitude envelope
A = np.ones(t.size)
A[1*t.size/5:2*t.size/5] = .2

# Time signal with noise
y = A*np.sin(2*np.pi*f*t)
y += (np.random.rand(N)*2-1)*.2

# Filter noise
filt = True
if filt :
    cutoff = int(N*f*1.5)
    yrms1 = (y**2).mean()**.5
    yF = np.fft.fft(y)
    yF[cutoff:]=0               # Filter
    y = np.real(np.fft.ifft(yF))
    yrms2 = (y**2).mean()**.5   # Keep signal energy
    y *= yrms1/yrms2

# Calculate envelope
envelope = np.abs(sig.hilbert(y));

# Plots
plt.figure(1, figsize=[10,7], dpi=100)
plt.clf()
plt.plot(A, label='Actual envelope')
plt.plot(y, label='Time signal')
plt.plot(envelope, '.', label='Calculated envelope')
plt.legend(loc='lower right')
plt.show()
plt.draw()
plt.pause(.1)

EDIT:
The demonstrates performance in a sliding window. Just replace the line with the envelope calculation above with this code. Note that the envelope at each location is only ever calculated using N_win surrounding points.

envelope = np.ones(N)*np.nan
N_win = 50      # needs at least twice the period or so
lag = 20        # needs to be >= period.
for i in range(N-N_win) :
    envelope[N_win+i-lag] = np.abs(sig.hilbert(y[i:i+N_win]))[-lag]

Note, that for my implementation of this to work well, you need to be looking ahead of where you are calculating the envelope by about 1 period. For real time processing, that means your 'envelope' signal will have to lag behind your main signal by at least one period. This is because of edge effects of the Hilbert transform. I don't think there is really any way around this.


A question about sampling rate of cosine signal

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


I would guess that the most prevalent difference between accents (of the same language) is how the vowels are pronounced. Humans distinguish vowels mainly by the frequency of the formant, i.e. regions of high energy in the harmonic spectrum.

A good way to display this is the "formant map" which plots the frequency region of the second formant against the frequency region of the first. See for example: http://www.phy.davidson.edu/fachome/dmb/digitalspeech/formantmap.gif

In this 2D map every map has a specific region. For different accent the centroids of this regions probably shift around. That is certainly true for different languages. An interesting approach would be to measure this map for, say, an Australian vs. an American English speaker and look for systematic differences.


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


Channel equalization affect on input signal

This appears to be a homework problem so I will avoid providing the complete solution here, but the following should suggest a path toward a solution:
Typically this problem is done when the channel is NOT known and results in a least squared approach given the effects of the added noise $v(n)$. Solving for the equalizer coefficients when the channel is a known quantity results in a similar process to solving for inverse convolution using a matrix equation approach as is linked here: Compensating Loudspeaker frequency response in an audio signal with a convolution equation of the form:
$$h[n] * w[n] = \delta[0]$$
Where $h[n]$ represents the non-causal channel response and $w[n]$ represents the non-causal equalizer coefficients. I use "non-causal" to mean each is centered about $n=0$ which simplifies the equation to remove any delay in the result, as in $\delta[0]$ vs $\delta[n-m]$. The result in turn for the causal solution simply changes the indexing back such that the first sample starts at $n=0$ rather than the middle sample.
In basic English this is saying knowing $h[n]$, solve for $w[n]$ that would result in an impulse response when the two systems are cascaded, since the output and input of a system whose impulse response is an impulse would be identical.
This would result in $y[n] = s[n]$ when the noise is not added, and since the noise is white, this would also be the solution that any least squares approach using actual signals $s[n]$ and $x[n]$ would converge to. As constructed I do not see this as a least squares problem but simply an inverse convolution problem. This is similar to the use of least-squares versus zero-forcing equalizers and has similar construction: the least-squares solution converges to the zero-forcing solution as SNR increases.
In this case we solve for $w[n]$ using inverse convolution.


Sub-Band Extraction Techniques for Complex (I/Q) Signal

I wanted to ask what are the techniques which can be used to extract this sub-band data?

That's a description of a filter bank.
Yes, the FFT can be used for such applications. You'll find that OFDM, which powers DVB-T, 4G/5G, WiFi, … (basically all high-speed wireless terrestrial links) does exactly that.
You'll also find that if you find the inherent sinc-shape of subcarriers in an FFT filterbank insufficient, but still need regularly spaced channels, polyphase filterbanks are the way to go.


What you're looking for is a color map. As its name implies, this process maps an input (single-component, or grayscale) value within a predetermined range to some other color based upon some mapping. The implementations I've seen typically use linear RGB interpolation between a number of control points specified throughout the allowed input range. The color map used in the image you showed is a pretty common one, starting with dark blue, then through cyan, green, yellow, orange, and ending with red. It is referred to as the jet colormap in MATLAB.

So, as an example, for an input intensity on the range $[0,1]$, you might have an RGB lookup table of:

0.0  -> (0, 0, 128)    (dark blue)
0.25 -> (0, 255, 0)    (green)
0.5  -> (255, 255, 0)  (yellow)
0.75 -> (255, 128, 0)  (orange)
1.0  -> (255, 0, 0)    (red)

Then, if you wanted to map the value $0.2$ to a color, then you would linearly interpolate between the first two datapoints as follows:

color = 1 / 5 * (0, 0, 128) / 0.25 + 4 / 5 * (0, 255, 0)
      = (0, 204, 26)

which would yield a greenish cyan, as it's close to the green control point.


Expected value of a Delta function times a function

If $h(t)$ is continuous at $a$, then $h(t)\delta(t-a) = h(a)\delta(t-a)$. Since the
$\operatorname{sinc}$ function is
continuous everywhere,
$$\begin{align}
Y(t) &= \sum_{n=-\infty}^{\infty} Z_n\delta(t-n\tau)h(t)\\
&= \sum_{n=-\infty}^{\infty} Z_n\delta(t-n\tau)h(n\tau)\\
&= \sum_{n=-\infty}^{\infty} Z_n\delta(t-n\tau)\operatorname{sinc}(n)\\
&= Z_0\delta(t)
\end{align}$$
since $\operatorname{sinc}(n) = \begin{cases}1,&n = 0,\\0,& n\neq 0.\end{cases}$
You (or your textbook's authors) probably meant to write something else.
You sure you are not missing a $\star$ indicating a convolution between
that $\delta(t-n\tau)$ and $h(t)$?


When calculating SNR, is noise included in signal?

The SNR is the ratio of $S$, the power in the signal $s(t)$, and $N$, the power in the noise $n(t)$. So, $$\text{SNR}_\text{dB} = 10\log_{10}\frac{S}{N}.$$

When using RMS values, we have $S = s_\text{RMS}^2/R$ and $N = n_\text{RMS}^2/R$. The resistances cancel out and then $$\text{SNR}_\text{dB} = 20\log_{10}\frac{s_\text{RMS}}{n_\text{RMS}}.$$


What information is lost, if any, in normalizing Seismic/Earthquake time-frames in an STFT to equal variance?

Although the specifics of your data will be the ultimate factor in how this sort of normalization will impact what you're trying to measure, I think it's useful to look at extreme cases to get a sense of what might go wrong.

For example, suppose the time samples of your STFT just happened to come out one-hot, with 1s in either the highest-frequency or lowest-frequency bins and 0s everywhere else :

STFT = [[1, 0, 0, 0, 0],
        [0, 0, 0, 0, 1],
        [1, 0, 0, 0, 0],
        ...
        ]

If you subtracted out the mean values from each time sample, then clearly you'd get some "bleeding" of the power information across frequencies. In particular, the middle-frequency bin would get some nonzero value, even though in your unnormalized STFT values, that bin never has any power.

What about just dividing each time sample by its standard deviation, like you're asking about ? In this case you'll just be changing the value of the nonzero frequency bin. This would certainly impact the mean value across time of each frequency bin, since all nonzero values in your spectrogram would now be 1+$\epsilon$ rather than just 1.

This is just one possible thought experiment that you could do to explore these effects. I'd say in general that the statistics of the frequency bins across time would definitely be affected by normalizing each time frame independently, but again that the specifics will depend heavily on the data that you're analyzing.


Is this an infinite response or a finite response?

A finite impulse response (FIR) filter has an impulse response of finite length. It can be implemented by a non-recursive structure, but there are also cases where the designer would choose a recursive structure for implementing an FIR filter (e.g., a recursively implemented moving average filter, as explained in this answer).

An infinite impulse response (IIR) filter has an impulse response of infinite extension. Consequently, it cannot be implemented by a non-recursive structure, but it requires a recursive structure.

With the above explanation in mind, I think it should be straightforward to decide which class the system in your question belongs to.


How to design digital filter from transfer function

Following OP comment requesting more details to Maximilian Matthé's answer:

Just like mentioned, the key to your question is the Bilinear Transform, where you substitute $s$ (or $j\omega$ in your case) with $\frac{2}{T}\frac{z-1}{z+1}$

This transform will have for effect to convert your transfer function from the continuous $s$ domain to the discrete $z$ domain.

I assume that you are fimiliar with the $z$ transform and that you know that a delay in the time domain is equivalent to multiplying by $z^{-1}$ in the z domain.  $Z\{f[x-n]\} = Z\{f[x]\}z^{-n}$

The implementation of a discrete filter (FIR or IIR) will happen in the time domain. Your software will implement something like

$$y[n] = b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]...+b_{N}x[n-N]-a_{1}y[n-1] - a_{2}y[n-2]...-a_{N}y[n-N]$$

In my previous equations :

$x$ is the input
$y$ is the output
$n$ is the index of the actual value in your buffer. $n-1$ being the previous value and so on.

If you bring that equation in the $z$ domain, you will have an equivalent like that.

$$Y = b_{0}X + b_{1}Xz^{-1}+b_{2}Xz^{-2}...+b_{N}Xz^{-N} -a_{1}Yz^{-1}-a_{2}Yz^{-2}...-a_{N}Yz^{-N}$$

After some algebra, you'll get the following form, which correspond to a transfer function (ratio of output on input)

$$\frac{Y}{X} = \frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}...+b_{N}z^{-N}}{1+a_{1}z^{-1}+a_{2}z^{-2}...+a_{N}z^{-N}}$$

This transfer function is what Maximilian Matthé answer refers to. In it's condensed form :

$$\frac{Y}{X}=\frac{\sum_{i=0}^{N}b_{i}z^{-i}}{1+\sum_{i=1}^{N}a_{i}z^{-i}}$$

Designing a discrete filter correspond to determining the coefficients of that transfer function. e.g. Finding the values of $a$ and $b$.

The value of your coefficients will change according to your sampling frequency, and you can confirm that considering that the Bilinear Transform introduce the sampling time $T$ in your equation.

You might want to consider the fact that the Bilinear transform is an approximation. If you look at Wikipedia's page, you will find this :

To compensate this approximation, you can substitute $\omega$ with

$$\omega_d = \frac{2}{T}\tan^{-1}(\frac{2\omega}{T})$$

This substitution is called frequency warping and will have a slight effect on the frequency response of your filter, but will most likely be negligible.

Finally, if you have access to matlab and the signal processing toolbox, you might want to have a look to bilinear() which will do all of this in a single command.


Finding the z-transform of $h[n] = a^n\cos(2\pi \frac{n}{F_s}f_0)$ for $n ≥ 0$ and zero for $n < 0$

The time domain signal (or impulse response)

$$h(n)=a^{n}\cos n\theta_0,\quad
\theta_0=2\pi\frac{f_0}{f_s},\; n\ge 0$$

is very common: it is a damped sinusoidal function (assuming $|a|<1$) which occurs frequently, because it is one possible response of a second order linear time-invariant system. So concerning your doubt, the cosine part is definitely meant to be part of the time domain signal. The poles and zeros can be found by rewriting $H(z)$:

$$\begin{align}H(z)&=\frac{1-az^{-1}\cos \theta_0}{1-2az^{-1}\cos\theta_0+a^2z^{-2}}\\&=
\frac{z(z-a\cos \theta_0)}{z^2-2az\cos\theta_0+a^2}\end{align}\tag{1}$$

From (1) it is easy to determine the zeros of $H(z)$:
$$z_{0,0}=0\quad z_{0,1}=a\cos\theta_0$$

For determing the poles, we write $H(z)$ as a partial fraction expansion:

$$H(z)=\frac{1}{2}\left [ \frac{1}{1-ae^{j\theta_0}z^{-1}}
+ \frac{1}{1-ae^{-j\theta_0}z^{-1}} \right ]\tag{2}$$

From (2) we see that the poles are given by
$$z_{\infty,0}=ae^{j\theta_0}\quad z_{\infty,1}=ae^{-j\theta_0}$$
We have complex conjugate poles because $h(n)$ is real-valued. Assuming that $h(n)$ is an impulse response, we can see from the poles that the system is stable if $|a|<1$ because then the poles are inside the unit circle of the complex plane.


Given a plot of both the magnitude $|H(\omega)|$ and its angle, How can you find the $H(\omega)$?

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

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

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

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

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

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


"Center" isn't the word that I'd use to describe what you want to do. Your signal does not have a zero-mean, so like @Dole suggested you should use a DC-Blocker (DC-Offset removing) filter to get the job done:

Here are your references:

https://ccrma.stanford.edu/~jos/filters/DC_Blocker.html
https://ccrma.stanford.edu/~jos/filters/DC_Blocker_Frequency_Response.html
https://ccrma.stanford.edu/~jos/filters/DC_Blocker_Software_Implementations.html

the code (from the above links) can be as simple as this, where x is the current sample and y is the filtered sample, while xm1 and ym1 are the x,y from the previous step (i.e. x[i-1], y[i-1]):

  y = x - xm1 + 0.995 * ym1;
  xm1 = x;
  ym1 = y;


When and how does one use a Wiener filter?

Perhaps a motivating radar/audio example would be an adaptive sidelobe canceller and adaptive noise cancellers. Rather than just show equations, let's walk through some text descriptions:

Let's say you are trying to record someone singing, but you're in a large auditorium. Your microphone is picking up the singer, but it's also picking up a lot of the reverberation. Now if we can characterize the reverberation somehow, we'll be able to remove it from our signal of the singer+reverberation. Enter another microphone (or a group of microphones), set up to just record the reverberation. With this setup, we can construct an adaptive filter that will remove the unwanted reverberation.

Another example would be some basic electronic protection for a radar system, where we want to remove jamming interference. In a radar system, you have a main channel which extracts the signal you're interested in. Sometimes, the case of some electromagnetic interference that originates from outside of the main antenna beam, i.e. from an antenna sidelobe, arises. This severely corrupts our main channel with noise, and may obscure target returns.

Intuitively, if we can somehow characterize just the interference, we can construct an adaptive filter that will remove the noise from the main channel. We can do this using single elements with isotopic (or nearly isotopic) responses. These auxiliary channels as they're called have very low gain compared to the entire antenna, and thus will receive only the interference and not a target which may be in the main beam.

The downside with wiener filters is that such a solution requires the noise to be stationary. Fortunately, we can use filter banks, time segmentation, and other signal processing techniques to make that possible.

Perhaps someone with a communication systems background could chime in and talk about adaptive equalizers, that's another application where you're trying to characterize some noisy channel, and you transmit a training signal, which the system knows about. The idea then is basically your first example: you know what you sent out, and you get some noisy version back. Using an adaptive filter, you can characterize that noise and now be able to remove it for anything else you send out (assuming the noise is perfectly stationary).

Hopefully that gives you some insight and motivation!


How do I filter out this horizon from my contour?

I have found a function provided by opencv2, threshold().

With this function, I can apply different combinations of threshold values so that this horizon can be filtered out.

The current values I am using is:

cv2.threshold(img, 100, 255, cv2.THRESH_BINARY)


Help implementing audio dynamic range compression

Here are some suggestions:

There are plenty of opensource implementations (Sox, Audacity, etc). Even if you don't understand them, you might be able to translate the code from C to javascript.
I'm not aware of a good explanation of the process online, but there are plenty of books on the subject:

Digital audio signal processing covers this topic and is well written. (As does DAFX, but DAFX is poorly organized and the coverage is less straightforward)
Digital Audio with Java also covers this topic and comes with working Java code that should be easy to translate to other languages, such as javascript. This book has many many flaws, but it is good for someone with no audio programming experience.

The principle is to create an envelope of the signal (controlled by attack and release), shape that envelope using using some transfer function (controlled by ratio, threshold and knee), and then apply that result back to the original signal. A makeup gain stage often follows.

@Deve's answer suggests some possible transfer functions.


Relation between FIR filter's transition width and phase delay

In this answer I gave two well-known heuristic design formulas relating the filter order of (linear-phase) FIR low-pass filters to the transition bandwidth. According to both formulas, the transition bandwidth is inversely proportional to the filter order. Since the phase delay of a linear-phase FIR filter equals $M/2$, where $M$ is the filter order, the phase delay $\tau_p$ and the transition bandwidth $\Delta\omega$ are also inversely proportional to each other:

$$\Delta\omega\propto \frac{1}{\tau_p}$$


How does the earphone cable in a mobile phone pick up the FM reception?

Your question is based on a misunderstanding. The RF signal is not directly sampled by the sound card.

Instead, the earphone cable just doubles as an antenna for a dedicated FM receiver; the soundcard hasn't got anything to do with FM reception.
In fact, it doesn't even "get" radio frequencies – they're filtered out before and fed to the FM receive chain.


Understanding MATLAB/Mathematica's bilinear interpolation

imresize gets you different results on Matlab and Octave. The one in Octave calls the interp2 function, which with 'bilinear' method will return you [100 75 50 25 0].

However, it seems like there is some extra processing in Matlab.  It performs affine transformation on the image. The interpolation is calculated on each pixel in the output image, then inverse map onto the indices of the pixels in the input image with the contributing weights. So I guess the mapped coordinates are not distributed uniformly (not exactly at 1.5 and 2.5 between [1,2], and [2,3] in your case).

Probably you can edit imresize_old to see more source code content than imresize.


Convolution is typically used in the context of filtering. You have an image, $I$, a filter kernel, $K$, and you convolve them together to get a filtered image, $J$:

$J = I \star K$

where $\star$ denotes convolution.

The nature of the filtering operation will be determined by the coefficients of the filter kernel, $K$.


How to estimate octave and size for visual features positioned at Harris corners

I'm not sure there is actually a good response to your precise question : the scale-space thing of SIFT and SURF was actually developed to automatically estimate the "good" relevant neighborhood size around a corner-like keypoint (which is what good features to track are).

Now, more positive answers would be:

build a database of keypoints and good matches (e.g. using square calibration patterns) and create an automated performance assessment on this database to find the correct size. This task can actually really be automated (see the work by Mikolajczyk and Schmid about point matching evaluation)
embed your good features in an image pyramid to also have some kind of scale associated with them. You can look for references on multi-scale FAST and Harris interest points, that do something very similar to this procedure.

To heuristically find the maximum block size, you can compute estimates of your images with a box blur (which is more or less what the blockSize operator does) and see when the corner disappears. Note however that more blur takes the corner away from its true location.

If you're really looking for some quick-and-dirty fix, try sizes between 5x5 and 11x11 (typical sizes used in stereo block matching). If you're looking for an intellectually satisfying criterion, then try to maximize the likelihood of good matching of two feature points under your noise level.


Given f(x) and its FFT F(u), I need to prove that df/dx = F(u)*2iu(pi)/n

This should have been a comment on the question but is too long to fit within
the 500-character limit.

The $f$ that is the subject of this question is a sequence or vector of $n$ numbers
$$[f_0, f_1, \ldots, f_{n-1}] ~~ \text{or} ~~ \left[f[0], f[1], \ldots, f[n-1]\right]
~~ \text{or} ~~ \left[f(0), f(1), \ldots, f(n-1)\right]$$
depending on choice of notation, and its Discrete Fourier Transform (DFT) $F$ (called the
FFT by the OP) is another sequence or vector of $n$ numbers. Thus,
it is not clear what the derivative $\frac{\mathrm d}{\mathrm dx}f(x) = f^\prime(x)$ means
since the argument of $f(\cdot)$ is an integer and not a continuous variable.  One
could attempt to make sense of what is intended by looking at the OP's approach to the
problem via (cyclic or periodic) convolution with the sequence $[1, -1]$ (actually should have been $[1,-1,0,0,\ldots, 0]$ which results in replacing $f_i$ by $f_{i+1}-f_i$.
Now, we are getting closer to the notion of differentiation since the derivative of
a continuous function is obtained by looking at the limit of $\frac{f(x+h)-f(x)}{h}$ as
$h \to 0$. Here we have the closest we can get since we are
using $\frac{f(i+1)-f(i)}{1} = f_{i+1}-f_i$, and of course $1$ cannot "approach" $0$.
Unfortunately, this does not give the desired result.  The $k$-th term of the
DFT of $[1,-1,0,0,\ldots, 0]$
is readily computed from the definition (no FFTs needed)
as $1 - e^{-j2\pi k/n}$, and so the DFT of the "derivative" sequence
$f_{i+1}-f_i$ has $k$-th term $F_k(1 - e^{-j2\pi k/n})$ which is not
what the OP wants to prove.

So, what does the OP want to prove? It is
$\frac{\mathrm df}{\mathrm dx} \mathbf{=} F(u)*j(2\pi u/n)$.
This needs clarification on two grounds. First, the left side
of this equality is a function of $x$, the right side a function
of $u$. So what does it mean to say the two are equal?
Second, what the right side means needs clarification. Does
the $*$ mean multiplication of the two complex numbers $F(u)$ and
$j(2\pi u/n)$, or is the $*$ denoting complex conjugation of $F(u)$
before the multiplication, or does $*$ denote convolution so
that $F(u)*j(2\pi u/n)$ means the $u$-th term of the (cyclic or periodic)
convolution of the sequences $[F_0,F_1,\ldots,F_{n-1}]$ and
$[1, j(2\pi/n),j(2\pi 2/n), \ldots, j(2\pi (n-1)/n)]$?
The closest match occurs if we assume that $=$ is a typographical
error and that what is meant is not the equality that the OP
claims is needed. Instead, suppose that what was really asked
for but mistranslated by the OP is that the $k$-th term of
the DFT of $\frac{\mathrm df}{\mathrm dx}$
by which is meant the sequence with $i$-th term $f_{i+1}-f_i$
is $F_ke^{j2\pi k/n}$ (another typo, the missing exponential!). As
noted above, this is not quite right: it should be
$F_k(1 - e^{-j2\pi k/n})$. Note, incidentally, that the inverse
DFT of the sequence with $k$-th term $F_ke^{j2\pi k/n}$ is a
cyclic shift of the sequence $f$; see, for example,
this
Wikipedia entry

But maybe as a non-expert in image processing I am missing something completely
and so I await clarification from the OP, or perhaps a more detailed
explanation from @JasonR of the method outlined in his comment on the OP's
question.


Order of operations when creating a spectrogram

To get a spectrogram using existing libraries, it is usually easier to do it for a whole auduio clip. For example in Python with librosa, it is a single call (assuming the audio is already loaded): librosa.stft.
A streaming implementation requires a bit more care. You will then need tohandling the window function, the FFT operation on each window, and the overlapped extraction of the windows. Not very complicated, but a bit harder.
The primary benefit would be lower-latency and lower memory requirements. If these are not key requirements, I would recommend the clip-wise solution.


Stationarity, discrete-translation operator, and the power spectral density matrix

Given the definition of the correlation matrix $\mathbf{R}_{\mathbf{x}}$ here, I am assuming that $\mathsf{E}[\mathbf{x}] = \mathbf{0}$. I do this because the correlation matrix is usually defined as $\mathsf{E}[(\mathbf{x} - \mathsf{E}[\mathbf{x}])(\mathbf{x} - \mathsf{E}[\mathbf{x}])^{\dagger}]$, where $\dagger$ indicates complex conjugate tranpose.

Note that since $\mathbf{R}_{\mathbf{x}}$ is a correlation matrix, it is Hermitian positive semi-definite, which means that it is unitarily diagonalizable and has all non-negative eigenvalues:
\begin{equation}
\mathbf{R}_{\mathbf{x}} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^{\dagger},
\end{equation}
where $\mathbf{U}$ is a unitary matrix and $\mathbf{\Lambda}$ is a diagonal matrix with non-negative numbers on its diagonal.

Let's consider the translation operator $\mathbf{T}$:
\begin{equation}
\mathbf{T} = \mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger},
\end{equation}
where $\mathbf{\Phi}^{\dagger}$ is the DFT matrix and $\mathbf{P}$ is diagonal. As I mentioned in a comment, $\mathbf{T}$ should be unitary, and that requires that the diagonal matrix $\mathbf{P}$ have complex exponentials on its diagonal. In fact, to make $\mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger}$ a translation operator, $\mathbf{P}$ should be
\begin{equation}
\mathbf{P} = \textrm{diag}(1,e^{2\pi i/M},\ldots,e^{2\pi i(M-1)/M}).
\end{equation}

I do not think this will have an impact on the proof, but it is important for applications.

Now let's assume that $\mathbf{R}_{\mathbf{T}\mathbf{x}} = \mathbf{R}_{\mathbf{x}}$:
\begin{eqnarray}
\mathsf{E}[(\mathbf{T}\mathbf{x})(\mathbf{T}\mathbf{x})^{\dagger}] &=& \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathsf{E}[\mathbf{x}\mathbf{x}^{\dagger}]\mathbf{T}^{\dagger} &=& \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathbf{R}_{\mathbf{x}}\mathbf{T}^{\dagger} &=& \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathbf{R}_{\mathbf{x}} &=& \mathbf{R}_{\mathbf{x}}\mathbf{T}
\end{eqnarray}
Since

$\mathbf{R}_{\mathbf{x}}$ and $\mathbf{T}$ are both diagonalizable
and they commute with each other,

they are simultaneously diagonalizable. This means that there is a single matrix that diagonalizes both.

We have been told from the beginning that $\mathbf{\Phi}$ diagonalizes $\mathbf{T}$, so now we know that $\mathbf{\Phi}$ diagonalizes $\mathbf{R}_{\mathbf{x}}$, too. This means that the unitary matrix $\mathbf{U}$ that diagonalizes $\mathbf{R}_{\mathbf{x}}$ must be $\mathbf{\Phi}$:
\begin{equation}
\mathbf{R}_{\mathbf{x}} = \mathbf{\Phi}\mathbf{\Lambda}\mathbf{\Phi}^{\dagger}.
\end{equation}
We have already established that $\mathbf{\Lambda}$ is diagonal with non-negative real numbers on its diagonal. The $\mathbf{S}_{\mathbf{x}}$ that we have sought, is $\mathbf{\Lambda}$.


Are there any pictures of a sample accurate spectrogram on the web or in books?

Producing accurate time-frequency analysis is an area of research. Some pointers are : Wigner-Ville transformation, re-allocation methods, wavelets.


In a PDM microphone, the sigma-delta convertor pushes the noise to frequency regions above the 0-20kHz spectrum. So if you'd take the FFT of the signal that comes directly from the PDM microphone, the straight zeros and ones, and only look at the bins that fall in the 0-20kHz region, you'd get the data that you need.
You can see this here, where a 16kHz sine wave has been converted to a PDM signal:

The noise floor on the left of the green 20kHz line is well below -70dB, but quickly increases to the right of it.
In the plot above, the PDM data rate is 768kHz. A common PDM microphone data rate will be more on the order of 3MHz. The problem with this is that you'll doing FFTs of a large amount of data and that you're throwing away a lot of unused data.
A moving average filter would indeed quickly reduce the number of data samples that go into the FFT. It will result in attenuation in the pass band that is too high for good audio quality, but it should be totally fine for an LED display.
I wrote about my own PDM microphone learnings in a series of blog posts.


Filter design with 8-tap latency on recursion

This is tricky. Any filter with 8 tab latency between coefficients can be represented in the z-domain as a rational function in $z^{-8}$ . That's basically the same as saying the impulse response has 7 zeros between each non zero tap. Any filter with this property is periodic in the frequency domain. You can pick any transfer function from -pi/8 to pi/8 but the regions from pi/8 to 3*pi/8,etc. will just be repetitions of that.

You may be able to do better with varying delays but this is a complicated design problem


Power Spectral Density of Finite Signals

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.


Spectral leakage from mathematical point of view

The reference is simple convolution.  So-called spectral leakage is due to multiplication with a window function.  Feed any longer signal to a finite length DFT or FFT, and a default rectangular window (if not some other shape of window) is usually applied to the longer signal.  That multiplication in the one domain results in (circular) convolution of the transforms in the other domain.  The results of this convolution are sometimes called "leakage".


Output of lowpass filter with damped sine wave input

I've finally figured this out thanks to a friend of mine. The thing is that $N(t)$ can be even or odd, and in both cases, the signal we get is easy to deal with.

We can calculate the desired expectation this way:

$$\mathbb{E}[Z(t)]=\mathbb{E}[h(t)*Y(t)]=\mathbb{E}[h(t)*Y(t)|N(t) =2n]\mathbb{P}[N(t)=2n]+\mathbb{E}[h(t)*Y(t)|N(t) =2n+1]\mathbb{P}[N(t)=2n+1]$$

But when $N(t)$ is even or odd, $Y(t)$ becomes simple:

$$Y(t) =
\left\{
 \begin{array}{ll}
  \cos^2(\omega_0t)  & \mbox{if } N(t)\mbox{ is even } \\
  -\cos(\omega_0t)\sin(\omega_0t)    & \mbox{if } N(t)\mbox{ is odd }
 \end{array}
\right.$$

So

$$\mathbb{E}[Z(t)]=\mathbb{P}[N(t)=2n]\cdot h(t)*\cos^2(\omega_0t)-\mathbb{P}[N(t)=2n+1]\cdot h(t)*(\cos(\omega_0t)\sin(\omega_0t))$$

And those two probabilities are easy to calculate, as well as the two convolutions left.


Estimating Damping Factor (Q) from Noisy Measurements

Our system has the impulse response (why we changed to imaginary instead of real? :) ):
$$h(t)=\mathbb{Im}(Ae(iwt)e(-bt))=Ae^{-bt}sin(wt)$$

With the following structure (ref.) (discretized):
$$H(s)=\frac{Kw}{(s+b)^2+w^2}$$

Under matlab, you only need to use procest:

procest(data,'P2U');

Which models the same structure under the form:
$$H(s)=\frac{K}{(1+2\zeta T_ps+T_p^2s^2)}$$


Confused over Stopband attenuation and stopband ripple

I understand that the stopband ripple can be given by the following:

  Stopband ripple = −20log10(δs) 40 = −20log10(δs) δs=0.01

Where does that come from? I have never heard the term "stopband ripple".

Stopband attenuation of 40 dB means that in the stop band any level of -40dB OR LESS is allowed. How much less than -40dB or where it hits -40dB isn't a controlled variable. -50dB, -60dB etc all are fine.


w = conv(u,v) returns the convolution of vectors u and v. If u and v
are vectors of polynomial coefficients, convolving them is equivalent
to multiplying the two polynomials. (from MATLAB documentation)

Both expressions are correct. In 2., whenever an addend of a convolution operation calls randn in a convolution computation procedure, the callee returns new values for both real and imaginary parts of h.


I guess image denoising and image super-resolution techniques could be used to enhance license plate images.
Here are a few relevant links:

Image denoising based - SNIDER: Single Noisy Image Denoising and Rectification for Improving License Plate Recognition
Image super-resolution based - Multiframe Superresolution of Vehicle License Plates Based on Distribution Estimation Approach


Berchin's FDLS arbitrary filter design algorithm

FDLS requires a causal frequency response. Your prototype frequency response has zero phase everywhere, which is most definitely not causal.

An IIR filter order of 50 is humongous. When FDLS has too many poles and zeroes available, it "tries" to cancel excess poles with excess zeroes. Unfortunately, due to numerical limitations, the cancellation is often pretty poor. For your notch filter, try 2nd or 4th order numerator and denominator, after making the response causal.

Though (number of measurements) = (number of coefficients) = [(numerator order) + (denominator order) + 1] is mathematically sufficient, I find that FDLS likes many more measurements than that. I recommend an absolute minimum of twice as many measurements as coefficients, and prefer 5x to 10x.

Greg Berchin


Filtering a signal to have a constant exponential moving average

I'll draw on my (limited) experience in audio mixing and mastering. Speech recorded through a microphone requires the same statistical properties (constant moving average and constant variance), and two straightforward processing stages ensure this.

An "exponential moving average" is another word for a 1-pole low-pass filter. This moving average can be brought close to zero by passing the signal through any high-pass filter, such as a "DC killer" or "remove subsonic rumble" filter that an audio processing package such as Audacity or Audition may provide. The stronger the high-pass characteristic, the closer the moving average will stick to the center line.

Second, you want to make the variance constant. That's a job for a level compressor. Find the RMS (square root of variance) of the signal over the next few milliseconds, and adjust the gain (amplification) of the signal over that period to bring the variance up or down to the variance you want. Audio processing packages provide this as well.

You mention numerical stability. The standard DSP literature describes ways to implement a robust high-pass filter. For instance, FIR filters tend to be very stable.


I would suggest you first demodulate your signal. Take the received signal y[t] = x[t]*cos(2*pi*1000*t) and multiply it by cos(2*pi*1000*t). Remember, $ cos(A)cos(B) = 1/2 (cos(A-B) + cos(A+B)) $. You'll be left with x[t]*0.5*(cos(0) + cos(2*pi*2000*t) = x[t]*0.5*1 + x[t]*0.5*cos(2*pi*2000*t). Filtering the remaining signal with a low pass filter at 1k Hz, leaving only x[t]. From there you can try and answer the rest of your questions.

In Matlab this will look like:

t = 0:1/Fs:END_TIME;  % Fs is your sampling rate. This makes a vector of individual time ticks
y_demod = y.*cos(2*pi*1000*t);   % This multiplies every element in Y with the cos(...)
                                 % ... element. Make sure they're the same size.
a = (1/Fs)*(Fs/4);               % A cutoff frequency of approximately 1/2 the sampling rate
x_est = filter([1 1-a],a,y_demod);

The simple 1-pole low pass filter was taken from here: https://stackoverflow.com/questions/1783633/matlab-apply-a-low-pass-or-high-pass-filter-to-an-array


Vector Modulators and Single Sideband Mixer (frequency translator) are the common and relatively simple analog component for rejecting the image signal.

The top block diagram in this post is the architecture of a vector modulator and single sideband mixer and explains in further detail how it provides image rejection using complex signalling.

Frequency shifting of a quadrature mixed signal

Normally we see these components with a real input where an internal 90° splitter immediately converts the input to a complex signal and then it is further processed. However there is no reason these can't also be used in identical fashion with the quadrature splitter bypassed in the case of a low IF IQ signal already exists.

This would be the favored approach for eliminating the image signal.

However to implement a generic complex filter in the analog domain, you would need  to do the following with four real filters as shown in the diagram below. This is related to performing full complex multiplication with real signals as given by:

$$(s_I+js_Q)(c_I-jc_Q)$$
$$= (s_Ic_I+s_Qc_Q) + j(s_Ic_Q-s_Qc_I)$$

Corresponding to the multiplication of a signal (s) with a coefficient such as in an FIR filter (c) showing that we have an I and Q input, multiplied by and I an Q coefficient, resulting in a I (real part given by $s_Ic_I+s_Qc_Q$) and Q (Imaginary part given by s_Ic_Q-s_Qc_I).

Equivalently complex analog filtering (the filtering of a complex input with a complex impulse response) is formed with four real filters as shown in the graphic below. Each filter has a real impulse response given by $h_i(t)$ or $h_q(t)$ corresponding the real and imaginary portion of the complex impulse response desired.

Also to add, I did come across this: https://ieeexplore.ieee.org/abstract/document/4151912

Note that SAW filters are an analog equivalent of an FIR filter ("all=zero" filter, with poles at infinity, formed by the summation of weighted delays)


How does a shift in time domain result in phase shift in frequency spectrum?

Consider first that in the time domain that a phase rotation or shift in phase or phase shift specifically is done by multiplying the waveform with $e^{j\phi}$ which rotates or shifts the complex phase for each sample in time by $\phi$. Realize that $e^{j\phi}$ is a complex phasor of magnitude $1$ and angle $\phi$. In contrast a time delay is given as $x(t-t_o)$ as the OP has indicated.

This is no different in the frequency domain where we have a complex waveform just the same but with a different independent variable (frequency instead of time). Similarly we can multiply the frequency domain waveform by the $e^{jt_o \omega}$ which for a given sample at frequency $\omega$ will cause a constant phase rotation of that sample.

In the frequency domain $e^{jt_o \omega} = cos(t_o \omega) + jsin(t_o \omega)$ is a constant phase shift of $t_o \omega$ for each sample given by $\omega$ and is not a scaling of frequency as the OP suspects since we are in the frequency domain (in the time domain a scaling of frequency would be of similar form given as $y(t) = cos((\omega t_o) t)$ but that is not what is occuring here since we have $Y(j\omega) = X(j\omega)e^{-j\omega t_o}$.

Consider a very simple case of $X(j\omega_1) = Ae^{j\phi_1}$, which means that the magnitude of frequency $\omega_1$ is $A$ and the phase of this frequency is $\phi_1$.  After the time delay given by $x(t-t_o)$ the result would be:

$$Ae^{j\phi_1} e^{-j\omega_1 t_o} = Ae^{j(\phi_1-\omega_1 t_o)}$$

Where $\omega_1 t_o$ represents a constant phase shift or rotation.


Implications of $X( -j\omega ) = X^*(j\omega)$

The fact is that if $x(t)$ is real and then $X(-j\omega) = X^*(j\omega)$ (it's easy to prove) and, because of duality, if $X(j\omega)$ is real, so also then $x(-t) = x^*(t)$.  in your case, you have both $x(t)$ and $X(j\omega)$ being real.  keep in mind that anything that is purely real is also equal to its complex conjugate.

so you can say:

"if $x(t)$ is real and even (that is $x(-t) = x(t)$) so also is $X(j\omega)$ real and even ($X(-j\omega) = X(j\omega)$)."


Harmonic Product Spectrum limitations in pitch detection

It could be that too few harmonic partials are present in the signal at these higher pitches. The HPS algorithm is pretty simple and relies on those upper harmonics to keep stacking up until the fundamental emerges from the background. Of course, we should wonder, what's your sampling rate? If it's 8000 hz, then there's only room for 3 harmonics of a 1000 hz pitch...


If you add two cosines you simply get the product of the sum and difference frequencies.
$$cos(x) + cos(y) = 2\cdot cos \left( \frac{x+y}{2} \right) \cdot cos \left( \frac{x-y}{2} \right)$$
If the frequencies are very close together, than the difference frequency is close to zero and that's what creates the "beat".
The exact definition of what constitutes a beat and what doesn't depends on your application. For example: in audio I would call any frequency difference of less than maybe 4 Hz or so a beat since it creates the typical beat sound.


How can one know the polynome equation based on the output of a system?

The right-hand polynomial is the Z-transform $X(z)$ of the discrete-time function $x[n]$, and is given by

$$X(z)=\sum_{n=0}^{\infty}x[n]z^{-n}$$

(unilateral or bilateral does not matter here, since $x[n]$ is causal anyway). So, if you now insert $x[n]$ you get

$$X(z)=\sum_{n=0}^{\infty}z^{-n}=\sum_{n=0}^\infty\left(z^{-1}\right)^n$$

You can use formula for the geometric series, as long as $|z|>1$ (i.e. $|z|>1$ is the region of convergence of the Z-transform):

$$X(z)=\sum_{n=0}^{\infty}\left(z^{-1}\right)^n=\frac{1}{1-z^{-1}}=\frac{z}{z-1}.$$

For the second case, i.e. $x[n]=a^nu[n]$, you can write the following:

$$X(z)=\sum_{n=0}^{\infty}a^nz^{-n}=\sum_{n=0}^{\infty}(az^{-1})^n$$

Then, again applying the geometric series formula you get (given that the series converges, i.e. $|az^{-1}|<1$)
$$X(z)=\frac{1}{1-az^{-1}}$$


Hi Nouali Ibrahim Yassine,

Without any additional information about the two signals, what you are trying to do is a blind source separation. This means that you will have to use a set of advanced statistical tools ( principal components analysis or independent component analysis for example, I would go more for the Independant component analysis) to try to figure out different signal patterns.

If you can have additional information about your signals (I guess you know your sensors, their number, their signal characteristics such as spectral distribution or some pre recorded patterns that you can use to correlate with the mixed signal, their inter-correlation) You may find a simpler solution to your problem.

In any case, here is some lecture about the methods that I stated. There may be other ones that are better suited for your problem.

Blind source seperation in general :
https://perso.univ-rennes1.fr/laurent.albera/alberasiteweb/bss.html
Principal component analysis :
http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/
Principal component analysis versus Independant component analysis :
https://www.quora.com/What-is-the-difference-between-PCA-and-ICA

Best regards,

Mourad


How we do frequency shaping and amplitude shaping in hearing aid

Merium:

To develop a better understanding of how frequency shaping is used in hearing aid processing, consider the following paper, "Physiological assessment of contrast-enhancing frequency shaping and multiband compression in hearing aids", Ian C. Bruce, Physiological Measurement, Volume 25, pages 945-956, 2004.  It is an older paper, but it provides a good overview of some ideas you may want to consider using in your project.  If you are a registered student at a college or university, you should have no trouble finding this article using the library's resources.

Due to the fact that common hearing impairment is associated with loss of high frequency sensitivity, I would suggest implementing a simple filter that provides gain for high frequencies and maintains the gain at lower frequencies.  You could do this using a multiband approach, which would give you more flexibility in setting the gain for each band, or a across single band that might be more readily realized in the DSKC6713 framework.  Have a look at this article to understand the different types of hearing loss.

As far as implementation goes, I would suggest using Matlab to develop and test your processing scheme prior to implementing it in the DSP C++ framework.

I hope this helps.

Michael.


Fourier Transform of a Pulse Train (Frequency Comb Simulation)

Let us first see the theory.  By pulse train I guess you mean a periodic function from minus infinity to plus infinity.  Your pulse train will have the Fourier series(FS) in stead of Fourier transform(FT).  Fortunately the FT of a single pulse and the FS of the pulse train are related.  Let the pulse train be periodic with P Hz.  Then the FS of the pulse train will be sampled version of FT with a sampling rate of P Hz.  i.e. except for a scaling factor, the FS values are the samples of the FT sampled at P Hz. And the single pulse FT is a sinc function.  Thus pulse train will have FS coefficients which are samples of the sinc function.


Two-GSMK channels demodulation: 3dB difference in BER

Yes that is correct, you are only getting half the signal: with cosine on both inputs half the signal moves to baseband and the other half moves to twice the frequency offset. I suspect the reason it is not working in the sine cosine case is the sign of that is such that you translated the entire signal upward to twice the frequency offset instead of to baseband. If that is the case you could invert the sine input and it should work without the 3 dB loss.


difference between mean filter and order statistic filter

Quick Answer

The main difference of this filters is how it perform the operations.

Mean Filter

Brief Description

Mean filtering is a spatial filter, and it's a simple, intuitive and easy to implement method of smoothing images, i.e. reducing the amount of intensity variation between one pixel and the next.

How It Works

The idea of mean filtering is simply to replace each pixel value in an image with the mean (`average') value of its neighbors, including itself.
Mean filtering is usually thought of as a convolution filter and it's based around a kernel, which represents the shape and size of the neighborhood to be sampled when calculating the mean. Often a 3×3 square kernel is used like this
 $$
   K =
  \left[ {\begin{array}{cc}
                \frac{1}{9} & \frac{1}{9} & \frac{1}{9}\\
                \frac{1}{9} & \frac{1}{9} & \frac{1}{9}\\
                \frac{1}{9} & \frac{1}{9} & \frac{1}{9}\\
              \end{array} } \right]
$$

Note: Check this link for more details in mean filtering.

Order Statistics Filter

Brief Description

This type of filter is based on estimators and  is based on "order", the sense of order is about some quatities like $\operatorname{min} $ (first order statistic), $\operatorname{max}$ (largest order statistic) and etc...

How It Works

Given $N$ observations $ X_{1}, X_{2}, X_{3}, \dots X_{N} $ of a random variable $X$, the order statistics are obtained by sorting the $\{X_{i}\}$ in ascending order. This produces $\{X_{(i)}\}$ satisfying:

$$X_{(1)} \leq X_{(2)} \leq X_{(3)} \dots \leq X_{(N)}$$

where $\{X_{i}\}$ are the order statistics of the N observations.
So, an Order Statistic Filter (OSF) is a estimator $ F(X_{1}, X_{2}, X_{3}, \dots X_{N})$.

Some common filters which fit the order statistic filter framework are:

The linear average, which has coefficients:
  $$ \alpha_{i} = \frac{1}{N} $$
The median filter, which has coefficients:
$$ \alpha_{i} = \left\{ \begin{array}{ll}
     1 & i = (N+1)/2\\
     0 & \text{otherwise} \end{array} \right.$$
The trimmed mean filter, which has coefficients:
$$ \alpha_{i} = \left\{ \begin{array}{ll}
     1/M & (N - M + 1)/2 \leq i \leq (N + M + 1)/2 \\
     0 & \text{otherwise} \end{array} \right.$$

Please check the references for more details in the Todd Veldhuizen web page.


Length of DFT defines whether I see the harmonic or the dominant frequency

Zero padding does not increase the frequency resolution— it just interpolates more samples in frequency but does not give you any more information. The frequency resolution in Hz of an non-windowed (rectangular window) Fourier Transform is 1/T where T is the length of your time domain waveform in seconds. You must have a longer data sample in order to get higher resolution or see lower and lower frequencies. An extreme case is that you need an infinite amount of time in order to observe a true DC signal.

Note that zero padding gives you more samples of the Discrete Time Fourier Transform (The DTFT is the result if the time axis was allowed to go from minus infinity to plus infinity; by zero padding you approximate that). The Discrete Fourier Transform in contrast has a time axis this only goes from 0 to N-1.


What does $I_{seg}$ in 3GPP TS 38.212 version 15.2.0 Release 15 represents?

$I_{\mbox{seg}} = 1$.

if payload $A>1013$ or $A\ge 360$ and $G>1077$. Divided two block segment
PUCCH $n_{\mbox{max}}=10$;
PDCCH $n_{\mbox{max}}=9$;


Channel encoding after digital modulation

The role of the modulator in your block diagram is to make an analog signal, while in classical communications, the data and the channel coder is discrete. So no, in the sense of your block diagram, that can't work.

Note that you stated trouble because your signal wasn't binary in your approach: Not all channel coders work on binary symbols. For example, Reed-Solomon Codes are used exactly because they work on symbols from a larger field than the binary numbers. (That brings advantages when dealing with burst errors.)

Non-binary LDPC codes exist, and there's other codes, too, that don't use binary numbers. Generally, very little in our modern communications is "forced" to be binary, it just happens so that $\mathbb F^2$ is a very easy-to-deal-with finite field.
From the top of my head, I can't think of a single channel coder type which couldn't exist as a base-3 instead of base-2 algorithm, for example. (Building a good and performant channel decoder might be hard, though.)

TL;DR:

Your modulator needs to stay between your digital processing, which includes the channel coding, and the analog channel, but the notion that "digital == binary" is plain wrong.

Wth the very widely used Reed-Solomon codes you have an excellent counterexample to that, and you might want to look into these codes.


increase of encoding rate increases bandwidth?

As per @JasonR's suggestion (but no longer succinct...., though
hopefully still correct...)

Error-control coding creates redundant symbols that are also transmitted
over the channel. A simple description of what happens is as follows.
The data symbols are divided into chunks of $k$ (data)
symbols and to each such chunk is appended a chunk of $n-k$
redundant symbols to create a codeword of $n$ symbols. The
rate of the code is $R = \frac{k}{n}$ since $n$ symbols must
be transmitted over the channel in order to get $k$ data symbols
across more reliably. The device that
creates and appends the redundant symbols is called an encoder.
Note that the encoder accepts $k$ symbols as input and spits
out $n > k$ symbols as output.  Now, it is possible for the encoder
to send the $n$ symbols to the transmitter for transmission using the
system and signaling design, modulation method that is already
set up. In this case, the channel signaling rate is not changed
at all and so the bandwidth remains the same. However,
the codeword, consisting of $n$ symbols, takes longer
to transmit than the $k$ data symbols would have using the
same transmitter and channel that we have available.
There is also the issue of what happens to the next chunk of
$k$ data symbols that are knocking at the encoder input while
the encoder is still busy sending the redundant symbols to
the transmitter.  So this scheme can be implemented if
we don't really care about how long it takes to transmit
the data (which is presumed to be a finite-length file or
audio track (say) and not a (potentially) infinite never-ending
stream of data all of which needs to be delivered to the
receiver, and are willing to have lots of buffers
between the encoder and the transmitter to save the codewords
being generated until the transmitter gets around to transmitting
them.  We will also need buffers at the receiver where the
redundant symbols will be stripped off, and the data
symbols delivered to the destination in chunks of $k$
symbols followed by $n-k$ intervals of dead silence: buffering
can be used to save the $k$ symbol chunks and patch them
into a continuous stream of data.

The alternative, which most people opt for, is to modify
the transmitter and receiver so that each channel symbol
takes less time (by a factor of $R$) to be transmitted,
so that the $n$ symbols of the codeword can be transmitted
in the same time that the $k$ data symbols would have been in
absence of error-control coding. The increase in the
signaling rate means that the transmitter output requires a
larger bandwidth (by a factor of $R^{-1}$) for transmission.
The channel thus must have larger bandwidth too, which can
create problems if wireless channels are being used.
Frequency bands are allotted by government agencies and if
your transmitter is already using the full spectrum allotted
to you, a Plan B and a complete re-design might be needed.
These systems with larger bandwidth also need buffers
at the transmitting end as well as the receiving end,
but only two relatively small ones of $k$ symbols at both
ends. The two buffers precede the encoder and are used
in ping-pong fashion. The data source is reading in
$k$ symbols into one buffer, say one symbol every $T$
seconds, and thus taking a total of $kT$ seconds to fill
the buffer. During this time, the encoder reads out
the $k$ data symbols from the other previously filled buffer
at a rate of one symbol every $RT$ seconds for a total
of $k(RT)$ seconds. The data symbols are also sent on to
the transmitter while the $n-k$ redundant symbols are created
within the encoder, which then spends the remaining
$kT-k(RT) = (n-k)(RT)$ seconds sending the $n-k$ redundant
symbols to the transmitter at a rate of one every $RT$
seconds.  At the end of a $kT$ second
period, the data source starts to fill the buffer just
emptied by the encoder, while the encoder starts to
accept data from the buffer just filled by the data
source. If you are still following all this, note
that the encoder and the transmitter are working
with a clock that is running at a higher rate
(by a factor of $R^{-1}$) than
the clock of the data source, while the clock
signals to the buffers are alternately at the
low rate (when connected to data source) and
at high rate (when connected to encoder).
Similar things occur in reverse between the
receiver output and the data destination
but I will skip the details.

A third alternative, as pointed out by Bryan,
is to change the modulation format
on the channel so that, for example, instead of transmitting one channel bit
with each channel symbol (as in BPSK), you transmit more channel bits per
symbol (e.g. use QPSK to send two bits per channel symbol). In this case,
there need not be any increase in bandwidth (which is primarily determined
by the symbol rate (rate in baud)). A common use of the BPSK/QPSK usage
is the use of rate-$\frac{1}{2}$ convolutional coding which would double
the bandwidth if BPSK were used, but by switching to QPSK, the bandwidth
remains the same even when coding is used. However, as also noted by Bryan,
the transmitter and receiver become more complicated.


For the Hanning window, the zero crossings occur at a much wider bandwidth (trading for lower side-lobe height) compared to a rectangular window (sharper main-lobe width but higher side-lobe). You need to be careful with the requirements regarding the choice of window given the trade off between side-lobe height and main-lobe width.


Adaptive thresholding based on number of pixels

The book Digital Image Processing by Gonzalez and Woods provides a good introduction into standard techniques for automatic binarization. The Otsu technique is a popular one, and it's good practice to implement it on your own at least once.

http://www.imageprocessingplace.com/

But above all: give it a try! Look at a histogram of grayscale values in your image and then ask yourself whether the distribution of grayscale values immediately suggests a technique. Even if you're not sure how to work with the histogram, you might have an idea of the solution that you'd like to implement.

It's not easy to learn image processing just from online sources in bits and pieces. If you plan to work in image processing, and if you can afford it, try to have your own copy of at least one image processing book for reference. The Gonzalez and Woods textbook is a good standard reference, and is used in a number of undergraduate courses.

You can use the OpenCV library, buy the O'Reilly book, read forum posts and questions about OpenCV, etc., but in my opinion the Gonzalez and Woods book will ease you into the subject more gracefully.

As always, if you are trying to solve a problem related to a specific image, please post the image. People who work in image processing for a living can often tell you a lot based just on an image; otherwise it can be difficult to provide more specific advice than just "read this chapter." A binarization technique that works well for one set of images can perform poorly for another set of images.


"but what are the other two plots (3rd and 4th) ?..."

Good questions.  This is the way I understand the two diagrams:

3) -

The Wavelet transform can be thought of as a bank of filters with progressively wider bandwidths (working from the bottom to the top of figure 1).  If you look at the output of the transform for a specific scale (e.g. scale 24) such as figure 3, you can think of this graph as a plot of your original signal after it has been passed through a bandpass filter where the characteristics of the filter are defined by the wavelet and the particular scale (24).  This plot is the magnitude of the signal in the particular passband at sequential increments of time.  If you reconstructed your signal from these coefficients you would end up with the signal resulting from the bandpass filter at that scale.

4)-

The continuous wavelet transform is particularly good at making irregularities in a signal "visible."  In this context, an irregularity is a discontinuity or a sharp change or "cusp" in your signal.  When sharp changes occur in a signal, the CWT will produce large coefficient values across many or all scales around the time of the irregularity.  The function that produced diagram 4 attempts to locate such irregularities by examining the peak values of the CWT.

As you can see from your data, this is not an exact science.  You're data doesn't really contain any "irregularities", but the algorithm picked up on some random changes and produced an output.

This type of processing is mostly useful in cases where the signal of interest is fairly regular but includes occasional sharp deviations.  Imagine what would happen if you had a voice signal and you lost 30 contiguous samples (VoIP packet loss) and then spliced the signal back together.  The resulting signal would be relatively smooth, but where the packet loss happened, the signal would have a sudden change in value.  The "local maxima" plot can reveal the location of this kind of disturbance in your data.


Good Heuristic Point Where Sending FFT to GPU is Computationally Worth It?

I'm afraid there won't ever be a rule of thumb!

The reasons are manifold, mainly, that both the systems you're considering and the problems you're trying to solve vary across a very large range.

Problem-based aspects

You say you want to do an FFT – but that's always but half of what you actually want to do!

Need the FFT converted to it's abs², then mapped to colors, then displayed on a screen? Do it in the GPU, it's right where it belongs; fosphor does that, at an easy 200MS/s on capable PC/GPU combinations:

In that case, the size of the FFT doesn't even matter. Your data is going to be processed further on by the GPU, so do the FFT there.

On the other hand, you might want to do something that depends on a lot of checks on individual elements in the FFT, in the CPU? Possibly just one FFT, and then not many more?

In that case, your theoretical throughput doesn't help you at all. Just waiting for the data to get out of your CPU cache, back into coherent memory, so it can be DMA'ed to the GPU, where you then initiate an FFT (possibly wasting an interrupt/context switch on the way), just to wait till it's done, the GPU DMA'ed the data back to your main memory and you got it into your CPU's cache won't pay, even for medium sized FFTs.

So: this whole "high-latency math accelerator business" really only ever pays off if you can do something sensible while you wait. If you can't, there's a huge latency penalty.

System-based aspects

Ok, not going into too much detail here, but:

DSP systems are either CPU or Memory-bandwidth limited
If your GPU operation helps with CPU limititation, but puts additional data movement load on the memory interface, while in fact the rest of your system is memory-bandwidth limited, you're hurting yourself.
Same applies the other way around: Maybe your algorithm (the FFT in your specific size of interest) is CPU limited, but your GPU acceleration leads to additional interrupts
What is an FFT size that your CPU can do very well? That is probably defined by the sizes of it's L1 and L2 caches. A Xeon number-crunching CPU will have dozens of Megabytes of those, whereas an ARM running in a Jetson NVidia SoC will not.
What is the FFT size your graphics card is good at? There's a humongous difference in numbers of parallel threads, their flexibility and memory bandwidth across cards.
What is a metric for "good", at all? Just some strange ratio of throughput and latency, but maybe also energy and leaving the right kind of ressources free for other work?
What is your CPU<->main memory interface? Is it a quad-channel DDR4 interface running at nearly 2GHz, or is it single-channel DDR?
What is your GPU<->GPU memory interface?
What is your GPU<->main memory interface?
How well does your CPU<->GPU communication work for your specific use case?
Is there a high load on e.g. the PCIe bus, for example, because the same PCIe switch has to handle the data flowing in and out of your high-rate system (eg. storage, but more likely, 10 Gigabit ethernet or video data)?

So, the answer probably won't be satisfactory, but it really is:

Somewhere above 64 bins, somewhere below 220 bins, for a single-precision FFT. It depends.


Difference in information between IQ Sampling vs. Real Sampling at $2F_s$

The bandwidths might be the same.  I agree with Jason’s assessment.

Is there a difference?

Yes.  Let’s say we construct both signals such that their frequency domain’s are equivalent from $0-F_s$.  For the real signal the frequency domain would be limited to $X(F_s-f) = X^{*}(F_s+f)$.  For the complex signal it would be $X(f) = X(F_s+f)$.

Can the IQ signal be recreated from the real signal?

Yes.  Again, we assume the frequency domains are equal from $0-F_s$.  First, we filter the real signal to remove the frequency content from $F_s-2F_s$, which would result in an interpolated version of the IQ signal. Then we downsample by throwing away every other sample, yielding the IQ signal.  There are practical considerations, but that would be one way to go about it.


Indeed, you can apply the same technique:

As you said, given a sine signal with time shift $t_1$, which is $x(t)=\sin(w(t-t_1))$, this becomes $x(t)=\sin(w t - wt_1)=\sin(wt+\phi)$ and $\phi=-wt_1$.

Now, if you encounter a signal $x(t)=\sin(\frac{\pi}{4} t - \frac{\pi}{4})$ you identify $w=\frac{\pi}{4}$ and $\phi=-\frac{\pi}{4}$. What remains is to solve $\phi=-w t_1$ for $t_1$ and you have it.


Analyzing a signal that contains frequency content at Fs/2 doesn't seem to work unless there is a phase shift

the discrete function $$x_q[n]=\sin(\pi n)$$ is always zero for all of the integers $n$.

the discrete function $$x_i[n]=\cos(\pi n)$$ is always $(-1)^n$ for integer $n$.

so this general sinusoid at Nyquist that has a phase term:

$$\begin{align}
x[n] &= A \cos(\pi n + \theta) \\
     &= A \big( \cos(\pi n) \cos(\theta) - \sin(\pi n) \sin(\theta) \big) \\
     &= \big(A\cos(\theta)\big) \cos(\pi n)  +  \big(-A\sin(\theta)\big) \sin(\pi n) \\
     &= \big(A\cos(\theta)\big) (-1)^n  +  \big(-A\sin(\theta)\big) \cdot 0 \\
     &= B (-1)^n \\
\end{align}$$

So now, how do you tell the difference between a sampled sinusoid that had amplitude $A$ and a phase angle $\theta$ and another sampled sinusoid that has amplitude $B \triangleq\big(A\cos(\theta)\big)$ and no phase shift from the cosine?


Optimizing Two Parallel IIR Filter Sweeps

Two (linear, time invariant) filters in parallel, $H_1(z)$ and $H_2(z)$, are equivalent to a single filter $H(z) = H_1(z) + H_2(z)$ (i.e. the sum of both transfer functions).

You don't give any information on how those filters are represented. Normally, they are represented in terms of the recursive ($a_k$) and non-recursive ($b_k$) coefficients (see Wikipedia for reference).

Now, if you are going to add both transfer functions, then you need to find the common denominator (which is the product of two polynomials, one from each filter), and more cumbersome yet is to find the new numerator polynomial coefficients.

If you are doing this automatically (since you adjust your filters periodically), I suggest you keep your algorithm as is.

Executing an IIR filter is not very CPU intensive since the order of the filters is normally not greater than 10 or so.

Note: if the filters were FIR, then you simply add the filter coefficients (there are no recursive coefficients (and no division of polynomials) to worry about). But FIR filters are normally much more CPU intensive than comparable IIR.


spectral subtraction for stationary noise

In almost all practical applications of spectral subtraction it's unrealistic to expect that the noise can be completely eliminated. Even if this were possible, the price to be paid would be severe distortion of the desired signal.

What you can try is tweak the parameters. The most obvious parameter to tune is the (over)-subtraction factor, which determines how much of the estimated noise is subtracted from the noisy spectrum. It is not uncommon to have over-subtraction factors of around $2.0$ or $2.5$. Of course you have to make sure that the resulting estimated speech (magnitude or power) spectrum remains non-negative.


How to design a continuously time-variant digital filter?

One approach would be to try to remove the frequency chirp from the observed data, thus translating all of the echoes to approximately baseband. I find this to be most straightforward by converting the observation to an analytic signal, then multiplying by a complex exponential whose instantaneous frequency is equal to $-1$ times the frequency chirp profile (while keeping its phase continuous). After dechirping the received data, you could then apply a lowpass filter to suppress any other sources that don't overlap in frequency with your chirp profile. If your followon analysis methods need to see the frequency ramp, you could reapply the chirp again by multiplying with another complex exponential.

The passband width of the lowpass filter defines how tightly around the transmitted tone that you reject other frequency components. The width of the passband would also need to be selected while taking the expected two-way propagation time of the transmitted signal; at time $t$, assuming a low-to-high frequency chirp, you might be transmitting frequency $f_t = f_c + \Delta f$, while the receiver is observing a delayed version of what you transmitted some time ago, for instance $f_r = f_c$. Your lowpass filter must have enough frequency coverage to cover the frequency slew of your chirp profile over the expected range of time delays. At the same time, though, you have an incentive to make the passband width as narrow as possible to reject other signal sources that are nearby your chirp profile in frequency, so as so often occurs in engineering, you have a tradeoff to examine.


Should scipy.signal.coherence be 1 for single input and output signals?

Your colleague is correct.

In

  Carter, G. Clifford. "Coherence and time delay estimation."
  Proceedings of the IEEE 75.2 (1987): 236-255.

a closed form expression for the pdf of the estimate $\mid \hat{\gamma^2}\mid$, the  $p(\mid \hat{\gamma^2}\mid |\; N,\mid \gamma^2 \mid)$, is given and the confidence pdf for $N=1$ is flat on $[0,1]$


Is Z-transform of $\sin(\omega_0n)$ same as that of $\sin(\omega_0n)u[n]$

The two sequences $\sin(\omega_0n)u[n]$ and $\sin(\omega_0n)$ are very different, so their transforms can't be identical. The $\mathcal{Z}$-transform of $x[n]=\sin(\omega_0n)u[n]$ can be computed as

$$X(z)=\sum_{n=0}^{\infty}\sin(\omega_0n)z^{-n}\tag{1}$$

which converges for $|z|>1$ to the expression given in the table you linked to. However, the $\mathcal{Z}$-transform of $x[n]=\sin(\omega_0n)$ would be

$$X(z)\overset{?}{=}\sum_{n=-\infty}^{\infty}\sin(\omega_0n)z^{-n}\tag{2}$$

which doesn't converge anywhere, and hence does not exist.

Consequently, the $\mathcal{Z}$-transform is no suitable tool for computing the response of a discrete-time LTI system to a sinusoidal input. However, note that the discrete-time Fourier transform (DTFT) of $x[n]=\sin(\omega_0n)$ does exist. It contains Dirac delta impulses:

$$\text{DTFT}\{\sin(\omega_0n)\}=-j\pi\left[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)\right]\tag{3}$$


According to my reading, LDPC is almost the best channel coding we can use for channel coding, but I have a question regarding that kind of coding.

Mentally, something being "the best" should always instantly raise a mental flag for you, saying "under which conditions, according to which measure".
It is right that iterative LDPC decoders can achieve maximum likelihood performance (i.e. being the best possible decoder), but only when

The code is large and
there's an infinite number of iterations.

While 2. is never fulfilled, there's often a number of iterations after which gains are small enough for people to simply stop and limit complexity.
The first condition, however, is often a pretty complicated to fulfill one:

will the LPDC will take the whole 2048 (bits) together?

That is a small LDPC code (128 bits would really be tiny, and I don't think I've seen any sensible OFDM application do that; least-rate IoT modes might be interested in doing that on the uplink, but that doesn't match the OFDM approach).
At a small size like 2048 bits, less complex codes and decoders might be comparable or even better suited for your use case and error model.  (PS: Have an error model before deciding on the code you use! There's an awesome website, http://pretty-good-codes.org/ , which is sadly offline now, which compares many codes with metrics.)
Try to put more bits into a single codeword if you want to harness the powers of LDPC codes. For example, in DVB-T2 (which is an OFDM system), the short codewords are 16000 bits long, the normal ones are 64800 bits long.

I mean it will take bits from the whole bits for coding and the rest for data OR it will take group by group, which means for example group of bits, and then code them for data and for parity or coding.

Neither. You take blocksize · rate (so, in your proposed system, 2048·1/2=1024, but really, use larger blocks and established LDPC codes) information bits and encode them as one. These are not systematic, usually, so you don't get separate redundancy bit. You get a code word that is blocksize bits long, and doesn't contain the original bits in any structured way, usually. (Systematic LDPCs are usually undesirable.)
You will need to use a decoder to get the original bits from the codeword.


Why lattice-structured all-pole IIR filter becomes stable, if the k-coefficients satisfies $|k_i| < 1$?

I'm referring to the formulas on page 9 of this document. The denominator polynomial of the transfer function at lattice stage number $m$ is given by

$$\begin{align}A_m(z)&=A_{m-1}(z)+k_mz^{-1}B_{m-1}(z)\\
\text{with}\quad B_m(z)&=z^{-m}A_m(1/z)\\
\text{and}\quad A_0(z)&=B_0(z)=1\end{align}\tag{1}$$

where $k_m$ is the $m^{th}$ reflection coefficient.

You can prove that $A_m(z)$ has all its zeros inside the unit circle if this is the case for $A_{m-1}(z)$ and if $|k_m|<1$. This follows from Rouché's Theorem, which says that if $F(z)$ and $G(z)$ are two complex-valued functions that are analytic in and on some closed contour $C$ (which in our case is the unit circle), then $F(z)$ and $F(z)+G(z)$ have the same number of zeros inside $C$ if $|G(z)|<|F(z)|$ is satisfied on $C$.

The proof goes as follows. Define $F(z)=z^mA_{m-1}(z)$ and $G(z)=k_mz^{m-1}B_{m-1}(z)$, which, according to $(1)$, means that $F(z)+G(z)=z^mA_{m}(z)$. The zeros of $F(z)$ are the $m-1$ zeros of $A_{m-1}(z)$ plus one zero at $z=0$. Consequently, if all zeros of $A_{m-1}(z)$ are inside the unit circle, then also all zeros of $F(z)$ are inside the unit circle. From the definition of $B_m(z)$ in $(1)$ it follows that $B_{m-1}(z)$ has the same magnitude as $A_{m-1}(z)$ on the unit circle (i.e., for $|z|=1$). This means that $|F(z)|<|G(z)|$ is satisfied on the unit circle if $|k_m|<1$. In this case $F(z)+G(z)$ and, consequently, $A_m(z)$ have all their zeros inside the unit circle. This is the case because $F(z)$ has all its $m$ zeros inside the unit circle, and according to Rouché's theorem $F(z)+G(z)$, which has exactly $m$ zeros, has the same number of zeros inside the unit circle as $F(z)$. So $F(z)+G(z)$ must have all its $m$ zeros inside the unit circle.

This proves that all denominator polynomials $A_m(z)$ in the iteration given by $(1)$ correspond to stable transfer functions if the iteration is started with $A_0(z)=B_0(z)=1$ and if $|k_m|<1$ is satisfied for all $m$.

I've learned the above proof based on Rouché's theorem from Jackson's excellent book Digital Filters and Signal Processing.
Another good explanation of lattice structures can be found in this book by Zelniker and Taylor.


Algorithm correcting chromatic aberration on hyperspectral images

Having worked with hyperspectral images from experimental setups, I know it's hard to fix these little details.

The best way my group has to fix alignment and chromatic shifts is to use calibration images: Take a picture of something simple for which you know how the result should look like, and use the same correction on every pictures you take hereafter.

A simple example I'm thinking would be to take a picture of a small white dot against a black background. Every misaligned color will create a colored dot at a different place. You then find the mathematical operation necessary to put all these colored dots back together, creating a single white one (it will probably be a translation with a color dependant magnitude).

IF this shift is always the same from picture to picture, you can implement it in your data reduction code. IF not, you will have to take a calibration picture before every picture you take.

If you're lucky, it will be a translation. If you're not, the transformation might also depend on where you are in the picture. In that case, use a grid of white dots.


How can I generate time-variant Rayleigh channel in MATLAB

Yes, you can have time variant channel in MATLAB either by using comm.RayleighChannel()  , you can read about it in help of matlab.

or the other way you use the same command you have mentioned but with each symbol the channel should be changed. so it will be like you will generate a channel which is Quasi-static. You above command will become:

h = 1/sqrt(2)*[randn(n_symbol,M+1) + 1i*randn(n_symbol,M+1)];

where n_symbol represents the number of symbols you have. then to do the convolution, must be use the loop of for.

for jj = 1:nSym

  new_signal(jj,:) = conv(h(jj,:),x(jj,:));   %% x is your signal

end

good luck


The impulse response and frequency response of a length N zero-phase FIR lowpass filter can be designed using the Remez exchange algorithm.

N = 11; % filter length - must be odd

b = [0 0.1 0.2 0.5]*2; % band edges

M = [1 1 0 0 ]; % desired band values

h = remez(N-1,b,M); % Remez multiple exchange design

The impulse response h is returned in linear-phase form, so it must be left-shifted (N-1)/2 = 5 samples to make it zero-phase.
Answer taken from https://www.researchgate.net/post/How_can_I_design_a_null_phase_FIR_filter.

I implemented this solution some time ago and still works very well in my code.


Which IIR filters approximate a Gaussian filter?

The Deriche and van Vliet filters are heuristics.  In both cases they choose the locations of poles and zeros to minimize either the RMS difference or the maximum difference of the filter's impulse response from a Gaussian.

Both filters are causal-anti-causal pairs.  So I think they have no phase error or group delay, but you need to be able to run them backwards on the data as well as forwards.  This makes them popular in image processing, but perhaps restricts their applicability elsewhere.

That they are heuristics is evidenced by the wealth of papers that tweak them.  For example a google search (while I was looking for the link to the Deriche paper) turned up this one that tries to solve the problem that the Deriche derivative-of-gaussian filter doesn't have an exactly 0 DC response.  There are also some interesting issues about initializing the boundary conditions correctly.

I have found the following overview a good resource:
Dave Hale, Recursive Gaussian ﬁlters, Colorado School of Mines Center for Wave Phenomena report CWP-546.


How can detect the Left Ventricle (LV) Cavity into the MR dicom file?

I'm looking at this DICOM now and it is very hard to tell where I am in the brain (since there is only one slice). Anyhow, tumours tend to be bright irregularities. The thing I think is a tumour certainly isn't in the LV.

You can just draw an ROI over it using something like FSLVIEW if you are planning to do it manually, or you can have some fun and use something that looks for irregularities across hemispheres in the edges as detected by an algorithm.

http://www.medicalimagingandgraphics.com/article/S0895-6111(11)00079-6/abstract

Hope this helps...


AR model order selection for half second EEG fragments

Do you want to select the model order adaptively based on the content of the input EEG, or choose the AR order as a design parameter for a system that will always be the same?. In the second case (not adaptive), maybe if you have some sample EEGs you could try a a lot of models from order = 0, to some number N (large enough but not so that running all the models takes forever). And Graph the Error between the EEG signal and white noise filtered by the model. You should arribe at a graph or Error vs order. At some value of N the error only marginally decreases. That is your model order.

EDIT

Matlab offers AIC and others.


Is there any computational method to prove whether a series is stationary or not?

1 Proof

You're misusing the word proof. Remember that a proof, even one led with stochastic methods, always leads to an absolute "If A, then B, no doubt" statement.

Since your $x$ is a realization of a stochastic process of which you don't know many properties, pretty much everything is possible – for example, a perfectly white process might have a realization that by sheer luck is extremely periodic. Now, if you have a finite observation with discrete values, which is exactly what computational implies, you can only make a statement that seeing such a realization is extremely unlikely, yet not impossible – which shows the difference between a proof and a strong suspicion based on good indication¹.

You can often calculate likelihoods that a specific observation stems from a given emitting quality (here: realization of a stochastic process being caused by a stationary process), and say that with a confidence of $x \%$ you can say that the process behind this finite observation is stationary.

Having a confidence that you need to incorporate into your statements validity is the crux of every descriptive science – and, having worked with medical scientists before, it happens very easily to claim you've proven e.g. effectiveness of a specific cure, but that is a different, non-mathematical, non-DSP meaning of proving things. Just because something worked 100 out of 100 times doesn't mean it's right. It might just imply some high probability it is, and understanding that has very interesting consequences for understanding scientific publications – for example, even for high-quality medical papers whose statements are true with a confidence of $0.98$, the probability that of 35 papers all are correct is less than $\frac12$. Social sciences often have to work with lesser confidences, due to many things only being observable on smaller sample sets.

2 Inference from a single Observation

Without making a very strong statement on the stochastic process², you can't say anything about different realizations of the same process by looking at a single realization.

In other words, having a single $x_k[n]$ will not allow you to make a statement on $x_l[n],\,l\ne k$, unless you claim that you can deduct how different realizations will behave just by looking at one limited one. That's a statement that I personally find just as strong (or actually, stronger) than stationarity.

3 Stationarity

This is a property of a stochastic process. If you just say "I want to show whether this is stationary", people will assume you mean strictly³ stationary.

This implies that the joint probability distribution⁴ is invariant to shifts in time. Well, your discrete observation started at a specific time, so no statement can be given (see 2).

Now, since you're looking at a digital signal (with finite observation duration/data points), this isn't as much of a setback as one would instantly presume: due to the observation duration being limited, we obviously can't make any statements about how things act when looking at things that happen at scales much larger than the observation.

Hence, here, the term of wide-sense⁵ stationarity (WSS) saves the day by reducing the statement from

  "this process' joint probability distribution is time-invariant"

to

  "when comparing a realization of the process with a time-shifted version of itself, the stochastic properties of that only depend on the time shift, not the absolute times."

Now, this second statement is based on autocorrelation – and that is easily estimated/computable.

In essence: if different realizations have the same autocorrelation function, the underlying process is wide-sense stationary.

Being limited observations, the estimated autocorrelation function will not be identical – but it can be similar.

By defining a metric for similarity and a threshold for that metric based on error probability, statements on the probability of that process being WSS can be made.

In the most intuitive case, similarity could be understood as the root difference of the squares of the autocorrelation estimates.

By definition, the power spectral density (PSD) of a process is the Fourier transform of the autocorrelation – and happily, thanks to the Wiener-Chintschin theorem, this is, for WSS processes, equivalent to the expectation value of the squared Fourier transform.

Hence, comparing the PSD estimate done by squaring the discrete Fourier transform (DFT) of observations of different realizations of the same process can be a way to describe how likely a process is stationary – especially if your underlying phenomenon is harmonic/periodic, doing the DFT is something that might be helpful, anyway, and with the FFT, calculating $\left(\text{DFT}(x_i)\right)^2$ is more computationally efficient than estimating autocorrelation based on convolution.

¹ This is really a side note: Watch 12 Angry Men, the 1957 Movie. Claiming something is proven is a big statement.

² ergodicity

³ strictly, or strongly stationary – I must admit, not being a native speaker, I'm not 100% sure there's no difference here; if someone could comment on that, I'd be grateful.

⁴ your signal is 1D, so this is equivalent to the probability density function being invariant

⁵ weak or wide-sense – same as ³


In general this is not possible, since you do not have full spectral information for each pixel, but only coordinates in some low-dimensional color space (i.e. RGB).

If your image data is in RGB, then the best you can do is to simply take the green color channel (which will be most responsive to 550 nm) and throw away the red and blue.

Depending on the actual response of your RGB channels, you might be able to better estimate the 550nm contribution via a linear combination of red and green (see Spectral Sensitivity).


How fast is the signal you are trying to analyze? The maximum theoretical frequency you can handle is going to be 100KHz (sampling rate / 2). Giving a "true" answer with 1 usec resolution is not possible with a 1MHz sampling rate. Practically this may not matter.

Without specific timing constraints it is hard to give a good answer. To get your signals to have the same sampling rate you can do linear interpolation between the points of your slow data or average your fast data.

I recommend averaging your fast data. If you are working in counts from an ADC summing 5 of the fast samples and dividing by 5 will be fast. Your compiler will likely use reciprocal multiplication to avoid the division operation.

Trying to line up the signals based on a single peak is likely to give incorrect results.
An alternative method is to maintain an auxiliary array with a selection of the N highest peaks stored. Do not simply get the N largest points; they will likely be clustered together and will not give a good result. Instead track the largest point and do not allow a new peak unless the signal has dropped below some threshold value (say 90% of the peak) to introduce hysteresis. The exact number of peaks should be chosen on based on the shape of your data and the desired level of confidence as well as memory and CPU time.

To match the peaks scan through the second data set maintaining a best fit offset. When you find a possible match to the first peak (a value within some tolerance) you can then compare the "future" peaks. Finding the offset with minimal error for the N peak points is not too onerous (sqrt((x_1 - xbar)^2 + ... (x_n - xbar)^2)).

This can be done in only one pass through each of the data sets with only the extra memory required being the auxiliary array that keeps the N peak points.


How do you apply a filter after DFT on an image?

Based on your C++ snipped, the correct way to do complex multiplication between samples of DFT sequences should be as follows:

$$ Z[k] =  X[k]\cdot Y[k] = (a + jb)(c+jd) = (ac-bd) + j (ad+bc) $$

Where the real part of the result is $ac-bd$, and the imaginary part is $ad+bc$.

Then you can take the inverse 2D-DFT.

Note that, for an alias free circular convolution implementation, your DFT length should be at least $100 + (5-1)/2 = 102 \times 102$ points. Then after the inverse DFT of $ 102 \times 102 $ points, discard the first $2$ rows and columns and reatin the remaining $100 \times 100$ part which are the samples of the filtered image.


Quantify the difference between two audio spectrograms

You could do a cross correlation between signals to even out the time offset.
But why go for spectrogram in the first place? Do you have acoustic conditions varying over time? If not, go for the spectrum, should do the same job but without the time alignement problem.


Calculating Cross Power Spectral density between two complex signals

Cross correlation is defined as $R_{xy}(\tau) = E\{x(t)y^*(t+\tau)\}$. You missed to take the conjugate of $x_2$. So $x_1x_2^* = (c+id)(a-ib) = ac+bd +i(ad-bc)$.

Cross PSD is the fourier transform of $R_{xy}(\tau)$, not the two equations you mentioned.