Does the 'FFTW' have a cost of k*N*log_2(N)?

You are trying to describe a complex phenomenon with a simple formula. That gets you only a rough approximation or an estimate but depending how much accuracy you need, you will have to dive deeper into the specifics. A few points

Any actual operation count will depend heavily on the hardware: Instruction set, pipeline, memory architecture, cache strategies, hardware accelerators, etc. In many cases the choice of the "best" algorithm is actually hardware dependent. Your best shot here is benchmarking and profiling on the target hardware.
The good old N/2*log2(N) formula basically describes the "number of butterfly operation" in a standard Cooley Tuckey FFT. So a 32-point FFT has 16*5 = 80 butterfly operations. There are 5 stages of 16 butterflies each. This is typically (but not always) the bulk of the work so it's a decent estimate for the overall complexity.
One butterfly has one complex multiply and two complex adds (or on
add and one subtract if you are picky about). The number of operations required will depend widely on the Hardware. An Analog Devices Sharc SIMD processor can do this with about 3 operations, an general purpose processor might take up to 20 operations.
There is a lot of other stuff that isn't part of the butterflies: bit-reversing, looping & branching, coefficient calculation storage and retrieval. These can factor heavily into your MIPS count as well and are even more hardware dependent.
The first two stages (for decimation in time) use only simple coefficients like +1, -1, +j, -j and don't need a full complex multiply. Most algorithms take advantage of this and so the first 2*log2(N) butterflies are a little cheaper. By how much depends again on the hardware.
Non power of two FFTs get more complicated. Let's look at a 96 point FFT. First you need 3 FFTs of 32 points each, so that is 3*5*16 = 240 butterflies. Then you need to assemble the result which takes 2 complex multiplies and 2 complex adds per sample, so it's another 96*2=192 extra complex multiplies and 192 complex adds. Total would be 432 mpy and 672 adds. A 128 FFT would have 448 mpy and 896 adds, so it may be better just to zero pad and use the larger power-of-two FFT as it greatly simplifies coefficient calculation and storage.
The more prime factors you have in the FFT length the more complicated and less efficient the algorithm gets.
In your specific case, I'd zero pad the 1920 to 2048 (unless you need true circular convolution, which I don't believe you do here). It's highly unlikely that any prime factorization (128*3*5) would do better. The 1080 is awkward, since it's just a tad over 1024. I'd probably zero pad to 1280 and do a 256*5 or zero pad to 1152 and do 128*9. You can also try the 8*27*5 and see if that's better.


Spatial and temporal frequency of a wave

The two frequencies you are referring to are the spatial frequency and temporal frequency of the wave, and you are correct in your reasoning on converting one to the other. The spatial frequency refers to how many complete periods the signal goes through for a given unit of distance (eg. cylcles/m) while the temporal frequency refers to how many complete periods the signal goes through for a given unit of time (eg. cyles/s or Hz).

The physical meaning of the two frequencies should be fairly easy to understand if you imagine the wave you are sampling is a simple sine wave. Then it is easy to see that for any constant point in time, your signal will be a frozen sin wave which is periodic in space, and thus has a "spatial frequency". Conversely, for any constant point in space, the value of the signal will oscillate over a regular period in time, and thus have a "temporal frequency" component. You can also imagine that the faster the wave was travelling through space, the faster the signal would change at any one point in space, giving rise to a higher temporal frequency and underlining the relationship between the two properties via the wave's speed.

Regarding using all the sample points to estimate the wave frequency, I believe you should be able to use the known speed to scale one of your dimensions to match the other, since $distance=velocity*time$. Then a 2D FFT of the matrix should give you information based on all of the available samples in whichever dimension you choose.

This article gives a brief description of the relationship between the two in the context of radar imaging and this article talks about the use of the 2D FFT for analyzing higher dimensional signals.


Converting a lowpass filter to a highpass filter. FIR filter-type 1

If you have a real-valued amplitude function $A_{LP}(\omega)$ of a low pass filter with unity gain in the passband, then it is easy to see that
$$A_{HP}(\omega)=1-A_{LP}(\omega)\tag{1}$$
is the (real-valued) amplitude function of a high pass filter.
A type I length $N$ linear phase low pass filter's frequency response is given by
$$H_{LP}(\omega)=A_{LP}(\omega)e^{-j\omega(N-1)/2}\tag{2}$$
Consequently, according to $(1)$, the desired high pass amplitude function is obtained by computing
$$\begin{align}H_{HP}(\omega)&=\big(1-A_{LP}(\omega)\big) e^{-j\omega(N-1)/2}\\&=e^{-j\omega(N-1)/2}-H_{LP}(\omega)\tag{3}\end{align}$$
The inverse discrete-time Fourier transform of $(3)$ is
$$h_{HP}[n]=\delta\left[n-\frac{N-1}{2}\right]-h_{LP}[n]\tag{4}$$
Note that for type I linear phase filters the filter length $N$ is odd, so $(N-1)/2$ is an integer shift. Eq. $(4)$ shows that apart from sign inversion we only modify the center bin of the low pass filter in order to obtain a type I linear phase high pass filter.


What is Wavelet packet transform (WPT) time complexity

I just found the answer in a paper. so thought it might be useful for someone.
If we decompose a signal with length of $N$ to $L$ level then time complexity is
$O(L\times Nlog_{2}N)$


Why is the dB plot of a signal usually non-positive?

In power spectra of signals, 0 dB is some agreed power level. It can for example represent the level of a full-scale (FS) sinusoid. In that case it would be rare to see values above 0 dBFS (sinusoid), as it would require peak cancellation as with $1.125\sin(x) + 0.125\sin(3x)$ to stay within a numerical amplitude range of $[-1, 1]$.  When a power spectrum is given in dB, the 0 dB reference should be specified. Otherwise one can't be sure what it means.

You saw the above 0 dB values not with a power spectrum of a signal but with a transfer function, or equivalently with a magnitude frequency response. In a magnitude frequency response, 0 dB means a magnitude or gain of 1, that there is neither attenuation or amplification. If you followed this standard convention and saw values above 0 dB, it means that the corresponding frequencies are amplified.


How do Anti-Aliasing-Filters in audio signal processing work?

I think you make this needlessly complicated. Let's do a specific example: let's say we want to go from 48 kHz up by 4/3 to 64 kHz and your signal bandwidth is 20 kHz.

In the first step, you would up-sample by a factor of 4 by inserting 3 zeros between each sample. Your new sample rate is 192 kHz. Your original spectrum is preserved and you get mirror spectra centered around 48 kHz & 96 kHz.

To eliminate this you need a lowpass filter with cutoff of 20 kHz at 192 kHz sample rate and sufficient attenuation at 24 kHz.

Finally you just down sample the result by factor of three by simply throwing away every second and third sample. This works without aliasing, since the lowpass filter has already eliminated all content that would alias.

Regarding your graphs: in general I would avoid drawing anything past the Nyqusit frequency. It's just periodic repetition and doesn't add any useful information and it can get confusing.


RF for Data Transmit and Receicve for Beginners

A possible combination is a copy of Digital Communications by Proakis & Salehi along with GNU Radio.

With SDR, you actually create / process the signals that carry your data over the RF medium.

Therefore, to an extent, "freshman level theory" is all you are stuck with if you have to write an end-to-end SDR system.

But, you could keep the book (or, a book) as a reference for points you don't understand when going through something like GNU Radio's PSK Demodulation resources for example.

I am not sure what your expectation is from SDR, but taking the comment (about SS7) as a starting point, it's not like the SDR has an API which you use to simply "send data". Although, I am sure that you could add that layer too once you have your basic DSP dataflow set up.

The beauty of SDR is that you take full control of every building block from transmission to reception. So, there is nothing stopping you from writing a Phase Locked Loop (which you are going to need) from first principles or use one of the existing ones (or here).

But you still have to have some awareness of what is the "loop bandwidth" for example, in the case of the PLL.

Hope this helps


Retrieve frequency value of an 2D frequency spectrum

You have a signal in a 2D-frequency domain. You know the maximum frequencies in both directions, I assume. Hence, you know which bandwidth is represented by one pixel and then... just count the pixels from the center to the point of interest along the axes and your done. Depending on your application, taking the Euclidean distance from the center to the peak is also an option, which does end in one instead of two frequencies.

If you would like to do a more sophisticated approach, you could fit a Gaussian distribution (or another function, if you have one in mind that is more appropriate in your case) to the peak. That way you could get a sub-pixel accuracy of the location of the peak.


My favorite "Rule of thumb" for the order of a low-pass FIR filter is the "fred harris rule of thumb":
$$
N=\frac{f_s}{\Delta f}\cdot\frac{\rm atten_{dB}}{22}
$$
where

$\Delta f$ is the transition band, in same units of $f_s$
$f_s$ is the sample rate of the filter
$\rm atten_{dB}$ is the target rejection in dB

For example if you have a transition band of 100 Hz in a system sampled at 1 kHz, and your rejection requirement is 50 dB in the stop band then the order can be approximated by:
$$
N = \frac{1\ \rm kHz}{100\ \rm Hz} \cdot\frac{50}{22} = 23\ \rm taps\tag{rounding up}
$$
Thank you Fred Harris!
Note another more detailed formula that takes into account passband ripple is Kaiser's formula thanks to James Kaiser from Bell Labs, which I included in my graphic below.
For most applications I have done, the Fred Harris approach has been just fine, as given a certain rejection, the resulting filters using traditional filter design algorithms such as Parks-McClellan and Remez have exceeded my passband ripple requirements upon meeting the rejection requirement.  (What I usually do is estimate the order, design the filter with that order, inspect the result and increase or decrease order from there to fine tune). The results from the estimates are just that: estimates, and can vary widely depending on overall design parameters and not to be assumed to be an exact solution.

For those familiar with filter design using window approaches, review of the boxcar or rectangular window (which is simple truncation) reveals why it takes $f_s/\Delta f$ taps (which is the same as $2\pi/\Delta\omega$ if the units for normalized frequency are radians/sample as is often done) to complete the transition band. See the images below that help explain this.
The top image below shows the expected Sinc in frequency for a rectangular window in time, in this case as a non-causal rectangular pulse centered on $t=0$. This is then repeated in discrete forms as a causal waveform starting at $t=0$, both with the Discrete Time Fourier Transform (DTFT) and Discrete Fourier Transform (DFT) where the difference is the samples in time extend to $\pm \infty$ for the DTFT resulting in a continuous waveform in the frequency domain. In both cases the result is an aliased Sinc function that is periodic over the interval $f= [0,f_s)$, with the key point that for $N$ samples in time of the rectangular function, the frequency response will have its first null at $f = 1/N$ (Where $f$ is the normalized frequency with 1 being the sampling rate).

This next image below shows the rectangular window approach to filter design (which I would never recommend but is informative). The first plot in the upper left hand corner shows the target frequency response for our filter as an ideal "brick wall" response. Please don't confuse that with the "boxcar window" (or "rectangular window") which is also a rectangular shape-- the window is in the time domain!
In order to realize such a filter, we would use the impulse response of the desired frequency response as the coefficients in our FIR filter (the coefficients of the filter are the impulse response---put an impulse in and out come all the coefficients!).  The impulse response for a rectangular frequency (brickwall) response is the inverse FT which is a Sinc function, in the time domain, shown in the lower left hand corner as the "Required Impulse Response". A Sinc function extends to plus and minus infinity, so to actually realize such a filter we would need an infinitely long filter and it would have infinitely long delay. Obviously we can't do that, so we truncate the coefficients to something realizable. The longer the filter, the closer we approximate the ideal brickwall response, but also the longer the delay will be (and the more resources we would need in terms of the filter construction; more taps).
Truncating the impulse response in the time domain is mathematically identical to multiplying by a rectangular window in the time domain. (Note that the impulse response is also delayed by half the duration of the window in order for the system to be causal).  Multiplying in the time domain is equivalent to convolution in the frequency domain. The frequency domain (FT) of the impulse response prior to truncation is our original desired brickwall frequency response. The frequency response for the rectangular window is a Sinc function in the frequency domain.
So when we truncate the desired impulse response (multiply in time by a rectangular window), we convolve the desired frequency response with a Sinc function, resulting in an approximation of our target frequency response as shown in the upper right-hand corner of the image below.

A key take-away for Sinc functions in general is the first null is $1/T$ where $T$ is the duration of the rectangular function.  For a sampled system, the first null would be at $2\pi/N$ where $N$ represents the number of samples for the duration of the rectangular function.  In the images a normalized radian frequency is used for the frequency axis- (if that confuses you just know the $2 \pi$ is the radian frequency for the sampling rate).  So in the process of the convolution, the sharp brickwall transition spreads out and in this case goes to 0 (our $\Delta \omega$) over a frequency of $2\pi/N$!  So here $$N = 2\pi/\Delta \omega$$ and of course the filter is poor with sidelobes etc.  Note this: This transition from the Sinc function is the sharpest available for a given number of taps; it has the best resolution in frequency but the poorest dynamic range (rejection). Other window typologies (Blackman, Blackman-harris, Kaiser (my favorite), etc) will significantly improve dynamic range but always at the expense of the transition.

So from the above, we see the origin of the $2\pi/\Delta\omega$ that is used in the approximation formulas, and we also see why there is an additional multiplication factor increasing the number of taps above this for typical filter designs; the rectangular window would give us the best transition possible with $N$ taps where $N = 2\pi/\Delta\omega$ but have very poor rejection. More taps are used to smooth the time transition further beyond the sharp transition of the rectangular window, providing greater rejection at the expense of transition bandwidth.


Why does a pole close to the unit circle result in an enhanced Q-factor?

Suppose some $H(q)=\frac{A(q)}{B(q)}$ where q is some complex variable and $A,B$ are functions of $q$.

Whether in the s or z planes, to evaluate the magnitude of $H(q)$ at some $q$, you evalaute and sum all distances from $q$ to the locations of the zeroes (i.e. the magnitude of $A$) and similarly for $B$ and create the fraction above for that specific $q$.

When does $H(q)$ becomes big? When either its nominator is big...OR its denominator is very small. In fact, if the denominator was ever to become exactly zero, then $H(q)$ would "shoot to infinity". It would acquire a very large value.

When does $H(q)$ becomes small? When either its nominator is small...OR its denominator is very big. The smallest value $H(q)$ can acquire is zero.

  ...why exactly the poles give such a change in Q-factor?

Given the location of the poles in your diagrams, it is useful to see what happens as your $q$ (or a particular $z$, in the case of the z-plane), approaches the bearing of the pole. In other words, as the $z$, changes frequency and walks around the z-plane, it approaches the location of the pole. EXACTLY AT the bearing of the pole, the distance to the zero is maximum. BUT! the distance to the pole is at its minimum. Therefore, the denominator becomes small and the response "shoots to infinity".

When the pole is closer to the center of the circle, its distance to the particular $z$ value is large. The closer the pole moves to the edge of the circle the smaller its distance will be with a particular $z$ and the higher the output of $H$ will be.

To explain the "peakiness" (or "...Q-factor...") you have to look at the rate by which the distance to the pole is diminishing as $z$ approaches the bearing of the pole.

  ...a pole closer to the unit circle will on one hand mean an "stronger" presence of a certain frequency it won't be decaying exponentially. But on the other hand a pole on the circle would mean a marginally stable system. In other words this means having a frequence which is perfectly present results in a +/- unstable system. Isn't that contradictory?

When examining the response of a system via the z-transform, the "presence of a frequency" means that the system resonates at that frequency. At extreme conditions, the system is so marginally stable that if you happen to excite that frequency even by a tiny little amount, it will go into oscillations that might take a very long amount of time to die out. In other words, the system will "ring".

A prime example of this are 1-pole resonators. These are bandpass filters with very high Q-factor. If you pass audio through such a resonator, you will hear the filter ringing as it effectively resonates at its own frequency.

  ...How should I modify my pole zero plot if I'd like to render undershoot in stead of overshoot like here?

If you are asking how to generate a notch filter, instead of a resonator (i.e. a filter that "cuts" one specific frequency instead of boosting it), then you can simply subtract the output of the 1-pole resonator from the rest of the signal OR, you could try to create a "well". A region of the z-plane where the conditions are such, for your $H$, that it dips to zero at a specific bearing to a particular $z$.

Have a ponder about how that looks like and if you cannot "see it", have a look at this link.

Hope this helps.


When Two Sampled Sinusoidal Are Orthogonal?

If you have two sinusoidal length $N$ sequences

$$x_1[n]=\sin\left(n\omega_1+\phi_1\right),\qquad n\in[0,N-1]\\
x_2[n]=\sin\left(n\omega_2+\phi_2\right),\qquad n\in[0,N-1]\tag{1}$$

then for orthogonality we require

$$\sum_{n=0}^{N-1}x_1[n]x_2[n]=0\tag{2}$$

For the given sequences $(1)$ we obtain

$$\begin{align}\sum_{n=0}^{N-1}x_1[n]x_2[n]=&\sum_{n=0}^{N-1}\sin\left(n\omega_1+\phi_1\right)\sin\left(n\omega_2+\phi_2\right)\\=&\frac12\sum_{n=0}^{N-1}\cos\left[n(\omega_1-\omega_2)+\phi_1-\phi_2\right]-\\&\frac12\sum_{n=0}^{N-1}\cos\left[n(\omega_1+\omega_2)+\phi_1+\phi_2\right]\tag{3}\end{align}$$

where I've used

$$\sin(\alpha)\sin(\beta)=\frac12\big[\cos(\alpha-\beta)-\cos(\alpha+\beta)\big]\tag{4}$$

The two terms on the right-hand-side of $(3)$ have the form

$$\sum_{n=0}^{N-1}\cos(n\omega+\phi)=\cos\left(\frac{N-1}{2}\omega+\phi\right)\frac{\sin(\omega N/2)}{\sin(\omega/2)}\tag{5}$$

For large $\omega$ (i.e., $\omega$ close to $\pi$), the magnitude of the term $(5)$ becomes small. Assuming that $\omega_1+\omega_2$ is large (but not greater than Nyquist), people tend to neglect the last term on the right-hand side of $(3)$.

From $(5)$, the first term on the right-hand side of $(3)$ is zero if

$$\frac{\sin((\omega_1-\omega_2) N/2)}{\sin((\omega_1-\omega_2)/2)}=0\tag{6}$$

is satisfied, i.e., for

$$\omega_1-\omega_2=\Delta\omega=\frac{2\pi k}{N},\quad k\in\mathbb{Z},\quad k\neq 0\tag{7}$$

With $\Delta \omega=2\pi\Delta f/f_s$ we obtain for the frequency difference

$$\Delta f = \frac{kf_s}{N},\quad k\in\mathbb{Z},\quad k\neq 0\tag{8}$$

Note that choosing the frequency difference according to $(8)$, the two sinusoids will be approximately orthogonal.
If we don't neglect the second term on the right-hand side of $(3)$ we can make them exactly orthogonal. However, in that case we cannot anymore freely choose one of the two frequencies, but the frequencies need to satisfy

$$f_1=\frac{(k_2-k_1)f_s}{2N},\quad f_2=\frac{(k_1+k_2)f_s}{2N},\\\quad k_1,k_2\in\mathbb{Z},\quad k_2>k_1\tag{9}$$


Downsampling and filtering (convolution)

I think you want to take the "central" part of the convolution. Convolution will give a signal that has length(c) = length(a) + length(b) - 1. If you have Matlab or Octave, there is an option when using the convolution function to just take the central part.


Why unit impulse function is used to find impulse response of an LTI system?

It is important to realize that due to the LTI property the response to an impulse $\delta(n)$ describes the system completely, because any input signal can be written as a weighted sum of shifted impulses:

$$x(n)=\sum_kx(k)\delta(n-k)\tag{1}$$

Let $\mathcal{T}\{\cdot\}$ denote the system operator, so the response to a shifted impulse is

$$\mathcal{T}\{\delta(n-k)\}=h(n-k)\tag{2}$$

From (1) and (2), due the linearity and the time- (shift-)invariance of the system, the response to any input signal $x(n)$ can be written as a weighted sum of responses to shifted impulses:

$$y(n)=\mathcal{T}\{x(n)\}=\sum_kx(k)\mathcal{T}\{\delta(n-k)\}=
\sum_kx(k)h(n-k)$$

This last equality is the convolution sum describing the input-output relation of a discrete LTI system.


You're not actually mixing two complex exponentials, because otherwise you'd get a complex exponential with the sum of the frequencies of the individual complex exponentials. That was already mentioned in a comment by Marcus Müller.
What you are doing is compute a signal with real part $\cos(\omega_0t)\cos(\omega_1t)$ and with imaginary part $\sin(\omega_0t)\sin(\omega_1t)$. If you lowpass filter that signal, only the components with the difference of the two frequencies remain, but these components are the same for the real and for the imaginary part, because
$$\cos(\omega_0t)\cos(\omega_1t)=\frac12\big[\cos[(\omega_0-\omega_1)t]+\cos[(\omega_0+\omega_1)t]\big]\tag{1}$$
and
$$\sin(\omega_0t)\sin(\omega_1t)=\frac12\big[\cos[(\omega_0-\omega_1)t]-\cos[(\omega_0+\omega_1)t]\big]\tag{2}$$
Clearly, the components with the difference frequency are the same.
Also, if you want to filter out a component with frequency $\omega_x$ (in your case the component with the sum frequency), it is not wise to use a filter with cut-off frequency $\omega_c=\omega_x$ because the attenuation at $\omega_c$ is usually quite modest, in your case $3$dB. Choose a frequency that is much closer to the lower frequency component to get a decent attenuation of the higher frequency component.


This has little to do with intuition. Tensors are rigorously defined mathematical objects, and in general simple arrays don't qualify as tensors. Specifically, signals are not tensors.

In the language of mathematics, and there the field of differential geometry, a tensor is a linear function with several arguments. It is therefore also called a multi-linear form. It is used to describe properties of manifolds at a single point.

A point on a (smooth) manifold comes with a tangent space that contains all possible "direction" at that point, where a direction is really just a tangent vector at that point. This tangent space at a single point is a vector space and it has the same dimension as the manifold.

The tangent space has a natural so called dual space, which contains the linear functions that map tangent space vectors to real numbers. This is the cotangent space, and it's also a vector space of the dimension of the manifold.

A tensor is now (roughly) a linear function that maps any number of dimensions between the tangent space, the cotangent space and the real numbers.

Tensors typically don't come as singles but in so called bundles, which are the disjoint union of the tensors of all points of the manifold. For example if you have a vector at every point of the manifold then the whole structure is called a vector bundle. Similarly, tensor bundles associate a tensor with every point on the manifold, usually in a smooth way.

In signal processing you can encounter tensors, but usually they are called differently. For example in video processing the image is your manifold (or rather, a function on a 2-dimensional manifold giving the brightness or color for each point) and the velocity field that describes the local motion is a rank 1 tensor field (or bundle) on that manifold.

Tensors also show up a lot in volume data processing, where tensors can describe local properties of the volume data. For example if you have a doppler ultrasonic tomography image, then the velocity data in each voxel is a tensor field as would be mechanical stress in tomographic material analysis.

For a 1-dimensional signals tensors often come in the form of derivative operators. For example if you have a signal $s(t)$, then the operator $\frac{\partial}{\partial t}$ is a tensor field of rank 1, or just a vector field. More generally, if you multiply that operator with a smooth function $g(t)$, then all possible vector fields on the manifold the signal lives on are of the form $g(t)\frac{\partial}{\partial t}$

To sum up, tensors are a description of differential properties of manifolds. If you want to understand them, you'll have to study differential geometry.


Why do musical instruments have different colour in phone?

Typical telephone codecs use a sampling frequency of 8kHz and a bandwidth of 7kHz.  That means that your signal will get cut off at frequencies higher than 3.5kHz.  That's ok for most speech but not sufficient for reproducing the timbre of most music instruments.

You'll likely fare reasonably well with instruments like recorders (which don't have a lot of overtones).  The worst combination is high-pitched and overtone-rich, like violin, oboe, soprano sax.  Also instruments with a full spectrum will sound different, like high hats, snares, gongs.


How sine wave is generated using lookup tables

Lookup tables are incredibly useful.

In audio, when you're doing a chorus or flange effect, or a synthesizer with a modulating envelope, you want to be able to cycle through a sine wave in a lookup table and not just keep track of some infinitely growing index n.

You'll want to use some mathematical program like MATLAB, Octave, or Mathematica to generate as many points as is practical (I wouldn't go with less than 512 or 1024 points) and store that in an array. This is your lookup table.

For a fractional index, what you'll end up doing is linearly interpolating around the nearest two indices at that fractional value. Let's say your index is 8/3. You've called up your index in a 1024 point table four times. What do you do? Now you're at index 32/3.

You'll need to find the lower index floor(32/3) and the higher one ceil(32/3). You'll need some conditionals in case you're at the end. I'll leave that exercise to you.

So now you've got your neighbor indices, you'll evaluate the interpolated value at that fractional index using the linear interpolation formula.

y = x*y2 + (1-x)*y1;

where x is the difference between your fractional index and the nearest index rounded down, y1 is the nearest value of the lookup table rounded down (you found the index earlier, now just use it to index into the lookup array) and y2 is the same thing for the nearest value rounded up.

You would track the index in a variable and increment it each loop of your program by 8/3 and do a wraparound when necessary (when you go beyond the length of the lookup table).

The trickiest thing is keeping track of the indices. The nearest value rounded up may be the 0th index.


Edge map based on a Haar Wavelet Transform

One single level of a standard separable 2-channel wavelet transform, denoted by $i$, uses a low-pass $l$ and a high-pass $h$  filters (followed by downsampling). Traditionally, one applies $l$ and $g$ on the rows of the image, putting the downsampled low-passed coefficients on a left-half, and the downsampled high-passed coefficients on a left-half. Then the resulting image is processed in the same way column-wise.

This results in a set of four separable 2D filters:

$l$ on rows, $l$ on columns (LL)
$l$ on rows, $h$ on columns (LH)
$h$ on rows, $l$ on columns (HL)
$h$ on rows, $h$ on columns (HH)

producing the arrangement below. This is iterated (levels) for on the image, then on the low-pass/low-pass images (LL)

Among those filters, three (LH, HL, HH) can be considered as edge detectors ($h$ is similar to a 2-point gradient) in different directions: more horizontal (HL), more vertical (LH), somehow diagonal (HH).

The corresponding coefficients, as least their combined energy, can be interpreted as a measure of edge strength, a map computed for each level.

Remember that for a continuous 2D field modeling an image $f(x,y)$, the image gradient

$$\nabla f =\left[\frac{\partial f(x,y)}{\partial x} ,\frac{\partial f(x,y)}{\partial y}  \right]^T$$

represents the local directional change in image intensities, whose norm is a measure of gradient "magnitude". You can compute several norms, the most common being the Euclidean, the Taxicab or the max norms.

The Haar filter $[1\,-1]$ is one possible discretization for the 1D gradient, the three 2D Haar wavelets compute an horizontal, a vertical and a somewhat diagonal gradient, hence their combined energy (Euclidean norm) is an estimate of the true gradient magnitude.


For an LTI system, why does the Fourier transform of the impulse response give the frequency response?

Let $h(t)$ denote the impulse response of an LTI system. Then, for any input $x(t)$,
the output is $$y(t) = \int_{-\infty}^\infty h(\tau)x(t-\tau)\,\mathrm d\tau.$$
In particular, the response to the input $x(t) = \exp(j2\pi ft)$ is
$$\begin{align}
y(t) &= \int_{-\infty}^\infty h(\tau)\exp(j2\pi f(t-\tau))\,\mathrm d\tau\\
&= \exp(j2\pi ft)\int_{-\infty}^\infty h(\tau)\exp(-j2\pi f\tau)\,\mathrm d\tau\\
&= H(f)\exp(j2\pi ft),\tag{1}
\end{align}$$
where $H(f)$ is the Fourier transform of the impulse response
$h(t)$.  In words, for an LTI system with impulse response $h(t)$, the
input $\exp(j2\pi ft)$ produces output $H(f)\exp(j2\pi ft)$.
This is precisely what we have as the definition of the
frequency response of an LTI system (call
this $FR(f)$ for now):

  for every frequency $f$, the response
  to $\exp(j2\pi ft)$ is $FR(f)\exp(j2\pi ft)$ which is just
  a (complex) constant times the input complex exponential.

But $(1)$ shows that $FR(f)$ is just $H(f)$, the Fourier
transform of the impulse response $h(t)$, which is what you
wanted to prove.


Steady state value of a complex convolution

You can't apply the final value theorem here. First of all, the correct form of the final value theorem for the $\mathcal{Z}$-transform is

$$\lim_{n\to\infty}y[n]=\lim_{z\to 1}(z-1)Y(z)\tag{1}$$

Eq. $(1)$ assumes that the poles of $(z-1)Y(z)$ are inside the unit circle, which is not the case in your example. In your case, the limit $(1)$ does not exist.

What you have to do is find the frequency response of the system, which is just the transfer function evaluated on the unit circle $z=e^{j\omega}$ (if it converges on the unit circle). Note that the region of convergence (ROC) you've found does not make sense because $|z|$ must be a real number. The requirement for the sum to converge is

$$\left|\frac{j}{2z}\right|<1\quad\Longrightarrow\quad |z|>\frac12\tag{2}$$

Since $H(z)$ converges on the unit circle $|z|=1$, the frequency response is given by

$$H(e^{j\omega})=\frac{1}{1-\frac{j}{2}e^{-j\omega}}\tag{3}$$

An LTI system's response to a sinusoidal input $x[n]=\cos(\omega_0n)$ is

$$y[n]=|H(e^{j\omega_0})|\,\cos[\omega_0n+\phi(\omega_0)]\tag{4}$$

where $\phi(\omega)$ is the system's phase response:

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{5}$$

So in the limit, after all transients have died out, the system's (steady-state) response to a switched sinusoidal input $x[n]=\cos(\omega_0n)u[n]$ is given by $(4)$. What is left now is to determine the magnitude and phase of $H(e^{j\omega})$ given by $(3)$ and plug these values into Eq. $(4)$.


How exactly does 8PSK work on each I and Q leg?

You are mapping a group of 3 bits (an 8PSK symbol) to one of 8 coordinates in the I-Q plane.  I'm not an FPGA guy, but you can do this with a "switch" statement or a lookup table.  I don't know your particular mapping, but a straightforward mapping from bits to I & Q values would be:

$$(0,0,0) \rightarrow (1, 0)$$
$$(0,0,1) \rightarrow (1/\sqrt{2},1/\sqrt{2})$$
$$(0,1,0) \rightarrow (0, 1)$$
$$(0,1,1) \rightarrow (-1/\sqrt{2},1/\sqrt{2})$$
$$(1,0,0) \rightarrow (-1, 0)$$
$$(1,0,1) \rightarrow (-1/\sqrt{2},-1/\sqrt{2})$$
$$(1,1,0) \rightarrow (0, -1)$$
$$(1,1,1) \rightarrow (1/\sqrt{2},-1/\sqrt{2})$$

Note that, in the above mapping, if you take entries 2, 4, 6, and 8, and ignore the least significant bit, you get a mapping for QPSK.


How is digital signal processing related to Computer Vision?

My attempt: You can treat an image as a discretization of the perspective projection of a continous 3D region, given a near and far plane. Using this discretization, filters used in signal processing can be used to modify the signal (a 2D array) for various purposes. These purposes include things like edge detection, image denoising, object detection, instance/semantic/general segmentation, and various other tasks. Notably 2-D wavelet transforms borne out of signal processing are used to transform image signals to another basis using function spaces spanned by a finite basis of orthogonal functions and perform filtering in this space, after which you invert the transform to achieve a modified image.


If sensor A has a defect, the clear answer is to only use sensor B.
A preferred solution to minimize noise would be to do a weighted average based on the quality of each sensor, when that can be actively characterized. This can be easily done for the case the OP has presented of taking a reading from a single rotating axis.
The optimum combining for the general case of independent signals in the presence of white noise is to do a weighted average of each signal using the signal to noise ratio, which is referred to as maximal ratio combining. In practical application we determine the extent that we can use these assumptions, and our ability to extract an SNR metric depends also on how stationary the signals are. With consideration to this we arrive at a reasonable approach for the OP specifically.
A standard average of the result of two sensors given as x and y would be done as:
$$r_{avg} = \frac{x + y}{2} \tag{1} \label{1}$$
Alternatively, with maximal ratio combining, the two sensors could be combined as follows:
$$r_{opt} = \frac{\sigma_y^2 x + \sigma_x^2 y}{\sigma^2_x  + \sigma^2_y  } \tag{2} \label{2}$$
Where:
$\sigma^2_x$, $\sigma^2_y$: variance of each sensor's measurement
It is intuitively clear From \ref{2} how this would be superior, since we are now doing a weighted average with favoritism toward the less noisy sensor in the average of the two. [Addressing FAT32's good comment to an earlier version of this post; in addition to the variance estimate $\sigma^2$, a moving average estimating the mean $\mu$ of each sensor can be tracked for further robustness against erroneous measurements or a completely faulty sensor.]
Notice that when we have no further information as to the quality of each sensor, \ref{2} becomes \ref{1} when we otherwise then set $\sigma^2_x = \sigma^2_y$. To the extent that the measurement results are stationary, we are able to easily extract the first and second order statistics (mean and standard deviation as $\mu$ and $\sigma^2$) directly from the measurements of each sensor to then be able to apply maximal ratio combining. In many cases stationarity can be reasonably assumed over a finite time interval, allowing in the OP's case to easily estimate the mean and standard deviation from each sensor over the maximum time duration allowed by the dynamics of the experiment, bandwidth of sensors etc, and then combining as in \ref{2}.
For a very robust (but detailed) approach, see How to interpret Allan Deviation plot for gyroscope? for further details on using the two-sample variance to dynamically determine the best averaging time for the mean and variance estimates. Thus it can be envisioned where the measurement system is continuously monitoring a moving windowed mean, variance and two-sample variance of the sensors to provide the maximal ratio solution, where the window duration is updated based on minimizing the two sample variance. That level of complexity is likely unnecessary, as it may be sufficient to choose a moving window that is known to be a stationary interval and simply compute the mean and variance over that duration to get the result defined in \ref{2}; but shows an interesting robust approach of what else could be done.
In summary, Maximal Ratio Combining is superior solution over simple averaging. A simpler, suboptimum but effective strategy is to select the sensor with the highest SNR and ignore the other, or doing the simple average and disabling a sensor upon a threshold noise condition (which should be included regardless, such as a stuck sensor or excessive offset condition, to Fat32's points in the comments). Worst case would be to arbitrarily just average both for the very reason the OP brings up when it is entirely feasible that one of the sensors could become defective. A simple average would only make sense when there is high confidence both sensors are without defect and we believe to have similar performance.
This is similar to the strategies used in diversity receivers where two independent channels are optimally combined:
https://en.wikipedia.org/wiki/Diversity_combining
Also to include from the comments the link Fat32 found which refers to this approach as the "Fraser-Potter fixed-interval smoother", from the sensor community while I was exposed to equivalent "Maximal Ratio Combining" from the wireless communications community.
https://en.wikipedia.org/wiki/Sensor_fusion
As a demonstration per my friend Fat's request, below shows the results from a simulation of two sensor readings with a background trend, offset, and slope differential, where the statistics of the signals are not known apriori (everything is extracted from the signals themselves).
The top plot represents the raw measurements from the two sensors, with one sensor significantly noisy. The bottom plot shows the comparative results of simply averaging the reading from each sensor to produce an average result after each update (standard average) and the maximal ratio combining average where each of the two samples being averaged is first weighted by a variance estimate derived from the data from each sensor. The estimators used a moving window to measure the variance of each sensor over the past 500 samples.

An error signal was generated comparing the final results to the avg of the underlying sensor signals without the added noise and is plotted below showing the clear performance improvement in overall noise of the final result.

I'll also include the results of the Allan Deviation (root of the 2-sample variance) showing the significance behind it. I chose to use 500 samples from a review of the ADEV confirming for us in this particular case that statistical computations up to 500 samples will improve the overall estimates, anything beyond that or less than that will be degraded (this includes our estimates of the variance of each sensor for the maximal ratio combining to get each result as well as how long we can average those results over time for the final sensor result). This isn't a necessary step but something I do personally check in any statistical analysis and here demonstrates its clear utility. Basically this is saying that even though there is a background slope, given the noise level itself we can assume stationarity in this case over a 500 sample time interval. This means we would improve our noise estimates by using up to the past 500 samples for each measurement to determine the statistics for the weighting. (And using as much as 2000 samples would provide as good of an answer as 100 samples!). Also most significantly, note the clear difference between the standard average and the maximal ratio combining, in that the latter is better than the best sensor alone!
The horizontal axis for the plot below is in log samples, showing the minimum of the more dominant noise at a 500 sample averaging time.
The 30 second primer in how to read the ADEV plot: The horizontal axis indicates an averaging duration (here in samples) while the vertical axis is an indication of the standard deviation of the mean of a random variable after being averaged over that duration. This relates to our confidence in any estimator that may be used on the underlying data. For example, when determining the average value over $N$ samples, if the samples are random white noise and stationary, the standard deviation will go down by the rate of $1/\sqrt{N}$ which is exactly what we see in the plot for all curves out to 100 samples (and out to 500 samples for the blue curve of the noisiest signal). After that it bottoms out and turns up due to the underlying drift in the measurement. This tells us for all practical purposes we can assume the noise on the data is white and stationary out to the duration in time where the ADEV plot is still going down by the rate of $1/\sqrt{N}$, and that we would minimize the uncertainty (noise) of any estimators that include a mean as part of their computation (rms, variance, etc) by using data over that duration: very simply it tells us the true extent as to how long in time we can average for and keep improving a result; average over a longer duration and results actually get worst!   (In this case we could simply remove the linear drift to them take advantage of an even longer averaging time - but the minimum would typically occur due to the signal statistics being no longer stationary due to more random effects that cannot be predictably removed.)

TO DO: The Kalman Filter similarly determines optimum weights for the averaging between the two samples. Show how in this case the Kalman would result in the same weighting solution (it is the optimum combining for this case) and why, and show under what conditions specifically the use of the Kalman would be warranted...


If you are asking about the grid resolution of the PSD then it is given by:

$$ \frac{{F}_{s}}{N} $$

Where $ {F}_{s} $ is the sampling rate and $ N $ is the number of samples.
Pay attention that you can increase the grid resolution by Padding the data with zeros.
Though it will increase the grid resolution it won't increase the resolution of the DFT (In the sense of resolution to see to close, in frequency, harmonic signals).


How to Match 2 Signals which Are a Shifted and Scaled Version of Each Other

Matlab has a function called xcorr() that computes the cross correlation between two vectors.  The cross correlation will be highest when the two signals overlap.  The following code will compute the shift in units of samples:

x = xcorr(s1,s2)
[val,idx] = max(x);
shift = length(s1) - idx

The scaling factor should be given by:

scale = max(s1)/max(s2)


Let there be two microphones, $x_1(t), x_2(t)$ capturing a source, $s(t)$, with some additive noise $n_{1,2}(t)$ that is uncorrelated with the source and with each other (i.e., its effect can be ignored while calculating the cross-correlation). Let $\tau_0$ (in seconds) be the time difference of arrival between the two microphones. The time domain equations of the microphones are
$$x_1(t) = s(t) + n_1(t) \\
        x_2(t) = \alpha_0 s(t - \tau_0) + n_2(t) $$
Now, taking the FFT and converting these to the frequency domain would give us,
$$
        X_1(\omega) = S(\omega) + N_1(\omega) \\
        X_2(\omega) = \alpha_0 S(\omega) \exp{(-j\omega \tau_0)} + N_2(t)
$$
The cross correlation function of the two microphone signals (or the cross-power spectrum in the frequency domain), can be written as
$$      \begin{aligned} R_{x_1,x_2}(l) &= \mathbb{E}[x_1(t) x_2(t-l)] \\
        R_{x_1, x_2}(\omega) &= X_1(\omega)X_2^*(\omega) \\
                             &= |S(\omega)|^2 \alpha_0 \exp{(j\omega \tau_0)}
        \end{aligned}
$$
The Generalized Cross-Correlation function with Phase Transform (GCC-PHAT) at lag $l$ is defined as
$$ \begin{aligned}
\tilde{R}_{x_1,x_2}(l) =& \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{R_{x_1,x_2}(\omega)}{|R_{x_1,x_2}(\omega)|}e^{j\omega l}d\omega \\
=&  \frac{1}{2\pi} \int_{-\pi}^{\pi} e^{j\omega(l + \tau_0)} d\omega \\
&= \delta(l - \tau_0)
\end{aligned} $$
where $\delta$ is the Dirac-delta function which is non-zero at $l = \tau_0$, which is the time-difference of arrival (TDOA). Typically, the lag is calculated in samples, so to convert that to seconds, you have to divide it by the sampling rate.
I hope this helps in understanding GCC-PHAT, which is useful tool for TDOA estimation, and hence, source localization.


Finding pitch from a wideband spectrogram

If your window length is shorter than the pitch period of a voiced utterance, the spectrogram will not be able to capture the fundamental frequency. This is the problem, yes?

But evidence of the fundamental frequency will still be available in the spectrogram because spoken voice has plenty of harmonic content. If you can identify the upper harmonics, showing up as horizontal stripes in the spectrogram, and find the size of the gap between them, this will be the fundamental pitch. Say your FFT resolution is very coarse under 200 hz, but you find a peak at 240 hz, 360 hz, and 480 hz. Well, this is a dead giveaway that the fundamental was 120 hz.

Moreover: The glottal impulse (vocal folds slamming together) is like a broadband impulse. You will get periodic vertical stripes along the spectrogram if the pitch is low enough, and your FFT window short enough.  The time difference between these stripes is the pitch period of the voice, so this is another clue you can use to determine the pitch even when the fundamental is ambiguous.


An Intuitive Way to Understand the Sharpen Convolution Filter in Image Processing

Consider this filter in 1D (for the sake of simplicity):

[-1 1 0]

This filter responds to pixel differences, i.e. edges (left pixel * (-1) + center pixel * 1 + right pixel * 0) which is (center pixel - left pixel).

Response of this filter is positive on rising edges and negative on falling edges.

The filter is not symmetric and hence it will "shift" the image to the left. To make this filter responding to "edge on the left" and "edge on the right" equally, we have to compensate it by similar filter:

[-1 1 0] + [0 1 -1] = [-1 2 -1]

Now this can produce "edge image". We will enhance the edges by adding them to the original image:

[-1 2 -1] + [0 1 0] = [-1 3 -1]

This kernel have sum of 1 and does not produce negative values. It is the desired sharpening kernel.

Summary: Sharpening =  edge detection + original image.


Unilateral Laplace Transform's Differentiation Property

You need to look at the derivation of that property. Integration by parts gives
$$\begin{align}\mathcal{L}\{f'(t)\}&=\int_{0^-}^{\infty}f'(t)e^{-st}dt\\&=f(t)e^{-st}\Big|_{0^-}^{\infty}+s\underbrace{\int_{0^-}^{\infty}f(t)e^{-st}dt}_{F(s)}\\&=\lim_{t\to\infty}f(t)e^{-st}-f(0^-)+sF(s)\tag{1}\end{align}$$
The first term in $(1)$ is only guaranteed to vanish for $t\to\infty$ if $\text{Re}\{s\}>\alpha$ for some value $\alpha$. So even if you define a unilateral Fourier transform, that term may not vanish in general.
However, if $\lim_{t\to\infty}f(t)=0$ holds, we could use the same property with a unilateral Fourier transform. It's just much more common to use the well-established (unilateral) Laplace transform in cases where non-zero initial conditions need to be taken into account.


FIR filter design for response inside unit circle

Note that an FIR filter's transfer function is completely determined by its behavior on the unit circle (or any other closed contour). So if you were to define a desired $\mathcal{Z}$-transform (instead of a desired frequency response on the unit circle), you could just evaluate this desired transfer function on the unit circle and use one of the many standard design procedures.

If it were convenient for you to specify a desired response on a circle in the $z$-plane with radius $r$, you could also just use a standard frequency domain design procedure, and obtain your final coefficients by computing

$$h[n]=r^{n}\tilde{h}[n]$$

where $\tilde{h}[n]$ are the coefficients obtained from the standard design procedure.

Of course you can also directly specify the desired zeros in the complex plane, as mentioned in Jason's comment. Note that you can't specify the poles because for a causal FIR filter they must lie in the origin of the $z$-plane.


How to prepare different input size in CNN

You can research "Spatial Pyramid Pooling (SPP)"
Btw, as far as I know, you don't need make same size for input for convolutional layer but FCN layer


find impulse response from step response

You can find the impulse response

Let's take the case of  a discrete system.

If $s[n]$ is the unit step response of the system, we can write

$$s[n]= u[n]\ast h[n]$$

where $h[n]$ is the impulse response of the system and $u[n]$ is the unit step function.

Now using commutative property you can write $$s[n]=h[n]\ast u[n]$$

Expanding convolution we get $$s[n] = \sum_{k=-\infty}^{\infty}h[k]u[n-k]$$

which can also be written as

$$s[n]=\sum_{k =-\infty}^{n} h[k]$$

So we can see that unit step response is like an accumulator of all value of impulse response from $-\infty$ to $n$. So now impulse response can be written as the first difference of step response.

$$h[n]=s[n] - s[n-1]$$


Maybe you are looking for something called "audio fingerprinting". This might employ DSP techniques, but in this form the question might be percieved as somewhat broad. There is no one correct algorithm of this kind, rather a bunch of good ideas that can be combined in different ways.


Capacity of cascade binary symmetric channels

A binary symmetric channel (BSC) can be characterized by its complemented probability $p$. Its well-known capacity is

$$C = 1 - H(p) = 1 - (-p\log(p) - (1-p)\log(1-p))$$

where $H(p)$ is binary entropy function:

A $L-$concatenated BSC, which is also a BSC characterized by $p_L$, can be visualized as in the figure below

The complemented probability $p_L$ is derived

\begin{align}
p_L &= p_{L-1} (1-p) + (1-p_{L-1}) p \\
&= p + (1-2p)p_{L-1}\\
\implies 1 - 2 p_L &= (1-2p)(1-2p_{L-1}) \\
\implies 1 - 2 p_L &= (1-2p)^L \\
\end{align}

Thus
$$p_L = \frac{1}{2}(1 - (1-2p)^L)$$

If $p=0$ then $p_L = 0$, $H(p_L) = 0$ and $C = 1$.

If $p=1$ then $p_L = 0$ or $p_L = 1$ depending on $L$ is pair or impair; $H(p_L) = 0$ and $C = 1$.

If $0 < p < 1$ then $\lim_{L \to \infty} p_L = 0.5$; $H(p_L) = 1$ and $C = 0$.

Conclusion: if the unit BSC is not certain $(p \neq 0, 1)$, the capacity of infinitely-concatenated BSC tends to $0$.


The requirement for being able to detect the "profile" of hand gestures is not a straightforward one to satisfy. There are issues of form and rotation / scale independency which are not "easy".

There are two main problems here:

Segmenting the hands
Figuring out the gesture.

To segment the hands, let's assume that the hands are imaged over a constant single colour background. This means that the overall hand region can be isolated with a "simple" threshold.

The second problem is one of recognising form. In other words, what are the "shape" differences between "rock", "paper" and "scissors"?

Therefore, whatever is used should be "able" to discriminate between "holding two fingers extended" versus "all fingers retracted in the shape of a fist", etc.

The "easy" way out is to use template matching. Very "simply", take pictures of the hand gestures for each hand and for various orientations and use them to "search" for their existence in an image using cross correlation (or other similar concept in some feature space).

This however, done in the "simple" way is very costly (computationaly) and also scale dependent. So, if the hands are near the camera they are larger, if they are further away they are smaller and therefore very difficult to "catch" from static images that were acquired at some standard distance.

So, an alternative would be to use Morphological Operators and specifically Skeletonisation.

The main idea here is to apply skeletonisation on the extracted hand regions, quite a few times to end up with a set of "lines" describing the shape you are trying to recognise. In this way, in one pass you effect the skeletonisation and in another pass you collect the (very few, compared to the original) pixel locations that compose the gesture. To then figure out the gesture itself, you work with features extracted from the pixel locations.

The added benefit of this type of a solution is scale invariance. The relative positions of the "finger segments" of the "scissors" gesture are the same whether you hold it near the camera or far away from the camera (Two lines forming a V ending up in a big blob). With careful selection of the features and search process you can even use the same algorithm and training dataset for the left and right hands separately (for example, using graph search).

Hope this helps.


To calculate the energy of a continuous signal you use the equation:  $$E_\infty = \int\limits_{-\infty}^{\infty} |x(t)|^2 dt$$   and for the power:   $$P_\infty = \lim\limits_{T\to\infty} \frac{1}{2T}\int\limits_{-T}^{T} |x(t)|^2 dt$$  If your signal has finite power ($0<P<\infty$) then it's a power signal. If the signal has finite energy then it's an energy signal (Notice that a signal can be neither type but can't be an energy and a power signal at the same time).


Consider some actual use cases for further insight:

Unmodulated Carrier: The peak to average ratio is 3 dB given the relationship between the rms and peak level of a sine wave.

Unfiltered Single-Carrier BPSK or QPSK: This would be the modulation of rectangular shaped symbols, so would have a 3 dB peak to average ratio just as in the case of an unmodulated carrier. Great for peak-average ratio but terrible for bandwidth constraint. Given we have to share bandwidth efficiently (in terrestrial communications) we don't have this luxury. GPS is an example where the signal can be essentially unfiltered and we have the benefit of power efficiency.

Root-Raised Cosine BPSK or QPSK: This is what is traditionally done  but the tighter the filtering (smaller $\alpha$ in the filter design) results in more overshoot of the symbol transitions and slower trajectories that may go through the origin resulting in a higher peak-average ratios.

pi/4 QPSK: This modulation changes phase by $\pi/4$ between symbols thus avoiding having the slower transitions from filtering from ever going through the origin (Meaning minimizing the lowest instantaneous power in the signal envelope, increasing the average relative to the peak). This modulation is therefore considered for lower power (battery operated) applications since we are able to drive the transmit power amplifier further into saturation meaning higher efficiency.

GMSK: Gaussian Minimum Shift Keying and other constant-envelope modulations have a 0 dB peak-average ratio and offer moderate bandwidth efficiency so are an excellent choice for power efficient (battery operated) applications.

OFDM: To the extent each bin in the FFT representation of an OFDM symbol is independent and in consideration of each other randomly generated, the resulting signal as the summation of multiple independent identically distributed random processes will be well approximated by a Gaussian distribution given the central limit theorem. A Gaussian distribution has a relatively high peak to average ratio, and we need to decide on a significance to establish what is the "peak" (usually on a distortion criteria from clipping anything beyond the established peak). In SC-OFDM we do have some control over the association of the signals in each bin, allowing for lower peak-average ratios than MC-OFDM where each bin can represent a different data source.


Coarse summary of various methods for detecting when a signal 'starts'

If both signals are audio signals from the same source (just different paths to get there) you can do

Calculate the envelope (various methods for that, depends on the signals and the accuracy you want)
Down sample
Do a running cross correlation (at the down sampled rate)


How many multiplications and additions does it take to compute an FFT of a signal?

If you are interested in more involved numbers of operations, you can check A modified split-radix FFT with fewer arithmetic operations, by Johnson and Frigo, IEEE Transactions on Signal Processing, 2007, authors of the "Fastest Fourier Transform in the West" (FFTW). For instance, they claim that the complexity of the classic radix-2 is about:
$$ 5N\log_2 N$$
while the split-radix yields:
$$4N \log_2 N − 6N + 8\,,$$
and they achieved the improvement:
$$\frac{34}{9} N \log_2 N − \frac{124 }{27 }N − 2 \log_2 N − \frac{2}{9}
(−1)^{\log_2 N} \log_2 N + \frac{16}{27}
 (−1)^{\log_2 N} + 8\,.$$
As said by others, you can see this number depends on implementations, and that the actual gain is subtle.


What is the meaning of positive and negative shift in signal correlation?

Your Eq(3), for positive shift, is correct. And your Eq(5), for negative shift, should also be the same, but is wrong for $\tau < 0$ case.

The lower limit of the integral in Eq(4) should be $-\tau$ but you took it $\tau$. The reason is the step function $u(t+\tau)$ will have the nonzero range for $t > -\tau$.

You will get the same result with the correct limit.


Transient response of system with single pole $0 \le p < 1$

If you have a single real-valued pole, you get a term $k\,p^{(n-n_0)}u[n-n_0],n\ge 0$, with some constant $k$ and some delay $n_0\ge 0$ in the system's impulse response ($u[n]$ is the unit step). Clearly, if $|p|\ge 1$ the transient will never settle. If $-1<p<0$, the transient will alternate in sign (i.e., "oscillate"), like e.g.,

$$\left(-\frac12\right)^nu[n]=\left\{1,-\frac12,\frac14,-\frac18,\dots\right\}\tag{1}$$

So for a transient that decays (i.e., for the system to be stable) and that does not oscillate you need $0\le p<1$.

Note that if you only have one real-valued pole $p$, stability implies $$-1<p<1\tag{2}$$ If you exclude systems with oscillating impulse responses (i.e., with a negative pole), you're left with the requirement $0\le p<1$. But, as you've noticed, this requirement is not only dictated by stability (which only requires $(2)$), but also by the requirement of having no oscillations.


My suggestion would be in following step.

If we know that outer shapes (rectangles) are at regular interval - first we should do some sort of segmentation and/or localization.
Remove portions of background which are of no use and each block is identified and processed separately.
Now if you are expecting lines only horizontally, you can traverse vertically and create scan lines from top-to-bottom. If the region has no dark separations - the scan lines will appear straight albeit with noise. If it crosses the dark line the intensity in the scan line will drop. This will tell the position of dark strip in the block.
A set of 4 - 5 scan lines (or more) will provide a good co-relating answer.


FFT error when number of samples growth - Matlab/octave

An FFT using symbolic math might be possible, but would be many orders of magnitude slower.  (I'm guessing at least 10,000 times slower, except for a set of exact equation input signals).  You would have to use a symbolic math package instead of Matlab.  (Perhaps Mathematica, Maxima, or Maple?)

Instead you might be able pre- and post-process certain inputs to produce known outputs.  For instance, you can check if the input is strictly symmetrical around input sample 0.  If so, set all the imaginary components of the FFT output to zero.  If strictly odd (anti-symmetric), set all the real result components to zero.  If strictly real, set all the negative frequency outputs to be exact complex conjugates of the positive frequency ones. etc.


Proving Nyquist Sampling Theorem for Strictly Band Limited Signals (Whittaker Shannon Interpolation Formula)

Approaching The Sampling Theorem as Inner Product Space
Preface
There are many ways to derive the Nyquist Shannon Sampling Theorem with the constraint on the sampling frequency being 2 times the Nyquist Frequency.
The classic derivation uses the summation of sampled series with Poisson Summation Formula.
Let's introduce different approach which is more similar to function analysis - Building an orthogonal space and using projection for analysis and synthesis.
Forming Orthonormal Basis
In this section we'll define an orthonormal base and derive the decomposition and composition process.
Definitions
First one should define the space of Band Limited Functions. The space of Band Limited Functions is defined by:
$$ \mathcal{B}_{ {W}_{s} } = \left\{ f \left( x \right) \mid F \left( \omega \right) = \mathscr{F} \left\{ f \left( x \right) \right\}, \; F \left( \omega \right) = 0, \; \forall \, \left| \omega \right| > {W}_{s} \right\} $$
In words, it means that for each function $ f \left( x \right) \in \mathcal{B}_{ {W}_{s} }$ its fourier transform $ F \left( \omega \right) = \mathscr{F} \left\{ f \left( x \right) \right\} $ vanishes for frequencies $ \left| \omega \right| > {W}_{s} $.
The inner product in this space is given by:
$$ \langle f \left( x \right), g \left( x \right) \rangle = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) g \left( x \right) dx, \; T = \frac{2 \pi}{ {W}_{s} } $$
One could easily show that this is a indeed an inner product space with valid inner product.
The main claim is the orthonormal basis of this space is given by:
$$ f \left( x \right) = \operatorname{sinc} \left( \frac{ x - n T }{T} \right) $$
Where $ \operatorname{sinc} \left( x \right) $ is the Normalized Sinc Function given by $ \operatorname{sinc} \left( x \right) = \frac{ \sin(\pi x) }{ \pi x } $.
The basis functions are parameterized by the parameter $ n $. Basically we have shifted and scaled function as the basis.
Proof of The Orthonormal Property
One must show the orthonormal property of the basis under the defined inner product:
$$ \begin{aligned}
\langle f \left( x \right), g \left( x \right) \rangle & = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) g \left( x \right) dx = \frac{1}{T} \int_{- \infty}^{\infty} \operatorname{sinc} \left( \frac{ x - n T }{T} \right) \operatorname{sinc} \left( \frac{ x - m T }{T} \right) dx && \text{} \\
& \overset{1}{=} \frac{1}{T} \int_{- \infty}^{\infty} \left( \frac{1}{2 \pi} \int_{- \infty}^{\infty} T \, \Pi \left( \frac{ \omega }{ {W}_{s} } \right) {e}^{-j \omega n T} {e}^{j \omega t} d\omega \right) \operatorname{sinc} \left( \frac{ x - m T }{T} \right) dx && \text{} \\
& \overset{2}{=} \frac{1}{T} \int_{-\infty}^{\infty} \frac{1}{2 \pi} T \, \Pi \left( \frac{\omega}{ {W}_{s} } \right) {e}^{-j \omega n T} \left( \int_{- \infty}^{\infty} {e}^{j \omega x} \operatorname{sinc} \left( \frac{x - m T}{T} \right) dx \right) d\omega && \text{} \\
& \overset{3}{=} \frac{1}{T} \int_{-\infty}^{\infty} \frac{1}{2 \pi} T \, \Pi \left( \frac{\omega}{ {W}_{s} } \right) {e}^{-j \omega n T} T \, \Pi \left( \frac{-\omega}{ {W}_{s} } \right) {e}^{j \omega m T} d\omega && \text{} \\
& \overset{4}{=} \frac{T}{2 \pi} \int_{ - \frac{ {W}_{s} }{2} }^{ \frac{ {W}_{s} }{2} } {e}^{j \omega \left( m - n \right) T} d\omega && \text{} \\
& \overset{5}{=} \begin{cases}
1 & \text{ if } m = n \\
0 & \text{ if } m \neq n
\end{cases}
\end{aligned} $$
Where:

Since $ \mathscr{F} \left\{ \operatorname{sinc} \left( \frac{ x - n T }{T} \right) \right\} = T \, \Pi \left( \frac{ \omega }{ {W}_{s} } \right) {e}^{-j \omega n T} $.
Changing order of integration for converging integrals.
Applying $ \mathscr{F} \left\{ \operatorname{sinc} \left( \frac{ x - m T }{T} \right) \right\} \left( - \omega\right) $.
Integration boundaries according to the Rect function (Multiplication).
Integration over a cycle or over a constant.

Since we proved the suggested base is indeed an orthonormal basis of the space the following holds:
$$ \forall f \left( x \right) \in \mathcal{B}_{{W}_{s}} , \; f \left( x \right) = \sum_{n = -\infty}^{n} \langle f \left( x \right), {g}_{n} \left( x \right) \rangle {g}_{n} \left( x \right) $$
Where $ {g}_{n} \left( x \right) = \operatorname{sinc} \left( \frac{ x - n T }{T} \right) $ and $ \langle f \left( x \right), {g}_{n} \left( x \right) \rangle $ is the projection of $ f \left( x \right) $ on $ {g}_{n} \left( x \right) $.
Projection Process
As written above, using the result of a projection of a function in the space onto the basis one could reconstruct it as:
$$ f \left( x \right) = \sum_{n = -\infty}^{n} \langle f \left( x \right), {g}_{n} \left( x \right) \rangle {g}_{n} \left( x \right) $$
The question is, what's is the projection of a general function in this space? Well, it turns out it can be shown in a sloed form way:
$$ \begin{aligned}
\langle f \left( x \right), {g}_{n} \left( x \right) \rangle & = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) {g}_{n} \left( x \right) dx = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) \operatorname{sinc} \left( \frac{ x - n T }{T} \right) dx && \text{} \\
& \overset{1}{=} \frac{1}{T} \int_{- \infty}^{\infty} \left( \frac{1}{2 \pi} \int_{-\infty}^{\infty} F \left( \omega \right) {e}^{j \omega x} d\omega \right) \operatorname{sinc} \left( \frac{ x - n T }{T} \right) dx && \text{} \\
& \overset{2}{=} \frac{1}{2 \pi T} \int_{- \infty}^{\infty} \left( \int_{- \infty}^{\infty} \operatorname{sinc} \left( \frac{ x - n T }{T} \right) {e}^{j \omega x} dx \right) F \left( \omega \right) d\omega && \text{} \\
& \overset{3}{=} \frac{1}{2 \pi T} \int_{- \infty}^{\infty} T \, \Pi \left( \frac{-\omega}{ {W}_{s} } \right) {e}^{j \omega n T} {F} \left( \omega \right) d\omega && \text{} \\
& \overset{4}{=} \frac{1}{2 \pi} \int_{ - \frac{ {W}_{s} }{2} }^{ \frac{ {W}_{s} }{2} } F \left( \omega \right) {e}^{j \omega n T} d\omega && \text{} \\
& \overset{5}{=} f \left( n T \right)
\end{aligned} $$
Where:

Since $ \mathscr{F} \left\{ f \left( x \right) \right\} = F \left( \omega \right) $.
Changing order of integration for converging integrals.
Applying $ \mathscr{F} \left\{ \operatorname{sinc} \left( \frac{ x - m T }{T} \right) \right\} \left( -\omega \right) $.
Integration boundaries according to the Rect function.
Applying Inverse Fourier Transform at $ x = n T $.

Wrapping it yields:
$$ f \left( x \right) = \sum_{n = -\infty}^{n} \langle f \left( x \right), {g}_{n} \left( x \right) \rangle {g}_{n} \left( x \right) = \sum_{n = - \infty}^{\infty} f \left( n T \right) \operatorname{sinc} \left( \frac{ x - n T }{T} \right) $$
Which is known as the Whittaker Shannon Interpolation Formula.
Conclusion
In the process above the analysis and synthesis of Band Limited Functions is shown using Orthonormal Basis. If one set $ T $ to be the Sampling Interval, usually denoted by $ {T}_{s} $ then the Sampling Frequency is given by $ {F}_{s} = \frac{1}{ {T}_{s} } $.
Now, since the function in frequency domain stretches in the range $ \left[ -\frac{{F}_{s}}{2}, \frac{{F}_{s}}{2} \right] $ we indeed have the known relation that the sampling frequency has to be twice the one sided support of the function in frequency.


Stability of open-loop transfer function from its Nyquist plot

As mentioned by Tendero in a comment, you should be interested in the stability of the closed-loop transfer function

$$Q(s)=\frac{H(s)}{1+H(s)}=\frac{s^3(s+1)}{s^4+s^3-1}\tag{1}$$

where I've assumed unity feedback. Since the stability of $(1)$ is determined by the zeros of the denominator $1+H(s)$, we're interested in the encirclement of the point $-1+j0$ (and not in the encirclement of the origin). In your example, the number of clockwise encirclements of the point $-1+j0$ happens to equal the number of encirclements of the origin, i.e., we have $N=1$. Since the given $H(s)$ has no poles in the right half-plane, the number of zeros in the right half-plane equals the number of clockwise encirclements, i.e., $Z=N=1$, which means that $Q(s)$ has one pole in the right half-plane, and, consequently, the closed-loop transfer function is unstable.

This is also easily verified by computing the poles of $Q(s)$, which are the roots of the polynomial $s^4+s^3-1$:

roots([1,1,0,0,-1])
ans =

  -1.38028 + 0.00000i
  -0.21945 + 0.91447i
  -0.21945 - 0.91447i
   0.81917 + 0.00000i

This shows that one of the four roots is in the right half-plane, as predicted by the Nyquist plot.


Linear phase FIR filter for impulse responses that don't appear symmetric

Your confusion is understandable.

If you consider the definition of linear phase FIR filter and the associated symmetry conditions on their impulse responses, then you can arrive the conclusion that the first two cases

$$ h_1[n] = [0,0,0,1,0] $$

and

$$ h_2[n] = [0,0,0,0,1,0,0,0,0,0,1] $$

are non-symmetric. However, as you use zeros and ones in those impulse responses, it can be seen that the following two new augmented impulse responses are equivalent to those non-symmetric looking ones:

$$ h_3[n] = [0,0,0,1,0,0,0] $$

and

$$ h_4[n] = [0,0,0,0,1,0,0,0,0,0,1,0,0,0,0] $$

And they are symmetric and linear phase according to the definition as well. This happens because of the zeros and use of ones in a specific way.

However, the following impulse response

$$h_5[n] = [0,0,1,0,0,0,0,1,0,1,0]$$

is not symmetric, is not linear phase and cannot be made linear phase by augmenting with zeros...


What is the theoretical upperbound of information that could be transmitted by a device per second?

There are many factors involved in understanding the theoretical limits to communication. What follows is just a brief introduction that only scratches the surface.

First, let's consider a simple scenario: there is no noise and no distortion of the signal being transmitted. We do allow for attenuation. Under these circumstances, you can transmit up to $R_p=2B$ pulses per second. $B$ is the bandwidth you have available. For example: if you have 10 kHz of bandwidth, you may transmit up to 20,000 pulses in one second. $R_p$ is called the Nyquist rate, after Harry Nyquist.

However, your question is about bits, not pulses. It turns out that you can transmit bits using pulses. Let's say you want to transmit the bit sequence $1001$. You could do this by transmitting pulses with amplitudes $+1, -1, -1, +1$. (The pulse shape does not matter, but they do need to "fit" in the bandwidth you have available.) If you do this, your bit rate $R_b$ is equal to your pulse rate, $R_b=R_p$.

You can transmit more than one bit per pulse, though. Consider allowing the pulses to have four amplitudes $+3,+1,-1,-3$. Then you could transmit the sequence $1001$ with only two pulses, of amplitudes, say, $+1,-3$. You're free to map bits to pulse amplitudes in any way you want. In this scenario, you have $R_b=2R_p$.

So, in theory, if there is no noise, then the pulse rate is limited to twice your bandwith, but the bit rate can be as large as you want, by allowing the pulses to take more and more different amplitudes.

Now consider a more realistic scenario, where noise is present. Now, the received signal is not just what you transmitted, but it includes noise too. Assume that you're transmitting with power equal to $S$ watts, and that the noise over your bandwidth $B$ has power $N$ watts. Now your bit rate is no longer as large as you want; it cannot be larger than

$$C=B\log_2(1+S/N).$$

$C$ is known as Shannon's channel capacity, after Claude Shannon. However, keep in mind that transmitting at, or near, capacity may involve a very complex coding system. In practice, a simple system will operate well below the channel capacity. Even very complex systems may not achieve the capacity.

What happens if you only have a simple system, and you attempt to communicate at a significant fraction of the capacity? In that case, you're likely to have errors in your communication. An error happens when the receiver mistakes a bit 1 for a bit 0 or vice versa. Depending on your application, a certain amount of errors can be tolerated. For example, in voice communications one error out of 1000 bits (a bit error rate (BER) of 0.001) is acceptable. Other systems are intolerant to errors: for example, you don't want your hard drive to corrupt your files, not even once a year.

In more complex situations the capacity may be much smaller and the systems that achieve high bit rates become much more complex. For example, in wireless communications the transmitted signal suffers not only from noise, but also from echos and Doppler shift.

In summary, as you see, there is no simple answer. Let me know if you need clarification, or more details.


Downsampling by non-integer factor Filter requirement

If the data is already bandlimited to below half the new sample rate, then a low pass filter may not be required.  If not, low pass filtering is required to prevent aliasing, which can corrupt your result.  Low pass filtering can be combined with Sinc kernel interpolation of all the new points needed at the new sample rate, e.g. only one step, no up/down pair of sample computations required.


Code Snippet to Convert ADC Samples to I/Q

Consider a real carrier at frequency $f_o$ with phase modulation given by
$$s(t) = 2\cos(2\pi f_o t + \phi(t))$$
Where for simplicity of explanation the amplitude was normalized to be $2$.
This real carrier is at positive and negative frequencies as given by Euler's formula, both of which contain the phase modulation (conjugate symmetric):
$$s(t) = 2\cos(2\pi f_o t + \phi(t)) = e^{j( 2\pi f_o t + \phi(t))}+  e^{-j (2\pi f_o t + \phi(t))}$$
If we multiply the signal by $e^{-j2\pi f_o t}$, which is by Euler's formula $\cos(2\pi f_o t) - j\sin(2\pi f_o t)$, it will move the positive frequency to baseband (DC) and the negative frequency to a negative frequency at twice the rate of the carrier (which we must filter out):
$$e^{-j2\pi f_o t} s(t) = e^{-j2\pi f_o t}e^{j( 2\pi f_o t + \phi(t))}+  e^{-j2\pi f_o t}e^{-j (2\pi f_o t + \phi(t))}$$
$$= e^{j \phi(t)} + e^{-j(4\pi f_ot + \phi(t))}$$
Thus the low pass filter is very wide, passing bandwidth $bw$ at baseband and filtering out the signal up at $-2f_o$.
This is the primary task in recovering the signal at baseband (here in this example the complex signal $e^{j\phi(t)}$ with $I = cos(\phi(t))$ and $Q= sin(\phi(t))$. Any further down-sampling to lower sampling rates after the filtering is just best practice since it is desirable in terms of cost, size, and power to run at the lowest sampling rate possible (minimizes resources needed).
In general, if a modulated signal does not have a symmetric (complex conjugate symmetric) spectrum about its carrier then a complex signal must be used to represent that signal at baseband (using two real signals as I and Q), since any real signal has a complex conjugate symmetric spectrum so one real signal alone would not be able to represent the asymmetric spectrum. The GPS C/A code waveform as the OP is using is however symmetric but if there is any residual small frequency offset (due to Doppler and local clock offsets) then there will be a small positive frequency spectrum offset mostly overlapping with the small negative frequency offset with only one real signal at baseband. By having a complex baseband we can represent that spectrum offset in frequency and therefore the receiver can properly remove that- so the complex baseband processing simplified the receivers ability to resolve frequency and phase offsets between the received signal and the receiver’s local reference.
Please see this additional post for further examples of translating to other frequencies with complex signals:
Frequency shifting of a quadrature mixed signal


Show That a 2D Linear Transform $ T \left( \cdot \right) $ Is Homogeneous

The idea here is that any Linear Operator is Homogeneous though not every Homogeneous Operator is Linear.

The classic proof indeed is to build the zero term as a sum (Difference) of 2 other elements in the domain. Then use Linearity to show it must be zero:

$$ T \left( 0, 0 \right) = T \left( {m}_{1} - {m}_{1}, {n}_{1} - {n}_{1} \right) = T \left( {m}_{1}, {n}_{1} \right) - T \left( {m}_{1}, {n}_{1} \right) = 0 $$

  Maybe the better question to ask would be to prove if a transformation is heterogeneous, is it linear?

As I wrote above, if an operator is heterogeneous is must not be linear.

Think of the simplest Linear Operator in 1D: $ T \left( x \right) = a x $. Clearly it is a linear function and indeed homogeneous.
Let's look on Affine Function: $ T \left( x \right) = a x + b $. Usually this is also called Linear Function while it is not unless $ b = 0 $. Why? It is not heterogeneous.

It is easy to see that $ T \left( {x}_{1} - {x}_{1} \right) \neq T \left( {x}_{1} \right) - T \left( {x}_{1} \right) $.


Find the frequency of a pole (Discrete Form)

Usually a discrete-time notch filter is implemented by an IIR (infinite impulse response) filter, which has a numerator as well as a denominator polynomial. The notch frequency - the frequency where the input signal should be completely suppressed - is determined by the numerator polynomial, i.e. by the zeros of the filter's transfer function, not (directly) by its poles. Of course, a practical design will also place the poles close to the notch-frequency, but the actual notch is implemented by a zero on the unit circle of the complex plane, i.e. by the numerator polynomial.

The numerator of a second order notch filter looks like this:

$$N(z)=z^2-2\cos(\omega_0)z+1\tag{1}$$

You can verify that this polynomial has two complex conjugate zeros at $z=e^{j\omega_0}$ and at $z=e^{-j\omega_0}$, so the notch frequency is given by the normalized angular frequency $\omega_0$.


Terminology: spectrum, spectrogram, spectrograph, sonogram, etc

It depends on context.

In signal processing, a spectrum (plural is spectra) shows the frequency content of an entire signal.  It's a 1-dimensional function of amplitude (vertical axis) vs frequency (horizontal axis):

Spectra are often shown with a logarithmic amplitude axis (such as dB), but this isn't necessary.

A machine that produces a spectrum is usually called a spectrum analyzer.  In other fields, the machine is called a spectrograph or spectrometer.

A spectrogram shows how the frequency content of a signal changes over time. It's a 2-dimensional function of amplitude (brightness or color) vs frequency (vertical axis) vs time (horizontal axis):

Sometimes this is called a sonogram.  The time and frequency axes are sometimes swapped.  If amplitude is shown as a 3D surface rather than using color, it's called a waterfall plot.

Confusingly, a machine that produces a spectrogram is also called a spectrograph, or spectrograph is used as a synonym for spectrogram.

Also the line is kind of blurred, because if you view the spectrum of a live signal on a spectrum analyzer, it's displaying the spectrum of small chunks of the signal and you're seeing how it changes over time, which is essentially the same thing as a spectrogram.

I think the important distinction is just the way they're displayed:  A spectrum is a 1D plot and a spectrogram is a 2D plot.


Trouble showing Time Invariance of recursive system

You've actually proved time-invariance already. It's just a matter of clean notation to see this. Let's use $y_1[n]$ to denote the response to a delayed input $x[n-k]$. The sequence $y_1[n]$ satisfies the following difference equation:

$$y_1[n]-4y_1[n-1]+4y_1[n-2]=20x[n-k]+10x[n-k-1]\tag{1}$$

The delayed response to the input $x[n]$ satisfies

$$y[n-k]-4y[n-k-1]+4y[n-k-2]=20x[n-k]+10x[n-k-1]\tag{2}$$

Now we have to check if $y_1[n]$ and $y[n-k]$ satisfy the same difference equation. Comparing $(1)$ and $(2)$ we see that this is indeed the case.


The concept of the bandwidth has several definitions.

For simple filters, one definition of bandwidth is the -3 dB frequency (for a lowpass filter) $f_c$ where the filter gain at $f_c$ is $H(f_c) = H(0)/\sqrt(2)$; i.e. $f_c$ is the half-power or cut-off frequency. This definition could be applied for a bandpass filters too but not for highpass filters. Note again that this definition of bandwidth requires that the filter frequency response be a continuous waveform, otherwise the the frequency will be undefined.

For signals, a 3dB bandwidth concept will not work and instead something like a 90 % energy (or power) definition is invoked. In this case the frequency spectrum of the signal is used to find the minimum frequency $f_b$ at which 90% of the total enegy of the (baseband, real) signal is contained.

In your case you have a signal with a rectangular frequency spectrum, which defies a 3 dB bandwidth, but enables the computation of a 90 % bandwidth definition.

However note that all these definitions make sense under typical waveform assumptions, but for your case neither -3dB nor 90% enegy definitions adequately define the signal's bandwidth.

You signal is bandlimited in the sense that the ractangular function is exactly zero out of its lobe and that frequency shall be taken as the bandwidth of this signal.


Difference Between Convolution and Multiplication of Frequency Response and Frequency Spectrum

Few notes:

Pay attention to the "Factors" in the DFT. Some use $ \frac{1}{2 \pi} $ some use others, pay attention you normalize accordingly.
Multiplication of the DFT is equivalent of Circular Convolution.
You should pad the signals accordingly to get the "Classic Convolution".


It depends on your band-width of interest. If your band-width of interest is very very small than the bandwidth of quantization-noise then you may say that noise is white.


At the end what has proven to be the best solution was onset detection based on either high frequency or energy content. Before it could work I had to use high-pass filter to cut out first 1 kHz, since it contained too much noise.

Once I had noise-only area I could use its profile to reduce noise from rest of the sample.

One library I found particularly useful was Aubio. It has a good set of examples and provides a lot of algorithms to choose from for onset detection.


Relationship between real and imaginary part of a real-valued and causal system

I suppose you've learned about the Hilbert transform relationship between the real part and the imaginary part of the Fourier transform of a real-valued and causal signal. So you can obtain the real part from the imaginary part using this relationship. Since the imaginary part is just composed of sinusoids, the real part can be obtained very easily (note that the Hilbert transform of a sine is a negative cosine).


Upsampling - What purpose does the interpolation filter have?

In step 1 you insert a zero-valued sample between each pair of successive original samples. This is a gain 0.5 operation. For example consider 0 Hz. You can calculate the signed 0 Hz amplitude or direct current (DC) offset simply as the mean of the sample values. Inserting the zero samples halves the mean, so the gain is 0.5, which you need to compensate for by a gain of 2.

Your 2x-oversampled samples represent values of the continuous band-limited signal exactly at the times of the original samples or exactly half-way between them. As you suggest, in step 2 you can use a discrete-time filter, which can be a finite-impulse-response (FIR) or an infinite-impulse-response (IIR) filter. If you would have an arbitrary real oversampling ratio such as 3.12897684124, you couldn't use a discrete-time filter, but would have to resort to continuous-time interpolation, which can be done using piece-wise polynomial splines.


Compressive Sensing Through MATLAB Codes

L1magic. This is the toolbox associated with the original paper.
CompSens. This looks like it's in C, but you could possibly call it with mex -- not sure.
Model-based compressive sensing toolbox.

  Most of the code is plain Matlab code
  Each folder in the package consists of a CS recovery algorithm based on a particular signal model, and a script that tests that recovery algorithm. The names of the scripts typically end with '_example.m'

Manopt. Possibly what you use to actually make the algorithms included in other toolboxes.

  Optimization on manifolds is a powerful paradigm to address nonlinear optimization problems.

Some examples from UMich.
Sparsity.

  This toolbox implements several algorithms to compute sparse expansion in redundant dictionaries and to solve inverse problems with sparse regularization (and also TV regularization).

But all of that, and more, is included in this list of toolboxes.

I've found the hard part is finding psuedocode -- that's where they actually describe the algorithm. Here are some examples of algorithms that included the psuedocode:

one paper that includes the actual algorithm (and the actual code available here in the function IST.
a beginners guide to compressed sensing.
Orthogonal matching pursuit
CS-KLT
Some more algorithms from Microsoft
MCMC for universal CS
An example
Compressive sensing of streams of pulses
GPU Accelerated Greedy Algorithms for Compressed Sensing
Adaptive sensing camera
Compressive sensing meets game theory


Fourier Transform/Series DFT/DFS textbook problem (simple?)

Your result

$$x[n]=\sum_{k=-9}^9a_ke^{j2\pi kn/N}\tag{1}$$

with $N=6$ is correct. If you compare $(1)$ to the IDFT (or DFS)

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi kn/N}\tag{2}$$

you'll notice that in $(2)$ each term $e^{j2\pi kn/N}$ occurs only once for each value of $k$ (for fixed $n$), whereas in $(1)$, due to the periodicity of the complex exponential, each term $e^{j2\pi kn/N}$ occurs several times. E.g., for $k=0$ in $(2)$, in $(1)$ you get contributions for $k=0$, for $k=6$, and for $k=-6$ (because $N=6$). Consequently, comparing $(1)$ and $(2)$ gives

$$X[0]=N(a_0+a_6+a_{-6})$$

For all other values of $k$, the procedure is completely analogous: search for indices $k$ in $(1)$ which have the same complex exponential term $e^{j2\pi kn/N}$. You basically have to keep adding and subtracting the number $N=6$ to the indices, as long you remain in the range $|k|\le 9$, so e.g. for $X[2]$ you get contributions for $k=2$, $k=2+6=8$, and $k=2-6=-4$, because they all have the same exponential term $e^{j4\pi n/N}$.

PS: The term $2\pi$ in your textbook result doesn't make sense to me, it should be $N$, as already pointed out in a comment.


Generating spectrum of a spectrum / identifying frequencies of ripples present in a spectrum

You are not doing anything "wrong" but sooner or later you will end up "re-inventing" the tool that is doing exactly that already. That tool is the Cepstrum.
You wish to estimate the frequency of ripples that show up in the frequency domain and you would like to do this via the Discrete Fourier Transform (DFT) and that is fine.
(Your assumption that if you take the DFT again it will send you back to the time domain is not entirely correct though)
Suppose your time domain signal is $x[n]$ and after the DFT it becomes $X[k]$.
Now, $X[k] \in \mathbb{C}$ but you don't care about the phase of the phase component (hope you understand why I wrote it this way). So, we take the magnitude of $X[k]$, let's call it $Q[k] = |X[k]|$. Notice here, $k$ still denotes frequency, capital letters still denote spectra.
$Q[k] \in \mathbb{R}$, so if you send it through the DFT, it does not care that this used to be complex or the result of another DFT or anything. It still sees a real "signal" and it decomposes it to a sum of sinusoids. If you $W = \mathcal{F}(Q)$, your $W$ will now be complex again, but this time it describes the fluctuations of $Q$.
Now, the thing is that because of the way the sums are put together within the DFT, the dynamic range of $Q$ will be all over the place. Specifically, the low frequencies (DC for example) can end up having HUGE sums and the high frequencies can end up having tiny little values. If you don't account for that, your $W$ will have a strong $\frac{1}{k}$ influence.
To account for this, instead of $W = \mathcal{F}(Q)$, you $W = \mathcal{F}(\log(Q))$. The effect of $\log()$ is to "compress" large values and boost small values. It is the same "trick" we do when we plot a spectrum in a logarithmic scale rather than linear scale.
But, if you do all this (which is a very reasonable thought progression), you have derived the cepstrum. Specifically, the Real Cepstrum.
Depending on your application's context, following up the cepstrum further might provide some more insights or uncover work that has already been undertaken in your field around what you are looking at.
Finally, I would also like to underline what Matt.L says, because you might realise that this "oscillation" is only described by a handful of harmonics and there might be a relationship that captures it that can be optimised by the data itself. This would be very useful, rather than trying to recover just one (dominant (?)) frequency.
Hope this helps.


Whether convolution followed by down sampling equals wavelet decomposition

Simply put, no!

  I know that wavelet decomposition equals convolution followed by down sampling by two

No again. Some linear wavelet decompositions can be turned into a cascade of convolutions followed by downsampling. But wavelets are much more general. Their discrete implementations involve intricate properties in filters, downsampling, and levels.

If you just filter a signal  by some filter, and downsample, in general, you get almost nothing but aliasing.


Why is there a negative in front of the phase response equation for this complex exponential?

For the magnitude, the denominator term will be $$1-0.9e^{-j\omega}=1-0.9(cos\omega - jsin\omega) = 1-0.9cos\omega+0.9jsin\omega$$

Hence the plus sign.

For the phase, if you multiply the numerator and denominator with the complex conjugate,

$$H(e^{j\omega}) = \frac {1-0.9cos\omega-0.9jsin\omega}{(1-0.9cos\omega+0.9jsin\omega)(1-0.9cos\omega-0.9jsin\omega)}$$

The denominator will become real.
Therefore,
$$\angle H(e^{j\omega}) = arctan(\frac {-0.9 sin\omega}{1-0.9cos\omega})$$
Taking the minus sign outside is how you get the negative in the phase.


What does double hat mean in Spectral theory?

It seems this notation is used for empirical transfer function estimates (see pg. 207, 8.4 Fourier Analysis). To wit, on pg. 79 of Modeling: identification and simulation of dynamical systems we have

  The double hat denotes that the ETFE is a very rough estimate of the transfer function.


Calculate bandwidth of a signal for Nyquist–Shannon sampling theorem

So apply the following principals to get to a solution:

Shifts in the time domain don't affect the bandwidth, so the $-4$ and the $+3$ can be ignored.
Reversal of the time domain axis doesn't affect the bandwidth, so you can ignore the negative sign on the $-2t$.
The $\mathrm{sinc}()$ function has a Fourier transform that is a rectangle function of a particular width.
Scaling the time axis (lengthening or contracting) of a function, scales the frequency axis in the opposite manner (contracting or lengthening) for the Fourier transform.
Multiplication of functions in the time domain, results in the convolution of those two functions' Fourier transforms in the frequency domain.
The convolution of two Fourier transforms in the frequency domain, will result in a Fourier transform whose bandwidth is the sum of the bandwidth of those initial two Fourier Transforms.

That should give you enough to come up with an answer.


Can an audio recording tell whether it is made in a large or small room?

There could be a number of approaches for this task, depending on the number of channels (microphones) in your input and structure of available data.

Given a labeled data set from a previously known set of rooms, you can train a neural network (NN) to classify in which one of them a new sample was recorded.

A second approach would be to estimate the parameters of the room impulse response (RIR). Certain methods could be under the limitation of a single microphone. A microphone array would supply more information like in this paper. Depending on the shape of the array, different methods can be utilized.

Also, do you have any other prior knowledge on the room? For example, can you get a clean (without speech) audio segment? Do you know where is the microphone and sound sources are located on the room? These priors can point you to methods that utilize this information and improve the estimation.


Finding zebra-like pattern in image (Detection of structured-light fringe centerline from photo)

I suggest the following steps:

Find a threshold to separate the foreground from the background.
For each blob in the binary image (one zebra stripe), for each x, find the weighted center (by pixel intensity) in y direction.
Possibly, smooth the y values, to remove noise.
Connect the (x,y) points by fitting some kind of curve. This article might help you. You can also fit a high-level polynomial, though it is worse, in my opinion.

Here is a Matlab code that shows steps 1,2 and 4. I skipped the automatic threshold selection. Instead I chose manual th=40:

These are the curves that are found by finding the weighted average per column:

These are the curves after fitting a polynomial:

Here is the code:

function Zebra()
    im = imread('http://i.stack.imgur.com/m0sy7.png');
    im = uint8(mean(im,3));

    th = 40;
    imBinary = im>th;
    imBinary = imclose(imBinary,strel('disk',2));
    % figure;imshow(imBinary);
    labels = logical(imBinary);
    props =regionprops(labels,im,'Image','Area','BoundingBox');

    figure(1);imshow(im .* uint8(imBinary));
    figure(2);imshow(im .* uint8(imBinary));

    for i=1:numel(props)
        %Ignore small ones
        if props(i).Area < 10
            continue
        end
        %Find weighted centroids
        boundingBox = props(i).BoundingBox;
        ul = boundingBox(1:2)+0.5;
        wh = boundingBox(3:4);
        clipped = im( ul(2): (ul(2)+wh(2)-1), ul(1): (ul(1)+wh(1)-1) );
        imClip = double(props(i).Image) .* double(clipped);
        rows = transpose( 1:size(imClip,1) );
        %Weighted calculation
        weightedRows  = sum(bsxfun(@times, imClip, rows),1) ./ sum(imClip,1);
        %Calculate x,y
        x = ( 1:numel(weightedRows) ) + ul(1) - 1;
        y = ( weightedRows ) + ul(2) - 1;
        figure(1);
        hold on;plot(x,y,'b','LineWidth',2);
        try %#ok<TRYNC>
            figure(2);
            [xo,yo] = FitCurveByPolynom(x,y);
            hold on;plot(xo,yo,'g','LineWidth',2);
        end
        linkaxes( cell2mat(get(get(0,'Children'),'Children')) )
    end
end

function [xo,yo] = FitCurveByPolynom(x,y)
   p = polyfit(x,y,15);
   yo = polyval(p,x);
   xo = x;
end


how to use image filters for spectrograms?

A spectrogram of a monochannel (real) audio file would generally be a complex 2D matrix, hence a sort of complex image.  So you can use 2D filters directly on the complex matrix in the workspace. However, filtering complex data can be complicated, due to the non-linear coupling between the phase and the magnitude: if you want to preserve the Hermitian properties, which is useful when you want to invert the spectrogram into a real (filtered) signal. So traditional algorithms like  enhancement, shrinkage, adaptive filtering, masking are possible.

Thus, no need to go through an intermediate binary format like:

  J = imread('example.jpg');

that will quantize or compress values and somehow loose the complex aspect.

Note that processing complex values requires some care. In many situations, people only want to analyze the magnitude of the spectrogram, which then become a mere positive real-valued 2D matrix. Then, all sorts of 2D processing are available. In that case, you should take care of preserving the positivity. Most linear filters with non-negative coefficients will do the job. Non-linear processing tools based on histograms, on ranking, may require a simplification of the floating-point values.


How can I increase the zero-padding on the DFT?

Your code does not implement the zero padding method. Below is the corrected version that does. Note that zero padding does not improve frequency resolution. It mereley provides interpolated samples of the original spectrum.

   clear all
clc

fs = 50; %50Hz
T = 1/fs;
N = 8;
t = (0:N-1)*T;

y = 5+cos(2*pi*12.5*t)+sin(2*pi*18.75*t);

M = 8;
K = M*N;       % DFT increased size
yk = [y zeros(1,(M-1)*N) ];  %zero padded input

X=zeros(1,K);   % high density

for k = 0:K-1
    for n = 0:K-1
        X(k+1) = X(k+1) + yk(n+1)*exp(-j*(2*pi/K)*k*n);
    end
end

x_mag = abs(X);

f_2 = fs*(0:(K-1))/K;
stem(f_2,x_mag,'-ok','linewidth',2,'MarkerFaceColor','black');