Confusion regarding sinusoidal synthesis formula,

Simple:

Euler's formula:

$$ e^{i\theta} = \cos(\theta) + i \sin(\theta) $$

and

$$ e^{-i\theta} = \cos(\theta) - i \sin(\theta) $$

Add them together:

$$ e^{i\theta} + e^{-i\theta} = 2 \cos(\theta) $$

You should see that in your equation.

From there it is usually stated:

$$ \cos(\theta) = \frac{ e^{i\theta} + e^{-i\theta}}{ 2 } $$

See my article

The Exponential Nature of the Complex Unit Circle

for a more expansive explanation.


First, please read this answer of mine for a detailed description of matched filters for real-valued signals. In particular, note that what I called the matched filter for a signal $x(t)$ is a(n LTI) filter with
impulse response $h(t) = x(-t)$ which is better described as the time-reversed signal rather than the "inverse" of the impulse response as you call it: to many people, inverse would suggest $[x(t)]^{-1}$ or $-x(t)$
depending on whether they were thinking of multiplicative or additive inverse.. Note also that convolving the channel output with the
time-reversed channel impulse response does not "cancel" the effect of the channel on the signal at the channel input. In general,
matched filtering distorts the signal being filtered tremendously because
the goal of matched filtering is to produce a peak response at the desired sampling time. This is quite different from the audio enthusiast's desire of
reproducing the input signal with the highest possible fidelity, and
thus using an "inverse" filter to remove all traces of the distortion
caused by the channel filter.  See
the last part of my answer cited above for some pictorial illustrations of how matched filtering distorts the signal so as create a peak response
at the sampling instant.

Second, there are no complex signal in real life; they exist only in the diseased imaginations of DSP engineers.  What we have in real life
is a real-valued bandpass signal
$$x(t) = x_i(t)\cos(2\pi f_0 t) - x_q(t)\sin(2\pi f_0 t)
= \Re\left(x_c(t) e^{j2\pi f_0t}\right)$$
where $x_i(t)$ and $x_q(t)$ are slowly varying signals compared to
the carrier signal of frequency $f_0$ Hz, and
$x_c(t) = x_i(t) + jx_q(t)$ is called the complex baseband equivalent
of the RF signal $x(t)$. The filter matched to $x(t)$ thus has
impulse response
\begin{align}
x(-t) &=  x_i(-t)\cos(-2\pi f_0 t) - x_q(-t)\sin(-2\pi f_0 t)\\
&= x_i(-t)\cos(2\pi f_0 t) - [-x_q(-t)]\sin(2\pi f_0 t)\\
&= \Re\left([x_i(-t) -j x_q(-t)]e^{j2\pi f_0t}\right)\\
&= \Re\left([x_c(-t)]^* e^{j2\pi f_0t}\right)
\end{align}
The complex baseband equivalent of the matched filter impulse
response is thus $[x_c(-t)]^*$ and thus when the DSP engineer
proceeds to match-filter the complex baseband equivalent input $x_c(t)$,
she must convolve $x_c(t)$ not with $x_c(-t)$ but rather with
$[x_c(-t)]^*$, cos doing the former instead of the latter
would result in her being in a state of sin.


Explanation relevant to the filter response ($\mathrm{sinc}$ interpolation) - Equation $(2)$ in my question

Eq$(1)$ doesn't look right. It should be:

$$x'(t) = \sum_{n=-\infty}^{+\infty} x(nT_s) \mathrm{sinc}(t-nT_s)$$

For your doubt, consider subsampling (lowering $F_s$) as equivalent to reconstructing $x(t)$ and then sampling it at a lower frequency.

Eq$(2)$ comes from the combination of $2$ lowpass filters:

The ideal reconstruction filter (from eq.$1$)
a new anti-alias filter previos to sampling at the lower $F_s$

These $2$ filters combined should result in eq$(2)$ (several lowpass filters in cascade is equivalent to one filter with the minimum bandwidth, and the product of the gains).


Does perceived ADC resolution, range and precision depend on fixed point format? What is the best way to deal with this for signal processing?

it may turn out to be very inconvenient, but you can do essentially anything in fixed-point that you can do in floating.  as long as you got "Add-with-carry" as an instruction and enough MIPS to do all of the multiword instructions you may need.  the 12 bits that come from your ADC don't care if it's fixed or floating and your ADC resolution is the same regardless.  can't promise that the quality of code coming out of a C compiler will be any good.


The reason is the decorrelation due to the varying phase shift between the input and output. See the plot below showing the phase of the OP's high pass filter, and how closely it follows the cancellation when plotted on a log log scale (the higher frequency components of the two signals are better aligned in time and therefore have higher correlation leading to better cancellation).

The correlation between signals vs phase is proportional to the sine of the phase (signals in quadrature are uncorrelated consistent $\sin(\pi/2)=0$. For small angles the sine of the angle is the angle in radians hence we see the close match between the response of the phase and the cancellation of your signal as in your plots.

Due to the small angle, the OP's cancellation signal, assuming perfect amplitude matching, will have a magnitude that is directly proportional to the phase itself in radians between the two signals. Consider the cancellation of two signals with identical magnitude but a difference in phase:

$$y_1= Ae^{j\phi_1}$$
$$y_2= Ae^{j\phi_2}$$

$$y_1-y_2 = Ae^{j\phi_1}- Ae^{j\phi_2} = Ae^{j\phi_1}(1-Ae^{j(\phi_2-\phi_1)})$$
$$= Ae^{j\phi_1}e^{j(\phi_2-\phi_1)/2}(e^{-j(\phi_2-\phi_1)/2}-e^{j(\phi_2-\phi_1)/2})$$
$$ = -Ae^{j\phi_1}e^{j(\phi_2-\phi_1)/2}2j\sin\big((\phi_2-\phi_1)/2\big)$$
$$ = -2Aje^{j(\phi_1+\phi_2)/2}\sin\big((\phi_2-\phi_1)/2\big)$$
$$ = 2Ae^{j(\phi_1+\phi_2-\pi)/2}\sin\big((\phi_2-\phi_1)/2\big)$$

For a small difference angle $\phi_2-\phi_1$, $\sin\big((\phi_2-\phi_1)/2\big) \approx (\phi_2-\phi_1)/2$, so the above cancellation becomes:

$$\approx A(\phi_2-\phi_1)e^{j(\phi_1+\phi_2-\pi)/2}$$

And thus for small angles the magnitude of the cancellation is directly proportional to the difference in phase. And the cancellation ratio as the magnitude of the cancelled signal compared to the magnitude of the input signal would be the difference in phase directly:

$$|y_1-y_2|/|y_1| =| \phi_2-\phi_1|$$

To achieve what the OP is trying to do, the phase response of the filter must be compensated for before doing the subtraction. This would be much easier if a linear phase filter was used instead of the Butterworth filter as the phase compensation of the input signal could then be accomplished with a simple single tap delay. Of course, to completely cancel the signal the amplitudes would need to be exactly matched in the region of cancellation as well.

In practical application even with close matching a finite delay will exist between the two identical signals prior to subtraction, resulting in a high pass response with a cutoff of $1/T$ given a finite delay $T$. With a fixed delay $T$ the lower frequency components much less than $1/T$ in frequency would have a higher cancellation, and the decorrelation will go proportional to $f$ as the frequency approaches $1/T$. Signals much higher than $1/T$ would be decorrelated so would result in a $3$ dB increase as expected for the the sum of two uncorrelated signals.


Have you tried just median filtering the data ? That's what I'd try as a first step. Winsorizing really doesn't work on a data set that has a changing mean.

It looks like python has scipy.signal.medilt to perform this.


peak limiting/audio compression formula needed

Two problems here: how to get a reliable estimate of the level, and how to compress the data.

Use robust statistics on the original (not peak-limited) data like median or quantiles instead of a running average to make your "typical level" detection robust to outliers.
$k \times \tanh(\frac{x}{k})$ works nice as a $C^\infty$ compression formula and it is actually what is happening in some audio circuits (using OTAs). To get adaptive compression that preserves the dynamics of the original signal and just remove transients, make k track the smoothed "average" level.

Blue: original signal
Green: 2 x median of the absolute value over a sliding window as a "typical level" detection
Red: tanh compression (formula given above with k equal to the level plotted in green)


At what stage do we compute the approximations and details while performing a DWT?

This expression is more a discretization of a continuous wavelet transform than an actual DWT (discrete wavelet transform), provided $\psi$ is a genuine wavelet. It only computes the wavelet coefficient $c_{j,k}$ associated to a specific shifted and dilated continuous wavelet $\psi_{j,k}(t)$. This yields a frame-like wavelet decomposition, if you picture all possible $c_{j,k}$. A frame is a kind of overcomplete basis satisfying some bounds; in other words, more projection vectors than the dimension of the space, still ensuring a certain stability.

It is not a discrete wavelet transform per se, since time is still continuous. This is a bit like Fourier series. Approximations and details for DWT are more classically obtained from a direct discrete scheme. From a discrete signal $x[n]$, at one single level, it is possible to build two downsampled sequences:

$a(n) = (h_0\ast x)\downarrow(2)[n]$ (approximation)
$d(n) = (h_1\ast x)\downarrow(2)[n]$ (details)

with $h_0$  and $h_1$ complementary (low-pass and high-pass) filters, such that one can recover $x$ from half-length sequences  $a(n)$ and  $b(n)$ only. Coefficients for $h_0$  and $h_1$  can be obtained from certain wavelet shapes only. Passing the signal in a filter before downsampling reduces aliasing artifacts, while preserving invertibility.

Approximations and details (discrete wavelets and wavelet packets) can be obtained by iterating the above at different levels.
Stationary wavelets are a version of DWT without subsampling (but generally the same wavelets). They often behave better in case of higher noise,for denoising, or model mismatch in adaptive filtering (some details in A Primal-Dual Proximal Algorithm for Sparse Template-Based Adaptive Filtering: Application to Seismic Multiple Removal).


Amplitude in frequency domain not as expected with discontinues time domain signal

If you are seeing half the magnitude (amplitude is the height of the signal, magnitude is the size of a complex value), that means you are using a 1/N normalized DFT (BTW, my preference).  The sum will be the same (since the other half is orthogonal), you are just dividing it by twice as large of a number.  You should also find that all the even numbered bins are zero and the odd bins have values.  This is because the even bin sinusoids split evenly on your half way mark and the odd ones are integer and a half and therefore no longer orthogonal on those regions with your signal.

Hope this helps,

Ced

Followup

Demo for MBaz:

The output is:

Full: -0.000000 + i * 0.000000
Half: 230.265249 + i * 97.354586
Half Magnitude = 250.000000
Half Phase     = 0.400000

The code is:

#include <math.h>
#include <stdio.h>

//===============================================
int main( int argCount, char* argValues[] )
{

        double theSignal[1000];

        double theReal[1000];
        double theImag[1000];

        double theAngleFactor = 2.0 * M_PI / 1000;

        double theFrequency = 50.0;  // Cycles per frame

//--- Build the Arrays

        double thePhi = 0.4;

        for( int n = 0; n < 1000; n++ )
        {
            double theAngle = (double) n
                            * theAngleFactor;

            double theArg = theFrequency * theAngle;

            if( n < 500 )
            {
                theSignal[n] = cos( theArg + thePhi );
            }
            else
            {
                theSignal[n] = -cos( theArg + thePhi );
            }

            theReal[n] =  cos( theArg );
            theImag[n] = -sin( theArg );

        }

//--- Calculate the Bin Value

        double theBinReal = 0.0;
        double theBinImag = 0.0;

        for( int n = 0; n < 1000; n++ )
        {
            theBinReal += theSignal[n] * theReal[n];
            theBinImag += theSignal[n] * theImag[n];
        }

//--- Calculate for Half the Interval

        double theHalfReal = 0.0;
        double theHalfImag = 0.0;

        for( int n = 0; n < 500; n++ )
        {
            theHalfReal += theSignal[n] * theReal[n];
            theHalfImag += theSignal[n] * theImag[n];
        }

        double theMag = sqrt( theHalfReal * theHalfReal
                            + theHalfImag * theHalfImag );

        double thePhase = atan2( theHalfImag, theHalfReal );

//--- Show the Results

        printf( "Full: %f + i * %f\n" ,
                theBinReal, theBinImag );

        printf( "Half: %f + i * %f\n" ,
                theHalfReal, theHalfImag );

        printf( "Half Magnitude = %f\n", theMag );
        printf( "Half Phase     = %f\n", thePhase );
}
//===============================================

Appendum

To elaborate on my original answer.  I changed the code so the second half is defined like this, e.g. 10Hz:

                theSignal[n] = cos( 2.0 * theArg + thePhi );

Here are the sums:

Full: 230.265249 + i * 97.354586
Half: 230.265249 + i * 97.354586
Half Magnitude = 250.000000
Half Phase     = 0.400000

I left the amplitude at 1 so the magnitude of the first half unnormalized DFT is 250 (=500/2) as expected.  The sums are the same.

Epilogue - The smoking gun

Here are the unnormalized DFT calculations for the signal as specified by the OP near the bins of interest.

Bin 50 is the 5Hz bin.

Bin 100 is the 10Hz bin.

Sample rate is 100Hz.  N is 1000 for 10 seconds.

For normalized DFTs:

If you divide the left half sums by 500 (N for the left half DFT), your magnitude will be 1.0.
if you divide the whole sum by 1000 (N for the whole), your magnitude will be 0.5.

Notice that the 5Hz signal doesn't have much impact on the right side, and the 10Hz signal doesn't have much impact on the left.  As I've explained previuosly, the large values for the odd bins are due to the basis vectors in those regions being a whole integer plus a half number of cycles.

        L E F T    H A L F      R I G H T  H A L F           W H O L E
Bin       Real       Imag          Real       Imag          Real       Imag
====  =====================   =======================   =====================
 45     66.907      0.000        -7.771     -0.000        59.136      0.000
 46      0.000      0.000         0.000      0.000         0.000      0.000
 47    109.279      0.000        -7.959      0.000       101.320      0.000
 48      0.000     -0.000        -0.000     -0.000        -0.000     -0.000
 49    321.420      0.000        -8.166     -0.000       313.254      0.000
 50      0.000   -500.000         0.000      0.000         0.000   -500.000
 51   -315.264     -0.000        -8.392     -0.000      -323.656     -0.000
 52      0.000      0.000         0.000      0.000         0.000      0.000
 53   -103.118     -0.000        -8.641      0.000      -111.759     -0.000
 54      0.000     -0.000        -0.000     -0.000         0.000     -0.000
 55    -60.736      0.000        -8.915     -0.000       -69.651      0.000

        L E F T    H A L F      R I G H T  H A L F           W H O L E
Bin       Real       Imag          Real       Imag          Real       Imag
====  =====================   =======================   =====================
 95     -4.985      0.000       -65.080      0.000       -70.065      0.000
 96     -0.000     -0.000        -0.000     -0.000        -0.000     -0.000
 97     -4.714     -0.000      -107.504      0.000      -112.218      0.000
 98     -0.000      0.000         0.000     -0.000         0.000     -0.000
 99     -4.467      0.000      -319.694      0.000      -324.161      0.000
100     -0.000     -0.000         0.000   -500.000         0.000   -500.000
101     -4.240     -0.000       316.942     -0.000       312.701     -0.000
102      0.000      0.000        -0.000      0.000        -0.000      0.000
103     -4.033     -0.000       104.751     -0.000       100.718     -0.000
104     -0.000      0.000        -0.000      0.000        -0.000      0.000
105     -3.841     -0.000        62.325     -0.000        58.484     -0.000


The bi linear transform is the transform from the Laplace Transform Domain to the Z Transform.

The Laplace Transform Domain is a regular plane.
This transform transforms vertical lines in the Laplace domain into circles in the Z Domain.
Hence the Fourier Vertical Line in Laplace Domain (The Y Vertical Lines) is transformed into the unit circle in the Z Transform.

I'm a little inaccurate here, since the whole left side of the plane is transformed into the inner part of the Z Plane and the right hand of the Laplace transform it transformed to the outer part of the Z Plane (Out of the unit circle).

This transform is used to derive Digital Filters from their Analog Laplace domain representation.


How to deal with audio with different length in Audio Classification

If the signals were of the same duration but obtained using different sampling rates, so that they have different sample lengths, you can resample the shorter sequence with a higher sampling rate to match that of the longer sequence.


Efficient parallelized polyphase decimation FIR filter

Assume you do a polyphase decomposition into $M$ shorter filters, because that's really the best you can do in exploiting the decimation properties of your system.
Your $M$ shorter filters are just that: FIR filters.
You can then go ahead and take  these filters and decompose it into $c=\left\lfloor\frac KM\right\rfloor$ filters each.
That way, you get $cM$ parallelly computable filters, $K-M < cM \le K$, where $cM=K$ exactly when your $K$ is an integer multiple of $M$.
Note that I'd be rather surprised if you could get a huuuge performance boost out of this, assuming this is about computation on say a general-purpose CPU. The memory / cache bandwidth (and CPU cycles) occupied by deinterleaving or dealing with non-contiguous samples can kill performance rather quickly.
What I've seen is that people instead go and implement their $M$ shorter filters (which still might be massive!) as FFT-based fast convolution filters. Some FFT implementations (FFTw, specifically) allow you to multi-thread your FFT.
If you're not writing software for machines that prefer linear memory accesses, but, say, FPGA gateware, things look different. Often, significant gains come from exploiting symmetry of filter coefficients, or zero-taps in Nyquist-M filters; make sure your restructuring doesn't break these advantages!


I somehow "proved" that given any LTI system, its transfer function has to be constant. What am I missing?

When the input to an LTI system is an eigenfunction, then:

$$ x(t) = e^{j\omega_0 t} => y(t) = H(\omega_0) e^{j\omega_0 t}$$

where $H(\omega_0)$ is the associated eigenvalue of the LTI system; the transfer function evaluated at the particular frequency:

$$ H(\omega_0) = H(s)|_{s = j\omega_0}. $$

Here not the transfer function $H(s)$, but just one value of it is found by just one eigenfunction.

But if you use the general expression for the eigenfunction, then the input is $x(t) = e^{j\omega t}$ and the output will be $$y(t) = H(j\omega) e^{j \omega t}. $$ Here you assume that $H(j\omega) = a$ ; a constant,  and therefore in the fraction expression you will find out that $H(s) = \frac{Y(s)}{X(s)} = a = H(j\omega)$. Here, the constant $a$ is actually the unknown variable $H(j\omega)$...

Note that the eigenvalue $H(\omega)$ will be a complex-constant for a given particular frequency $\omega_0$. But as the given frequency changes, the constant (in general) will change too. It's just like an ordinary function $f(x)$ treated as a single value for a particular $x = x_0$ but a rather a varying mapping for general $x$ values...


Note that in this case you can see that the system is causal only from the given implementation. It's important to understand that you can't see it from the difference equation (if no initial conditions are given), and in general you can't see it from the transfer function either (if no region of convergence is given). The only case for which the expression of the transfer function is sufficient to conclude that the corresponding system is non-causal is if the transfer function is not proper, i.e., if the degree of the numerator (as a function of $z$, not $z^{-1}$!) is greater than the degree of the denominator. Also take a look at this related answer.
Let's look at an example:
$$y[n]=ay[n-1]+bx[n]\tag{1}$$
Does $(1)$ describe a causal system? Well, it can describe a causal system, but it can also be read as
$$y[n-1]=\frac{1}{a}\big(y[n]-bx[n]\big),\qquad a\neq 0\tag{2}$$
which can clearly describe a non-causal system if interpreted as "at time instant $n-1$, compute $y[n-1]$ from $y[n]$ and $x[n]$ according to the equation".
The corresponding transfer function
$$H(z)=\frac{b}{1-az^{-1}}\tag{3}$$
just conveys the same information as the difference equation $(1)$. Consequently, the expression $(3)$ doesn't reveal the system's causality either. Only if we know the region of convergence (ROC) of $(3)$ can we judge whether the system is causal or not. If the ROC is $|z|>a$, the corresponding system is causal, and and if the ROC is $|z|<a$, the system is non-causal.


Dependence of correlation peak height on GPS C/A signal's SNR, noise floor level (dBFS), correlation period length and ADC resolution

Even if your amplified signal matched the quantization noise floor of the ADC (representing a 3 dB degradation) you would still achieve a significantly larger and clearer peak from what you are showing.

The search bin size for 1 C/A code duration of 500 Hz is a good choice for the largest step (since the correlation versus frequency offset would follow a Sinc function response with the first null at $1/T$ where here $T$ is 1 ms, so the correlation would go to zero with a 1 KHz offset from true frequency.

With a 1 C/A code duration you can expect a processing gain of 30 dB ($10Log_{10}(1023)$) and the unspread GPS signal itself can be anywhere from -23 dB to -15 dB below the noise floor typically within the 1 MHz bandwidth of the main lobe of the GPS Signal. (You can measure your post correlation SNR by to see if it is in the ballpark of these levels plus the processing gain).

As for gain, the quantization noise of the ADC can be approximated by being 6 dB/bit + 1.76 dB below the full scale level of a sine wave, and that quantization noise is the total noise power spread evenly across the Nyquist bandwidth. I typically amplify the signal so it is at least 10 dB above this noise level at the input to the ADC, resulting in a 0.4 dB noise figure from the quantization noise. The 0.4 dB quantity comes from the power sum of the two noises (amplified thermal noise as the received signal and quantization noise): $10Log_{10}(1+10^{d/10})$ where d is the delta in dB that the received signal is above the quantization noise.

So if you ADC starts clipping with a sine wave that is at 10 dBm (for example), then for a 7 bit ADC, the quantization noise will be at 6dB x 7 bits + 1.76 = 43.8 dB lower or -33.8 dBm, spread across +/-1.023 MHz in your case with a sampling rate of 2.046 MHz. So within the 1 MHz main lobe of the GPS signal, the quantization noise will be 3 dB lower (since 1 MHz is half of this ADC bandwdith), so that puts the quantization noise in band at -36.8 dBm.

The received signal assuming a front-end noise figure of 3 dB is -174 dBm/Hz (thermal noise floor) + 60 dB (1 MHz bandwidth) + 3 dB (noise figure) = -111 dBm (these are typical numbers, your mileage may vary).

So if we want to place the recieved signal to be 10 dB above the ADC quantization noise floor in order to minimize further SNR degradation while maintain as high of a dynamic range as possible (the typical trade space and why we don't want to apply too much gain--we get little in SNR improvement but lose our ability to handle stronger signals or interference) we would need to have the received signal that starts at -111 dBm amplified to be -36.8 dBm (quantization noise floor) + 10 dB = -26.8 dBm.

With all those assumptions, your recommended gain is 84 dB. This hopefully shows you how to determine it for your specific case.


What acoustic properties of sounds enable us to distinguish the owner of the voice?

Some possible perceivable differences might include:

differences in pitch frequency (bass vs. soprano voice, etc.)
differences in overtone series (formants, head resonances, accent)
variations in modulation of pitch over time (vibrato, tremolo, accent)
differences in voiced to unvoiced ratios (raspiness?)
variations in typical durations of each vowel or consonant
variations in cadence modulation (accent?)
differences in choice of phonemes for each word or phrase (accent?)


Relation between the DTFT and the spectrum of a sampled signal

If you have a continuous-time signal $x(t)$, then the two signals you're talking about are

$$\begin{align}
x_c(t) &=x(t)\cdot\sum_{n=-\infty}^{\infty}\delta(t-nT) \\
 &=\sum_{n=-\infty}^{\infty}x(t)\delta(t-nT) \\
 &=\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT) \\ \tag{1}
\end{align}$$

and you define

$$x_d[n] \triangleq x(nT)\tag{2}$$

The first signal given by $(1)$ is technically a continuous-time signal, even though it is only non-zero at discrete times $t=nT$. The reason why it is considered a continuous-time signal is because it can and must be transformed using the continuous-time Fourier transform (CTFT). So $(1)$ is the continuous-time representation of a sampled signal. Eq. $(2)$ is the discrete-time representation of the same signal. Here the sampled signal is represented as a sequence of numbers. You can't apply the CTFT to $(2)$, but you must use the discrete-time Fourier transform (DTFT).

The nice thing is now that the CTFT of $x_c(t)$ given by $(1)$ and the DTFT of $x_d[n]$ given by $(2)$ are identical. So if the CTFT

$$\begin{align}
X_c(j\Omega) &= \int_{-\infty}^{\infty}x_c(t)e^{-j\Omega t}dt \\
   &= \int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x(nT)\delta(t-nT) e^{-j\Omega t}dt \\
   &= \sum_{n=-\infty}^{\infty}x(nT) \int_{-\infty}^{\infty}\delta(t-nT) e^{-j\Omega t}dt \\
   &= \sum_{n=-\infty}^{\infty}x(nT)  e^{-j\Omega nT} \\
   &= \sum_{n=-\infty}^{\infty}x_d[n]  e^{-j n(\Omega T)} \\
\end{align}$$

and the DTFT:

$$X_d(e^{j\omega})=\sum_{n=-\infty}^{\infty}x_d[n]e^{-jn\omega}$$

we have:

$$X_d(e^{j\omega})=X_c\left(\tfrac{j\omega}{T}\right)\tag{3}$$

In sum, the signals $(1)$ and $(2)$ are just two different representations of the same signal, and their spectra (one defined by the CTFT, the other defined by the DTFT) are identical.


You can investigate that if the size of the DCT block (which is the whole image in here) is $x\times y$, then

c=d(1,1)/sqrt(x*y);

That is, the DC element of the DCT matrix (i.e. d(1,1)) is actually the mean of the image scaled by $\sqrt{xy}$.

So you are basically subtracting a scaled version of the DC component form all coefficients, which makes not much sense at the moment, except if you provide more context.


How to Design Nonlinear Phase FIR Filter?

Although not a comprehensive list, the following are filter algorithms available in Matlab for designing non-linear phase FIR filters:

firlpnorm: Nonlinear phase equiripple or least squares, uses Lp norm: for equiripple filters p = infinity, for least squares p =2.

cfirmpm: complex and nonlinear phase equiripple filters.

See the Matlab documentation for full details on how to design non-linear phase using these algorithms.


How to segment phone call audio into silence/non silence?

There are a bunch of parameters that you can look at:

Overall energy
Short term spectrum: Speech has a fairly distinctive "pink-like" spectrum and noise (which is happening during the non-speech parts) tends to be white if it's electrically dominated or "red" (i.e. low frequency heavy) if it's acoustic background noise or microphone noise
Amplitude statistics. Most noise signals have a Gaussian distribution, speech is closer to a Laplace distribution

I think a combination of these three should give a fairly robust detection scheme.


Are There Common Values of Standard Deviation for Gaussian Noise of an Image?

You can easily have a look on the values of the STD on Image Denoising Papers:

The range of 1-15 is considered low.
The range 15-30 is considered medium.
The range 30-50 (Even above) is considered high.

The above values is for images in the range [0, 255].


Find Fourier series $f(t) = e^{jx t}$ , $−\pi < t < \pi$

There is no need for you to find Fourier series just compare with equation $x(t)=\sum_{-\infty}^{+\infty}C_n*e^{j*w_n*t}$ we will know $C_n,w_n.$


Can a Goertzel filter response be mapped to a Direct Form II Biquad IIR Filter

try it, you'll like it.

the form in the Jocobson & Lyons paper is Direct Form II, with the coefficients explicit.  the form in the CMSIS source is a transpose of a Direct Form I, i believe, and the one in the Lyons & Bell paper is a Direct Form II (not transposed, but set $a_0=1$ because they should not have had that $a_0$ gain block in the path.).  the $b_0, b_1, b_2$ coefs come from the numerator and the $a_1, a_2$ come from the denominator and need to have their signs changed (since the scaled feedback signals are added and not subtracted in the feedback).

are you floating-point or fixed-point?  if the latter, the untransposed Direct Form I will serve you better.


Since $x(t)$ is a random process, its Fourier transform generally doesn't exist. In its original form, the process $x(t)$ doesn't even have a power spectrum in the conventional sense because it is cyclostationary (and not stationary). It can be made stationary by adding a random phase $\theta$ that is uniformly distributed on $[0,T]$ and independent of $a_n$:

$$x'(t)=x(t+\theta)=\sum_na_nh(t+\theta-nT)\tag{1}$$

This process is stationary and its power spectrum can be computed, which is usually done by first determining its autocorrelation function, and then applying the Fourier transform to the autocorrelation function. The corresponding derivation can be found in most textbooks on digital communications:

$$S_{x'}(f)=\frac{1}{T}|H(f)|^2S_a(f)\tag{2}$$

where $H(f)$ is the Fourier transform of the pulse $h(t)$, and $S_a(f)$ is the power spectrum of the discrete-time random process $a_n$.

If we can assume that the symbols $a_n$ are uncorrelated, the power spectrum $S_a(f)$ is constant, $S_a(f)=P_a$ (where $P_a$ is the power of $a_n$), and $(2)$ becomes

$$S_{x'}(f)=\frac{1}{T}|H(f)|^2P_a\tag{3}$$

Finally, the power (not energy) of $x'(t)$ is the integral over its power spectrum (note that $h(t)$ is band-limited):

$$P_{x'}=\int_{-\infty}^{\infty}S_{x'}(f)df=\frac{1}{T}P_a\int_{-\infty}^{\infty}|H(f)|^2df=\frac{1}{T}P_aE_h\tag{4}$$

where $E_h$ is the energy of the impulse $h(t)$.


Hilbert Transformer processing block by block

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

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

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

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

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

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

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

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

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

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


I have two explanations with the intent of providing additional intuitive insight beyond the concise mathematical explanation; first, an explanation from the Fourier Series Expansion along with the idea of repeating the truncated time domain sequence which shows the implied discontinuities that result, requiring more frequency components to reconstruct than would really exist if the waveform was not truncated. And second, an explanation from looking at the DFT as a bank of not-so-great filters.

First Explanation: Fourier Series Expansion  and Time Periodicty

The Fourier Transform of a time limited sequence is identical to the Fourier Transform of a sequence that is periodic for all time.

This is true for the Fourier Transform as well as the Discrete Fourier Transform and is known as the Periodicity Property of the DFT:

Periodicity Property: Given the N-point DFT vector X[k], with an Inverse DFT samples x[n], with both k and n ranging from 0 to N-1; if n lies outside of the range of 0,1,2...,N-1, then

$x[n] = x[modulo(n,N)]$

Likewise, given a N-point time sequence x[n], with a DFT X[k], if k lies outside of the range of 0, 1,2,..., N-1, then

$X[k] = X[modulo(k,N)]$

To help gain an intuitive understanding, a key take-way is that anything that is sampled in one domain, becomes periodic in the other domain. Similarly, anything that is periodic in one domain, becomes sampled (discrete-valued) in the other domain. Here specifically, "sampled" means that the signal will only exist as non-zero values at discrete locations in the domain (stream of impulses).

Sampling in One Domain -> Periodicity in the Other Domain: Consider the spectrum of a 3 Hz analog cosine wave before and after sampling with an A/D converter. The digital spectrum can be viewed as periodic; a cylinder view of the digital spectrum is also a valid view to explain the periodicity, but I find this extension to the analog frequency domain (by having the frequency extend to +/- infinity) helps some gain an intuitive view of the signal processing involved.

Periodicity in one Domain -> Sampling in the other Domain: Consider the Fourier Series Expansion as a simple example demonstrating this property.  The Fourier Series Expansion is done over a finite time domain interval from 0 to T. When decomposed into separate frequency components, the only frequencies used are DC, the fundamental frequency 1/T, and integer multiples of 1/T (harmonics). In effect, since the frequencies can only exist at multiples of 1/T (and DC),  the frequency domain has been sampled.

Further, if we reconstruct the time domain waveform, by summing together the individual frequency components, we can also view the implied periodicity in the time domain if we allow the frequency components to extend beyond the interval 0 to T. It is because of this periodicity that frequency components cannot exist at any other frequencies than multiples of 1/T (because of the opposite condition: if they did exist, then they would not begin and end consistently in the 0 to T time interval, and therefore periodicity cannot exist).

Understanding the above will hopefully help provide an intuitive explanation of spectral leakage.  So now I will repeat a main point:

The Fourier Transform of a time limited sequence is identical to the Fourier Transform of a sequence that is periodic for all time.

Spectral Leakage with the "Fourier Series Expansion" View:

Consider two sinusoidal waveforms in the time interval from 0 to T, the first with an integer number of cycles over the time interval, and the second case with a non-integer number of cycles.

Clearly in Case 1 we can see with the repetition view that even repeating our pure sinusoid remains a pure sinusoid, but in Case 2 our sinusoid now suffers from abrupt transitions, and using the Fourier Series Expansion view of reconstruction, would requred several frequency components to reconstruct such a waveform in the time domain.

Second Explanation: Filter Bank View of the DFT

Another intuitive explanation for spectral leakage (and helps significantly with understanding the DFT in general) is what I call the filter bank view of the DFT.
To see this, consider a simple 4 pt DFT as shown in the figure below and observe for each bin, we are effectively rotating the signal, and then passing the rotated values through a 4 tap unity gain FIR filter. For the first bin, which corresponds to DC, there is no rotation, so we just sum the four samples, and for the other bins we progressively rotate at higher frequencies as we move through the DFT bins:

(Side note- If we performed a streaming DFT, where we calculated a new 4 point DFT on a 4 point sequence as we scanned through a waveform, it would exactly be such a filter bank, but regardless if we do that or not, this view gives great insight into spectral leakage in addition to the usual convolving of a sinc function in the frequency domain that the mathematical explanation reveals)

Now consider the frequency response for each equivalent FIR filter, using the coefficients given in the DFT (For example use freqz([coeff]) in Matlab or Python) as shown in the figure below. (Note that the filter coefficients are in reverse order from that shown in the DFT matrix. This is because the coefficients of the filter are reversed in time; the earliest samples in time are at the far end of the FIR filter state register, so reversed in time):

Here is the main point: Since each Filter in the construction of the DFT is fundamentally a unity gain FIR filter, the shape of this filter in frequency approaches a sinc function as the length of the DFT gets longer (and is an aliased sinc function for small N). So we will call these sinc filters and note that a sinc filter has relatively high sidelobes and has an envelope that rolls off very slowly with frequency (at 1/f for a pure sinc). With the phase rotators in the DFT, we just move the main lobe of this sinc filter to each bin of interest, but the side lobes that exist for each bin allow frequencies at other locations to cause energy to appear in this bin. The amount of leakage is completely predicted by these filters.

Using this view, consider an input single-toned signal with a frequency that is somewhere between two frequency bins as shown in the figure below (an input that exists exactly on any bin will have an integer number of cycles in the time domain interval and therefore NO spectral leakage as we showed earlier). The top filter shows the amplitude at this frequency location that will "leak" to the first bin. The second filter shows the amplitude (slightly higher), the third bin (that our frequency in closest to) will have the highest response and the fourth bin will be lower.

SUMMARY

I have presented two explanations with the intent of providing additional intuitive insight beyond the concise mathematical explanation of multiplication by a rectangular window in the time domain is convolution in the frequency domain (and the leakage we see therefore is the result of a sinc function convolving in frequency with our waveform of interest that we truncated in time); first, an explanation from the Fourier Series Expansion along with the idea of repeating the truncated time domain sequence which shows the implied discontinuities that result, requiring more frequency components to reconstruct than would really exist if the waveform was not truncated. And second, an explanation from looking at the DFT as a bank of filters, and poor filters at that (specifically unity gain filters that approach a sinc function frequency response as N increases). Such filters provide the best frequency resolution but poorest dynamic range.


long decaying tails which can make them robust

Why would filter long decaying tails make a signal robust

can make them robust, not makes them robust!
You need a long filter length for narrow bandwidth and steep transition width.
So, with a long filter, you can make very narrow subcarrier channels.
Those being narrow, they don't "see" multipath, they only see a flat channel.


If you know the number of sinusoidal components, parametric estimators such as MUSIC may be suitable. ESPRIT might be another.  Both seem similar to you proposed optimization/minimization search.

But parametric estimators are reported to produce bogus decompositions if given the wrong number of components.  Possible reason: For a decomposition into an arbitrary number of arbitrary frequency sinusoids, there can be a near infinite number of solutions that might be arbitrarily close (e.g. you can get a bunch of cosines to converge to one sine wave of yet another frequency.)


Why dont sin and phase shifted cosine overlap?

You're confusing degrees and radians.

$$\sin(x)=\cos(x-\pi/2)$$

$\pi/2$ radians corresponds to 90 degrees.


Spectrum, FFT magnitude, and pwelch are pretty much the same thing, i.e. a metric of energy density as a function of frequency. They all can be configured in different ways: window length, window shape, overlap, FFT size, lin/log representation, smoothing, different averaging methods, etc.
All these parameters need to be chosen carefully depending on the requirements of your application: how fast do you need to track, what frequency resolution do you need,  what's the  noise spectrum look like etc.
Choosing the parameters properly is actually more important than the core algorithm (which is more or less the same anyway).


Why does Pre-warping work with a notch filter but not a low pass filter?

The bilinear transform compresses the analog frequency in the s-domain associated with  $f=\infty$ to the point $z=-1$ in the z-domain (which is where the frequency is $f_s/2$ where $f_s$ is the sampling rate. Consider what the magnitude of your notch filter is at $f=\infty$: It has no effect on the magnitude or phase of your signal; there is no filtering impact there whatsoever.
Compare in contrast to this your low pass filter (and how much of your "desired" shape is already beyond the Nyquist boundary which no pre-warping would be able to reconstruct). As you approach the Nyquist boundary the ability to compensate with pre-warping becomes increasingly difficult (you will see the same with the notch if it is brought close enough such that there is continued attenuation beyond the boundary).
The degree of this frequency warping and how it pronounced is depicted in the chart example that I pasted below (not the OP's sampling rates) which shows which frequency on the vertical axis for a given sampling rate of 20 rad/sec (125.6 Hz) would map to the digital frequency domain extending from $\pi$ to $\pi$ rad/sample.
With prewarping, especially at high frequencies, it is very difficult to match an analog shape over a wider frequency range given this mapping is according to the following transcendental equation (would require Taylor series expansion out to an increasingly higher order as we near the Nyquist boundary to match effectively, but is easily matched with the first or second terms only at the lower frequencies where the slope becomes approximately linear):
$$\omega_a = \frac{2}{T}tan(\omega_d T/2)$$
(as referenced here:  https://en.wikipedia.org/wiki/Bilinear_transform)
But for a single point, or narrow range as typical for a notch filter, this repositioning of a single notch zero would be trivial to do and with high accuracy for the repositioned notch location (but distortion right at the Nyquist boundary if the target filter has residual insertion loss beyond Nyquist).


I've found a proper definition thanks to the following article:
http://www.rehab.research.va.gov/jour/2013/509/jrrd-2012-12-0233.html

This is what I've edited Wikipedia with:
$$f_{sma}= {1 \over {T}} {\int_{0}^{T} |x(t)-a_x|+|y(t)-a_y|+|z(t)-a_z| dt}$$

(It wasn't because of text-search but images-search. Algorithms are sometimes surprising)


How to apply Prewitt Edge Detection algorithm on a given image?

I was able to answer my question after doing more research. The key is to slide the 3x3 kernel and perform Gx and Gy (gradient for the x and y directions) and calculate the output image values. For example, to calculate the output image top left value -

[x - -]
[- - -]
[- - -]

we have to calculate both gradients (x and y)->

Gx

[-1 0 1]
[-1 0 1]
[-1 0 1]

we take the first chunk of the original image (e.g. sliding the kernel over the top values of the input image)

[1 2 4]
[5 2 5]
[5 4 2]

Gx = 4+5+2 - 1-5-5 = 0

Gy = 5+4+2 - 1-2-4 = 4
We take the absolute values of both and sum them so |4| + |0| = 4

And this is the first top value of the output image.

[4 - -]
[- - -]
[- - -]

The same way we calculate the rest. I hope that will help anyone else facing the same problem.


Edit:
I've come to realize that my definition below of "Frequency Resolution" is completely wrong (as well as OP's question). Frequency resolution is how similar the window function's magnitude in frequency space is to the Dirac delta function. This is because the product of the window and the signal in the time domain becomes convolution in the frequency domain (and a convolution with the Dirac delta function is a sampling which would give perfect frequency resolution) The fatter the mainlobe (quantified by it's variance), and the higher the sidelobes, the worse the frequency resolution. Additionally, Time resolution can be quantified a the variance of the window function in the time domain.

Frequency Resolution is not Bin Resolution/Width. In the graph below notice that the lobes dont get closer (frequency resolution) even though bin width is decreasing.

Credit: Dan Boschen
Frequency resolution is rather a property of the Fourier transform of the rectangular function (i.e. the sinc function).

We must window functions to work with Fourier transforms (even when working theoretically). As a consequence we are always working with $f(t)w(t)$ rather than the function $f(t)$ itself (here $w(t)$ is a rectangular function). By the Convolution theorem the Fourier transform of a windowed function is always a convolution of $\hat{f}$ with $\hat{w}=$ sinc. Notably when $f$ is sinusoidal, $\hat{f}$ will be a Dirac delta function and the convolution will just be a sampling of a sinc function. Thus we periodically lose frequencies completely when windowing, the periodicity of this loss is the frequency resolution.
Since, on windowed functions, the DTFT is a periodic approximation of the CTFT it also acquires these properties.
The confusion arises because when we dont pad zeros to the DFT (i.e. only sample $f(t)w(t)$ where $w(t)=1$), the bin width is equal to the Frequency resolution.

However, we may also pad zeros (i.e. also sample $f(t)w(t)$ where $w(t)=0$) and this results in the DTF better interpolating the DTFT of $f(t)w(t)$. Confer with the first graph.

To see why the Fourier transform of the rectangular function is the a sinc funcion watch this video and consider the winding of the sinusoidal functions (its quite involved though)

To answer OP's example the bin resolution is $$\frac{F_s}{N} = \frac{2000}{2000}=1$$ where $F_s=2000$Hz is the sampling rate, and $N$ the DFT size.
The frequency resolution is what the bin resolution would be if we just sampled in the window (no zero padding)
$$\frac{F_s}{M} = \frac{1}{T}=2$$ where $M$ is the number of samples in the window, $T$ is the duration of the sample, and $F_s=M / T$.


Although it's not entirely clear to what you try to do (maybe there are better ways to accomplish your goal), I'll try to give some hints.

There exist "time stretching" techniques that can slow down a sample without changing its pitch, but they are not trivial to implement yourself. Some audio editors have plugins that support this.

As a poor man's implementation, if the sound has a single pitch, you can try something like the following (I've successfully used this in the past to extend a long, sung, single note):

First I cut the sound in three parts:

Sound -> | Begin | Middle | End |

Always cut at zero crossings.

Then, if you make "Middle" short enough, like around 0.1 seconds or so, you can repeat it a number of times times to stretch your sound a bit:

| Begin | Middle | Middle | Middle | Middle | Middle | Middle | End |


Comparison of Linear Convolution and N Point DFT

The result {4,1,2,3} is the circular convolution of {1,2,3,4} and {0,1,0,0} which you correctly get by taking the inverse DFT of the product of the DFTs of the two sequences.

We can check this by doing the circular convolution the long way via matrix multiplication as follows:

$x(n) = [1,2,3,4]$

$y(n) = [0,1,0,0]$

To solve for the circular convolution $x(n) (*) y(n)$, (where I use $(*)$ to denote circular connvolution) we can do that with multiplication as $Xy$, where $X$ is a matrix formed by repeating $x(n)$ in each column with a circular rotation of shift n with n = 0,1,2,3 and $y$ is $y(n)$ as a column vector:

$$ Xy =\begin{bmatrix}
1 & 4 & 3 &2 \\
2 & 1 & 4 & 3 \\
3 & 2 & 1 & 4 \\
4 & 3 & 2 & 1 \\
\end{bmatrix} \begin{bmatrix}0 \\ 1 \\ 0 \\ 0 \end{bmatrix} = [4,1,2,3]$$

Confirming the relatinship:

$$x[n](*)y[n] = IFFT(FFT(a)FFT(b))$$

Where the ($*$) operator here is specifically circular convolution. And the multiplication of the two FFT results is the element by element product.

Also to add to this: In comparison, if you instead conjugate one of the FFT results prior to multiplication as follows, the result is a circular cross-correlation of the two sequences!

$$XCORR = IFFT(FFT(a)FFT(b)^*)$$


Using Principal Component Analysis for corner feature extraction

Neural networks can cope with not linearly separable classes, but they're not magic either. They still assume smooth decision boundaries, and I don't think your data would have those.

And before you try to feed the same data to an SVM, RVM, random forest, or any other classifier: I don't think the corner coordinates are good features for any off-the-shelf classification algorithm. The typical input vector in machine learning is something like {age,weight,height,eye color}, where for every input vector, index 2 means "weight". In your input vector, index 2 can sometimes be a corner somewhere near the front wheel of the car, and in another case it's a corner near the roof.

Have you tried simpler shape descriptors, like Fourier descriptors or maybe Hu moments?


The ScaleValues property is easy to incorporate. Multiply the entire input signal by the first scale value, then the output of each section by subsequent values.


Adding Noise in Frequency Domain vs Time Domain

I usually like to add noise with my own code.
Adding Whit Noise is really easy task given the signal as you have.

All you need is to calculate your signal second moment at the frequency and add noise to the frequency bins such that the second moment of the noise creates your desired SNR.

Since the DFT is unitary transform, adding white noise at frequency domain is equivalent to adding noise at time domain.

There 2 options here, are you looking for adding White Noise (Namely spread all over your frequency) or Band Limited (Colored) noise?

For white noise just add Noise to each bin in the frequency domain such that the noise Second Moment (The Means squared multiplied by the number of bins) is as your desired SNR.

For colored noise do the same but only for the relevant frequency bins (Where the second moment is calculated only on the number of "Active" bins).


generator matrix coefficient in convolutional code

It is a systematic  rate-$\frac 12$ convolutional with constraint length $2$.


Is there a difference between Image Alignment and Stereo Rectification for stereo correspondence or matching?

There two ways to look at this problem.

In simple terms, image rectification warps both images onto a common coordinate frame by typically estimating the transformation using the epipolar geometry.
Image alignment finds the transformation from one image to the other. It doesn't guarantee any constrains on the epipolar geometry and only one single image is warped.

Image rectification might benefit from image alignment (of local patches) to obtain the correspondences. Once these are found, the epipolar lines can be forced to be parallel.


FIR filter design by the Fourier transform method

Some papers provide rule-of-thumb formula for the filter length (and coefficient quantization), for instance:

Accurate estimation of minimum filter length for optimum FIR digital filters, 2000, Ichige et al.
On computational complexity in digital filters, 1981, Bellanger
FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation, 1975, Rabiner et al.

The one I generally  use in lectures is borrowed from the second paper, for a low-pass filter with the following design:

Then, the estimated filter length $N$ (apparently $2M+1$ for you) is, with sampling frequency $f_s$:

$$ N = \frac{2 f_s}{3(f_2-f_1)}\log_{10} \left(\frac{1}{10\delta_1\delta_2}\right)$$


Audio signals on line level outputs are AC signals with no DC bias, so they run positive and negative around the 0V level. Many devices support 2VRMS output for full scale signal, so that is about 5.66 Vpp.


Contourlet transform for image compression

I do not know of any book solely devoted to contourlets. There is a chapter in Beyond Wavelets, 2003 (G. Welland, Ed.) by Do and Vetterli, named Contourlets. You also have a 100-page paper (almost a book) where contourlets are put in the perspective of Multidimensional Filter Banks and Multiscale Geometric Representations (2012).

At the risk of mildening your teacher's assertions, good compression is the result of a transformation with nice qualities and an efficient coding scheme. The "best" transform can fail to properly compress an image if you do not have a neat coding scheme associated with it. Contourlets, for instance, are slightly redundant. And its compaction capabilities are thus limited by the redundancy. To my knowledge, redundant transforms have not proved yet they are much better than critically-sampled ones for compression (except perhaps at large compression ratios).
Even if you are not talking about research papers, let me mention a couple of different contourlet versions. The CRISP-contourlet, which is non-redundant, and The Nonsubsampled Contourlet Transform, fully redundant.

When trying to understand a transformation, especially with wavelet-like decompositions, it is often useful to look at the interplay between the math and the sampling schemes.

In addition to Minh N. Do thesis, I am trying to find a copy of A. L. da Cunha PhD.


How to add a particular noise on a signal?

Please see this post that details most of the noise considerations for an ADC

What are advantages of having higher sampling rate of a signal?

In particular the relationship for the noise floor given a full scale sine wave:

$$SNR = 6.02 \text{ dB/bit} + 1.76 \text{ dB}$$

Using your numbers:

1.5Vpp sinewave. If we assume this sinewave is full scale (any higher voltage would start to add clipping noise), then this is the level of 0 dBFS. In particular the variance of this or the standard deviation is used in the dB ratio computation according to:

$$SNR = 20Log_{10}(\sigma_s/\sigma_n)$$

or

$$SNR = 10Log_{10}(\sigma_s^2/\sigma_n^2)$$

Since the noise itself is given in units of dBFS, there is nothing to do further with the 1.5V level. This level is important to determine if the signal is indeed at full scale, and if we wanted to determine the noise level in units of volts.

120 dBFs/Hz with a 200 MHz Sampling Rate:

The quantization noise can be reasonably assumed to be a uniformly distributed white noise (see link above). So the total noise in 200 MHz of bandwidth would be:

$$-120 \text{ dBFS} + 10Log_{10}(200e6) = -37.0 \text{ dB}$$

(Note that this is the double-sided noise in the spectrum extending from $-f_s/2$ to $+f_s/2$)

From the formula above we see that this would be 6.16 bits of precision. A 6 bit converter with the sinewave right at full scale would be very close the the added noise that was given in the example.

To "add" such noise to the signal, scale the signal to fit peak-peak with a 6 bit system and round to closest integer to create the resulting quantization noise.

So the +/-1.5V peak to peak sine wave would then be scaled to $0$ to $2^6$, or $-2^5$ to $+2^5-1$ if you want to represent a signed signal.  The noise will occur from either rounding or truncating to integer values. (There is no noise difference in rounding or truncation; truncation simply introduces an offset).

As far as the FFT, each bin in the FFT if not windowed has a bandwidth 1 bin wide. So the FFT noise floor will go down according to the number of bins. This is explained further in this post:

How large FFTs can pull signals out of the noise floor?


Significance of modular arithmetic in DFT?

This answer to a differently worded question answers your question also.

And this answer is relevant.

And this answer and this answer.

I confess that I take great exception to anyone who denies that the Discrete Fourier Transform (and inverse DFT) fundamentally periodically extend the finite set of data presented to it.  I'm a nazi about it.  The DFT assumes that the $N$ samples of $x[n]$ passed to it are one period of a periodic sequence where $x[n+N]=x[n]$ for all integer $n$. And the inverse DFT does the same in that the $N$ samples of $X[k]$ passed to it are one period of a periodic sequence where $X[k+N]=X[k]$ for all integer $k$.

The modulo arithmetic done to the indices of $x[n]$ and $X[k]$ simply guarantee this periodicity while also guaranteeing that neither index, $n$ or $k$, go outside of the bounds of the finite sized arrays in memory representing $x[n]$ and $X[k]$.  That is

$$ 0 \le n = \operatorname{mod}_N(n+mN) < N $$
and
$$ 0 \le k = \operatorname{mod}_N(k+mN) < N $$

for any integer $m$.


Oversampling here would be extremely easy, so do it!

The benefits are of a noise nature (which is probably why you're doing the whole lock-in amplifier dance in the first place); you get a part of your oversampling factor as SNR gain.


Prior to upsampling, you have a white signal meaning every single frequency in the Nyquist bandwidth from $-\pi$ to $\pi$ is represented. This is a requirement to obtain an impulse (because the Fourier transform of an impulse is a white spectrum). The reason that white signals are often used as inputs for purposes of system identification is that they excite all frequencies of the system. This is important when trying to determine the system's response.

After upsampling and interpolating by $N$, you obtain a sequence with a spectrum that is more or less white in the bandwidth from $-\frac{\pi}{2N}$ to $\frac{\pi}{2N}$. The rest of the band (if your interpolation filter is working correctly) will be nearly zero or null. That means that these frequencies are no longer contributing to your coefficient calculation. In turn, this means that you cannot obtain a delta function since it requires every frequency present to be represented properly. The closest approximation to a delta function that you can obtain is a sinc pulse. This is why you are observing this result.

Summary: You cannot estimate the response at frequencies that are not contained in your input signal. Your input signal and desired signal are no longer white (because they are bandlimited by the interpolation filter). Therefore, the adaptive filter does not estimate those frequencies. This leads to the sinc response you are observing which is the closest band-limited response to the true system response (a delta function).


I have worked with this kind of a positioning system before (NOT underwater, but still sound based), and I hope I can help. Your first option is probably not the best idea, even in a controlled environment there are many factors which can change your signal's amplitude and therefore making it unreliable information.

I think your second option is a better way to go. It could be even better if you can embed some time information to your transmitted signal as well (a timestamp like GPS). With a timestamp containing information about when the signal was fired, you can calculate the time difference without ping-back, which would save you some distance error. Or you would have two methods to determine the position and you can perhaps use some sort of filtering depending on a set of initial results (to decide which one works better).

Another method that I would like to suggest is an array configuration. If you can make sure that your transmitters will transmit at ABSOLUTELY at the same time (with respect to each other), you can use the phase differences of the received signals in your boat to determine the range and the respective angles(with some ambiguity). So again you wouldn't need to ping back your signal. But in this configuration the absolute positions of your transmitters must remain the same, and the positioning of your receivers in your boat must remain the same as well.

When we used the ping method, our biggest problem was echo. I don't know exactly how echo works underwater but if it exists, I would suggest you to be careful about it. In our project we had only used the ping without any extras that I suggested to you, and I must admit it wasn't very reliable.

Good luck with your project, I hope you do better than us, and I would like to hear about the results.


Is there any free software that can recognize and classify sounds?

The audio feature extraction libraries given above could be a good start, with some caveats:

The features you need for this task are plain MFCCs, so most of those libraries are totally overkill since they also extract features which are more meaningful in a musical context.
Some of those libraries have been designed for offline analysis (processing a pre-recorded file rather than an ongoing stream of audio buffers). Your application in home-automation suggests that you want something capable of producing on-the-fly results as ambient audio is recorded.
None of those libraries actually address the very problem of recognizing sounds ; they just extract audio features upon which classification can be performed, so you need to post-process their output with machine learning software - using a model previously trained on examples of each class of signals you want to discriminate.

You'll end up having less plumbing to do if you use a speech recognition toolbox like HTK or Sphinx. HTK is generic enough for your task - it offers you both feature extraction and recognition/classification with HMM (a very common machine learning technique not only used for speech recognition but also for audio categorization in video analysis), and has extensions for on the fly recognition.

From what I remember, the MERL system you mention, and pretty much every work they did on audio for multimedia document analysis, is based on HMMs on the MPEG-7 audio descriptors.


How to Remove Double Lines Detected Along the Edges by Edge Detector?

I ran the following code:
clear();
close('all');
mI = im2double(imread('bCfdb.png')); %<! Loading the image
vBlurStd = [0, 0.1, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2];

mII = sum(cat(3, 0.299, 0.587, 0.114) .* mI, 3); %<! Y (Luminosity like channel)

hFigure = figure('Position', [100, 100, 1200, 900]);
hTiledChartLayout = tiledlayout(3, 3);

kk = 0;
for ii = 1:3
    for jj = 1:3
        kk = kk + 1;
        if(kk == 1)
            mE = edge(mII, 'canny');
        else
            mE = edge(imgaussfilt(mII, vBlurStd(kk)), 'canny');
        end
        nexttile();
        imshow(mE);
    end
end

set(hTiledChartLayout, 'TileSpacing', 'tight', 'Padding', 'tight');

This is the result:

What Does the Script?

We set an array of different blur radius to evaluate the algorithm at - vBlurStd.
We load the image and convert it into Double in the range [0, 1] - mI = im2double(imread('bCfdb.png'));
We extract the Y channel from the YCbCr color model - mII = sum(cat(3, 0.299, 0.587, 0.114) .* mI, 3);.
The algorithm is basically blurring the image with Gaussian Blur and then applying Canny Edge Detector - edge(imgaussfilt(mII, vBlurStd(kk)), 'canny');.

The images are left to right, top down by their blur radius.


You compute the energy (before the log operation) as

$$\max\left\{\sum_i s^2(i),\epsilon\right\}\tag{1}$$

with some small constant $\epsilon$ (they suggest $\epsilon=2e-22$ which gives $\ln\epsilon\approx -50$) to prevent the energy to become too small. The reason for this is that the logarithm of zero is $-\infty$ and gives you numerical trouble. So you always want to take the logarithm of a strictly positive number, never of zero. As far as I understand, the choice of $\epsilon=2e-22$ corresponding to a log value of $-50$ is rather arbitrary.


First you have to remember a Laplace transform property:

$$ e^{a t} u(t) \longleftrightarrow \frac{1}{s-a} ~~~,~~~ \mathcal{Re}\{s\} > \mathcal{Re}\{a\}   \tag{1} $$

$$ -e^{a t} u(-t) \longleftrightarrow \frac{1}{s-a} ~~~,~~~ \mathcal{Re}\{s\} < \mathcal{Re}\{a\}   \tag{2} $$

This property states that a given Laplace transform $H(s)$ can correspond to multiple inverse transforms, depending on the region of convergence ROC.

Therefore given $H(s) = \frac{1}{s-a}$ , you can find two possible inverses as $x(t) = e^{at} u(t)$ or $x(t) = -e^{at} u(-t)$ depending on whether the ROC is to the left or right of the pole location. Note that one of them is causal and the complementary is anti-causal.

Now, given your transfer function

$$ H(s) = \frac{s-1}{(s+1)(s-2)} = \frac{2/3}{s+1} + \frac{1/3}{s-2} ,  \tag{3} $$

it has two poles at $s = -1$ and $s=2$. There will be three possible ROC's with three different inverses :

$$
\begin{align}
\mathcal{Re}\{s\} < \mathcal{Re}\{-1\} & \implies h(t) = -\frac{2}{3} e^{-t}u(-t) - \frac{1}{3} e^{2t} u(-t)  \tag{4} \\
\mathcal{Re}\{-1\} < \mathcal{Re}\{s\} < \mathcal{Re}\{2\}  &\implies h(t) = \frac{2}{3} e^{-t}u(t) - \frac{1}{3} e^{2t} u(-t)  \tag{5} \\
\mathcal{Re}\{s\} > \mathcal{Re}\{2\} & \implies h(t) = \frac{2}{3} e^{-t}u(t) + \frac{1}{3} e^{2t} u(t)  \tag{6}
\end{align}
$$

The impulse response in (4) is anti-causal and unstable.

The impulse response in (5) is non-causal and stable.

The impulse response in (6) is causal and unstable.


Developing Lead Compensator to Decrease Settle Time

I deleted a previous attempt using a approach with the single lead compensator and provide this alternative solution below.
The root locus shows all closed loop pole locations, which follow the root locus trajectories as the loop gain is adjusted. From that we can choose closed loop pole locations that meet our requirements for settling time, loop bandwidth, stability etc. In general the dominant poles will set the loop characteristics such that we can ignore all "weaker" poles. Dominant poles are the ones closes to the imaginary axis (in the left half plane specifically for stability). This is because the impulse response decays faster as the pole is further from to the imaginary axis (consider that a real pole at $s = -x$ has an impulse response given by $e^{-xt}$), so if one pole is significantly closer, it will be continuing to decay long after the further pole has mostly diminished.
Here is the root locus for the OP's plant prior to any compensation. All open loop poles will terminate in the open loop zeros as the loop gain is increased, and the open loop zeros never change, becoming the closed loop zeros. Since there are more finite poles than finite zeros (as a proper system), there is an additional zero at infinity. So we see the two dominant open loop poles at $s=-0.22$ start out going up and down in the positive and negative imaginary direction, circle around and meet in the right half plane where one turns and terminates on the zero at $s=0.012$ while the other terminates in the zero at $s = +18$) (off the plot to the right). The other pole at $s=-45$ terminates in the zero shown at $s=-1.05$ while the remaining pole at $s=-1000$ terminates in the zero at infinity going in the direction of the negative real axis.

Here is the resulting root locus after using a compensator that uses two zeros at $s=-.22$ to cancel the two dominant poles and then a pole at $s= -1.05$ to cancel that low frequency zero and one more pole at $s= -2$ to terminate in the zeros in the right half plane (that can't be cancelled).  With a loop gain of 0.1164 the damping factor is 0.93 which provides for fastest rise time while keeping a reasonable overshoot (the distance from the real axis is the natural frequency while the distance from the imaginary axis is the decay time as previously described; by increasing the natural frequency we introduce some ringing into the loop but it will increase the rise/fall time by getting to the target level quicker--if we overdue that then we get excessive ringing/overshoot which is not desirable).

With the resulting comparative step response (this has a very large negative excursion but will be fully settled in less than the 3 seconds desired by the OP:

This solution using pole-zero cancellation will be very sensitive to the accuracy of the cancellation; any small error will result in a residual transient that will take a long time to settle.


How to find the color channel shifting of an RGB image that was supposed to be greyscale?

It seems like conceptually you're nearly there:

Take one of the three colour channels as reference (for example: red)
Match the green against the red
Match the blue against the red

The problem then is the matching with sub-pixel accuracy. To get
meaningfull numbers is a bit tricky, since the interpolation errors
any method will make have an impact on the accuracy.

There are two useable image registration methods I can think of:

Lucas-Kanade image registration. Using linear interpolation might
not give you accurate enough results, so consider bicubic or other methods.
Neil Dodgson has one made a nice overview. It's important that different
sub-pixel shifts of the interpolation kernel have a similar frequency transfer.
For the cubic family, the approximating b-spline is far better than cattmull-rom, in this case.
First upscale the image, then do pixel-precise image registration using your favorite
technique. (cross-correlation in the Fourier domain should do).
This only works if upscaling is done carefully. Bilinear or bicubic will most likely not give you enough accuracy. I can think of three ways:

a. Yen-interpolation. See equation 11 of his paper. Really slow, but optimal. Although you seem to have some alias in the signal so the 'bandlimited' assumption Yen makes might not hold.

b. Calculate FFT of the image, zero pad the high frequencies, inverse FFT.

c. Non-linear upscaling. Since the edges are quite sharp, the image is not properly bandlimited, that might be the main limitation of the previous methods. Edge Dependent Directional Interpolation might be better in this case.

Once you have the sub-pixel translation between the two, the correction to the green and blue channel is already solved, no matter which method you choose.


Power Spectral Density computation and units

I don't know if I get your question, but here is something you can look into.

Fourier transform is basically converting a time domain function to frequency domain. That being said, Fourier transform basic formula gives you the amplitudes of various components of a given signal in frequency domain. Note the word amplitude. That's what Fourier transform provides you, amplitudes. These amplitudes are functions of, of course, frequency. By exponential Fourier transform formula, we get these amplitudes as complex quantities. That might be your $g(t) \iff G(f)$, where $g(t)$ gives the amplitudes of the signal vs. time, and $G(f)$ gives the amplitudes of the different components of the signal vs. their individual frequencies.

Now you are interested in PSD. That's power. So now you don't want amplitudes vs. frequency, you want power vs. frequency, right? Power is proportional to the square of amplitude. So how do you convert amplitude to power? You square the amplitudes.

But, amplitudes occur as complex quantities in exponential formula, remember? So if you multiply one complex quantity with its conjugate, you get the square of its amplitude. That's why you might use

$$S_g(f)=\lvert G(f)\rvert^2=G(f)G(-f)$$

Here, $G(-f)$ is the conjugate of the complex quantity $G(f)$. At $f = 0$, you of course get a real quantity. In that case, you have,

$$G(f)G(-f)=G(f)G(f)=\lvert G(f)\rvert^2$$

Hope that helps.


As shown in other answers, one way to solve such problems is using the $\mathcal{Z}$-transform. However, your approach of finding a general and a particular solution also works, but you've got a few things wrong. Since this looks like homework, I'll just try to give you a few essential hints to enable you to solve the problem yourself (and to actually learn something).

You seem to have found the characteristic polynomial and its solutions:

$$\lambda^2-\lambda+\frac14\lambda=\left(\lambda-\frac12\right)^2=0\tag{1}$$

The characteristic polynomial has a double root at $\lambda=\frac12$, which must be taken into account when forming the general solution. That's where you've gone wrong. You treated the double root as two distinct single roots, which doesn't work. The solution of the homogeneous difference equation ($x[n]=0$) has the form

$$y_h[n]=C_1\left(\frac12\right)^n+C_2n\left(\frac12\right)^n,\qquad n\ge 0\tag{2}$$

A particular solution (taking into account the excitation $x[n]=u[n]$) can be found by inspection:

$$y_p[n]=C_3\tag{3}$$

First determine $C_3$ by simply plugging the particular solution $(3)$ into the difference equation. The final solution is then given by

$$y[n]=y_h[n]+y_p[n]\tag{4}$$

The constants $C_1$ and $C_2$ must now be chosen such that the given initial conditions are satisfied.


When i have already found beamforming,do i still need to use MRC?i have some problems for them

You are doing transmit beamforming and choosing your vector $\mathbf{f}_A$ so that you beamforming in the "best" direction. So you transmit the signal: $\mathbf{z}=\sqrt{P} \mathbf{f}_A x$, where $\mathbf{z}$ is a length $N_T$ vector. Now you receive the signal: $\mathbf{y} = \sqrt{P}\mathbf{H}\mathbf{f}_A x + \mathbf{n}$, which is a length $N_R$ vector.

The post-MRC SNR should still be a vector. You will have a SNR for each transmit stream, so $N_T$ SNR values. To compute the SNR of the $k^{\text{th}}$ transmit stream after MRC you take the channel gains from transmit antenna $k$ to each of the $N_R$ receive antennas, $\mathbf{h}_k$, and compute its output power: $\gamma_k = P ||\mathbf{h}_k^H \mathbf{H} \mathbf{f}_A||^2$. Or, you can also equivalently do it all at once by: $||\mathbf{H}^H \mathbf{H} \mathbf{f}_A||^2$.

Your formula is actually computing the pre-MRC SNR, whereas your friends was showing how to get the post-MRC SNR. Hope this helps!


Halving a floating point time series and it's effect on noise levels

...why does the graph show an increase in noise throughout the data?

It does not. The graph is a depiction. You should focus on the numbers themselves.

Your silence segment includes at last three relatively loud spikes.

Spikes are sharp and brief transitions of the amplitude in the time domain. The bandwidth of a spike is inversely proportional to its duration. That is, the shorter the pulse, the more sinusoids it takes to synthesize it. A quick way to show this is via the Discrete Fourier Transform of the sinc pulse.

Spikes and more generally impulse-like waveforms, manifest themselves in a spectrogram as vertical "line" structures. There are three, relatively clear vertical "lines" in your spectrogram towards the end of the first plot.

The point is that the spike is a sharp transition that stands out from its surrounding samples (in the time domain) and causes relatively large sums in the frequency domain.

This is what is picked up by your normalisation by np.max and is made worse by your choice of y_axis='linear'.

If you try to depict this sequence of numbers x=[8,3,6,4,7,6,120,5,3,6,5,3,2] then the visualisation algorithm will assign the colour at the top of the color map to 120. Effectively, you would get one really bright value which would "swallow" all the rest. If you were to depict the same sequence of numbers but instead of 120 you substituted with 6.5, then you would observe a much more "colourful" output.

This is what is happening in this case. When you plot the whole file, you include the pulses that "skew" your visualisation (but not the underlying numbers) but this is not happening when you simply use the first half of the recording.

To deal with this, you can either use a simple logarithmic transform like log(1+x) or simply switch your scale to logarithmic (i.e. Decibel).

Both will have the same effect: Large values will remain large but small values will be "boosted". This does not distort the number. It simply changes the way the number is presented.

This looks like the beginning of a vinyl recording which includes a bit of hiss along with the occasional pop, because of depositions on the record (?)

Hope this helps.


trying to understand fast fourier transform in matlab

Alternate simplistic answer to (2):

One often has to cut off a signal at a beginning and an end to fit it into any finite length FFT.  The sharp edges potentially produced by these two cuts (a rectangular window) in a longer time domain signal can result in various high and low frequency artifacts in the FFT result (especially for any signals or signal components that are not exactly integer periodic within the FFT's length).  A Hamming window cuts a longer signal into the FFT's length using a softer rounder (raised cosine) edge, thus producing far less high frequency artifacts in the FFT result.

Alternate detailed answer to (1):

Line 2: A standard FFT returns a complex result.  Many uses do not require a complex number result (with both real and imaginary components), so Line 2 uses the abs() function to return the magnitude value of each FFT result bin, which are strictly real valued (no imaginary stuff).

Line 3: An FFT of a strictly real input returns a result where the upper (or negative) half is identical in magnitude to the lower half (except for the DC and Fs/2 or Nyquist bin).  So line 3 just throws away all that redundant information (it's the same as what's kept, why repeat it?).

Line 4: However, if you throw away half the FFT result, the total energy gets reduced by half.  So Line 4 multiplies by 2 to fix that (everything except the DC and Fs/2 bins), and get the correct total energy.


What is the name of this transform similar to the Radon transform?

This is a simple Cartesian to Polar coordinates transformation. Instead of expressing the position of a pixel using $(x,y)$ you are expressing it using $(r, θ)$ where $r$ is distance from the centre $(x_0,y_0)$ and $θ$ is angle.

The images are exactly the same, wherever they are definable, but the "unrolling" results from still using $(x,y)$ to depict $(r,θ)$.

For more information, please see this and this.


Creating a noise model with a fixed PSNR

From wikipedia:

$$
    PSNR = 20 \cdot \log_{10}(MAX_I) - 10 \cdot \log_{10}(MSE)
$$

So, you could:

To identify the max value of your image: $MAX_I$
Generate noise samples for $\sigma = 1$
Solve for $MSE$ in the equation above.
Scale all you generated noise samples by the same factor to achieve that $MSE$ value.


Calculating the system output using frequency response

Yes, you are correct.  The system input is rather strange, though, since $cos(6\pi n + \frac{\pi}{6})$ is equal to $cos(\frac{\pi}{6})$.


What is localizability in Computer Vision?

In the robot navigation problem, the localization problem refers to the real time estimation of its position and orientation under various backgrounds. This is usually achieved by some natural landmark selection (laser points, camera views, etc.), and the features in the image (corners, tiny lines with different orientations, etc.). So the localizability reflects how easy, robust and consistent the robot's position and orientation can be estimated with the choice of feature points.

Hugues Sert et al's Mobile robot algebraic localizability includes more strict definition with details, and its relation to observability.


DFT for something similar to convolution

You might be looking for the convolution property of the DFT. If $z[n] = x[n] * y[n]$, where $*$ represents circular convolution, then $Z[k] = X[k] Y[k]$, where $Z[k]$ is the DFT of $z[n]$ and so on.

This property is used in fast convolution techniques like overlap-save to implement convolution in $O(N \log N)$ time instead of the $O(n^2)$ time required by direct implementation of the convolution sum.


Is what I'm after even possible?

yes.

Let's start with the simplification your camera was a perfect pinhole camera.

For starters, imagine you pick an arbitrary point object in your scene. You then get a camera matrix $\mathbf C$ that allows your to predict the 2D position $\mathbf y$ in your image of your 3D-positioned object at $\mathbf x$:

$$
\mathbf y = \mathbf {Cx}, \quad \mathbf y= \begin{pmatrix} y_1\\y_2\\1 \end{pmatrix}\in\mathbb R^3,\,\mathbf x= \begin{pmatrix} x_1\\x_2\\x_3\\1 \end{pmatrix}\in\mathbb R^4,\,\mathbf C \in \mathbb R^{3\times 4}
$$

Notice that $\mathbf C$ isn't reversible, ie. there's multiple $\mathbf x$ that explain a $\mathbf y$.

Now, translating your camera by a movement vector $\mathbf m$ is the same as shifting the scenery by $\mathbf{-m}$, which leads to a $\tilde{\mathbf{x}}= \mathbf x - \mathbf m$.

Correspondingly, we get a

$$\tilde{\mathbf{y}} = \mathbf C\left(\mathbf x - \mathbf m\right)\text.$$

Now, since $\mathbf C$ is a matrix, i.e. a linear operator, this leads to

$$\begin{align}
\tilde{\mathbf{y}} &= \mathbf C\mathbf x - \mathbf C \mathbf m\\
 &= \mathbf y -\mathbf C \mathbf m\quad \iff\\
\mathbf C\mathbf m &= \mathbf y - \tilde{\mathbf y}\text.
\end{align}$$

As mentioned above, $\mathbf C$ isn't reversible, which means that from a single photo you can't tell the 3D coordinate of a point without additional info. However, if you can identify three  objects at different points $\mathbf y_{a,b,c}$ and $\tilde{\mathbf {y}}_{a,b,c}$, you should be able to elimate the unknowns from the system of equations.

Also note that your real camera isn't a perfect pinhole camera, and that might be used to actually get a less underdefined system, so you might be able to mathematically derive a unique solution with less than three points.


Difference Equations, why M<=N for causality?

You're right, it's simply not true that $M\le N$ is necessary for causality. The difference equation in your question can always be implemented by a causal system.

However, note that the difference equation is not uniquely related to a causal system. A simple example is

$$y[n]=ay[n-1]+bx[n]\tag{1}$$

This can obviously be implemented by a causal system. But for $a\neq 0$ you can also rewrite (1) as

$$y[n-1]=\frac{1}{a}(y[n]-bx[n])\tag{2}$$

which suggests an anti-causal system, even though (1) and (2) are completely equivalent. This difference is reflected by the transfer functions ($\mathcal{Z}$-transforms) of the corresponding systems. For a causal system the region of convergence (ROC) is outside a circle enclosing all poles, whereas for an anti-causal system the ROC is inside a circle outside of which all poles are located.


Unable to understand the derivation of the update equation for LMS

so with an LMS filter, we have a time-variant $N$-tap FIR filter:

$$ y[n] = \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] $$

$x[n]$ is the input signal, $y[n]$ is the FIR output, and $h_n[k]$ are the FIR tap coefficients at the time of sample $n$.

with an LMS filter, we also have another input called the desired signal: $d[n]$.  we want our LMS filter to adapt so that the output $y[n]$ follows the desired signal $d[n]$.  from those two signals, we define an LMS error signal:

$$ e[n] \triangleq y[n] - d[n] $$

and this Least Mean Square adaptive filter will adjust its tap coefficients, $h_n[k]$, in such a manner as to eventually minimize the square of this error signal $e[n]$.

so the error signal is:

$$\begin{align}
 e[n] &= y[n] - d[n] \\
      &= \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \\
\end{align}$$

and the square of this error signal (which we want to lessen in an LMS filter) is:

$$ \big(e[n]\big)^2 = \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right)^2 $$

now, one thing Bernard Widrow must have been thinking is that, since the inputs $x[n]$ and $d[n]$ are given, the only controls we can adjust to lessen this square error $e^2[n]$ are the $N$ tap coefficients $h_n[k]$ for $0 \le k \le N-1$.  with $x[n]$ and $d[n]$ fixed, then $e^2[n]$ is a function, $f(\cdot)$, of the $N$ tap coefficients $h_n[k]$.  so at time $n$ this function of $N$ variables looks like

$$\begin{align}
 e^2[n] & \triangleq f\big(h_n[0], h_n[1], h_n[2], ... h_n[N-1] \big) \\
        &= \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right)^2 \\
\end{align}$$

at time $n$, the derivative of $f(\cdot)$ with respect to the $m$-th tap coefficient is

$$\begin{align}
 \frac{\partial}{\partial h_n[m]} \, e^2[n] &= \frac{\partial}{\partial h_n[m]} \, f\big(h_n[0], h_n[1], h_n[2], ... h_n[N-1] \big) \\
    &= \frac{\partial}{\partial h_n[m]} \, \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right)^2 \\
    &= 2 \, \left( \sum\limits_{k=0}^{N-1} h_n[k] \, x[n-k] - d[n] \right) \cdot x[n-m] \\
    &= 2 \, e[n] \cdot x[n-m] \\
\end{align}$$

now from this we can set all of these derivatives to zero and set up $N$ equations and $N$ unknowns and try to solve the mess.  i am sure Widrow did this.  but it doesn't lead to a very cheap (in terms of computational cost) algorithm.

but for each $m$, where $0 \le m \le N-1$, the quantity

$$ \frac{\partial}{\partial h_n[m]} \, f\big(h_n[0], h_n[1], h_n[2], ... h_n[N-1] \big) = 2 e[n] x[n-m] $$

is one component of the vector which is the gradient of $f(\cdot)$.

for the moment, try to visualize this as a function of two variables: $f\big(h_n[0], h_n[1]\big)$.  this is a surface above the point in a plane with coordinates $(h_n[0], h_n[1])$.  the gradient

$$\begin{align}
 \left(\tfrac{\partial f(\cdot)}{\partial h_n[0]}, \ \tfrac{\partial f(\cdot)}{\partial h_n[1]}\right) &= \Big(2 e[n] x[n-0], \ 2 e[n] x[n-1] \Big) \\
&= 2 \, e[n] \Big(x[n], \ x[n-1] \Big) \\
\end{align} $$

is the 2-dim vector that points in the direction of the steepest increasing slope of this surface.  so if $\tfrac{\mu}{2}$ is a tiny positive number, moving a tiny bit in the direction of the gradient should increase this function a tiny bit.

$$ f\big(h_n[0], \ h_n[1]\big) < f\big(h_n[0] + 2 e[n] x[n] \tfrac{\mu}{2}, \ h_n[1] + 2 e[n] x[n-1] \tfrac{\mu}{2} \big) $$

but $f(\cdot)$ is the square of error and we want to decrease it, so we move in the opposite direction of the gradient:

$$ f\big(h_n[0], \ h_n[1]\big) > f\big(h_n[0] - 2 e[n] x[n] \tfrac{\mu}{2}, \ h_n[1] - 2 e[n] x[n-1] \tfrac{\mu}{2} \big) $$

that's with two taps ($N=2$).  extending the idea of this "surface" to more dimensions than two, we expect

$$ f\big(h_n[0]-e[n]x[n]\mu, ... ,h_n[k]-e[n]x[n-k]\mu, ... \big) < f\big(h_n[0], ... ,h_n[k], ... \big) $$

so, for a small stepsize $\mu$, the adaptation algorithm (or coefficient "update equation") is

$$ h_{n+1}[k] \ \leftarrow \ h_n[k] \ - \ (\mu \, e[n]) \, x[n-k] $$

that is the simple unnormalized LMS adaptive filter.  you want $\mu$, also known as the adaptation gain, to be large enough so that the LMS filter adapts at an acceptably quick rate, but small enough so that it doesn't overshoot in the adaptation and become unstable.

because the stepsize is really $(\mu \, e[n]) x[n-k]$ and both $e[n]$ and $x[n-k]$ are expected to be proportional in magnitude to the mean $\big|x[n]\big|$, this means louder signals will cause the adaptation to happen faster than quieter signals.  so from this, the Normalized LMS adaptive filter is defined to divide that term by some measure of the mean of $\big|x[n]\big|^2$.

so the Normalized LMS adaptation algorithm (or coefficient update) is

$$ h_{n+1}[k] \ \leftarrow \ h_n[k] \ - \ \frac{\mu \, e[n]}{\overline{x^2[n]} + \epsilon} \, x[n-k] $$

where $\overline{x^2[n]}$ is some measure of mean of $x^2[n]$ for the samples most recent to time $n$.  you can simply square $x[n]$ and run that through a simple first-order IIR filter to get a measure of the mean:

$$ \overline{x^2[n]} \ \leftarrow \ p \, \overline{x^2[n-1]} \ + \ (1-p) x^2[n] $$

where $p$ is the pole of the one-pole IIR filter and $0 < (1-p) \ll 1$.  the quantity $\epsilon > 0$ is a tiny amount to add to the denominator so that we never divide by zero, even if the input $x[n]$ goes to zero.


Some problem with the expression of QAM basis and signal

Base choices are ambiguous. The $-$ is there because the author of the book you're reading decided it was better that way.

You'll probably find the cause for that in the quadrature mixer graphic that they use, where if the phase shift from $I$ to $Q$ is $+\frac\pi2$, since $\cos(x+\frac\pi2) = -\sin(x)$.


The power of a discrete-time signal $x[n]$ is given by

$$P_x=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{n=-N}^{N}|x[n]|^2$$

which is identical to the first formula in your question. The second formula is wrong and does not even make any sense due to the infinite limits of the sum.

I'd advise you to use text books or high-quality internet resources such as MIT OpenCourseWare.


Derivative with respect to complex conjugate

That's a trick which you will also find in a DSP context, that's why I choose to provide an answer here. It is related to the Wirtinger derivative, and you can find more details about it in this answer over at math.SE.

In practice this trick is often used to compute the extremum (minimum or maximum) of a real-valued function depending on a complex variable (or vector/matrix). In order to find the extremum, you formally take the derivative with respect to the complex conjugate of the variable of interest, set this derivative equal to zero, and from this equation derive the optimum value of the (possibly vector-/matrix-valued) variable.

As a simple example, take the minimization (with respect to the vector $x$) of the following error quantity:

$$\epsilon=||Ax-b||_2^2\tag{1}$$

where the vectors $x$ and $b$, and the matrix $A$ are assumed to be complex-valued. The error can be written as

$$\epsilon=x^HA^HAx-x^HA^Hb-b^HAx+b^Hb\tag{2}$$

where $^H$ denotes the Hermitian (conjugate) transpose. Now you take the derivative of $(2)$ with respect to $x^H$, formally treating this vector as being independent of $x$. The result is

$$\frac{\partial\epsilon}{\partial x^H}=A^HAx-A^Hb\tag{3}$$

Setting this derivative equal to zero and assuming that $A^HA$ is invertible, the optimal value of $x$ is given by

$$x_o=(A^HA)^{-1}A^Hb\tag{4}$$

where the matrix $(A^HA)^{-1}A^H$ on the right hand side of $(4)$ is the Moore-Penrose inverse.


How to find DC bias of periodic signal from samples that include a non-integer number of periods?

Low-pass filter with DC gain = 1 (or 0 dB).

i s'pose, if you want to be sophisticated and want a rock-solid DC that might vary very little, then the output of your LPF can go into a median filter.


identification of sounds within audio signals

Audio Fingerprinting can be used to match near identical recordings of sounds with a bit of additive noise. It will not work for similar sounds, cover versions etc.

Without knowing your application, I'd recommend giving Dan Ellis' audio fingerprinting code a try. You'll find MATLAB and Python versions here:

http://labrosa.ee.columbia.edu/matlab/audfprint/


Are discrete systems defined by LCCDE always LTI?

What you have is a so-called behavior which truncates the signal space $\mathcal{X}\oplus\mathcal{Y}$ of all possible signal pairs $(x, y)$ to the ones that satisfies this equality.

This does not specify an input-output mapping and hence a dependence and hence causality does not apply here. They simultaneously satisfy this. However regardless of causality, it is LTI since you can take arbitrarily and superposition applies.

Just like you can integrate any differential equation sufficiently many times and end up with a integral equation, you can also forward time shift these equations then if one of them depends on future inputs and that is your choice of dependent variable then you can speculate about causality.

This is actually the tool to allow for nondifferentiable input signals and other discontinuities and the results are investigated under the name of weak solutions


Discrete time sine wave generation near nyquist

There's nothing wrong about this – it's absolutely what you should be seeing.
The sine wave is intact, if you applied a reconstruction filter – there's no amplitude fluctuation of the wave, as soon as you limit its bandwidth to the nyquist bandwidth, the analog signal "has" to follow the full sines.


Doppler Shift of Phase Coded Pulse Compression Waveform

So this is something we run into with the radar community on a pretty regular basis.  For LFMs, the solution is a bit trivial, but for coded waveforms it becomes a real issue.  For us radar folk, we're dealing with targets that are at escape velocity (right around 11 km/s typically). I actually have some pending intellectual property being released via my company on this subject, so I'll try to discuss the solution without going too far into it.

Depending on the direction of the velocity vector, the Doppler distortion on an LFM will result in either an elongation of the waveform, or a "squishing" of the signal.  It's important to note here that for the LFM case, this pretty much just results in an effective chirp slope change, so I like to think of it as the start/stop frequencies have remained the same, but the effective pulse duration it took to get there has been distorted.

The solution for an LFM is simply just to account for the effective slope change that you expect. For example, if you're expecting an elongation by a factor of 2, "squish" your waveform by a factor of 2 when you transmit the waveform; we call this process pre-distortion, and the goal is to equalize the distortion effect from the target such that at your pulse compression filter, you get maximum pulse compression gain.  This of course relies on some a priori knowledge of the target velocity, which isn't totally unreasonable given how most of these systems perform tracking operations.

Coded waveforms present an issue though, we're now non-linear, so the solution isn't as easy on the surface.  If you chunk through the algebra a bit using generic phasor notation, you'll start to see some terms pop out.  It turns out the result is fairly similar: the chip widths change.  How the chip widths change, I'll leave that for you to figure out with your algebra.

The important thing to note here though is that, at least in my experience, you do not punch holes/zeros in your signal like you've asked.  The coded waveform chips will now either be shorter or longer, so you may have added or lost some phase information, but if you can account for that on transmit just like we did with the LFM, you can negate the distortion via pre-distortion of the waveform.  In short: the chips will change their duration, but there will not be gaps since all of the chips change.

Now as for any Doppler shift that's on the phase coded waveform, that can easily be taken care of if you have some expectation of target velocity: you simply just shift the pulse compression filter in frequency (or the received waveform itself, whichever is easiest for you to implement in your system) so that the Doppler shift is "matched" to your pulse compression filter.  Again, this requires some knowledge of target velocity.  If you didn't have knowledge of target velocity, you could do a multi-hypothesis approach.

The reason we typically don't have to worry about this for LFMs is evident in their ambiguity function: they are incredibly Doppler tolerant, i.e. Doppler shifts just cause slight range migration and decrease in pulse compression gain, but not signal loss.  With coded waveforms, their ambiguity function typically looks like a spike (often called "thumb tack ambiguity"), which results in complete signal loss with even the slightest Doppler shift.  Therefore, it's pretty important that we get that Doppler hypothesis correct!

This method is in open literature, and I'm pretty sure some of the Mark Richards signal processing books touch on it. This paper from MITRE touches upon a multi-hypothesis technique.  A lot of that paper is poorly written and sometimes incorrect in my opinion, but they have a good explanation of the multi-hypothesis approach, and a nice block diagram for it towards the end.


Analog QAM as an example:

https://en.wikipedia.org/wiki/Quadrature_amplitude_modulation#Analog_QAM

same idea for other IQ modulations schemes.


DWTs turn a signal into a list of coefficients, whose combined sizes are equal to or greater than the signal length. Coefficients are grouped in packs called subbands. Each subband gathers coefficients resulting from some band-pass filtering, followed by a subsampling.
They are not exactly narrow-band signals, but projections onto embedded subspaces. The "narrowband" signal for each subband could be recovered by inverse transforming those coefficients only.


Not ALL DACs (D/A converter) do this, but a traditional DAC operates on the digital samples as a zero-order hold; each sample is held in the analog at a constant value until the next sample. This process is identical mathematically to convolving the discrete samples (as weighted impulses) with a rectangular function of width T where T is the sample time.

The Fourier Transform of a rectangular function of width T is a Sinc function with first null at 1/T, where 1/T is the sampling rate.

Convolving in the time domain is equivalent to multiplying in the frequency domain, so the spectrum of your discrete signal is multiplied by the Sinc function in the process of analog reconstruction in the DAC due to the zero order hold.

There are other DAC implementations that reduce this effect such as interpolating DACs and "Return to Zero" (R2Z) DACs. An R2Z DAC forces the output to return to zero prior to the next sample, reducing T and therefore pushing out in frequency where the first null occurs (therefore reducing passband rolloff).