A reverb effect usually includes an added, delayed version of the original signal.

Calculate an autocorrelation. If there is significant amount of power at delays larger than what you'd expect for a normal-sized room (say, 1/5 of a second), that'd be a good indication.


just did a Google image search and these diagrams have popped up.

these are radix-2 butterflies with expanded detail (showing where the "twiddle factor", Q, goes) for various DIF and DIT FFTs.

a radix-2 FFT (the simplest in concept) uses only radix-2 butterfly operations, each with 2 complex numbers going and and 2 complex numbers coming out.

these four radix-2 butterflies combined

make a radix-4 butterfly:

the image on the right just above would be a short-hand notation for a radix-4 butterfly.  if there were 8 signals going into the dot and 8 coming out, it would be a radix-8 butterfly.

let $p \in \{1, 2, ...\}$ and $P \triangleq 2^p$.  essentially a radix-$P$ FFT is an FFT that uses radix-$P$ butterflies and a radix-$P$ butterfly is, in and of itself, a complete $P$-point FFT with some phase adjustment on the internal twiddle factors.

a split-radix FFT has butterflies of different sizes.  here is the best example i can think of:

a 256-point FFT can be done entirely using radix-4 butterflies because 256 is a power of 4: $256=4^4$.  there would be 4 FFT passes and all of the passes would have radix-4 butterflies.  a 256-point FFT can also be done with a radix-2 FFT and it would be 8 passes, all with radix-2 butterflies, instead of 4 passes.

now suppose you wanted a 512-point FFT.  if you started with a radix-2 FFT, all you need to do is double the size of the buffer and add one more pass because $512=2^9$.  so there is no split radix if you wanted to do that.  turns out that the radix-4 FFT is a little more efficient than the radix-2 so you might want to replace 8 contiguous passes of the radix-2 FFT with 4 passes of a radix-4 FFT and one final (or one initial) pass using 2-point butterflies.  that is a split-radix FFT.


How to scale measured QPSK points to determine magnitude error?

Even though have AGC etc in your system, it is not ensured that the received constellation points will lie exactly at the points that you would expect from the TX constellation.

There can be residual offset, gain and phase rotation (at least). In the following , there's a simple method to estimate and remove offset and gain. Essentially, you calculate the mean of the RX constellation, which is your estimated offset. Then, after offset removal, you measure the energy of the constellation and scale the constellation accordingly, such that the constellation has same energy as tx constellation energy:

N = 10000 # number of TX symbols
QAM = np.array([1+1j, 1-1j, -1+1j, -1-1j])
txEnergy = np.var(QAM)

tx = QAM[np.random.randint(4, size=(N,))]

chanGain = 2 # channel gain (combined of AGC, channel and all other processing)
offset = 0.3-0.5j # channel offset
sigma2 = 0.3  # noise variance

noise = np.sqrt(sigma2/2) * (np.random.randn(N) + 1j*np.random.randn(N))

rx = chanGain * tx + offset + noise  # emulate effects of channel + AGC + ...
# RX is the received signal after demodulation

# estimate mean and remove it from the data
rx_mean = np.mean(rx)
rx_mean_correct = rx - rx_mean

# estimate the energy of the constellation
rx_energy = np.var(rx_mean_correct)

print ("Estimated vs real offset: %s <-> %s" % (rx_mean, offset))
print ("Estimated vs real energy: %s <-> %s" % (rx_energy, chanGain**2*txEnergy))

# scale rx constellation to have same energy as tx constellation
corrected = rx_mean_correct / np.sqrt(rx_energy) * np.sqrt(txEnergy)

# calculate constellation MSE
MSE = np.var(rx_mean_correct - tx)
print ("RX MSE: %f" % MSE)

# plot the stuff
plt.plot(rx.real, rx.imag, 'x', label='Received')
plt.plot(corrected.real, corrected.imag, 'x', label='corrected')
plt.plot(tx.real, tx.imag, 'x', label='tx')
plt.legend();

Program output

Estimated vs real offset: (0.300844347234-0.50994604822j) <-> (0.3-0.5j)
Estimated vs real energy: 8.29630318416 <-> 8.0
RX MSE: 2.297161

In order to remove the phase rotation, you need to have some information about the transmitted symbols. I.e. after the received constellation has been processed as above, you need to compare the RX constellation with the known TX points to get the rotation between them. Then, you rotate the RX constellation back according to the estimated phase rotation.


For a continuous-time signal $x(t)$ with Fourier Transform $X(f)$, the
Fourier transform of $x(t-t_0)$ is $e^{-j2\pi ft_0}X(f)$, that is, the value
of $X(3)$, say, gets changed to $e^{-j2\pi 3t_0}X(3)$ which has the same
magnitude as $X(3)$ but a different phase.

For a discrete-time continuous-amplitude signal $x[n]$ with DFT $X[k]$
of length $N$, you have essentially discarded every piece of information
conveyed to you by $x[-1], x[-2], \ldots $ and by $x[N], x[N+1], x[N+2], \ldots$
and have retained only $x[0], x[1], \ldots, x[N-1]$. If you phase-shift
the DFT appropriately, then yes, you will get a circular shift of the $x[n]$
values. The result will be the linear shift of the $x[n]$ values only
in the fortuitous circumstance that the sequence $x[n]$ is itself periodic
with period $N$, and so the untrained observer such as myself cannot tell
if that symbol that just shifted into my DFT wndow is $x[-1]$ coming from
the left (linear shift) or $x[N-1]$ from the right (circular shift)
because by golly they look exactly the same to me.  Savvier denizens
of this forum will give you more elaborate answers explaining how you
can tell the difference between the two cases.


Non-negative matrix factorization for audio separation - why does it work?

Recall that the columns of $\mathbf{W}$ can be thought of as "basis" vectors (or elements of a dictionary - the building blocks of any signal) and elements in each column in $\mathbf{H}$ give the corresponding weights (that vary over time). This allows us to decompose the spectrogram based on not just frequency components but also temporal onset information; so it does more than what a simple bandpass filter or comb filter would do.

Non-negative matrix factorization of a spectrogram will not magically separate a piano and a singer singing in the same pitch. However it gives a useful approximation in terms of a sum of weighted basis vectors which (hopefully) can be split up into contributions from different sources because it's unlikely for each source to occupy exactly the same frequency bin at exactly the same time instant.

More concretely, let $\mathbf{V}$ be $M\times N$, $\mathbf{W}$ be $M \times K$ and $\mathbf{H}$ be $K \times N$. So we have $M$ frequency bins, $N$ time samples and $K$ decomposed components. If $\mathbf{w}_i$ are the columns of $\mathbf{W}$ and $\mathbf{h}_i$ are the rows of $\mathbf{H}$ we can write:
$$
\mathbf{V} \approx \sum_{i=1}^K \mathbf{w}_i \mathbf{h}_i^T.
$$
If we know there are only two sources in the recording (piano and singer) we can try to split up their contributions by choosing subsets of the columns of $\mathbf{W}$ and corresponding subset of the rows of $\mathbf{H}$. So if $S \subset \{1,\ldots,K\}$ then the part contributed by the piano can be denoted by:
$$
\mathbf{V}_{\mbox{piano}} = \sum_{i \in S} \mathbf{w}_i \mathbf{h}_i^T
$$
and the part contributed by the singer is:
$$
\mathbf{V}_{\mbox{singer}} = \sum_{i \in \{1,\ldots,K\}\setminus S} \mathbf{w}_i \mathbf{h}_i^T
$$

In reality, we will likely end up with a decomposition which never achieves this separation exactly. That is, there will be $\mathbf{w}_i$'s that have contributions from both the singer and the piano making it difficult to separate the two.

Here's a Python notebook showing this procedure for a mixture of drums and guitar: http://nbviewer.jupyter.org/gist/ingle/93de575aac6a4c7fe9ee5f3d5adab98f (Or if that doesn't work, here: https://gist.github.com/ingle/93de575aac6a4c7fe9ee5f3d5adab98f)

Note that the NMF algorithm only generates a decomposition. It cannot select subsets of $\{\mathbf{w}_i\}$ corresponding to each source. In the Python example, there's a manual step of figuring out which $\mathbf{w}_i$'s sound most like the guitar v/s drum. One can perhaps automate this step by noting that the drum's $\mathbf{w}_i$ vectors have more stuff in lower frequencies.

Analyzing each column (frame) of the spectrogram using a bank of bandpass filters is another way of decomposing the spectrogram. However, note that the decomposition generated by NMF is low rank i.e. parsimonious. In the Python example, it was much easier to manually select subsets of 16 columns of $\mathbf{W}$ corresponding to the two sources. With a bank of bandpass filters we would have had to turn many more knobs (# filters, locations of pass-bands for each frame) and the number of parameter combinations can grow very quickly.

References:

https://ccrma.stanford.edu/~njb/teaching/sstutorial/part2.pdf

http://musicinformationretrieval.com/nmf.html


ADC sampling rate and signal bandwidth for FHSS

So for frequency hopping spread spectrum (FHSS),when the receiver is unaware of the transmitter's hopping pattern, does that mean the ADC sample rate needs to be higher than twice the entire hopping bandwidth?

Yes.

If I use a 2Msps (4x data BW) ADC it looks like I will still recover my data; sure it under-samples some of MHz modulated frequency but it does capture entirely my data.

No.  See above.  It'll alias whatever is in the entire 24MHz band into the 1MHz that it can resolve, or you'll precede it with an anti-alias filter and lose $\frac{23}{24}$ of your signal.
On the bright side, if you're taking 24MHz worth of bandwidth and only using 1MHz of it, that's close to $14\mathrm{dB}$ of processing gain, so while you need an ADC that's 24 times faster, it can be 4 bits narrower.


Are there any good texts on EEG analysis geared towards programmers

Your application is not something that can be done very easily. The project would be well-suited for a PhD topic to be done over 4-5 years.

If you really want to work towards it, learn

1) Signal processing. Knowing just one technique (like the DWT) will not get you anywhere. Get a thick book on the basics.

2) Machine learning (supervised and unsupervised). Your application will need classifiers.

3) Basics of EEG. Again, get a thick book. Focus on continuous-data analysis.

4) Start reading publications in your field. Try sciencedirect.com, ieeexplore.org, google scholar, etc.

5) Offer to work as a research/teaching assistant with neuroscientists at your local university.


Is that to be expected when doing FFT analysis of any signal?

No. Either your calculation or the implementation of your analysis is off.

Secondly, I thought that the fundamental frequencies in my signal should have the highest amplitudes from all other peaks.

Why did you think that? There are many signals where the fundamental is lower than the harmonics or even entirely absent.

Finally, I know that the other peaks represent the rest of the harmonics of my signal

These are not all harmonics of your expected frequencies. You have also peaks at 0.25Hz, 0.5Hz. Could also be noise spurs, intermodulation, etc.

what do they exactly say about my signal?

Without knowing what that signal actual is, it's impossible to even  guess.


Digital Gaussian White noise signal generation in C++

If you are generating sample of white noise (Normal distribution value) with frequency of 48 kHz, then you already have this in frequency range of 0 - 24 kHz (at least you should if generator is truly random). If you use same samples with different sampling frequency, then they will also be a white noise.

Now if you want to make it a band-limited white noise, then you must use a digital filter. Choice is up to you, but for sake of simplicity I would recommend FIR filter. You must use it to filter in real time your incoming samples. How to do it? Well, that's totally up to you, you can even use short buffer to perform filtering by linear convolution for each generated sample but it won't be very efficient.


Relation between curve fitting and filtering

Generally speaking, I would not call a polynomial, in general,  a low-pass approximation to a signal. Per se, talking practically, standard definitions of the Fourier transform do not allow polynomials on $]-\infty,\infty[$ to have a useful spectrum. Generally speaking, polynomials (except constants) tend to infinity at fast rate, hence cannot be really limited in frequencies.

So it is difficult to say a polynomial filters out certain frequencies. However, locally, on an interval, they can pass nicely through data points, so they can smooth data, to some extend. Typical examples are sensor trends or drifts, when they are monotonous, and peaks or bumps.

You can find a lot of works on polynomial trend filtering, or polynomial peak filters. When the cost function for those fits is a squared error, you end up with a linear system in the polynomial coefficients, which turns into a weighted average of signal samples, that can be interpreted as coefficients of a linear filter. A famous example is the Savitsky-Golay filter, which was recently analyzed in frequency in 2011 by R. Shafer: What Is a Savitzky-Golay Filter?. This family of filters is also interesting as they provide smoothed derivatives.

But you can also fit a polynomial with other loss functions, like the $\ell_1$-norm, or with additional penalties. They are currently active topic in signal processing (I am working on some extension to Savitsky-Golay filters). Here are a few references, to start from:

On digital smoothing filters: A brief review of closed form solutions and two new filter approaches, 1986
Polynomial-Based Interpolation Filters---Part I: Filter Synthesis, 2007
$\ell_1$ Trend Filtering, 2009
Symmetrizing Smoothing Filters, 2013


1) Yes, your solution is correct.

2) If you assume that all a priori probabilities are equal, then the boundaries for AWGN is always the middle points between the possible values of X.  In this case, then, the decision boundaries are at 75 and 175.

I believe that this rule (decision boundary at the middle points) can be generalised to applying to any noise probability distribution that is symmetric and monotonically decreases as the distance from zero increases.


If a K-NN ($ k $ Nearest Neighbors) Algorithm Performs Very Well for Low $ k $, Can Something Be Inferred About the Data Set?

If the result is consistent with a large test set than it means your training data is dense and well define the degrees of freedom of the problem.

If the training set was small it means there is a good separation between your classes in space.


Contour Integral and Residue Theory for Inverse $z$-Transform

It's a consequence of Cauchy's Residue Theorem, which can be used to evaluate contour integrals of analytic functions. Take a look at this answer for an example.


detection of periodicities in n-dimensional signals

Would these be limited to detecting periodicities that are aligned with one of the basis vectors of the n-D space?

Well, if you say "sampling is sufficient", then it follows that the whole space is the span of the base vectors (in fact, you usually have them orthogonal).

From that follows that a non-base-vector aligned oscillation can be represented by a linear combination of oscillations that are aligned with base vectors.

From that follows that yes, with your requirement that sampling frequency is sufficient to begin with, the usual methods of 1D frequency analysis work.

Note that you can't generally detect all periodicities with just a sufficiently fast sampling – if something is periodic with a period  longer than you observation window, you have no information about that. So, observation length must be sufficient, too, unless your signal model restricts both highest and lowest frequency.


What is the physical meaning of the convolution of two signals?

There's not particularly any "physical" meaning to the convolution operation. The main use of convolution in engineering is in describing the output of a linear, time-invariant (LTI) system. The input-output behavior of an LTI system can be characterized via its impulse response, and the output of an LTI system for any input signal $x(t)$ can be expressed as the convolution of the input signal with the system's impulse response.

Namely, if the signal $x(t)$ is applied to an LTI system with impulse response $h(t)$, then the output signal is:

$$
y(t) = x(t) * h(t) = \int_{-\infty}^{\infty}x(\tau)h(t - \tau)d\tau
$$

Like I said, there's not much of a physical interpretation, but you can think of a convolution qualitatively as "smearing" the energy present in $x(t)$ out in time in some way, dependent upon the shape of the impulse response $h(t)$. At an engineering level (rigorous mathematicians wouldn't approve), you can get some insight by looking more closely at the structure of the integrand itself. You can think of the output $y(t)$ as the sum of an infinite number of copies of the impulse response, each shifted by a slightly different time delay ($\tau$) and scaled according to the value of the input signal at the value of $t$ that corresponds to the delay: $x(\tau)$.

This sort of interpretation is similar to taking discrete-time convolution (discussed in Atul Ingle's answer) to a limit of an infinitesimally-short sample period, which again isn't fully mathematically sound, but makes for a decently intuitive way to visualize the action for a continuous-time system.


Your reasoning seems reasonable and correct to me. Assuming that the filter's impulse response $h(t)$ is real, then its magnitude response $\left|H(\omega)\right|$ will be symmetric about $\omega=0$, so the negative-frequency terms will be scaled by the same values as those on the positive side of the frequency axis. I don't see any reason why you would need to multiply or divide by 2 anywhere.

The only thing that might differ if you see a solution somewhere else might be the scaling factor at the front of the expression. The typical transform using angular frequency would be:

$$
\mathcal{F}\{\cos(ax)\} = \pi\left(\delta(\omega - a) + \delta(\omega + a)\right)
$$

which assumes a Fourier transform definition of:

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

It looks like you have an extra factor of $0.5$ at the front of your expression that you might have gotten from the transform. Note that different implementations or expressions of Fourier transforms often have arbitrary scaling factors included; what's most important is that if you apply the forward and inverse transforms in sequence, you get the original signal back:

$$
\mathcal{F}^{-1}\left\{\mathcal{F}\{x(t)\}\right\} = \mathcal{F}^{-1}\left\{X(\omega)\right\} = x(t)
$$

or at least to some known scaling factor.


Is there any way to align temporally two signals in spectrogram?

You can find a large number of  "model-based echo suppression" or "adaptive acoustic echo cancellation" (for the keywords). They often combine a sparsifying transformation (wavelets, time-frequency decompositions) and some form of fitting or adaptation of the model to the time and frequency localization where the echo dwells. Underlying hypotheszs are:

the model somehow looks like the echo,
the similarity will increase in the transformed domain,
this representation will help the morphing, or make is more efficient (faster, etc.)

For voice, you find a look of literature, overviews, etc. You may also look at neighboring domains: vibration-based non-destructive testing or monitoring, ultrasounds. I would mention seimic processing, where echoes are sometimes named multiples (waves bounding several times between layers). Here are two examples I have been involved in. A first one quite fast which uses 1-tap adaptive filters (hint: a complex filter), the second  more involved with optimization constraints. They use time-scale wavelets, but could be readily adapted to short-term Fourier transforms or multi-input/multi-output filter banks.

Adaptive multiple subtraction with wavelet-based complex unary Wiener filters, 2012  (preprint)
A Primal-Dual Proximal Algorithm for Sparse Template-Based Adaptive Filtering: Application to Seismic Multiple Removal, 2014 (preprint)


Inverted values - digital resolution issue?

This looks like a "signed number" problem.

An unsigned binary number with n bits has a range from 0 to 2^n -1

But a Signed number weights the MSB as -2^(n-1) so the range is -2^(n-1) to 2^(n-1) -1

Consider an unsigned signal going from 0 to 255 being viewed as though it was a 2's complement number. Notice how it would look as though it flips from +127 to -128 and continue up to -1 (255 corresponds to -1).

Is that what you're seeing?

I notice that you're max digital value is 65536. But the maximum value a 16 bit unsigned number can hold is 65535 .. 2^16 -1


Yes for a given camera with a fixed focus lens, sensor size and sensor to lens distance the optical magnification;
$$ M = |{ \frac{ h_{obj} }{h_{img} } }|$$
which is object vertical length divided by its image length on the sensor (a number less than unity, hence a reduction,  unless a macro image is taken with  a macro lens wihich has a true magnification) will be a constant for most of the proper object distance range.


Communication between two USRP : mismatched spectra between transmission and reception

A few (hopefully useful) comments and ideas:

The HackRF One is not a USRP.
If you're receiving with an NI 2921 USRP, you should be using UHD to interface with it, not osmocom.
Use a standalone spectrum analyzer, first to double-check that nobody else is transmiting at that frequency, and then to study your transmitted signal.
Transmit a single sine wave and verify that both the spectrum analyzer and the USRP can see it. Verify that you're transmitting it with enough power.
If that works, try with two sine waves separated by a few kHz.
It is normal that the received signals are at a slightly wrong frequency. Each device has its own idea of what 1 Hz is; this is why receivers need to run a carrier synchronization algorithm.
Your devices may not support a sampling frequency of exactly 200 kHz and may be using the closest available frequency. Check any warnings on the console. This may result in impaired performance.
In general, pay attention to any output on the console: you may be overrunning or underrunning your flowgraphs, or something else may be wrong.


What is a "Hampel Filter" and how does it work?

The documentation clearly describes its function: "For each sample of x, the function computes the median of a window composed of the sample and its six surrounding samples, three per side. It also estimates the standard deviation of each sample about its window median using the median absolute deviation. If a sample differs from the median by more than three standard deviations, it is replaced with the median."

The $3\sigma$ rule is enforced via the estimator $\sigma \approx 1.4826\ \operatorname{MAD}$.

The most costly operation must be the median of seven elements, possibly achieved by sorting, which can be performed efficiently in a sliding window.

Possibly relevant reference: "Hampel F. R., ”The influence curve and its role in robust estimation,” Journal
of the American Statistical Association, 69, 382–393, 1974."


In principle what you describe is reasonable. Here's an example in Python that implements your idea:

Fs = 44100
T = 100
t = np.arange(0, T, 1/Fs)
N = 1024
t = t[:N*int(len(t)/N)]

f = 0.7*t*t+2000

signal = np.sin(2*np.pi*f*t)  # generate some signal for our analysis

chunks = signal.reshape((-1, N))   # create chunks of the signal of length N
C = np.fft.fft(chunks, axis=1)     # perform FFT of each chunk
C = C[:, :N//2] # keep only positive frequencies

Spec = 10*np.log10(abs(C*C))      # draw the spectrum
plt.imshow(Spec.T);

I can only recommend to try out your DSP algorithms not in Delphi but in some scripting language like Python or Matlab, because you are much more flexible in your analysis and can try out new ideas much more quickly.

To improve the spectrum image, you can use a different window function (i.e. no rectangular window as I use here) and overlap the different chunks. Then, it would be closer to the short-time Fourier transform as noted by MBaz.


If you just want to demodulate, then it is not necessary to recover I and Q components. What you have is most likely the QAM symbols. For example, for $16$-QAM they are $16$ different integer values that are assigned to different constellation points. You need to know the type of assignment (i.e. coding) that is used in the QAM constellation at the modulator.

At the $16$-QAM modulator, every $\log_2 16=4$ data bits are mapped to a certain constellation point (which ultimately is identified by an integer). For instance, $0101$ can represent the symbol $5$, and so on. Gray coding is usually used for this purpose, but it is not always the case. So you need to find out the exact constellation mapping and then, the demodulation process is nothing more than the inverse of that mapping. So when there is a $5$ your data bits are $0101$.


What is the purpose of Balance / Fader in a digital Audio amplifier?

In a stereo amplifier the balance controls allows to set different gains for the left and right channel. That would be useful if the listener is not located right on on the center line between the speakers or if the acoustic loading the two speakers is very different. The balance control can be used to reposition the image in front of the listener.
A "fader" is more common a car. In addition to being able to balance between left and right a fader can also balance between front and rear as well. This control is mostly outdated: in a well tuned car, a fader will do more damage than good and many cars these days already have "front" "driver" "rear" etc. modes which are individually optimized for the use case.

In a signal flow what is the most apt position to place balance and fader(for example after before a compressor unit or at the starting of the signal flow or at the end of signal flow) and why ?

That depends a lot on what other functions the amplifiers has. It's mostly "in the middle", after the "flavoring" controls (different modes, tone controls, dynamic loudness, dynamic range adjustments, ...) but before the speaker management (cross over, EQ, limiter, bass manager, etc.).


How to understand relationships between ellipse and second moment matrix of Harris corner detector?

Consider a 2-by-2 diagonal matrix $M$ with diagonal elements $a>0$ and $b>0$. Then the equation

$$x^TMx=const\tag{1}$$

(with $x$ a column vector with elements $u$ and $v$) becomes

$$au^2 + bv^2 = const$$

which is simply the equation of an ellipse (or a circle if $a=b$). If $M$ is not diagonal, equation (1) still describes an ellipse if $M$ is positive definite (i.e. it has positive eigenvalues), but then its major and minor axes do not coincide with the Cartesian axes.

Also have a look at this.


You can implement discrete wavelet transforms on output1 and output2. Use LL from output1 and LH, HL, and HH to form a new block, then implement inverse wavelet transform.


how to find the inverse response of a system

Dan's answer -- to compute $H(s)$ as normal, and then compute $1/H(s)$ -- is equivalent to your suggestion of swapping the $x$ and $y$ (or doing it in one step by solving for $H_{yx}(s) = \frac{X(s)}{Y(s)}$).
In the Laplace domain it's justified by noting that
$$\frac{H(s)}{H(s)} = H(s)H_{yx}(s) = 1$$
In theory this means that a system followed by its inverse system has an output equal to its input.  In practice you need to make sure that $H_{xy}$ is realizable, and you have to recognize that if you're actually chaining real systems, then your $H(s)$ is just a model that may not capture all of your real-world system's real-world dynamics.  So the actual end-to-end equation is something like
$$\hat X(s) = \left(H(s) X(s) + N(s)\right) \hat H_{yx}(s)$$
where $\hat X$ denotes that it's just an estimate, $N(s)$ is whatever noise is injected in your measurement of $Y$ and your implementation of $\hat H_{yx}$, and $\hat H_{yx}$ emphasizes that it's a guess at your system model, not the real thing.
It's also useful to note that doing this in the Laplace domain immediately lets you know whether what you're trying is practical.  If $H(s) = 1 / (\tau_0 s + 1)$, for instance, you'd immediately see that $1 / H(s)$ is, strictly, unrealizable (because you can't have a true naked differentiator).
However, as long as noise is negligible you could approach it arbitrarily closely by setting $\hat H_{xy}(s) = \frac{\tau_0 s + 1}{\tau_1 s + 1}$, which basically says you're approximating a derivative as $\frac{s}{\tau_1 s + 1}$ -- then you could use my end-end equation above, along with your knowledge of the nature of $X(s)$ and an estimate of the noise you're going to inject at $N(s)$ to estimate your end-end system fidelity.


How is the DTFT of a periodic, sampled signal linked to the DFT?

I think that your misunderstanding comes from the fact that you assume that $f[n]$ is periodic. If $f[n]$ were periodic then its DTFT would be a Dirac comb. What is happening here is that you want to find the $N$-periodic sequence whose DFT coefficients equal equidistant samples of the DTFT of $f[n]$. And it turns out that that sequence is given by
$$\tilde{f}[n]=\sum_{k=-\infty}^{\infty}f[n-kN]\tag{1}$$
which is - as expected - an aliased version of the original sequence. Remember that sampling in one domain corresponds to aliasing in the other domain.
A derivation of $(1)$ is given in this answer.


How to Simulate AWGN (Additive White Gaussian Noise) in Communication Systems for Specific Bandwidth

To create Band Limited AWGN all you need is randn in MATLAB.

The question lies only in how to set its Standard Deviation.
To illustrate that, let's assume our AWGN generator has PSD which equals to $ {N}_{0} $.
Namely we have AWGN with zero mean and variance equals to $ \delta(0) {N}_{0} $.

Assuming we have limited bandwidth channel, hence an ideal LPF is applied.
Assume its cutoff frequency is $ {F}_{LPF} = \frac{{F}_{S}}{2} $.

Hence the Variance of the band limited AWGN is (Multiplication in frequency) the integral over its PSD multiplied by the Norm of the LPF (The factor 2 is for integration over the range $ -{F}_{LPF}:{F}_{LPF} $):

$$ {Var}_{BandLimitedAWGN} = 2 {F}_{LPF} {N}_{0} = {F}_{S} {N}_{0} $$

Now, generate in MATLAB, using randn noise with the corresponding STD (By the data of your simulation).

Good Luck!

Some Remarks

The Variance of the noise is independent of the signal (At least in the classic model).
The Variance of the noise is a function only of the analog channel and the Analog to Digital converter.
The classic model assumes that if the signal is sampled at $ {F}_{s} $ an ideal LPF with a cut off frequency of $ {F}_{LPF} = \frac{{F}_{S}}{2} $.
As I said, this is the frequency which sets the variance of the Band Limited noise.
The input signal (Which is the transmitted signal + noise) may have any bandwidth it might have, after the LPF its bandwidth is limited.
In order to minimize the energy of the noise in the system the LPF band width and the sampling rate should be as low as possible (Namely, the bandwidth of the signal in interest).
Though if the next step is "Matched Filter" the SNR will be maximized for any finite energy white noise (Or colored if the "Colorization" is known and the Matched Filter is accordingly updated).


Upsampling (by $M$) and downsampling (by $N$) are interchangeable if and only if $M$ and $N$ are co-prime (i.e., they have no common factors other than $1$). This is the case in your example because $2$ and $3$ are co-prime.

Think about your specific example in the time domain. Introduce a signal $u[n]$ between the upsampler and the downsampler. In the first case you have

$$u[n]=\begin{cases}x[n/3],& n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}$$

and

$$y[n]=u[4n]=\begin{cases}x[4n/3],& 4n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{1}$$

In the second case you have

$$u[n]=x[4n]\tag{2}$$

and

$$y[n]=\begin{cases}u[n/3],& n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{3}$$

With $(2)$ you get

$$y[n]=\begin{cases}x[4n/3],& 4n=3k,\quad k\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{4}$$

which is identical with $(1)$.


My preamble bits not enough for locking the timing recovery for a gaussian channel

Timing recovery is implemented as a control loop (measure the timing error and adjust the sample offset) and therefore would have an acquisition time depending on the control loop parameters (loop BW). A reasonable estimate for acquisition time assuming the start of acquisition is in the linear capture range of the loop is the first order 10% to 90% rise time estimate of $t_r = 0.35/BW$ where $t_r$ is the estimated rise time from 10% to 90% and $BW$ is the loop BW in Hz. This is accurate for a first order loop and a reasonable estimate whenever there is a dominant pole as that would set the settling time.

If the Initial state is not within the linear operating range of the loop or if there is lower SNR, this can take much longer. The best approach is to model this in a time domain simulation and observe the timing error versus time with step response testing to evaluate convergence time/ behavior and validate that the loop BW is as predicted from a more detailed control loop model from which these parameters can be computed and help identify issues such as capture range and non-linear behavior. A simulation can also help identify issues related the SNR as well as noise enhancement in the loop. For this purpose the preamble would be treated as infinitely long until it is certain the loop can work in the desired length in all conditions (will this loop operate with carrier offsets or will that be addressed before timing recovery, etc)


How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period?

This is basically what @hooman suggests: fit a parabola to the three points near the peak of the sample cross-correlation of the data.

Using the formula for $p$ here:
$$
p = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma}
$$
where $\alpha,\beta,$ and $\gamma$ are the values of the sample cross-correlation just before the peak, at the peak, and just after the peak respectively.

I've assumed that the signal is a white noise sequence and that a simple FIR filter can be used to simulate a fractional delay.

The plot below shows the FIR filter in the top plot, and 100 estimates of the fractional delay using this technique.

R Code Below

#35112
library(signal)

fractional_delay_filter <- function(delay, length)
{
  eps <- 10^-10
  t <- seq(0,length-1)
  t[t == 0 ] <- eps
  filter <- sinc(t-delay,length)
  return(filter)
}

# https://ccrma.stanford.edu/~jos/parshl/Peak_Detection_Steps_3.html
peak_location  <- function(ykm1,yk,ykp1)
{
  return(0.5*(ykm1 - ykp1)/(ykm1 - 2*yk + ykp1))
}

estimate_delay <- function(x,xd)
{
  cc <- ccf(xd,x)
  idx <- which.max(cc$acf)
  return(peak_location(cc$acf[idx-1], cc$acf[idx], cc$acf[idx+1]) + cc$lag[idx])
}

true_delay <- 10.5
d5p1 <- fractional_delay_filter(true_delay,32)
fz <- freqz(d5p1)
#plot(fz)
gd <- grpdelay(d5p1)
#plot(gd)

T <- 1000
Nres <- 100

results <- rep(0,Nres)
for (n in seq(1,Nres))
{
  x <- rnorm(T,1)
  xd <- filter(d5p1, 1, x)

  #plot(x, type='l')
  #lines(xd[1:T],col='red')å

  results[n] <- estimate_delay(x,xd[1:T])
}

par(mfrow=c(2,1))
plot(d5p1, type='l', col='blue')
title('Fractional filter')

plot(results,col='blue')
lines(c(1,Nres), c(true_delay,true_delay), col='green')
title('Estimates and true delay')

The sinc function is from the signal library (which you've said privately you have trouble installing in R). It's defined as:

sinc <- function(omega, N)
{
  eps <- 10^-56
  return( sin(omega*N/2)/(omega/2 + eps) )
}

So changing that sinc function to:

sinc_pjk <- function(x)
{
  return(sin(pi * x) / (pi * x))
}

and the test delay to 10.1 samples, Olli's version is less biased:

<>


The large spike at ~9MHz could be some form of interference from electronics near the receiver or internal to the receiver (maybe from using a poorly conditioned power source).  My best guess for the spike at DC is a DC bias from quantization.  If the digitizer truncates instead rounds, the resulting digital signal will have a non-zero mean and therefore a DC offset.  We had this happen in an FPGA downcoverter implementation at a previous job.  Once they changed the FPGA to round, the DC offset went away.  If you have ever used an Ettus USRP radio, you see the same DC offset in many of those boxes.


"I welcome the addition of a spectrogram showing many channels at one
  particular time interval."

If you don't want to include it, or don't have enough space, it is always worthwhile to put them in a supplementary file, attached to the journal, or separately on the Internet. I suggest the reviewer is looking for an onset phenomenon, and it would be easy to generate a few crops in  time interval $[1500-2500]$ side-by-side, potentially for only a tenth of your data (about 50 signals).


The benefits of a fractionally spaced equalizer

The more samples you have per symbol the less timing error you have since delay in samples is converted to delay in time by dividing by the sampling rate. I beleive the far more dominant reason for a fractionally spaced equalizer is the aliasing you mention as can be seen in the frequency response of the equalizer. For example, 2 samples per symbol can capture all the phase and amplitude distortion represented over the channel bandwidth uniquely as depicted in the graphic below. (This is assuming the signal was sufficiently filtered prior to being decimated to the fractional rate used).


How to recover signal repeated throughout other signal?

This is a source separation problem, which, in general is quite difficult to solve. One way to attack it is to use non-negative matrix factorization (NMF) with short-time Fourier analysis. Please see this question and a sample script here.

The overarching idea is to try to isolate your component sounds A, B, C,... using a time-frequency decomposition (or, perhaps, some other kind of subspace analysis using a different transform, like wavelet transform).

If the components don't always all play at exactly the same time, or have different onsets, you can hope to separate them temporally. If they don't all have identical frequency content, you have some hope of separating them in frequency. Or, a combination of both, which is what NMF of the STFT achieves.


How can I calculate the Number of bins per octave?

I do not know how to determine the ideal number of bins. However, I can shed some light on why in Librosa's CQT function, the default values are num_bins_per_octave=12 and n_bins=84. (This is my first answer so any feedback would be helpful. And because this does not directly answer the question, should I perhaps post this as a comment?)
The Librosa CQT function you mention can be found in the Librosa documentation here.
First some music theory
Western music consists of 12 identifiable types of pitches. Each of these 12 pitches is considered a new note. The 12 pitches are as follows:
C
C# / Db
D
D# / Eb
E
F
F# / Gb
G
G# / Ab
A
A# / Bb
B
The notes with slashes can be called by either name, depending upon which key we’re in. (It doesn't matter what a key is right now.)
These 12 pitches repeat in the same order throughout the complete span of human hearing. So after B, there would be another C, then another C#, another D, and so on and so forth, giving us multiple C pitches, C# pitches, D pitches, etc.
A specific group of these 12 pitches is defined by something called an octave. An octave is the distance between one note (like C#) and the next note bearing its same name (the next C# that’s either higher or lower).
The reason why an octave is divided into 12 pitches rather than, let's say 11 or 13 pitches, can be found here.
Librosa's default num_bins_per_octave is 12 specifically because there are 12 pitches in an octave.
The piano
A great visual of multiple consecutive octaves can be seen in this 76-key piano here:  Notice that the set of 12 unique pitches repeats over and over again on this keyboard as you move left to right (or vice versa), counting both the black and white keys.
A normal piano has 88-keys, 7 octaves plus four extra notes. Looking at the 7 octaves, each octave has 12 pitches giving us a total of 7*12=84 pitches on this part of the piano. This is where the 84 from in the n_bins=84 parameter of the Librosa CQT method comes from.
There are multiple reasons why the standard piano has 88 keys instead of 84 that you can read about here.
Hope this helps.


From the description I'll take you have a filter with 3 coefficients. (and the formulae are missing the symbol i). With a sample rate of 100 Hz you would require a filter of length 33 to make it cover the complete wave form of one cycle of 3 Hz signal. Just looking one tenth of a waveform a filter would have to do magic to attenuate (just) a 3Hz signal but leave 1Hz signals intact.

If you'll have octave or matlab, you might want to see the frequency response of your filter:

 `a=exp(i*pi*2/100);freqz([1,a*conj(a)],1)`

Here the 3dB cutoff is approximately at Fs/4 (one quarter of the sampling frequenzy).


I learnt to use HTK based on this very basic tutorial

http://www.info2.uqam.ca/~boukadoum_m/DIC9315/Notes/Markov/HTK_basic_tutorial.pdf

It details how to make a yes/no recogniser. It should get you off to a good start at least.  There is also some HTK code on the main HTK site that allows you to perform real time speech recognition.

To improve robustness, you might also consider training the recogniser on speech tokens that are mixed with the kind of noises that you'd expect to hear around your receiver. You could place a recording device near to the traffic light and record a day's worth of background material, then randomly select samples windows to add to your speech corpus at various levels.


You are not analyzing a sine wave, you are analyzing a windowed sine wave.

$$
y_1(t) = \sin(2\pi f_0t /f_s) w_{\rm rect}(n)
$$
where $f_0 = 510,000$, $f_s = 4,000,000$, and
$$
w_{\rm rect}(n) = \left\{  \begin{array}{cl}
1 & \mbox{where } 0 \le t \le (2^8 - 1)/fs\\
0 & \mbox{otherwise}
 \end{array}\right.
$$

One way to get closer to what you expect is to choose the duration so that your sine wave is "continuous" at each end.  That would mean setting

points = 40000; //2^8

for example.

If I do that, then I get

See this question and its answers for addressing a similar problem.

The reason I chose 40,000 points is to ensure that there is no discontinuity at the ends of the data. This is what happens when you plot [y2 y2] for 40,000 points vs 256 points as in your example.

Notice the discontinuity in the 256 case (second plot).

I agree that choosing a larger number of points is "interpolating the problem away".  There is no real way around that, however.  There will only be 255 zero points in the FFT of a signal of length 256, so if you're expecting to see zeroes everywhere except $f_0$, you're going to be disappointed.

Let's look at the analytic expression for the DFT of a windowed sine wave:
$$
X(\omega) = \sum_{n=0}^{N-1}\sin(2\pi \hat{f}_0 n)e^{-j2\pi nk/N} \\
= \sum_{n=0}^{N-1} \frac{1}{2}\left[e^{+j2\pi \hat{f}_0 n} - e^{-j2\pi \hat{f}_0 n} \right]e^{-j2\pi n\omega/N}\\
= \sum_{n=0}^{N-1} \left [ e^{+j2\pi n ( \hat{f}_0 - \omega )}  - e^{-j2\pi n ( \hat{f}_0 + \omega )}  \right]\\
= \frac{1 - e^{+j2\pi ( \hat{f}_0 - \omega ) N} }{1 - e^{+j2\pi ( \hat{f}_0 - \omega ) } } + \frac{1 - e^{-j2\pi ( \hat{f}_0 + \omega ) N} }{1 - e^{-j2\pi ( \hat{f}_0 + \omega ) } }\\
= \frac{e^{+j\pi ( \hat{f}_0 - \omega ) N}}{e^{+j\pi ( \hat{f}_0 - \omega ) }}
\frac{(e^{-j\pi ( \hat{f}_0 - \omega ) N} - e^{+j\pi ( \hat{f}_0 - \omega ) N})}
{e^{-j\pi ( \hat{f}_0 - \omega )} - e^{+j\pi ( \hat{f}_0 - \omega )} }
+
\frac{e^{-j\pi ( \hat{f}_0 + \omega ) N}}{e^{-j\pi ( \hat{f}_0 + \omega ) }}
\frac{(e^{+j\pi ( \hat{f}_0 + \omega ) N} - e^{-j\pi ( \hat{f}_0 + \omega ) N})}
{e^{+j\pi ( \hat{f}_0 + \omega )} - e^{-j\pi ( \hat{f}_0 + \omega )} }\\
= e^{+j\pi ( \hat{f}_0 - \omega ) (N-1)} \frac{\sin(\pi ( \hat{f}_0 - \omega ) N/2)}{\sin(\pi ( \hat{f}_0 - \omega )/2)}
+
e^{-j\pi ( \hat{f}_0 + \omega ) (N-1)} \frac{\sin(\pi ( \hat{f}_0 + \omega ) N/2)}{\sin(\pi ( \hat{f}_0 + \omega )/2)}
$$
which is where the sinc-like shape comes from.


How to modify zeros and poles in a delta-sigma modulator loop?

i dunno, i might have posted this before.  but here is a 2nd-order sigma-delta quantizer that operates not as an oversampled rate.  so you can hear the quantization noise, but you can also hear the music beneath it.  and the output is binary: $y[n] = \pm A$.
your professor is mistaken if he/she thinks they can model this without the $z^{-1}$ in the feedback path nor leave out the gain of the comparator.
%
%
%
%
%
%
%   simulated 1 bit sigma-delta converter:
%
%
%            x(n)-y(n-1)    w(n)                v(n)                 ( mean(y^2) = A^2 )
%
%   x ---->(+)--->[1/(z-1)]--->(+)--->[1/(z-1)]--->[Quantizer]----.---> y = +/- A = quantized value
%           ^                   ^                                 |
%           |                   |                                 |
%           |                   '----[-fbg]<----.                 |
%           |                                   |                 |
%           '------[-1]<------------------------'------[1/z]<-----'
%
%
%
%
%
%   "linearized" model:
%                                                          .---- q = quantization noise  ( mean(q) = 0 )
%                                                          |
%                                                          |
%            x - y/z        w                   v          |         ( mean(y^2) = G^2*mean(v^2) + mean(q^2) )
%                                                          v
%   x ---->(+)--->[1/(z-1)]--->(+)--->[1/(z-1)]--->[G]--->(+)-----.---> y = G*v + q
%           ^                   ^                                 |
%           |                   |                                 |
%           |                   '----[-fbg]<----.                 |
%           |                                   |                 |
%           '------[-1]<------------------------'------[1/z]<-----'
%
%
%
%
%
%
%
%
%           W = 1/(z-1)*(X - Y/z)
%
%
%           V = 1/(z-1)*(W - fbg*Y/z)
%
%             = (X - Y/z - fbg*Y*(z-1)/z)/(z-1)^2
%
%             = (X*z - Y*(1+fbg*(z-1))) / (z*(z-1)^2)
%
%
%           Y = G*V + Q = G*(X*z - Y*(1+fbg*(z-1)))/(z*(z-1)^2) + Q
%
%             = G*X/(z-1)^2 - G*Y*(1+fbg*(z-1))/(z*(z-1)^2) + Q
%
%
%           Y + G*Y*(1-fbg + fbg*z)/(z*(z-1)^2) = G*X/(z-1)^2 + Q
%
%
%           Y = (G*X/(z-1)^2 + Q)/(1 + G*(1-fbg + fbg*z)/(z*(z-1)^2))
%
%             = (G*X/(z-1)^2 + Q)*(z*(z-1)^2)/((z*(z-1)^2) + G*(1-fbg + fbg*z))
%
%             = z*(G*X + Q*(z-1)^2)/(z^3 - 2*z^2 + (G*fbg+1)*z + G*(1-fbg))
%
%             = z*(G*X + Q*(z-1)^2)/(z*(z-1)^2 + G*fbg*z + G*(1-fbga))
%
%
%    as z -> 1  (DC)
%
%           Y  ->  z*X/(fbg*z + (1-fbg)) =  X/(fbg + (1-fbg)/z)  -->  X
%
%

if ~exist('mean_vv', 'var')
    linearized_model = 0                % run this with 0 the first time to define G and mean(q^2)
end

if ~exist('A', 'var')
    A = 1.0                             % comparator output magnitude
end

if ~exist('fbg', 'var')
    fbg = 2.0                           % feedback gain to internal integrator
end

%
%   if there is an input soundfile specified, use it.  else, create a sin wave
%

if exist('inputFile', 'var')

    [inputBuffer, Fs] = audioread(inputFile);

    fileSize = length(inputBuffer);

    numSamples = 2.^(ceil(log2(fileSize(1))));  % round up to nearest power of 2

    x = zeros(numSamples, 1);                   % zero pad if necessary

    x(1:fileSize) = inputBuffer(:,1);           % if multi-channel, use left channel only

    clear inputBuffer;                          % free this memory
    clear fileSize;

    t = linspace(0.0, (numSamples-1)/Fs, numSamples);   % time

else

    if ~exist('numSamples', 'var')
        numSamples = 65536                              % number of samples in simulation
    end

    if ~exist('Fs', 'var')
        Fs = 44100                                      % (oversampled) sample rate
    end

    if ~exist('f0', 'var')
        f0 = 261.6255653                                % input freq (middle C)
    end

    if ~exist('Amplitude', 'var')
        Amplitude = 0.25                                % input amplitude
    end

    t = linspace(0.0, (numSamples-1)/Fs, numSamples);   % time
    x = Amplitude*cos(2*pi*f0*t);                       % the input

end

sound(x, Fs);                                   % listen to input sound
pause;

y = zeros(1, numSamples);                       % the output (created and initialized for speed later)

if linearized_model
                                                % artificial quantization noise for linearized model
                                                % mean(q) = 0, var(q) = mean(q^2) = mean(y^2) - G^2*mean(v^2)
                                                % does not have to be uniform or triangle p.d.f.
    q = sqrt(6.0*(A^2 - G^2*mean_vv))*( rand(1, numSamples) - rand(1, numSamples) );
else
    q = zeros(1, numSamples);
end

sum_yv = 0.0;
sum_vv = 0.0;

w = 0;
v = 0;
for n = 1:numSamples

    if linearized_model

        y(n) = G*v + q(n);                      % here the comparator is modelled as a little gain with additive noise

    else

        if (v >= 0)                             % the comparator
            y(n) = +A;
        else
            y(n) = -A;
        end

        q(n) = y(n) - (sum_vv+1e-20)/(sum_yv+1e-20)*v;

    end

    sum_yv = sum_yv + y(n)*v;                   % collect some statistics on v
    sum_vv = sum_vv +    v*v;

    v = v + w  - fbg*y(n);                      % second integrator
    w = w + x(n) - y(n);                        % first integrator

end

if ~linearized_model                            % don't recalculate this if using the linearized model
    mean_yv = sum_yv/numSamples;
    mean_vv = sum_vv/numSamples;
    G = mean_yv/mean_vv;                        % the apparent comparator gain (assuming stationary input)
end

%
%
%
%     Y = ((G*z)*X + (z^3 - 2*z^2 + z)*Q) / (z^3 - 2*z^2 + (G*a+1)*z + G*(1-a))
%
%
%
Hx = freqz([0  0 G 0], [1 -2 G*fbg+1 G*(1-fbg)], numSamples/2);
Hq = freqz([1 -2 1 0], [1 -2 G*fbg+1 G*(1-fbg)], numSamples/2);

plot(t, y, 'b');
sound(y, Fs);                                   % this could sound pretty bad
pause;

Y = fft(fftshift(y .* kaiser(numSamples, 5.0)'));
Q = fft(fftshift(q .* kaiser(numSamples, 5.0)'));

f = linspace(0.0, (numSamples/2-1)/numSamples*Fs, numSamples/2);

plot(f, 20*log10(abs(Y(1:numSamples/2)) + 1e-10), 'b');
hold on;
plot(f, 20*log10(abs(Q(1:numSamples/2)) + 1e-10), 'r');
plot(f, 20*log10(abs(Hq) + 1e-10), 'g');
axis([0 Fs/2 -50 100]);
hold off;
pause;

semilogx(f(2:numSamples/2), 20*log10(abs(Y(2:numSamples/2)) + 1e-10), 'b');
hold on;
semilogx(f(2:numSamples/2), 20*log10(abs(Q(2:numSamples/2)) + 1e-10), 'r');
semilogx(f(2:numSamples/2), 20*log10(abs(Hq(2:numSamples/2)) + 1e-10), 'g');
axis([Fs/numSamples Fs/2 -50 100]);
hold off;
pause;

semilogx(f(2:numSamples/2), 20*log10(abs(Y(2:numSamples/2)) + 1e-10), 'b');
hold on;
semilogx(f(2:numSamples/2), 20*log10(abs(Hq(2:numSamples/2)) + 1e-10), 'r');
semilogx(f(2:numSamples/2), 20*log10(abs(Hx(2:numSamples/2)) + 1e-10), 'g');
axis([Fs/numSamples Fs/2 -50 110]);
hold off;


In sentences in normal human speech, words are usually not separated by silent frames or thresholds.  Humans guess (perhaps some form of expectation driven pattern matching) at word separations, and are sometimes wrong.


Downsampling and Gaussian Filtering in the Context of Scale Space Pyramids

You're correct, it has to do with the Cut Off frequency of the Gaussian Blur Filter in its Frequency Domain.

In order to see it, just apply a DFT (Using MATLAB it can be achieved by fft / fft2) and look on the absolute value.
Look for the -3dB point and you'll see.

There is also an intuitive explanation on the original article which say that blurring using Gaussian Kernel with $ \sigma $ means removing any detail which its size is smaller than $ \sqrt{ \sigma } $.


methods of computing fixed point atan2 on FPGA

You can use logarithms to get rid of the division. For $(x, y)$ in the first quadrant:

$$z = \log_2(y)-\log_2(x)\\
\text{atan2}(y, x) = \text{atan}(y/x) = \text{atan}(2^z)$$

Figure 1. Plot of $\text{atan}(2^z)$

You would need to approximate $\text{atan}(2^z)$ in range $-30 < z < 30$ to get your required accuracy of 1E-9. You can take advantage of the symmetry $\text{atan}(2^{-z}) = \frac{\pi}{2}-\text{atan}(2^z)$ or alternatively ensure that $(x, y)$ is in a known octant. To approximate $\log_2(a)$:

$$b = \text{floor}(\log_2(a))\\
c = \frac{a}{2^b}\\
\log_2(a) = b + \log_2(c)$$

$b$ can be calculated by finding the location of the most significant non-zero bit. $c$ can be calculated by a bit shift. You would need to approximate $\log_2(c)$ in range $1 \le c < 2$.

Figure 2. Plot of $\log_2(c)$

For your accuracy requirements, linear interpolation and uniform sampling, $2^{14} + 1 = 16385$ samples of $\log_2(c)$ and $30\times 2^{12} + 1 = 122881$ samples of $\text{atan}(2^z)$ for $0 < z < 30$ should suffice. The latter table is pretty large. With it, the error due to interpolation depends greatly on $z$:

Figure 3. $\text{atan}(2^z)$ approximation largest absolute error for different ranges of $z$ (horizontal axis) for different numbers of samples (32 to 8192) per unit interval of $z$. The largest absolute error for $0 \le z < 1$ (omitted) is slightly less than for $\text{floor}(\log_2(z)) = 0$.

The $\text{atan}(2^z)$ table can be split into multiple subtables that correspond to $0 \le z < 1$ and different $\text{floor}(\log_2(z))$ with $z \ge 1$, which is easy to calculate. The table lengths can be chosen as guided by Fig. 3. The within-subtable index can be calculated by a simple bit string manipulation. For your accuracy requirements the $\text{atan}(2^z)$ subtables will have a total of 29217 samples if you extend the range of $z$ to $0 \le z < 32$ for simplicity.

For later reference, here is the clunky Python script that I used to calculate the approximation errors:

from numpy import *
from math import *
N = 10
M = 20
x = array(range(N + 1))/double(N) + 1
y = empty(N + 1, double)
for i in range(N + 1):
    y[i] = log(x[i], 2)

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y[i] + (y[i + 1] - y[i])*j/M
        if N*M < 1000:
            print str((i*M + j)/double(N*M) + 1) + ' ' + str(a)
        b = log((i*M + j)/double(N*M) + 1, 2)
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

y2 = empty(N + 1, double)
for i in range(1, N):
    y2[i] = -1.0/16.0*y[i-1] + 9.0/8.0*y[i] - 1.0/16.0*y[i+1]

y2[0] = -1.0/16.0*log(-1.0/N + 1, 2) + 9.0/8.0*y[0] - 1.0/16.0*y[1]
y2[N] = -1.0/16.0*y[N-1] + 9.0/8.0*y[N] - 1.0/16.0*log((N+1.0)/N + 1, 2)

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y2[i] + (y2[i + 1] - y2[i])*j/M
        b = log((i*M + j)/double(N*M) + 1, 2)
        if N*M < 1000:
            print a
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

y2[0] = 15.0/16.0*y[0] + 1.0/8.0*y[1] - 1.0/16.0*y[2]
y2[N] = -1.0/16.0*y[N - 2] + 1.0/8.0*y[N - 1] + 15.0/16.0*y[N]

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y2[i] + (y2[i + 1] - y2[i])*j/M
        b = log((i*M + j)/double(N*M) + 1, 2)
        if N*M < 1000:
            print str(a) + ' ' + str(b)
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

P = 32
NN = 13
M = 8
for k in range(NN):
    N = 2**k
    x = array(range(N*P + 1))/double(N)
    y = empty((N*P + 1, NN), double)
    maxErr = zeros(P)
    for i in range(N*P + 1):
        y[i] = atan(2**x[i])

    for i in range(N*P):
        for j in range(M):
            a = y[i] + (y[i + 1] - y[i])*j/M
            b = atan(2**((i*M + j)/double(N*M)))
            err = abs(a - b)
            if (i*M + j > 0 and err > maxErr[int(i/N)]):
                maxErr[int(i/N)] = err

    print N
    for i in range(P):
        print str(i) + " " + str(maxErr[i])

The local maximum error from approximating a function $f(x)$ by linearly interpolating $\hat{f}(x)$ from samples of $f(x)$, taken by uniform sampling with sampling interval $\Delta x$, can be approximated analytically by:

$$\widehat{f}(x) - f(x) \approx (\Delta x)^2\lim_{\Delta x\rightarrow 0}\frac{\frac{f(x) + f(x + \Delta x)}{2} - f(x + \frac{\Delta x}{2})}{(\Delta x)^2} = \frac{(\Delta x)^2 f''(x)}{8},$$

where $f''(x)$ is the second derivative of $f(x)$ and $x$ is at a local maximum of the absolute error. With the above we get the approximations:

$$\widehat{\text{atan}}(2^z) - \text{atan}(2^z) \approx \frac{(\Delta z)^2 2^z(1 - 4^z)\ln(2)^2}{8(4^z + 1)^2},\\
\widehat{\log_2}(a) - \log_2(a) \approx \frac{-(\Delta a)^2}{8 a^2\ln(2)}.$$

Because the functions are concave and the samples match the function, the error is always to one direction. The local maximum absolute error could be halved if the sign of the error was made to alternate back and forth once every sampling interval. With linear interpolation, close to optimal results can be achieved by prefiltering each table by:

$$y[k] = \cases{\begin{array}{rrrrrl}&&b_0x[k]&\negthickspace\negthickspace\negthickspace+ b_1x[k+1]&\negthickspace\negthickspace\negthickspace+ b_2x[k+2]&\text{if } k = 0,\\
&c_1x[k-1]&\negthickspace\negthickspace\negthickspace+ c_0x[k]&\negthickspace\negthickspace\negthickspace+ c_1x[k+1]&&\text{if }0 < k < N,\\
b_2x[k-2]&\negthickspace\negthickspace\negthickspace+ b_1x[k-1]&\negthickspace\negthickspace\negthickspace+ b_0x[k]&&&\text{if } k = N,
\end{array}}$$

where $x$ and $y$ are the original and the filtered table both spanning $0 \le k \le N$ and the weights are $c_0 = \frac{9}{8}, c_1 = -\frac{1}{16}, b_0 = \frac{15}{16}, b_1 = \frac{1}{8}, b_2 = -\frac{1}{16}$. The end conditioning (first and last row in the above equation) reduces error at the ends of the table compared to using samples of the function outside of the table, because the first and the last sample need not be adjusted to reduce the error from interpolation between it and a sample just outside the table. Subtables with different sampling intervals should be prefiltered separately. The values of the weights $c_0, c_1$ were found by minimizing sequentially for increasing exponent $N$ the maximum absolute value  of the approximate error:

$$(\Delta x)^N\lim_{\Delta x\rightarrow 0}\frac{\left(c_1f(x - \Delta x) + c_0f(x) + c_1f(x + \Delta x)\right)(1-a) + \left(c_1f(x) + c_0f(x + \Delta x) + c_1f(x + 2 \Delta x)\right)a - f(x + a\Delta x)}{(\Delta x)^N}
=\left\{\begin{array}{ll}(c_0 + 2c_1 - 1)f(x) &\text{if } N = 0, \bigg| c_1 = \frac{1 - c_0}{2}\\
0&\text{if }N = 1,\\
\frac{1+a-a^2-c_0}{2}(\Delta x)^2 f''(x)&\text{if }N=2, \bigg|c_0 = \frac{9}{8}\end{array}\right.$$

for inter-sample interpolation positions $0 \le a < 1$, with a concave or convex function $f(x)$ (for example $f(x) = e^x$). With those weights solved, the values of the end conditioning weights $b_0, b_1, b_2$ were found by minimizing similarly the maximum absolute value of:

$$(\Delta x)^N\lim_{\Delta x\rightarrow 0}\frac{\left(b_0f(x) + b_1f(x + \Delta x) + b_2f(x + 2 \Delta x)\right)(1-a) + \left(c_1f(x) + c_0f(x + \Delta x) + c_1f(x + 2 \Delta x)\right)a - f(x + a\Delta x)}{(\Delta x)^N}
=\left\{\begin{array}{ll}\left(b_0 + b_1 + b_2 - 1 + a(1 - b_0 - b_1 - b_2)\right)f(x) &\text{if } N = 0, \bigg| b_2 = 1 - b_0 - b_1\\
(a-1)(2b_0+b_1-2)\Delta x f'(x)&\text{if }N = 1,\bigg|b_1=2-2b_0\\
\left(-\frac{1}{2}a^2 + \left(\frac{23}{16} - b_0\right)a + b_0 - 1\right)(\Delta x)^2f''(x)&\text{if }N=2, \bigg|b_0 = \frac{15}{16}\end{array}\right.$$

for $0 \le a < 1$. Use of the prefilter about halves the approximation error and is easier to do than full optimization of the tables.

Figure 4. Approximation error of $\log_2(a)$ from 11 samples, with and without prefilter and with and without end conditioning. Without end conditioning the prefilter has access to values of the function just outside of the table.

This article probably presents a very similar algorithm: R. Gutierrez, V. Torres, and J. Valls, “FPGA-implementation of atan(Y/X) based on logarithmic transformation and LUT-based techniques,” Journal of Systems Architecture, vol. 56, 2010. The abstract says their implementation beats previous CORDIC-based algorithms in speed and LUT-based algorithms in footprint size.


Achievable rate for MIMO channel with estimated channel matrix

I found the following relation which answers my question

$$R=\log\det \left(\boldsymbol{C}_{MMSE}^{-1}\right)$$

where $\boldsymbol{C}_{MMSE}$ is the covariance matrix of the error in the data $\boldsymbol{e}=\boldsymbol{x}-\hat{\boldsymbol{x}}$

I found the answer in the paper "An Iteratively Weighted MMSE Approach to Distributed Sum-Utility Maximization for a MIMO Interfering Broadcast Channel" but it is said to be well known.

Therefore, we need to assume some method for estimation the data using the estimated channel, for example, match-filter or MMSE estimator, and then calculate its error (as a function of the real and estimated channels).


Convolution : Anechoic Speech signal and Impulse response

I have no access to your audio files so I've downloaded:

IR from here (mono/r1_omni.wav) - it's a really long one
Anechoic recording from here (operatic-voice/mono/singing.wav)

Resampled voice signals:

Final convolved signal:

As for your questions:
1.
As you did the plot of IR in logarithmic scale it's clearly visible that towards its end there is another peak at 20s. Once convolving with speech it will yield another repetition of the speech signal. If you look at the first plot, your audio starts around 1s mark. As you are convolving the IR with the input signal, it will repeat it after 20s, and this is what you see in your convolved output around 21s.
Think about how an IR would look like if you want to produce an echo after 20s - one strong peak at 0s and another, smaller at 20s. I am guessing that this "peak" at the end is, in fact, the build-up before the main peak and should actually precede the main one (you can tell it by digital zeros and low energy samples in-between), or potentially a result of harmonic distortion that was moved to the end. Definitely, someone must've extracted/saved this IR in such a way.
Nonetheless, you can easily ignore anything past the 10s mark as it doesn't carry any extra information during convolution. TBH, I would ignore the part after 5s, where it roughly meets the noise floor. Otherwise, it's a waste of computational resources to convolve with full length.
Possibly you can automate the cropping. It looks like the IR file is stored in 16-bit audio, and when the noise-floor ends (around the 10s mark), you get sample values equal to 1. If you look when they start in the vector, you can find your cut-off point automatically, should you have more IR's.
2.
First of all your division by fsIR has no effect as it's normalized by the maximum value in the line that follows. Anyway, most likely your IR is reverberant and generally increases the energy. Think of it as a filter that has a lot of frequency boost. Here is a comparison for resampled anechoic recording and unscaled auralised signal.


Actually, it does not really matter, if a or b are histograms or time domain sequences. However, I suppose you want to convolve to histograms to get the probability distribution of a random variable that is the sum of two variables. Then, its PDF is the convolution of the PDF of the original variables. But, considering the convolution of two time-domain sequences might make it easier to understand.
The convolution theorem states fft(ab)=fft(a)*fft(b), where * denotes the N-point circular convolution and a and b are of length N points. Your statement of the convolution theorem fft(ab) = fft(a) . fft(b) is not correct.
As I said, the convolution theorem in discrete domain corresponds to circular convolution. Most probably, you want to have a linear convolution. Then, you need to zero-pad a and b before applying the convolution theorem.
If the original bins range from 0-3 (left-side of the interval), then the bin edges after convolution would range from 0-6, because the (linear) convolution has a length of (N+N-1) samples, where N is the length of the sequence (which is 4 in your case). You can also understand this from the fact that the convolution of the histogram is the PDF of the sum of both random variables: The minimum sum would be 0, the maximum sum would be 3+3=6.

For more information on circular convolution, you might have a look at an article I wrote or the Wikipedia entry.


You can figure out the filter's transfer function by rewriting the body of the preprocessor macro a bit. Let FILTER_TC(FC, FS) (which really just defines the cutoff frequency of the filter) be equal to k in the below expression.

#define FILTER_LP(FILTER_OUT, FILTER_IN, FC, FS) \
    FILTER_OUT = FILTER_OUT + ((signed long) k * ((signed long) FILTER_IN - (signed long) FILTER_OUT)) >> 16)

In the above, I just expanded the += operator to be a bit less terse. By inspection, we can write the discrete-time difference equation:

$$
y[n] = y[n-1] + k(x[n] - y[n-1]) \\
y[n] = (1-k)y[n-1] + k x[n]
$$

where $y[n]$ and $x[n]$ are equal to FILTER_OUT and FILTER_IN at time index $n$, respectively. This maps to a $z$-domain transfer function of:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac{k}{1 - (1-k)z^{-1}}
$$

This is a lowpass filter, although it is not a Butterworth filter. Its frequency response is plotted below for k = 0.25.

In the above frequency response plot, the blue trace shows the magnitude response and the green trace corresponds to the phase response. It shows that the example you gave isn't a particularly good lowpass filter, even for a first-order filter. You could do better with respect to high-frequency attenuation by using a correct Butterworth structure, for example.

The bit manipulations and explicit casts to signed long types are used because the filter is implemented using fixed-point arithmetic. This is a scheme by which you use integers to represent real numbers using a fixed quantization level. Fixed-point implementations of DSP algorithms were almost ubiquitous in the past, but floating-point processing is becoming more and more common. I would expect to find fixed-point implementations primarily in low-power or embedded contexts today.


Why are linear phase filters called so, if they provide Constant delay instead of linear

Your premise is wrong. Linear phase filters offer linear phase.

Phase is not the same as time delay – at $f_1$, time delay $t$ will lead to a phase of $\varphi_1=-\omega_1\cdot t$, whereas at $f_2$, the phase will be $\varphi_2=-\omega_2\cdot t$ ($\omega$ is $2\pi f$); as you can immediately see, phase is a linear function of frequency.

(It's always good to remind oneself that frequency is just the derivative of phase over time, $\omega=\frac{d\varphi}{dt}$, and thus phase is just an integral of frequency over time.)


This signal looks pretty good, so that a simple sign decision gives you "low" or "high" periods.
Then, you just start at the beginning, and count samples until you find your start bit.
From the start bit, you just sample your signal in 1 bit period, 2 bit periods, ... and get your transmit data stream. Remember, UART is by definition unsynchronized.


Your code is correct and the results are consistent. You may be surprised by them due to some 'hidden features'.

First,  conv2 returns by default the full convolution, such that the result is the size of the image plus a border of half the size of the kernel (that is, the total size if the size of the image plus that of the kernel). When you interpret your results, be aware of it!

Second, the results represent coefficients that are stronger for a higher correlation between your kernel and your local image patch: as expected you also extract the borders of the image. See in particular your leftmost result showing strong vertical line.

Last, imagesc scales by default the scale between the highest to the lowest coefficient. That's why in the leftmost result, you mainly see the border.

There are different options to conv2described in help conv2 which allow to control this behavior.

Be aware also that there many different definitions of kernels for detecting edges, such as log-Gabors

.

If you are interested in a full implementation (in python) you may have a look at:  https://pythonhosted.org/LogGabor/ (shameless self-plug ).


What are the limitations of FIR filters?

Note that this is not about general limitations of FIR filters, but about the special case of linear-phase FIR filters. If you understand why a type-II linear-phase FIR filter has a zero at $z=-1$, then the limitations of the other types should be obvious too.

It's always about zeros at either $z=1$ (DC) or $z=-1$ (Nyquist). Given the transfer function

$$H(z)=\sum_nh[n]z^{-n}\tag{1}$$

it is easy to see that

$$H(1)=\sum_nh[n]\tag{2}$$

and

$$H(-1)=\sum_n(-1)^nh[n]\tag{3}$$

Now look at examples of the four linear-phase FIR filter types. An example of an impulse response of a type-I filter is

h = [1 -1 2 3 2 -1 1]

Using $(2)$ and $(3)$ can you say anything about zeros at $z=1$ or $z=-1$?

An example of a type-II filter is

h = [1 -1 2 2 -1 1]

Clearly, applying Eq. $(3)$ results in a zero at $z=-1$, regardless of the actual values of the impulse response.

The following two filters are type-III and type-IV filters, respectively:

h = [1 -1 2 0 -2 1 -1]

h = [1 -1 2 -2 1 -1]

Now apply Eqs $(2)$ and $(3)$ and see what you get.

This answer contains more information on the four types of linear-phase FIR filters.


Intuitive meaning of complex time domain signal representation (quadrature sampling)

One common use of the IQ data from an SDR is to feed them to an FFT.  The upper and lower halves of an FFT result indicated different relationship between the (real) cosine and (imaginary) sine components of the complex FFT result. The direction of rotation of each complex component tells the rest of the signal chain whether the complex component is an upper or lower sideband of the assumed center tuned frequency of the IQ SDR data stream.

The quadrature component in the IQ input is what allows the two halves of an FFT result to not just be redundant complex conjugates of each other (which they would be given a strictly real input, no Q or zero Q).  Thus allowing separating the upper and lower sidebands as independent spectrum.

Also see: https://electronics.stackexchange.com/questions/39796/can-somebody-explain-what-iq-quadrature-means-in-terms-of-sdr/99617#99617


The half-band lowpass filter, has the following impulse response :
$$ h[n] = \frac{\sin(\frac{\pi}{2} n) }{ \pi n} $$

And for $x[n] = \delta[n]$ ,  $w[n] = h[n]$.

Then the output $y[n]$ will be :

$$y[n] = w[n] \big( 1 + (-1)^n \big) $$

which can be simplified by expanding the parenthesis as:

$$
\begin{align}
y[n] &= \frac{\sin(\frac{\pi}{2} n) }{ \pi n} \cdot \big( e^{j 2 \pi n} + e^{j \pi n} \big) \\ \\
y[n] &= \frac{\sin(\frac{\pi}{2} n) }{ \pi n} \cdot  e^{j \frac{3 \pi}{2} n} \cdot  \big( e^{j \frac{\pi}{2} n} + e^{-j \frac{\pi}{2} n} \big)\\ \\
y[n] &= \frac{ \sin(\frac{\pi}{2} n) }{ \pi n} \cdot 2 \cdot  \cos( \frac{\pi}{2} n) \cdot  e^{j \frac{3 \pi}{2}  n} \\ \\
y[n] &= \frac{\sin(\pi n) }{ \pi n}  \cdot  e^{j \frac{3 \pi}{2} n}
\end{align}
$$
Which can be shown to be :

$$y[n] =  \begin{cases}{ 1 ~~~~, n = 0 \\ 0 ~~~~ , n \neq 0  } \end{cases} $$

Therefore $y[n] = \delta[n]$. The system is a zero-phase all pass filter. And finally, the required quantity is
$$ M = \sum_{n=-\infty}^{\infty} y[n] = \sum_{n=-\infty}^{\infty} \delta[n] = 1 $$


Do we assume modulating and carrier signal uncorrelated in AM modulation?

The total power of an amplitude modulated signal is
$$\begin{align}\overline{s^2_{AM}(t)}&=\overline{\big(A+m(t)\big)^2\cos^2(2\pi f_ct)}\\&=\frac12\overline{\big(A+m(t)\big)^2}+\frac12\overline{\big(A+m(t)\big)^2\cos(4\pi f_ct)}\tag{1}\end{align}$$
The second term on the right-hand side of $(1)$ is zero if $m(t)$ is a lowpass signal, and if $f_c$ is sufficiently large (i.e., greater than the highest frequency in $m(t)$).
Consequently, we have
$$\begin{align}\overline{s^2_{AM}(t)}&=\frac12\overline{\big(A+m(t)\big)^2}\\&=\frac12\big(A^2+2A\cdot\overline{m(t)}+\overline{m^2(t)}\big)\tag{2}\end{align}$$
The term $\overline{m(t)}$ is the DC value of the message signal, which is assumed to be zero. With this assumption we finally arrive at
$$\overline{s^2_{AM}(t)}=\frac12\big(A^2+\overline{m^2(t)}\big)\tag{3}$$
In sum, the assumptions for which $(3)$ is valid are that $m(t)$ is a lowpass signal, that the carrier frequency is sufficiently high, and that $m(t)$ has a DC value equal to zero. It is usually straightforward to satisfy all of these assumptions to a very high degree of accuracy.


Manipulating RGB color by copying channel value

I missed the conversion part.
I converted the result image from BGR into RGB using this command :

cvtColor(image, image, COLOR_BGR2RGB);

Then I exported the result image into image file.

imwrite("image.png", image);

Turned out the export command will automatically convert the BGR color-space into the RGB color-space. This is why the color result is off because I did 2x conversion.
So, I don't need the conversion command.


Does windowing affect Parseval's theorem?

Parseval's theorem will hold, but take into account that your signal in the time domain will no longer be $x[n]$. Namely, if you have that

$$\sum_{n=0}^{N-1} \Big| x[n] \Big|^2  = \frac{1}{N}  \sum_{k=0}^{N-1} \Big| X[k] \Big|^2$$

then, if you window the signal $x[n]$ with a window $w[n]$, your signal will now be $\hat{x}[n]=x[n]w[n]$, and the theorem will now return the following equality:

$$\sum_{n=0}^{N-1} \Big| \hat{x}[n] \Big|^2  = \frac{1}{N}  \sum_{k=0}^{N-1} \Big| \hat{X}[k] \Big|^2$$

where $\hat{X}[k]$ is the DFT of $\hat{x}[n]$ (and obviously differs from $X[k]$).


What's the point of defining the signal over the whole time domain?

Bottom line: it makes the math easy.  In addition, there are certainly instances where we care that a signal has existed for some time in the past.

But how does it even make sense to write a negative value for time,

Think of $t=0$ as "now", $t < 0$ as "before now", and $t > 0$ as "later".  Does that make sense?
By extension, $t = 0$ is just any arbitrary moment in time; if you don't like it, shift everything over.

shouldn't any "real world" signal start at some instant

Aside from the big bang, no, why should there be some instant where a signal "starts"?
Note that the concept of signals that exist for all time comes from Fourier analysis.  If you have some reason to start considering signals at an instant, there is Laplace analysis, where time pretty much starts at $t = 0$.

does negative time do anything other than make us refrain from writing too many piece-wise defined functions?

Fourier analysis is a good way to analyze signals that are "just there", like radio signals or speech.  If your analysis revolves around a system that's in steady state, processing a signal that has had pretty much the same statistics that it has now and will have in the future, then treating that signal as existing for all time is both sensible, and makes the math easy.


Why does FEC seem to make BER worse at low Eb/N0?

This is a behaviour that is commonly outlined in some textbooks and tutorials on FEC, but usually in the form of an observation. For instance, Turbo Codes are sometimes characterised by their Threshold $E_b/N_o$, which can be considered as the minimum SNR per bit from where the BER starts to decrease. I have not seen a formal explanation for this property, but it can be deduced at least for some types of codes I assume. In this Wikipedia article, for instance, it is explained as a consequence of the risk-polling effect that comes from the noise averaging capability of FEC. Let me complement this, also by an intuitive approach.

Consider for instance the use of a $(n,k)$ binary block code with hard decision decoding. The maximum number of correctable bit errors for this code is commonly denoted by $t$.
If the number of bit errors within each transmitted block is higher than $t$, then the receiver will not be able to recover the original codeword. Moreover, the received block will be at roughly the same distance of multiple codewords within the set of size $2^k$. The decoder then picks one of those codewords with high uncertainty, which results in a BER close to 0.5.

Simply put, a few bit errors - though above the error correction capability $t$ of the code - may result in a complete block error, which explains why FEC underperforms the uncoded system for a certain SNR range.

This is why in systems that operate at low SNR or when errors appear in bursts, it is a good idea to use interleaving in order to spread the bit errors along several blocks and keep the number of bit errors per block below $t$.


Why can we have non-integer frequency bins in FFT?

Frequencies in DFT, $\omega = \frac{2\pi}{N}k, \ k = 0,1,2,...,N-1 $, only depends on the Length $N$ of DFT, and nothing else.

$$X[k] = \sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N} nk}, \qquad k = 0,1,2,...,N-1$$

It is very important to understand this expression as the projection of time-domain finite length sequence $x[n], n=0,1,2,3,...,N-1$, of length $N$, onto orthogonal Discrete Fourier Basis Vectors $W_k^{N} = e^{j\frac{2\pi}{N}nk}, k=0,1,2,...,N-1$. And, the DFT Coefficients $X[k]$ are the projection coefficients actually on the $k^{th}$  basis vector. So, with this idea in mind, we can say $x[n]$ is actually a vector $\vec{x}$ in an $N$-Dimensional vector space. And, this $N$-D vector space has a set on $N$ Fourier Basis which reveals the frequency information of that $\vec{x}$. We know that projection of a vector onto another depends on dot-product (inner-product), $\langle\vec{x},\vec{y} \rangle$. The projection coefficient of a vector $\vec{x}$ onto a vector $\vec{y}$ :

$$P_{\vec{x}} = \frac{\langle\vec{x}, \vec{y}\rangle}{\langle\vec{y}, \vec{y}\rangle} = \frac{\sum^{N-1}_{n=0}x[n]y^{*}[n]}{\sum^{N-1}_{n=0}y[n]y^{*}[n]}$$
And, if we are projecting a vector $\vec{x}$ onto a unit length vector $\hat{y}$, then the Projection coefficient is exactly the dot-product, because, the denominator of the above expression will become 1.

$$P_{\vec{x}} = \langle\vec{x}, \hat{y}\rangle$$
That is exactly what the DFT Coefficient is. Projection Coefficient of time-domain sequence $\vec{x}$ onto $\vec{W_k^N}$
$$X[k] = \langle\vec{x}, \vec{W_{k}^{N}}\rangle = \sum^{N-1}_{n=0}x[n]\,W_k^{N}[n] = \sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N} nk}$$

Now, notice that the $N$ Fourier Basis vectors are each at a digital frequency $\omega = \frac{2\pi}{N}k, k=0,1,2,...,N-1$. So, the DFT/FFT frequency bins just depend on the DFT Length and nothing else. DFT by itself is a perfectly fine operation on an $N$-dimensional vector. It is a linear transformation, just a change of basis from canonical basis set to Fourier Basis set, so as to reveal frequency information in that signal $\vec{x}$.

It reveals frequency information because projection(inner-product) of one vector onto another actually is a measure of similarity between the two vectors. So, each of the $X[k]$ is telling us how much similarity is between the time-domain sequence $\vec{x}$ and $k^{th}$ Fourier Basis, and, $k^{th}$ Fourier Basis is a sequence of digital frequency $\omega = \frac{2\pi}{N}k$.

As I said, DFT by itself is perfectly fine operation on any vector of finite length. You do not have to link these digital frequencies with sampling frequency. But if you want to, then, it is enough to know that when we sample a signal at sampling frequency $f_s$, then the maximum frequency that can be represented in the sampled signal is $f_s/2$. And the maximum digital frequency is $\pi$. Hence, $f_s/2$ maps to $\pi$. And therefore the digital frequency of $k^{th}$ Fourier Basis, $\omega = \frac{2\pi}{N}k$ maps to $f = \frac{f_s}{N}k$.


Designing appropriate filter by choosing values for $A$ and $B$

Bottom Line:
$$A <1$$
$$B =2\cos(\omega_n)$$
Where $\omega_n$ is the normalized angular frequency of the desired notch location (in this case for the OP with a sampling rate of 8KHz and notch at 400 Hz this would be $\omega_n = 2\pi400/8000 = \pi/10$), resulting in $B \approx 1.902$ and $A$ is the frequency notch bandwidth parameter; the closer $A$ is to 1 the tighter the bandwidth of the notch.
Details:
The straight-forward approach to eliminating narrow band noise at 400 Hz would be a classic 2nd order notch filter. This is detailed at this post:
Transfer function of second order notch filter
This approach is done by simply placing a zero at the frequency of interest, and by adding a pole close to the zero, but inside the unit circle of course to be stable, we can adjust the bandwidth of the notch by the proximity of this pole to that zero, meaning $|p|<1$, but close to 1.
Using real coefficients, meaning complex conjugate zeros and poles as detailed in the linked post, this results in the transfer function:
$$ H(z) = \frac{1+a}{2}\frac{z^2-2z\cos(\omega_n)+1}{(z^2-2az\cos(\omega_n)+a^2)} $$
Which as a difference equation is exactly in the form that the OP has used and therefore given this is the same filter, answers the OP's question. This is is derived from the z-transform as follows (we can ignore the gain scaling parameter $(1+a)/2$ which for most cases is close to 1, and express in decreasing powers of z (and assuming I didn't make simple algebra errors!):
$$ \frac{Y(z)}{X(z)} = \frac{1-2z^{-1}\cos(\omega_n)+z^{-2}}{(1-2az^{-1}\cos(\omega_n)+a^2z^{-2})} $$
$$ Y(z) (1-2az^{-1}\cos(\omega_n)+a^2z^{-2})= X(z)(1-2z^{-1}\cos(\omega_n)+z^{-2})$$
$$Y(z)- Y(z)2a\cos(\omega_n)z^{-1}+Y(z)a^2z^{-2} = X(z)-X(z)2\cos(\omega_n) z^{-1} +X(z)z^{-2} $$
And then with the inverse Z transform we get the desired coefficients that the OP wanted:
$$y[n] - 2a\cos(\omega_n)y[n-1] + a^2y[n-2] = x[n]-2\cos(\omega_n)x[n-1] + x[n-2] $$
With comparison to OP's expression:
$$y[n] - ABy[n-1] +A^2y[n-2] = x[n] - Bx[n-1] +x[n-2]$$
Results in the following:
$$A = a$$
$$B =2\cos(\omega_n)$$
Where $\omega_n$ is the normalized angular frequency of the desired notch (for the OP's sampling rate and notch frequency this would be $\omega_n = 2\pi400/8000 = \pi/10$) and $a$ is the notch bandwidth parameter ($a<1$, the closer a is to 1 the tighter the notch, refer to the linked post for complete details).


Discussion of simulation results of a signal sequence: amplitude

Phase offset and time delay (time offset) are NOT the same and I think this may be a primary source of confusion (Given the OP's other recent posts dealing only with time offset but here the formula used introduces a carrier phase offset). Hopefully the following diagrams help distinguish the two.
No Time or Phase Error
Below is the constellation diagram and eye diagram for the (raised cosine) QPSK waveform in the receiver after the matched filter, with no phase or time offset and showing the decision locations in red:

The eye diagram for the real component of the signal is shown below and should be clear that the imaginary component would look identical:
The blue circles are samples with 4 samples per symbol, one of which is at the ideal sampling location (at sample number 3 and 7 specifically; this eye diagram spans over two symbols).

TIME OFFSET
Next a small time offset (timing error) of about 1/3 of a sample is introduced to show the effect of time offset alone.

PHASE OFFSET
With the time offset back to 0, a phase offset of $\pi/10$ is introduced showing it's effect on the constellation and eye diagram.

Conclusion
Phase Offset and Time Offset are two different parameters that each independently need to be corrected (one can exist without the other). Further, phase offset is typically changing with time, which specifically then when changing at a linear rate is a static carrier offset since frequency is the derivative of phase with time, and is corrected with a carrier recovery loop (while time offset is corrected with a timing recovery loop).


Correctly scaling FFT of different lengths

Currently, I find my impulse response by performing an FFT of size 32k
  on my entire measurement signal, and my reference signal. I divide
  each by the scaling factor 1/sqrt(M)

Scaling here is unnecessary since you are going to divide the two spectra anyway. Any scale factor will cancel out.

  I then divide the two signals, and perform an inverse FFT, and then scale by dividing by sqrt(M).

Scaling here is wrong. Dividing the spectra creates a transfer function and the impulse response is simply the inverse FFT of the transfer function, regardless of FFT size (unless it creates time domain aliasing)

  I then once again perform an FFT on each window, and again scale with the 1/sqrt(M).

Wrong again. Don't scale FFTs for Transfer/Function<->Impulse response relationships

  And then get my magnitudes using: i = sqrt (re/im) mag = 20 * log10(i/M)

This makes no sense to me. You seem to divide the real part by the imaginary part. That could potentially be part of a phase calculation but has nothing to do with the magnitude. If you have $H(\omega) = a+ j \cdot b$ then your level in dB is simply $L(\omega) = 10 \cdot log_{10}(a^2+b^2)$

  I am pretty certain that the issue is to do with having not used any scaling factors for my windows used.

This has nothing to do with the windows. Stop scaling anything at all, and this should just work.


With slight modification you might want to use RootSift:

http://www.robots.ox.ac.uk/~vgg/publications/2012/Arandjelovic12/arandjelovic12.pdf

Also the other steps in the paper will guide to improve the recall rate.

Cheers,


If noise is your signal, what is your noise?

There is a derivation over on the math.SE site that might be useful.

Provided your data is Gaussian with variance $\sigma^2$.

If you use the unbiased estimator of the population variance:
$$
S^2_{N-1} = \frac{1}{N-1} \sum_{n=0}^{N-1} (x[n]- \bar{x})^2
$$
where $N$ is the number of points in your batch and $\bar{x}$ is the sample mean, then the variance of this estimate will be:
$$
\frac{2 \sigma^4}{N-1}
$$


confusion related to finding inverse Z-Transform using Complex Integral Method

The closed contour $C$ must lie inside the region of convergence, so for the ROC $|z|<2$ you have no poles inside $C$ for $n\ge 0$, hence $x[n]=0$ for $n\ge 0$.


If you know (or better: understand) the basic $\mathcal{Z}$-transform relation

$$x[n-k]\Longleftrightarrow z^{-k}X(z)$$

then from the given time domain equation you immediately get

$$Y(z)=X(z)+2X(z)z^{-1}=X(z)(1+2z^{-1})$$

Since the transfer function $H(z)$ is the ratio of $Y(z)$ and $X(z)$, you get

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


You can just write $(2)$ in $y_t$, so it becomes

$$y_t = \begin{bmatrix} b_1-a_1 && -a_2\end{bmatrix} \begin{bmatrix} x_{a,t}\\ x_{b,t}\end{bmatrix} +\begin{bmatrix} 1 \end{bmatrix}u_t $$

Actually, this representation is called controllable canonical form. However, your system is not strictly proper, so you need to split it as

$$G(z)=\frac{z^2 + b_1 z}{z^2 + a_1 z + a_2} = 1 + \frac{(b_1 - a_1) z - a_2}{z^2 + a_1 z + a_2}$$

See how the coefficients appear?


Noise reduction using multiple recordings of the same signal

Is it possible to reconstruct the original pure signal?

No, that is information-theoretical impossible. Also, that signal doesn't exist, probably, to begin with ;)

However, you can definitely increase the the SNR simply by averaging; that becomes pretty obvious when you consider the signal of interest to be correlated within your recording, whereas your noise is claimed to be independent.

If you sum up multiple recordings, the correlated signal's amplitudes will increase linearly with the number of observations, whereas your noises' standard variance only with the square number of that.

  I have two (audio) recordings of one analog signal, I used two separate devices (smart-phones) recording from two different locations.

Find a signal model for your analog signal – it's not just "any signal" (that would be completely indistinguishable from noise), but it's probably very specific. Ask a new question – "I have these and that signal, they follow my signal model (maybe you need help formulating one, another good question), how do find a low-variance estimator for the parameters for that signal model".


How to evaluate an impulse/response model in MATLAB

One way of evaluating the model response to arbitrary data is the lsim command. As you don't provide more information with respect to your data and/or model an example below (taken straight from the MATLAB help). See help lsim or doc lsim for more information.

% create model
%
num = [1 4];
den = [1 20 5];
G = idtf(num,den);

% evaluate response
%
t = 0:0.01:5;
u = sin(t);
y = lsim(G,u,t);

% plot
%
figure;
plot(t,y)

In addition, if desired, you can evaluate the model impulse, step and/or frequency response using the impulse, step and bode commands.


Does a linear phase FIR filter shifted in frequency preserve linear phase?

The generalized linear phase FIR filter has the following frequency response:

$$H(\omega) = A(\omega)~e^{j (\alpha \omega + \beta)}$$ for some constant $\alpha$ and $\beta$ and $A(\omega)$ being real. The group delay of this filter is independent of frequency: $$\tau = -\frac{d (\alpha \omega + \beta) } {d \omega} = - \alpha $$

Now if you linearly shift this filters frequency response then you get the new filter's frequency response as :

$$H_2(\omega) = H(\omega-\omega_0) = A(\omega-\omega_0)~e^{j(\alpha (\omega-\omega_0) + \beta)} = B(\omega)~e^{ j ( \alpha \omega - \alpha \omega_0 + \beta}) $$

which also shows the generalized linear phase property. The new group delay is:

$$\tau  = -\frac{d (\alpha \omega - \alpha \omega_0 + \beta) } {d \omega} = - \alpha $$

the same as the previous filter. And is also independent of frequency $\omega$.

Recall to page 295 of Discrete-Time Signal Processing 2E by A.Oppenheim for a more detailed discussion of the generalized linear phase systems.


Correct way to quantize accelerometer data with microsecond precision?

useless ramblings first

As an SDR person: 12h · 200 S/s = 12·3600s·200S/s = 8.64 MS is not a "large volume of data" (we usually deal with a multiple of that per second), so rest assured, most of your computation should probably be pretty fast to complete, which is good news!

Relating your Jitter to your sampling period

Now, with 200 Hz, your sample period is 5000 µs, so a 80 µs variance (something doesn't add up with the units here, it should be s², not s) would be an error of ca 1.6%, or -18dB. If you assume sufficient band-limitation, rule of thumb says that the resulting phase noise would be in the same order of magnitude, ie. -18dB.

It's generally a good idea to define first what noise level are you aiming for; otherwise you'll be optimizing something that is already good enough.

Of course, you could properly resample, with much less distortion, but frankly, accelerometer data doesn't sound like you necessarily have great SNR in the first place (but I might be totally wrong about that!).

Interpolator considerations

What you need, if you actually want to generate a perfectly periodically sampled signal, is an interpolator.

What confuses me is that especially jerk is a derivative measure, i.e. you'd be finding an interpolation of your signal more or less anyway; so, where's the problem in using a standard interpolation method (scipy comes with some) to get perfectly periodic samples, should you not choose to implement your own interpolating differentiator (hint: splines are of polynomial shape and that makes differentiation especially easy, so you can directly calculate a derivative from your data/time points, and be done with it, without first calculating a periodically-sampled signal representation).

Bandwidth considerations

If you know that your signal of interest is contained in frequencies significantly below 100 Hz, then consider this:

Assume your signal is a sum of individually delayed, scaled sines of different frequencies.

The point of maximum derivative, i.e. maximum steepness, is the point where our sampling time offset $\Delta t$ hurts the most.

A single sine has a maximum derivative of 1 (proof: the derivative of a sine is a cosine, which has a max of 1).

If each of your $N$ sines $s_i$ is scaled with an amplitude factor $a_i$ and has frequency $f_i$ (and phase==delay $\varphi_i$, which we do not care about right now), then its formula is

$$s_i = a_i \sin(2\pi f_i t + \varphi_i)$$

and its time derivative

$$\dot s_i = a_i\cdot 2\pi f_i \cdot \cos(2\pi f_i t + \varphi_i)$$

has a maximum of

$$\max \dot s_i = 2\pi a_i f_i = (\max s_i) \cdot 2\pi f_i \text.$$

In a linear approximation (and that's a prett good approximation for a range of less than 1.6%/2 = 0.8% of a cycle of a cosine around it's zero crossing; the non-scaled cosine is pretty much a 45° line there), the (defined-by-me) relative error (relation of error magnitude to amplitude of sine) you get is

$$\begin{align}
\max |e_i| &:= \max \left\lvert\frac{\dot s_i \cdot \Delta t}{a_i}\right\rvert\\
& = 2\pi f_i \Delta t\text,
\end{align}$$

which means that the lower the sine's frequency, the less you're sensitive to jitter.

Let's look at the case where you only care about frequency up to 50 Hz, or 5 Hz:

$$
\begin{align}
\max |e_i| &= 2\pi f_i \Delta t&f_i = 50\,\text{Hz}\\
&= 2\pi\, 5\cdot10^1\, \cdot\, 8\cdot10^{-5}\\
&= 2\pi\, 4\cdot10^{-3}\\
&\approx 2.5 \%\\[4em]
\max |e_l| &= 2\pi f_l \Delta t&f_l = 5\,\text{Hz}\\
&\approx 0.25 \%
\end{align}
$$

So, the lower your signal of interest's frequency is, the less you need to actually care about the jitter!


Downsampling audio for use in Machine Learning

Is this necessary if I only intend to use the data in the Neural Network toolkit provided by the repo I linked?

Yes.

Whether or not you are downsampling (instead of just decimating) has nothing to do with classification performance but rather, it is to preserve (as much as possible) the information  contained in the signal.

When changing the sample rate of a signal, you need to account for aliasing. Otherwise, the signal will be distorted in such a way that it would be impossible to reconstruct the original signal from the samples.

  If so, is there an industry standard method for changing sample rate?

Yes.

The "right" way to do it when downsampling is to first apply an anti-aliasing filter and then decimate the signal but when upsampling, you first upsample and then apply interpolation (which can also be expressed as a filter). Various platforms provide functions to do just that (e.g. Python, MATLAB).

  I realise that it probably depends on what features are being used.

To some extent. Perhaps you could simply leave everything unchanged, provided that you adjust certain points in the algorithm that are specifically hardwired for a sampling frequency that is different than the one you are using. However, the bandwidth for speech is limited. If your classifier is also trained on what is NOT speech, having more bandwidth means that you can discriminate more things. If the majority of your dataset is pure speech then having more bandwidth would not make a huge difference.

  Can I use my 48khz data - downsampled with no filtering - with the expectation that the model will work for prediction on real 16khz data?

No.

If you downsample without filtering you are very likely to produce distorted signals.

  And then for future readers sake, how about the other way? Say I had an 8khz file, could I increase the sample rate to 16khz without filtering?

No.

When you increase the sample rate, you have a more serious problem. You now have to fill "gaps" for which you do not have data. The filtering (or interpolation) stage fills in those gaps but does not add any new information. When you downsample, you are throwing away information. Upsampling cannot replace it.

Perhaps it would be beneficial for you to look up some basics of DSP such as sampling, sampling and reconstruction, aliasing and anti-aliasing filters. Also, to further understand what your features look like, how do they work and how to choose and use them, you might want to look up more information on the Short Time Fourier Transform and Multirate Signal Processing (more generally).

Hope this helps.


What should be the cut off frequency of a Filter operating on an OFDM signal?

You don't need that filter. Your previous questions still indicate you shouldn't be using a filter. OFDM is a filterbank in itself, and applying a filter beforehand will increase delay spread, and thus you get ISI¹.

Sure, with a filter, you'll save on the noise power you'd get in carrier sidelobes, but that should typically be lower than what you lose in ISI.

So, no, you don't apply a different linear filter. You simply ignore the FFT bins that don't contain signal.

¹ you'll get ISI under the assumption that your OFDM system is designed with the minimum amount of subcarriers needed to compensate the maximum expectable delay spread on the channel. That makes a lot of sense because using significantly more subcarriers has no advantage, but makes your system depend on a longer coherency time. If you convolve the channel output with your receive filters, you get a "convolved channel-filter", with a length that is thus longer than your channel, and thus worse than your OFDM system can deal with.


Optimal trade-off between oversampling and filter length

1-d MMSE derivative filter for uniform spectral prior

With some simplifying assumptions about the signal distribution, a filter with least square frequency response error is the Bayesian minimum mean square error (MMSE) filter for a uniform prior of the frequency spectrum. For such 1-d filters, an oversampling factor $\beta \approx 2\ldots2.5$ seems optimal, as it is in the sweet spot of the cost-benefit curves (Fig. 1).

Figure 1. Mean square frequency response error of least square optimal derivative filters as function of computational complexity as measured by $\beta N$ which is proportional to the total number of multiply-accumulates (MACs), for various oversampling factors $\beta$ for which the filters have been optimized. Also shown are traces for different numbers $N$ of non-zero filter coefficients. No advantage was taken from impulse response anti-symmetry or simple coefficients such as $1$.

More data is needed before general conclusions can be drawn.