How do you simultaneously undersample and oversample?

Undersampling is also known as "bandpass sampling" and "IF sampling".  Consider the sampling of an 11 Hz sine wave and a 1 Hz sine wave, both sampled by a 10 Hz sampling clock as in the figure below:
The samples that are produced (the red circles) for both waveforms are identical. Thus the 11 Hz "Intermediate Freuency" or IF has been downconverted in the sampling process to 1 Hz. We can use this technique to undersample a bandpass waveform, and thereby down-convert it at the same time (putting our waveform of interest in the first Nyquist zone: ±$F_s/2$ where $F_s$ is the sampling freqquency).  A key requirement that satisfies Nyquist is that the sampling rate used must exceed twice the bandwidth of the signal (note I said bandwidth, not highest frequency!).
That is "Undersampling" as most often described. To oversample a waveform, we sample at a rate significantly higher than the Nyquist criteria, with one motivation to increase the SNR as limited by quantization noise (by sampling hgiher, the same quantization noise is spread over a wider digital frequency, which we can subsequently filter down to our bandwidth of interest and thereby eliminate a significant amount of the noise).  Therefore to both undersample and oversample a waveform, we choose a sampling frequency that is significantly larger than the signal bandwidth (thus oversample), and then use a multiple of this sampling clock to down-convert the waveform of interest positioned at a higher IF frequency (undersample).
For example, consider a waveform that occupies 10 KHz of BW that resides at 20.25 MHz. If we sample this with a 10 MSps sampling clock, we would be both undersampling and oversampling at the same time. (The undersampling creates a digital waveform with 10 KHz of BW at 250 KHz, which is oversampled at 10 MSps).


The reason one uses a complex signal at baseband is so that one can distinguish between positive and negative frequency...otherwise, when the signal is mixed to a higher frequency, it will have two identical sidebands, taking up twice the spectrum.  Without knowing what you're modulating, I can only offer general advice...take your real baseband signal taking up bandwidth 0..Fmax, and multiply it by exp(2*piiFmax/2*t) (note the 2s cancel out, but the form is more recognizable).  This will heterodyne it down to a complex signal occupying -Fmax/2..Fmax/2, and that can be passed to the RF chip.


Methods for aligning FFT bins with an FSK signal

Before I address your question, I'll suggest that using a frequency domain technique (FFTs) to demodulate inherently time domain data (a time series of sequential data symbols) won't be very practical.  Symbol timing recovery is inherently a time domain process. ISI encountered in common modulations like GFSK/GMSK means you really have more than 2 frequencies for 2FSK symbols at the receiver.

Now on to your question of how to perform frequency alignment.

You requested algorithms for what is normally called Carrier Frequency Offset (CFO) estimation and correction.

In communications receivers, there are sometimes two stages of CFO estimation and correction: coarse and fine.  For 2-FSK burst communications, a single stage of coarse CFO correction is usually all one would need for non-coherent FSK demodulation.

CFO correction is generally performed as a (cyclic) digital spectrum shift by multiplication with a time-varying complex exponential in the time domain:

$$x_c[n] = x[n]\space e^{j\pi \frac{-f_{offset}}{F_s/2}n}$$

where $f_{offset}$ is the estimated CFO, and $F_s$ is the sample rate.

Although for FSK, if one has already computed the baseband waveform from the FSK-modulated symbols (i.e. the output of an FM discriminator), a simple level shift of the basedband waveform, to have it's peaks and troughs centered about 0, is equivalent to CFO correction.

CFO estimation algorithms can be divided up into two broad classes: Open-loop feed-forward algorithms and closed-loop feedback algorithms.  Closed-loop algorithms work a sample at a time, are usually used for fine frequency offset estimation, and usually involve some arrangement of PLLs.  Open-loop algorithms work on a block of samples at a time and are used for coarse frequency offset estimation.

You probably just want to use an open loop instantaneous frequency estimate for computing your CFO estimate.  The instantaneous frequency estimate of an FM/FSK waveform at any given sample is given by

$$y[n] = \dfrac{F_s/2}{\pi}\mathrm{arg}\left(x[n]x^*[n-1]\right)$$

Normally to generate an CFO estimate, one computes the instantaneous frequency over several symbols and looks at the average peak and trough levels to estimate how much the waveform level needs to be shifted for the waveform to be centered about 0.

In your proposed case of having all the samples for one symbol, assuming you don't have a terribly large CFO, computing the instantaneous frequency for that one symbol will tell you what 2FSK symbol you have without performing the FFT.


How to recover frequency and damping coefficient from superposition of damped oscillators?

If you are absolutely sure that you have a signal that is the superposition of damped oscillators then you can simply track their evolution in time.

Fs = 2000; % Sampling frequency (Hz)
T = 4;     % Time duration (seconds)
t = 0:(1./Fs):(T-(1./Fs)); % Time vector (in seconds)
p = 2.0.*pi.*t; % Phase vector
S = 1.0*sin(440.*p).*exp(-t) + 0.25 .* sin(800.*p).*exp(-t*10); % Signal

This signal looks like:

And via a Short Time Fourier Transform (STFT) (surf(stft(S))) it looks like this:

Granted, if you have beating frequencies, the task is getting a bit more challenging. For instance:

S = 1.0*sin(440.*p).*exp(-t) + 0.25 .* sin(442.*p).*exp(-t*2) + sin(800.*p).*exp(-t*10);

Would show up as:

But once you see oscillations on the damped decay, you know you have beating and you can recover how close the beating pair is by measuring the beating frequency and how quickly it fades (to derive its decay).

In a similar way (i.e. using the STFT) you can also recover the phases of the sinusoids.

Hope this helps.

UPDATE:

The STFT of the signal that was shared via pastebin is capable of distinguishing three components in this signal. This is as follows:

%Given a signal in s and that the signal package is installed (in Octave)
[S,F,T] = specgram(s);
surf(abs(S));

This returns:


You need to generate early reflections with a few taps of delays (= convolution with the sum of a handful of diracs) ; and the "tail" with what is usually implemented with a network of all-pass (AP) and comb filters.

The first part is trivial to implement but difficult to get to sound right. It might help to look at the positioning of peaks at the head of pre-recorded impulse responses to get an idea of what kind of responses are "natural". It's hard to make it easily parametrizable by a few settings, though you could get away with it by offering a bunch of presets for this part to get an "hybrid" reverb à la Virsyn Reflect.

The tail is where algorithmic reverbs shine (no pun intended), since it is easy to parametrize. The grandfather of algorithmic reverbs is Schroeder's algorithm. Note that it only generates a "tail", not early reflections - you could augment it with a few taps of delays to give it a head. It does not sound very good (very "grainy") but it's a decent start - with that you have something to mess with and it helps understanding the impact of each parameter. Many praised algorithmic reverbs, especially those from the 80s (Lexicon, Eventide, Publison) are just finely tuned topologies of all-pass filters and comb filters. I suspect this involved a lot of trial and error by trained ears who got very good at figuring out how changes in parameters or topology would sound, rather than any scientific approach. Here is another interesting read - showing Keith Barr's favorite reverb topology. The original midiverb might not be your definition of "high quality algorithmic reverb", but it sounded sweet, and it puts thing in perspective to know it was made with a "dsp" only capable of MACs with coefficients of 0.5...

Valhalla room is my favorite software algorithmic reverb - take some time to read its developer's blog to find inspiration about the kind of things that help getting it right.


I want to know if there is a practical way to find out if multipath propagation is significant for audio in a given environment?

Never done this but it might work. Record the audio simultaneously using multiple microphones located in the neighborhood of your main recording position. If the cross-correlations between signals from different pairs of microphones show peaks at different time lags than the autocorrelation, then there must be either multiple sources or multiple paths. Be careful not to confuse multiple paths with reflections from the microphone bodies.


Does orthogonal beamforming mean $\mathbf f_a \cdot \mathbf f_b =0$?

Orthogonality might be considered in literature as background knowledge. The definition is somehow a converse version. If $\cdot $ denotes an inner product (or a dot or scalar product), and whenever
 $\mathbf  f_a \cdot \mathbf f_b = 0$, then $f_a$ and  $f_b$ are said to be orthogonal.

Note that, if $\mathbf  f_a$ is zero, it is orthogonal to everything.


In more standard DSP terms, you have the following filter:

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

where $x[n]$ and $y[n]$ are the input and output signals at time $n$ respectively.

The transfer function (which you didn't ask for) is:

$$
H(z) = \frac{1-a}{1 - az^{-1}}
$$

so here is your single pole, at $z=a$ in the complex plane. This filter is also known as exponential smoothing, exponential moving average (EMA), or exponentially weighted moving average (EWMA).

The infinite impulse response is $h[n] = (1-a) u[n] a^n$. In layman terms, when the input signal is 0 except for $x[0]=1$, the output signal is an exponential $(1-a) \times a^n$ starting at $n=0$.

What you want is the step response (i.e. what happens if the input signal is a constant $K$ starting at time $n=0$).

In this case, the output signal is the convolution $h \star Ku$ of the impulse response and the step signal. This is (for time $n \ge 0$):

$$
y[n] = \sum_{k=0}^{k=n} K \times h[k] = K(1-a) \sum_{k=0}^{k=n} a^k
= K(1-a)\frac{1-a^{n+1}}{1-a}
= K (1-a^{n+1})
$$

As the time $n$ grows, $a^{n+1}$ vanishes, and the step response grows monotonically to its limit value $y = K$, which is the value of the input signal.

This is what you get in your simulation.


Mulitplication in Frequency Domain as a Similarity Measure of Audio Clips

The fft product alone is not the convolution, but the frequency domain of the convolution. To complete the operation the OP must also take the inverse FFT to get a circular convolution result.
$$CONV = \text{ifft}(\text{fft}(a) \text{fft}(b))$$
However, similarity would be determined using correlation not convolution. To do this, simply complex conjugate one of the FFT results as follows:
$$XCORR = \text{ifft}(\text{fft}(a) \text{fft}^*(b))$$
The above is the cross-correlation function (using circular convolution). The result is the correlation of $a$ and $b$ at repeated circular shifts in time.


How does sensor fusion work? I want to understand the math/physics/algorithms

Sensors provide analog or digital insights into a reality that is difficult to grasp. Just like your senses do.

Sensors are designed to seize a special portion of a measurable phenomenon: electrical, chemical, physical...  Just like your senses do. The usual five ones are sight, hearing, touch, smell and taste, but there are more to human skills.

Understanding a physical phenomenon via signal processing is like trying to tell the main ingredients from the  recipe of the meal you ordered in a restaurant last evening. You do not know what happened in the kitchen, yet you have chosen what to order, and have seen, tasted, smelled, even touched the plate.

Each of your senses and understanding provided you with hints or knowledge, but you can never be sure, as some cooks master sense deception. What you see as an artifial green marmelade made of strawberries can taste like apple jelly.

Based on you cooking expertise, the combination of all your senses can bring you closer to the the actual recipe, provided that:

you can use your senses often enough: if you take only one bite every 10 days, it is unlikely you get the menu. This is data sampling;
each sense is sensitive enough. This is sensor sensitivity;
your sense list is close to complete, and the coverage is sufficient for your purpose. If you cannot taste sugar anymore, you won't be able to cook some meals for people who still can taste it. This is measurement span;
you are able infer, or model some of the processes you are unaware of. This is modeling.

Sensor fusion is the art and science of combining sensory data, knowledge and models from disparate sources so that the resulting information has more validity or less uncertainty than individual  sources.

Examples: The act of averaging $N$ sensor data for the same deterministic signal, with different stochastic realizations of a noise term, is the most basic sensor fusion operation. With the model of independent Gaussian noises of deviation $\sigma$, the average theoretically yields a $\frac{\sigma}{\sqrt{N}}$ deviation(less uncertain). The combination of three channels (Red, Blue, Green) offers a color image (more valid than a monochannel gray-scale).

The first example combines information from a single sensor model. The second one from sensors operating in the same domain (electromagnetic waves), yet in different portions of the spectrum. In general, fusion operates on different sensors with different rates, ranges, domains and mostly units.

The problem is highly dependent on the phenomenon you are looking at, the available sensors, and the information you are looking for.

The physics tell you the potential information you can get from your sensors. The mathematics can model how they are related or complementary, or what is irrelevant information (noise). The algorithms will combine the previous knowledge as optimally as possible, in terms of precision, accuracy or speed.

The topic is related to the realms of Sensor fusion, Data fusion or Information integration, with a short overview in Principles and Techniques for Sensor Data Fusion. Many more books are available, as in Best book for learning sensor fusion, specifically regarding IMU and GPS integration.

Regarding your actual problem, a first step would consist in understanding What are the differences between a gyroscope, accelerometer and magnetometer? which could help you push techniques a little further. And achieve the goal of fusion: using sensor differences, sum their data, in the most clever way.


How do you change pitch with wavelet transforms in MATLAB?

I haven't implemented a pitch shifter using wavelets, but I am currently learning about them. From what I understand scales are the inverse of frequency, therefore to get your desired result you'd have to divide the scales instead of multiplying.

Furthermore, according to 'Applications of the Continuous Wavelet Transform in the Processing of Musical Signals', you also have to update the phase as adjusting the scales isn't enough. The algorithm outlined in the paper is:

semitone = 12;
pitchFactor = 2^(semitone/12);

wt = cwtft(input);

mag = abs(wt.cfs);
phase = angle(wt.cfs);
unwrapped = unwrap(phase);

wt.cfs = mag.*exp(1i*unwrapped*pitchFactor);
wt.scales = wt.scales./pitchFactor;

result = icwtft(wt);

Sadly, that's all I can offer until I learn more.


How to, efficiently, identify telephone number from the key pressed data audio file?

If you interested in computational efficiency then performing FFTs is not the best solution. Here's something you might try: apply the samples of your .wav file to a parallel bank of 8 Goertzel filters. (One filter for each of the 8 possible tone frequencies.) For example, below is the Goertzel filter for detecting the 770 Hz tone.

Monitor (measure) the output powers of the 8 filters every, say, 25 milliseconds. If the sample rate is 8000 Hz, then 25 msec means measure the 8 filter output power values every 200 input samples. So, in the figure below, apply the first 200 input samples to the filter and compute a single v(n) value. Then, and only then, compute a single complex-valued y(m) sample and single real-valued P770(m) sample. For the first 200 input samples do the same for the other 7 filters. At the end of the first 200 input samples see if two the eight Pxxx(m) power sample values exceeded some threshold value. If so, then you know a phone button key was pressed. Using some sort of logic, decode the two high-level Pxxx(m) power values to determine which key was pressed and store that single 'key number'. If no key was pushed during those first 200 input samples set the stored 'key number' value to Nan in Matlab.

Repeat the above processing for every block of 200 samples in your .wav file. If your .wav file had 16,000 samples then you should end up with a sequence of 80 'key number' values. Plot those 'key number' values versus the time index m.


You are correct. There is a trick you can use to code this process very easily. Assuming Matlab syntax, define a rectangular impulse response with M ones:

ir = ones(1,M);

Then, filter the upsampled train of impulses x:

zoh = conv(x,ir);

and you're done!


Restore real signal from its complex representation if sampled under Nyquist frequency

You simply can't.
The analytic signal contains the exact same amount of information as the original real-valued signal, and if either is undersampled, you get the loss of information that undersampling gets you.
That's it. Unless you've got any secondary information about the signal that you forgot to mention, there's nothing you can do.


How to mix 8bit/8kHz unsigned linear PCM data in AppsScript?

Assuming that abyte[] and bbyte[] are two arrays with samples valued from 0-255, this is the solution:

combined = [];
for(i = 0; i < abyte.length; i++)
{
    sample =  2*(abyte[i] + bbyte[i]) - ((abyte[i]*bbyte[i])/128) - 256
    combined[i] = Math.round(sample);
}

I have a full example with both wav files extracted into arrays (I did this using matlab).  Also I verified that it plays back fine using matlab.

https://jsfiddle.net/4hsfo96g/


How to deal with overlap in STFT when considering the highest energy among multiple channels

The obvious approach to this would be a crossfade between channels. This leads to problems, just as obvious, as you would get nasty comb filter effects at 2000Hz and above. So this is a no-go.

A solution just as simple would be to use just one channel, but apply a gain at each window, so that the energy is aligned with the one of the loudest channel at that point in time, i.e. the gain would be
$$a_i = \frac{\max\limits_{c}(E_i^c)}{E_i^1}$$
for $c$ as channel index, $i$ as time index and assuming you use channel $1$.
This way, you avoid any phase problems whatsoever and still have a signal, that suffices your energy requirement. Also, there is no non-linear processing and thus, no artifacts resulting from it, which would most certainly occur, if you follow some approach with time aligning the signals from different channels to avoid comb filter effects.

EDIT: I know that the overall SNR will not be improved by this approach, but I think the advantages are worth it, since SNR differences should be quite small, as distance between microphones is also small.


What distance metric can I use for comparing image features like elongation and solidity of each image?

If I understood you correctly, each contour is described by a 2-element feature vector $f = [e, s]$, where $e$ is elongation and $s$ is solidity.

In that case, you might want to try the Mahalanobis distance, which is defined as follows:
$$d(f_1, f_2) = \sqrt{(f_2 - f_1)C^{-1}(f_2 - f_1)}$$
where $f_1$ and $f_2$ are the feature vectors that you are comparing, and $C$ is the covariance matrix of your data set.


The following is all assuming we are interested in Finite Impulse Response (FIR) filter design.

The cutoff frequency is not really the determining factor in the filter length. The 3 main things that determine the filter length are: the amount of passband ripple, the amount of stopband attenuation and the relative size of the transition band, i.e., the size of the band going from the passband to the stopband. Cutoff frequency will affect the design if it is really close to 0 Hz or half your sampling frequency.

The smaller the size of the transition band the longer your filter is going to be. To specifically answer your questions.

The sinc filter design is largely by trial an error to achieve your desired specifications. You play with the length of filter, the type of window used, and the filter cut-off to try to meet your specifications. So it may not be the most suitable in terms of trying to meet your spec. That said, it will give you an idea of the size of filter that would be required. Other design techniques are likely to require filters of similar lengths.

Regarding the transition width and sampling frequency, as I mentioned above it is actually the opposite. The smaller the transition width, the longer your filter is going to need to be.

One alternative, is to actually downsample your signal to reduce your sampling frequency. Since you are attenuating all signals above 70 Hz, why not just use decimation? The advantages of this technique are:

Reduced computational load
The filters used can be much smaller because they can have larger transition widths.
If really needed, you can upsample to get the data back to the original sampling rate using the same filters you used for the downsampling.

The main disadvantage is that your overall design is a bit more complicated. A good reference for this is "Multirate Signal Processing" by Crochiere and Rabiner.

If you are looking for the smallest filter lengths, then using the Parks–McClellan algorithm based on the Remez exchange technique has the shortest filter lengths for meeting passband ripple and stopband attenuation levels.


Sparse recovery, Restricted Isometry Property for ILL-POSED problems

The purpose of compression sensing is to sample a signal below Nyquist given that the signal is sparse in a certain transformation. So it only makes sense when  $M < N$. Othwersie if $M>N$ or $M=N$ then what you have is conventional sampling anyways. An important property of RIP matrices is that it is guaranteed to only change the length of any vector  "very little" as long as the vector  is at least $K$-sparse (has at most $K$ non-zero coefficients. By definition its is the following: $$(1-\delta_S)||x||_2^2 \le ||A x||_2^2 \le (1+\delta_S)||x||_2^2$$ where $x$ is the sparse vector and $\delta_s$ is an arbitrary constant close to zero. So yes, $M$ should be less than $N$ for compressive sensing to make sense and infact all the known matrices that follow RIP are known to exist when $M = O(K(log(N/K))$ strictly. Finding these matrices is not trivial and is a NP Hard problem.

So you do not need to consider RIP for A when $M>=N$, because we are now not talking in compressive sensing terms. The $l_2$ terms should now be the objective because convergence can be achieved in the $l_2$ convergence.

If $M>=N$ and $x$ is not sparse then the regularization could still work but $\lambda$ should be adjusted accordingly with apriori knowledge of $x$ not being sparse. Although I dont see any point of the $l_1$ being part of the objective.


Why does speeding up audio files with VLC Player not match with doing the same thing with RPG Maker or Audacity?

i don't have any of those products, but evidently the VLC Player does what we call "time-compression" when you speed up.  and likely "time-expansion" when you slow down.  together we call those operations "time-scaling".  time-scaling and sample-rate-conversion, together, make a pitch-shifter.  so a time-scaler can change tempo without changing pitch, a sample-rate converter can change tempe and pitch together, and a pitch shifter can can pitch without changing tempo.

the other products, RPG and Audacity, likely just use the same interpolation we use for sample rate conversion to do the speeding up.   essentially the same as speeding up playback of analog tape.


Fourier Slice Theorem - Reconstruction Fourier Space

Fourier slice theorem states that Fourier Transform of your projections are equal to slices of 2D Fourier transform.

When you sample the projections over discrete angles the FT of projections become samples of 2D Fourier transform.

So you have to use your obtained samples to interpolate the remaining points.

in your code you put the samples in their position but the remaining points remain zero.


You got it mostly backwards, but otherwise OK :)
This can be pretty directly answered by writing down what the Cooley-Tukey FFT's twiddle factor $W_N$ is:
$$W_N=e^{-i\frac{2\pi}{N}}$$
and that's it. For example, in your picture, in Stage 0, you're multiplying with $W_2^0=\left(e^{-i\frac{2\pi}{N}}\right)^0=e^{-i\frac{2\pi0}{N}}=e^0=1$.

So k is the index number of the iteration

No, $k$ is the exponent of your twiddle factors. It goes from 0 to half the size of your transform (in each stage).

From the image is the first stage N = 8 (since there are 8 butterflies) or is N = 2

Stage 0 has subtransforms of size N=2, as you can see. (Also, just write down the 2-DFT's formula and compare.)


Depending on the actual recordings, the algorithm complexity could range from dead easy to really complex...

I'll take the studio recording case first, so I can assume :
- (Almost) no noise coming from outside (cars, trucks, bus...)
- Nobody slamming the door in the middle of the recording
- Voice samples are recorded at optimal level independently of who is speaking (some ppl speak very low, others loud and clear)

If that's the case, then you could implement the logic of what sound engineers call a noise gate : below a certain audio power threshold, it's silence, otherwise it's not.
You'd have to calculate mean power over a very small window and check it against the threshold. Probably take safety margins around "speeched" areas so you don't cut the start/end of the words.

Determining the correct threshold is a manual process involving trials and errors.

Noise gates sometimes use 2 different thresholds : one to determine speech has started and another to determine speech has stopped. The second one being lower than the first one so it doesn't stop in the middle of a word.

Determining threshold(s) for different files recorded at different levels (because of different speakers) might become tricky. One way to help is to first normalize the volumes : i.e. check for the max value and apply gain in order to bring this value to the max authorized one in the range you're working with (might be -1.0 to 1.0 or 0.0 to 1.0 or signed/unsigned integer values).

Another thing to define is the minimum length of "speeched" areas and the minimum length of silent areas.

If there's a bit more background noise (and I assume the noise volume doesn't change too much over time), you could use the property of white noise when calculating autocorrelation factors : it's mean value goes close to zero

Here is an example of "brute force" autocorrelation

int n = vocalSample.length;
for (int j = 0; j < n; j++) {
    for (int i = 0; i < n; i++) {
        correllation[j] += vocalSample[i] * vocalSample[(n + i - j) % n];
    }
}

then calculate mean value of the array.

Do this for very small arrays (e.g. 1ms at 22KHz -> 22 values), otherwise it will take ages.

I actually used this technique with some success on Ted talks to remove silence in order to perform speaker identification.

If you have hours and hours of recording, you might consider multithreading your calculations :-)
In my calculations (remove noise, normalize, extract vocal print), I crunched 10Gb in about 8 minutes using 8 worker threads.

State-of-the-art vad is something beyond what could be explained here, here is a link to an extremely interesting IEEE article on the subject : http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6403507&punumber%3D97


If I understand your question correctly you want this (index | binary | bit rev. | bit rev. index):

$$
0\quad        000\quad      000\quad 0\\
1\quad        001\quad      100\quad 4\\
2\quad        010\quad      010\quad 2\\
3\quad        011\quad      110\quad 6\\
4\quad        100\quad      001\quad 1\\
5\quad        101\quad      101\quad 5\\
6\quad        110\quad      011\quad 3\\
7\quad        111\quad      111\quad 7\\
$$


What are the Even and Odd Components of a 2D Signal

The even part of a 2D image $x[n,m]$ is defined as:

$$x_e[n,m] = \frac{x[n,m] + x[-n,-m] }{2} $$

and the odd part is
$$x_o[n,m] = \frac{x[n,m] - x[-n,-m] }{2} $$.


Why Does ASK Modulation Create Fourier Sidebands?

If you did a continuous on off keying of a 10101010... pattern, then you would see sidebands as described since this is simply an up-conversion of the Fourier Transform of a 50% duty cycle square wave (moved to any carrier frequency). However if the data pattern for this case of a rectangular on-off keying was random, the resulting spectrum would be continuous and under the envelope of a Sinc function (Sinc squared as a power spectrum). The first case itself is also under the envelope of the same Sinc function in frequency, but since the repetition rate of the underlying pattern was not random (in this case the single pulse on off cycle) and repeated at a specific repetition rate, only integer harmonics of that repetition rate can be in the spectrum, with all the even harmonics falling on the nulls of the Sinc envelope (and thus a square wave with 50% duty cycle only has odd order harmonics).
This is demonstrated below, showing a single rectangular pulse and the magnitude (in dB) of its Fourier Transform. If this was a carrier frequency turned on and off, the same spectrum would result just at that particular carrier (DC is just another carrier frequency). When we repeat that pulse in time at a periodic rate, the spectrum can only have non-zero content at integer multiples of that periodic rate (just as implied with the Fourier Series Expansion), but will still be within the envelope of the base pulse (in this case a Sinc function).  So in the second case it is repeated at a 50% duty cycle, such that the even harmonics are in the nulls of the Sinc (the waveform repeats at $1/(2T)$). The third case is the same pulse duration repeating at a 25% duty cycle, so the 4th harmonic lands on the first null (with the waveform repeating at $1/(4T)$).

For example, if we generated a long pseudo-random sequence of such pulses with a pattern that repeats once per second, with each pulse having a duration of 1 millisecond, we would see spectral lines every 1 Hz under the envelope of a Sinc function with it's nulls at 1 KHz spacing. Thus we see in the extreme case of a random pattern that never repeats, that the spectrum will again be continuous such as in the first case of a single pulse (with the primary difference that it will be a power spectrum rather than an energy spectrum since if the pulse continues for all time, it will have infinite energy but finite power).


The multiplication of the two noisy signal gives

$$(x_1+n_1)(x_2+n_2)=x_1x_2+x_1n_2+x_2n_1+n_1n_2=x+n\tag{1}$$

with the desired signal

$$x=x_1x_2\tag{1}$$

and the noise part

$$n=x_1n_2+x_2n_1+n_1n_2\tag{2}$$

Assuming all signals are independent of each other and have zero mean, we get for the signal power

$$\sigma^2_{x}=\sigma_{x_1}^2\sigma_{x_2}^2\tag{3}$$

and for the noise power

$$\sigma_n^2=\sigma_{x_1}^2\sigma^2_{n_2}+\sigma^2_{x_2}\sigma^2_{n_1}+\sigma^2_{n_1}\sigma^2_{n_2}\tag{4}$$

For the total SNR you get

$$\text{SNR}=\frac{\sigma_x^2}{\sigma_n^2}=\frac{\sigma_{x_1}^2\sigma_{x_2}^2}{\sigma_{x_1}^2\sigma^2_{n_2}+\sigma^2_{x_2}\sigma^2_{n_1}+\sigma^2_{n_1}\sigma^2_{n_2}}\tag{5}$$

With $\text{SNR}_1=\sigma_{x_1}^2/\sigma_{n_1}^2$ and $\text{SNR}_2=\sigma_{x_2}^2/\sigma_{n_2}^2$ this can be rewritten as

$$\text{SNR}=\frac{\text{SNR}_1\text{SNR}_2}{\text{SNR}_1+\text{SNR}_2+1}\tag{6}$$


Relation between Bandwidth and Baud-Rate for 8-PSK

The final statement is correct, for PSK with proper pulse shaping the baud rate and the bandwidth are the same (the bandwidth will typically be 20-30% higher than the symbol rate, but read on). The baud rate is the symbol rate or the rate of change from one symbol to the next in the modulated constellation.

Your earlier figures suggest unfiltered waveforms and "bandwidth" definition to be the first main lobe of the resulting Sinc function spectral shape. Specifically, with a random data pattern, the spectrum for an unfiltered PSK waveform will be a sinc function with a null to null bandwidth from the first null that is twice the baud rate.  (and power spectral density will be sinc squared). Perhaps this is the source of your confusion.

Once filtered (typically with raised cosine or root-raised cosine filters or other similar Nyquist filters) the bandwidth can approach the symbol rate, but will be slightly wider given the parameters of the Nyquist filter. (For an infinitely long filter the bandwidth can be exactly the symbol rate but that is not practical).

This also applies for all the same reasons for QAM which is shown in the figure below (I happen to have figures for QAM that help demonstrate this), and when you see and understand the additional explanation below you will know how it can be applied to determine the bandwidth of other modulated waveform. The figure below shows a QAM waveform versus time (the real portion of the waveform) as an unfiltered waveform in blue versus what it would look like with raised cosine pulse shaping in red, and the resulting spectrum for each case is shown directly below that. This waveform was modulated with a baud rate of 1 MS/s.

The unfiltered view is clear from observing the resultant modulated signal as a stream of impulses that are the convolved with a rectangular pulse that is the width of your symbol period. Again see the example related to QAM below, which equally applies to PSK. If the waveform was not pulse-shaped, the stream of impulses representing the modulation symbols would be convolved in time with a rectangular function of width $T$ (where $T$ is the symbol period) to result in the blue waveform above. Note for PSK this stream of impulses is complex values all with same magnitude by varying phase from symbol to symbol. The stream of impulses if indeed completely random and independent from symbol to symbol would have a white spectral density. Since the impulses are convolved with a rectangular function in time, its spectrum is multiplied by the Fourier transform in frequency, which for a rectangular function is a Sinc with a first null at $1/T$! Thus it should be very clear that the spectrum for a unfiltered PSK or QAM modulated waveform will be a Sinc function with a main lobe of twice the symbol rate. Pulse shaping reduces the bandwidth further.

For an immediate insight into why the bandwidth can reduce to half the main lobe at most, consider a modulation case where the waveform alternates on every symbol between two symbols. The result would be a square wave at a rate of half the symbol rate. This would represent the highest frequency data pattern possible (such as for BPSK the data pattern of 101010 vs 11001100), and filtering down to half the main lobe would allow the fundamental of this square wave to still pass, as a sine wave yet still containing the information (101010) for subsequent demodulation. Filtering any lower than half the main lobe and this particular pattern (representing the highest frequency possible) would be completely removed.


Is time varying channels same as circulary symmetric complex gaussian?

From the literature i understood that if the channel is non time varying then we can consider it as circularly symmetric complex gaussian distributed.

That's plain wrong. Some non-time-variant channels are $\sim\mathcal{CN}$, some aren't.

  So can we consider the same if the channel is time varying ?

No.


Natural, Forced and Total System Response - Time domain Analysis of DT

First of all, many textbooks use the following definition (but see this answer for a different one):

zero-state response (ZSR) = forced response

zero-input response (ZIR) = natural response

The zero-state response is the response of a system to an input signal given that the initial state of the system is zero. The zero-input response is the output of the system due to a non-zero initial state. For a discrete-time system, the initial state is given by the values stored in memory (i.e. in the delay elements). Since we're dealing with linear systems, the superposition principle can be applied, and the total response of the system is given by ZSR+ZIR. So if the initial state is zero, the total response equals the ZSR, and if there's no input signal, the output is given by the ZIR, which always decays if the system is stable.

In analog circuit theory you have exactly the same. There the initial state is the energy stored in capacitors and inductors. Your car engine analogy is less fortunate because it doesn't really work. You could come up with a mechanical analogy though as long as the system is linear.

It might be a good idea to formulate a new question concerning recursive systems (after having searched the internet and this site). But you're right that a system is recursive if not only the input but also past output values determine the current output. The reason why this is done is the fact that in this way certain system properties can be implemented much more efficiently than by only using the current and past input values. The problem with recursive systems is that they can be unstable, and that quantization effects (in signals and filter coefficients) can more easily result in numerical problems than for non-recursive systems.

PS: If you teach yourself I strongly recommend MIT's open courseware. You should probably start with 6.003 Signals and Systems.


Creating a digital filter, from Laplace to $\mathcal Z$-transform (zero order hold) to code?

The example I looked at used a tustin or bilinear conversion not a zero order hold(the default for matlabs "c2d" command). So this is more an answer to what i wanted to do rather than the question that i asked above.

I solved the following (converting the s domain function into code) by taking the s domain function.
$$\frac{as+b}{cs^2+ds}$$

and putting this into matlab (matlab command "g=tf([a b],[c d 0])"). Then performed the bilenear conversion with the matlab command "c2d(g,Ts,'tustin')" where g was my transfer funtion and Ts my sampling rate. This produced the output

$$\frac{ez^2+fz+g}{iz^2+jz+k}$$

The a and b coeficients can then be taken from this equation such that(if $i!=1$ the equation needs to be multiplied through by the inverse of "i"):
$b0=e$ $b1=f$  $b2=g$
 $a0=i$  $a1=j$ $a2=k$

this can then be converted to code by setting the initial states for simplicity let $$z0=z1=z2=0$$

then set up a loop that repeats the following algorithm

$$output=z0*b0+z1*b1+z2*b2$$
$$z2=z1$$
$$z1=z0$$
$$z0=input-a1*z1-a2*z2$$

For anyone else that got lost like me, this is known as an IIR filter and googling IIR filter design helped sooo much.


Is a DC signal 0 frequency or extremely high?

In most mathematical treatments of DSP, a "DC" signal acts like a frequency component with a frequency and phase of zero (0 Hz, cos(0)).


struggling to understand why Fourier basis is orthogonal

Just use the formula for the geometric series (I use $l=h-k\neq mN$):

$$\sum_{n=0}^{N-1}e^{-j\frac{2\pi}{N}nl}=\frac{1-e^{-j\frac{2\pi}{N}Nl}}{1-e^{-j\frac{2\pi}{N}l}}=\frac{1-e^{-j2\pi l}}{1-e^{-j\frac{2\pi}{N}l}}=\frac{1-1}{1-e^{-j\frac{2\pi}{N}l}}=0,\quad l\neq mN$$


DSP Benchmarks for ARM Cortex-M4 processor

You're best off looking one level up, at datasheets or published benchmarks for actual devices that are available. One notable manufacturer of Cortex-M4 parts is Freescale. The Cortex-M4 is just a processor core design that is licensed by silicon manufacturers as the basis for their microprocessors. Each manufacturer designs their own peripherals and memory architecture and stitches them together with the core design. Therefore, what you might get from a core-level benchmark is the number of cycles required to execute your algorithm implemented in assembly.

What this doesn't capture are implementation-specific issues that can affect performance, such as memory latency, interrupt response latency, etc. These quantities will vary from manufacturer to manufacturer (and in some cases, across devices in a single family). Therefore, you really need to take the entire microprocessor into account if you're looking to do something that is pushing the envelope of what the chip might be capable of.

Even better: you can often find low-cost development kits for ARM processors (examples for Freescale's Cortex-M4 family). With that in hand, you can get a very good feel for how performant your algorithm will be once implemented in-system.


How can I reduce noise errors when detecting DTMF with the Goertzel algorithm?

One technique to reject noise is to use the ratio of the Goertzel filter output power to the total signal power over the same 40 mS sliding time windows, as well as ratio of the absolute Goertzel power to some threshold, as part of your decision criteria.

You can adjust your buffer size to match this minimum DTMF signal time (2048 may be too long if the sample rate is well below 44.1kHz).  Using a non-rectangular window may require using a longer Goertzel filter than 40 mS.  Adjusting the Goertzel filter length to be close to an integer number of periods of the filter's center frequency will also help reject DC.

Added:
Using a sliding Goertzel (and sliding RMS power) will reduce the effects of framing.


What are the differences between fft and dct?

First, you should compare the DCT to the DFT, not the FFT. FFT is just a family of algorithms to approximate the DFT.

Also:

  [...]
  In particular, a DCT is a Fourier-related transform similar to the discrete Fourier transform (DFT), but using only real numbers. DCTs are equivalent to DFTs of roughly twice the length, operating on real data with even symmetry (since the Fourier transform of a real and even function is real and even), where in some variants the input and/or output data are shifted by half a sample. There are eight standard DCT variants, of which four are common.

Source: https://en.wikipedia.org/wiki/Discrete_cosine_transform


What unit is the z axis of a power spectrogram of audio?

The short answer is that without a calibrated data acquisition chain from the original sound as it was picked up from at a microphone and converted to a digital file, it doesn’t have any units.   Its just numbers that fit into 16 bits (or 24 or 8) of dynamic range.

The number you have is probably proportional to something in nature if the original sound was natural but something coming out of Garage Band is synthetic.


Is there such thing as harmony independent timbre characteristic of a sound?

well, since timbre is a multi-dimensional perceptual phenomenon, there is no single metric.  so lotsa different metrics.

possibly the most common metric is obtained by high-pass filtering the note, measuring the average power, and normalizing that against the power of the unfiltered note and also normalizing against the pitch.  and it yields a "brightness" measure.  the HPF is a digital differentiator.

let your input note be $x[n]$.

differentiator output:
$$ x_1[n] = x[n] - x[n-1] $$

"sliding" average-power of HPF

$$ x_2[n] = \sum\limits_{m=0}^{M-1} (x_1[n-m])^2 \, w[m]$$

$x_2[n] \ge 0$ and $w[m]$ is a window function.

Normalized against the note power (so that louder notes having the same timbre do not get a larger measure of brightness:

$$ x_3[n] = \frac{x_2[n]}{\sum\limits_{m=0}^{M-1} (x[n-m])^2 \, w[m] + \epsilon} $$

$x_3[n] \ge 0$ and you can multiply $x[n]$ by nearly any value you want and get virtually the same $x_3[n]$.  $\epsilon>0$ is necessary to add to the denominator to preclude division by zero.  but $\epsilon \ll |x[n]|^2$ for most non-zero $x[n]$.

now this is normalized against different note loudness, but is not normalized against the pitch of the note.  higher pitched notes will result in higher values of $x_3[n]$, even if the waveform is otherwise identical.

to normalize against the pitch of the note, you will need a pitch detector or a fundamental frequency detector.  that's another whole problem and i have discussed it here on SE.


What is the smallest sampling points that can recover the signal?

This depends if you want to do sampling or just model identification. Model identification means "figuring out the numbers in your model from actually measured data". If the model is good and the data is not noisy you only need a few points. The minimum is equal to the number of unknown model parameters. Let's say you don't know v0, cl, epsilon and t0. In this case you need at least four points to identify those. Your model isn't linear so there is some math involved, however it's not overly complicated and very workable.
If the measurements are noisy it would be good to take more points and do a least square error fit. The number of points depends on the amount of noise and precision required for the results.

You can also sample the sequence but you need to consider sampling theorem. Your time sequence represents a first order lowpass filter. The cutoff frequency is determined only through your model parameter "cl" and is simply cl divided by 2 pi. At this frequency the spectrum is at -3 dB, at ten times this frequency it's at -20dB and at a 100 times this frequency it's -40 dB etc. Since it never gets all the way to zero you have to live with some amount of aliasing. Again it boils down to how much error you can/want to tolerate.


What does it mean to combine signals coherently?

It means that they are in phase---think of the signals as waves. Adding noncoherent signals results in cancellation, or fading.


Hopefully you'll get a bunch of answers here from the very general to the super specific. I'll put in my two cents here.
The recommendations I would make are from the field of radar and communication systems. These systems tend to exercise almost all aspects of signal processing:

Signal generation and mixing
Sampling, decimating/upsampling
Signal conditioning
General filter design and implementation
Detection theory
Target/observable tracking
Advanced techniques such as eigen-decomposition and super-resolution

Unfortunately books in the radar and communications fields are not free, but in my opinion are worth their asking price (as far as texts are concerned).
Bassem Mahafza: Radar System Analysis and Design using MATLAB, Chapman & Hall
This book is an awesome treatment of many of the aspects of a radar/comms system. As the title suggests, buying this book gives you access to all of the MATLAB functions provided by the author. As you work through the theory of the book, many sections have accompanying examples that you can play with allowing you to both see implementation details and visuals.
Merril I. Skolnik: Radar Handbook, Third Edition, McGraw-Hill
This book is a very large collection of different authors' work. It covers almost everything you can imagine about signal processing in a radar system. Theory, practicalities, and implementation specifics are all covered here. This book is to be used as a reference for applications and does not really contain examples of how to code certain things.
Harry L. Van Trees: Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory, Wiley & Sons
This one is dense and I would consider to be the more hardcore of the three. This text exposes to you a bunch of different techniques in signal processing when you have access to an antenna array, which would not be possible otherwise. Applications include radar, astronomy, communications, seismology, medical diagnosis, etc.
I could list many more but these are the three that come to mind that are a good mix of theory, applications of that theory, and examples.


It depends on the application and on what basis you decide to look at. If there are only a few targets and little to no clutter, then a radar image can be considered sparse in the image domain i.e. it has only a few Fourier components. A possible example of this is ground to air or air to air radars.

Synthetic Aperture Radar (SAR) imaging of the ground - imaging point-like targets on a specular reflecting surfaces (road, highways, airport runways) may also be considered sparse in the image domain. Once you have speckle though (e.g. high res imaging of fields, trees etc) then the image is not sparse in the target domain and sparse reconstruction will have a hard time with this type of data. The spare reconstruction may be able to separate point like targets from the clutter because it reconstructs the targets and not the speckle background.

Sparse reconstruction has also been applied to Inverse-SAR imaging of ships at sea. The ship tends to behave as a collection of corner or other simple reflectors - the reflection off the ocean waves tends to be away from the radar. You do run into cases where the ocean surface is rough and you get significant returns from it.

Yet another approach in SAR imaging to assume that the gradient of the image is sparse. This tends to behave like an edge detector.

So, it really depends on the application of the radar and on what basis you decide to use.

The CoSeRa workshop is held almost every year and focuses on compressive sensing. Many of the papers are focused on radar applications (SAR, moving target indication, tomography). The papers from the last few years are available through the IEEE.


Confusion regarding pdf of circularly symmetric complex gaussian rv

Let me try to establish the relation between the univariate PDF for a real Gaussian and the univariate PDF for a complex proper (i.e. circular symmetric) Gaussian.

You know that
$p_x(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp(\frac{-x^2}{2\sigma^2})$
is the PDF of a real-valued Gaussian with variance $\sigma^2$. We write $x\sim\mathcal{N}(0,\sigma^2)$ to denote that $x$ is a random variable that follows a real-valued Gaussian with zero mean and variance $\sigma^2$.

Now, let $z=a+jb$ be a circular symmetric random variable with real part $a$ and imaginary part $b$. In a circularly symmetric Gaussian random variable, the real and imaginary part are i.i.d., i.e. $a\sim\mathcal{N}(0,\sigma^2)$ and $b\sim\mathcal{N}(0,\sigma^2)$. Since $a,b$ are independent, their joint PDF is the product of their PDFs.

$$p_z(z=a+jb)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-a^2}{2\sigma^2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-b^2}{2\sigma^2}\right)$$

which gives the PDF of the complex variable $z$ at the value $a+jb$. We can reformulate this to
$$\begin{align}p_z(z=a+jb)&=\frac{1}{2\pi\sigma^2}\exp\left(\frac{-(a^2+b^2)}{2\sigma^2}\right)\\&=\frac{1}{2\pi\sigma^2}\exp\left(\frac{-z^*z}{2\sigma^2}\right)\\&=\frac{1}{2\pi\sigma^2}\exp\left(\frac{-\|z\|^2}{2\sigma^2}\right)\end{align}$$

We also write for this case that $z$ follows a circular Gaussian distribution and denote this by $z\sim\mathcal{CN}(0,2\sigma^2)$. Why is it suddenly $2\sigma^2$ (i.e. the double variance compared to the uni-variate, real case? Well, because $E[\|z\|^2]=2\sigma^2$ since it consists of a real and an imaginary part (which are independent and hence their variance adds up together).

So, to sum up:

if $x\sim\mathcal{N}(0,\sigma^2)$, then $p_x(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{x^2}{2\sigma^2}\right)$
if $z\sim\mathcal{CN}(0,\sigma^2)$, then $p_z(z)=\frac{1}{\pi\sigma^2}\exp\left(-\frac{z^*z}{\sigma^2}\right)$

Does this answer your confusions?


How to realize Poles and zeros at infinity??especially through transfer function?

It's actually quite straightforward: positive powers of $s$ (or, in discrete-time, $z$), correspond to poles at infinity. Negative powers give you zeros at infinity.

Let's look at some examples. In continuous time, an ideal differentiator has the transfer function

$$H(s)=s\tag{1}$$

Clearly, $\lim_{s\to\infty}H(s)=\infty$, hence you have a pole at infinity (and a zero at $s=0$).

The transfer function of an ideal integrator is

$$H(s)=\frac{1}{s}\tag{2}$$

Here we have $\lim_{s\to\infty}H(s)=0$, i.e., a zero at infinity (and a pole at $s=0$).

In discrete-time, a delay of one sample corresponds to

$$H(z)=z^{-1}\tag{3}$$

and you get $\lim_{z\to\infty}H(z)=0$.

A (non-causal) advance by one sample is represented by the transfer function

$$H(z)=z\tag{4}$$

In this case you have a pole at infinity: $\lim_{z\to\infty}H(z)=\infty$

Note that in continuous time as well as in discrete time, systems with poles at infinity cannot be causal and stable.


help to calculate phase shift in certain frequency(identified from FFT)

Phase shift in time domain translates to frequency shift in the Fourier transform.

$FT of exp(j2πf0t)x(t) = X(f−f0).$

If there is a phase shift in time, the Fourier transform that you take is already shifted in frequency. Find the shift in frequency f0 and you should get the phase change.


If you only have 15 hours I would built it on a well known
easy to understand algorithm like Viola and Jones, which uses weak binary classifiers and combines them to detect faces. It is relativily fast algorithm that runs in realtime.

You can get more information on the algorithm in the following stack exchange Post


Is there a process for deriving special cases of convolution?

Well, indeed there are special cases for convolutions but yours is quite straightforward. You have to consider two cases: $t<0$ and $t>0$.

If you sketch the convolution process for these two cases, by sliding your time shifted and time reversed input on the same axis as your impulse response, you will end up with the following:

Can you handle the integrals shown?


No, in a conventional sense of a "periodic signal" phrase, but, if you permit me to delve into a math subtlety, differentiation can turn an aperiodic waveform to a periodic one:
$$
\frac{\mathrm{d}}{\mathrm{d}t}(a\cdot t + b\cdot \cos(\omega\cdot t)) = a-b\cdot\omega\cdot\sin(\omega\cdot t)
$$
Notwithstanding a dubious usefulness of this excursion.


Number of samples in a discrete signals after sampling a continuous periodic signal should be integer?

There is nothing unordinary about it. Signals can have a period that is not an integer multiple of samples, and that is normal. For example, if you record a 440 Hz sine wave with a sampling rate of 48 kHz, it is not an integer multiple. It takes 11 full sine wave periods to hit an integer amount of 1200 samples.


Implementation of Fourier Domain Denoising with Hard Threshold

Thresholding in the Fourier domain is an archaic method sometimes called spectral subtraction, used for background noise removal in speech. Bad results in image processing can be due to several factors:

Misinterpretation: you say "The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients". Not quite. You keep "high amplitude coefficients", no matter low or high frequency.
Test data: the Shepp-Logan model is known to be quite specific. Bad results on it do not imply bad results on real data, yet:
Transform choice: Fourier is very good for close-to-stationnary data. Standard images with smooth parts inside contours are not well represented in the Fourier domain, at least not concentrated enough for acurate data compaction and reconstruction with few (high) coeffcients. Derivatives, or multiscale directional derivatives like wavelets can help better.
Implementation: Fourier implementation via FFT suffers boundary artifacts
coefficients, not talking about the 1D implementation you seem to use (prefer fft2 here)
Coefficient selection: hard and scalar thresholding typically yields more important visual artifacts than "smoother" thresholding and vector or group thresholding.
Formulation: you can generally get better results when you do not use a "parameter-tuned algorithm", but prefer a variational formulation based on your priors, including noise statistics.

So you have plenty of room to improve your results. Two references for multivariate image denoising with automated threshold (for Gaussian noise):

A nonlinear Stein based estimator for multichannel image denoising, 2008
Wavelet transform for the denoising of multivariate images, 2010


two dimensional integration of a trigonometric function

Since, the Taylor series expansion for $e^z$ about $z=0$ is
$$e^z = \sum_{n=0}^{\infty}\dfrac{z^n}{n!}$$
then, ignoring the question of convergence, you can say
$$\begin{align*} \int_a^b \int_c^d e^{A\sin(x)\cos(y-B)}dx dy &= \int_a^b \int_c^d  \sum_{n=0}^{\infty}\dfrac{(A\cos(y-B))^n}{n!}\sin^n x \; dx \, dy\\
\\
&= \int_a^b   \sum_{n=0}^{\infty}\dfrac{(A\cos(y-B))^n}{n!}\int_c^d\sin^n x \; dx \, dy\\
\\
&=   \sum_{n=0}^{\infty}\dfrac{A^n}{n!}\int_a^b \cos^n(y-B)\int_c^d\sin^n x \; dx \, dy\\
\\
&=   \sum_{n=0}^{\infty}\dfrac{A^n}{n!}\left(\int_{a-B}^{b-B} \cos^n u\;du\right)\left(\int_c^d\sin^n x \; dx \right)\\
\\
&=  \left(u\big|_{a-B}^{b-B}\right)\left(x\big|_c^d\right) + A\left(\sin u\big|_{a-B}^{b-B}\right)\left(-\cos x\big|_c^d\right) \\
&\quad + \sum_{n=2}^{\infty}\dfrac{A^n}{n!}\left(\int_{a-B}^{b-B} \cos^n u\;du\right)\left(\int_c^d\sin^n x \; dx \right)\\
\end{align*}$$
Using integration by parts, the reduction formulas for the integrals of $\sin^nx$ and $\cos^nu$ work out to be
$$\int \sin^nx \; dx = -\dfrac{1}{n}\sin^{n-1} x \cos x +\dfrac{n-1}{n}\int \sin^{n-2} x \;dx$$
and
$$\int \cos^nu \;du = \dfrac{1}{n}\cos^{n-1} u \sin u +\dfrac{n-1}{n}\int \cos^{n-2} u\; du$$
For $|A|\gg 0$, you could use the Taylor series expansion of $e^z$ about $z=A$ for a series that converges more rapidly for that value of $A$.


How is time or sampling rate accounted for in this model?

Units of s represent the Laplace transform while units of z represent the z transform. It is often much easier to transform a time domain (units of t) signal to Laplace given that it translates integro-differential equations into simple algebra. The z transform does the same thing, but takes advantage of the repetition from sampling to make the transform even easier.

"s" is any complex variable (with real and imaginary values), usually expressed $s=\sigma + j\omega$. When you limit s to only be imaginary ($s=j\omega$), the result is the Fourier Transform, which is a single-valued function versus frequency. Thus for a system, $H(s)$ would be the general transfer function of s and $H(j\omega)$ is the frequency response.

Similary "z" is any complex variable, but due to the mapping from analog to digital in the construction of the z -transform, what is the frequency axis in the s-plane (the $j\omega$ axis) becomes the unit circle in the z plane ($e^{j\omega}$), where $\omega$ here is typically normalized radian frequency (see What is normalized frequency)  with a unique range of $0$ to $2\pi$. Thus a digital system would have a general transfer function given as $H(z)$ while the frequency response would be $H(e^{j\omega}$). Conveniently, the z transform of one sample delay is simply $z^{-1}$, which hints at the great utility of using the z-transform for discrete time systems rather than the Laplace transform (which you could certainly still do instead if you want to take a mathematical beating!).

(Note: The mapping from $j\omega$ to the unit circle applies to many but not all mapping techniques: see z-Transform Methods: Definition vs. Integration Rule)

That said, to model the complete system digitally as a discrete system, you would translate the blocks that are functions of s to be functions of z. It is the process of translating, or mapping, from s to z where the sampling rate comes into play. There are several mapping techniques to choose from and a few popular ones are the matched-z transform, the method of impulse invariance and the bilinear transform. (For more on that see How/why are the $\mathcal Z$-transform and unit delays related? )   A simple example of this is the mapping of integration and differentiation operations from s to z using the matched-z transform (which by the way is identical to the method of impulse invariance for an all pole system):

Integration in time is equivalent to multiplying by $\frac{1}{s}$ in Laplace (there we see the conversion to simple algebra!). This is mapped to the following in the z domain:
$$\frac{T}{1-z^{-1}}$$

Differentiation in time is equivalent to multiplying by $s$ in Laplace (along with subtracting the initial condition which I didn't include).  This is mapped to the following in the z-domain:
$$\frac{1-z^{-1}}{T}$$

where $T$ represents the time of one sample period.

Notice the factor of $T$ which is the sample period. The differentiation example is very easy to see what is going on: $1-z^{-1}$ is a simple difference of two successive samples:

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

Which is a (not the only) discrete approximation of a continuous time domain derivative. Consider taking the derivative of any time domain function: as we sample at higher and higher rates, the difference between two successive samples will get smaller, thus by dividing by the sample duration for the case of the derivative which is also getting smaller, we normalize the result to be independent of the sampling rate.

In your case specifically, this would occur in the mapping from s to z for the continuous time transfer function given by $H_b(s)$.


For digital signals, you represent the $s(n) = \delta(n)$ dirac impulse by

$$s[n] = \delta_\text{discrete}[n] = \begin{cases}
0& \forall n \ne 0\text,\\
1 & n = 0
\end{cases}\tag{*}$$

which is obvious if you think about the dirac impulse of having the function (in continuous math) to give you

$$\int\limits_{-\infty}^\infty f(x)\delta_\text{continuous}(x) dx = f(x) $$

and in the discrete case, your integral collapses to the sum of

$$\sum\limits_{-\infty}^\infty f[x]\delta_\text{discrete}[x] dx \overset!= f[x] \tag{$\dagger$}$$

The only function that fulfills $(\dagger)$ is $\delta$ from $(\text *)$.


Given a linear phase FIR filter $h[n]$ of length $L$, if you use the conv or filter function to produce the output samples $y[n]$ for a given input $x[n]$, then there will be a constant group delay $N$ associated with the type of FIR filter being used. So this delay of $N$ samples must be taken into account when evaluating the output result in applications where data phase is critical such as in the case of using a Hilbert transform to produce an analytic signal...

Assuming that you have a filter with odd length $L$ which can be expressed as
$$ L = 2 N + 1$$ then the associated group delay will be the integer $N$; i.e., the information in the output data will be shifted $N$ sample to the right, hence you may select the output after the Nth sample as such:
$$ y_i[n] = y_d[n-N] $$.

In a computer environment such as MATLAB/Octave this translates into matrix indices such that:

xi = yd(N+1:end-N)

where $y_d$ was the matrix that holds the result of the raw convolution output;i.e.,

yd = conv(x,h)

if you have used filter function to produce the raw data than the slicing
should be done as:

xi = yd(N+1:end)

where $yd$ now is the matrix that holds the result of the raw filter output;i.e.,

yd = filter(h,1,x)

When the filter length $L$ is an even number, there happens a half-sample delay of information associated with the raw-data output of the computation and you should consider that as well. Therefore It's better to use an odd length FIR filter unless you have convincing reasons to do the otherwise.


If your window is a rectangular window, then there is no mapping out since splitting the stream into frames is effectivley applying that window. Once you split the stream into frames, if your window is not rectangular, you would then multiply the window against the frame.

y[i] = x[i] * w[i], for i = 0 to N-1, then take the FFT of y.


How to compute magnitude and phase response from transfer function in Z-domain?

Use euler identity to expand the exponential to cos and sin. Then rationalize the denominator, and the last step is to group it to real part and imaginary part. To compute the magnitude just use the pythagoras theorem.

For an example I'll use simple transfer function like: $$H(z) = \frac {z}{z-0.5}$$
Then substitue the $z$:$$H(e^{jw})=\frac{e^{jw}}{e^{jw}-0.5}$$
Expand to cos and sin:
$$H(e^{jw})=\frac{cos(w)+jsin(w)}{cos(w)+jsin(w)-0.5}$$
Rationalize the denominator, group into real and imaginary part:
$$\begin{eqnarray*}H(e^{jw})&=&\frac{cos(w)+jsin(w)}{cos(w)+jsin(w)-0.5}\cdot\frac{cos(w)-jsin(w)-0.5}{cos(w)-jsin(w)-0.5}\\
&=&\frac {e^{jw}\cdot (e^{-jw}-0.5)}{(e^{jw}-0.5)\cdot (e^{-jw}-0.5)}\\
&=&\frac {e^{jw-jw}-0.5e^{jw}}{e^{jw-jw}-0.5e^{jw}-0.5e^{-jw}+0.25}\\
&=&\frac {1-0.5e^{jw}}{1.25-0.5cos(w)-j0.5sin(w)-0.5cos(w)+j0.5sin(w)}\\
&=&\frac {1-0.5cos(w)-j0.5sin(w)}{1.25-cos(w)}\\
&=&\frac {1-0.5cos(w)}{1.25-cos(w)} - j\frac {0.5sin(w)}{1.25-cos(w)}\\
\end{eqnarray*}$$
Compute the magnitude:
$$ \begin{eqnarray*} \left|H(e^{jw})\right| &=& \sqrt{\Re(H(e^{jw}))^2+\Im(H(e^{jw}))^2}\\
&=& \frac {\sqrt{(1-0.5cos(w))^2+(0.5sin(w))^2}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1-cos(w)+0.25cos^2(w)+0.25sin^2(w)}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1-cos(w)+0.25(cos^2(w)+sin^2(w))}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1-cos(w)+0.25}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1.25-cos(w)}} {1.25-cos(w)}\\
\end{eqnarray*}$$
And you're done!


Look at the magnitude and angle of the signal relative to your carrier and you should see the same thing you would see if you looked at your angle at baseband using the I and Q samples.

It is difficult to see anything with the carrier there (and I would never simulate with an actual carrier when I can do everything at baseband, just in case you put it at the carrier for simulation purposes). To inspect this waveform, I suggest rotating it to baseband by multiplying it by $e^{-j2\pi f_c t}$, low pass filter the high frequency component assuming your input was real, and then plot the complex signal on an I Q complex plot. (Real on one axis and Imaginary on the other). You should expect to see the signal moving through your 16QAM constellation points with a trajectory that depended on the pulse shaping filter used. Alternatively you could also look at the real and imaginary separately in an eye diagram.

Here are some plots I have showing that may help illustrate the different cases:

On the left is an eye diagram and the associated constellation diagram for QPSK, with two cases for the IQ (real-imaginary) constellation; one where it was derived from an incorrect sampling point (timing error) and the other where each symbol was properly sampled, as indicated by the vertical bars in the eye diagram. Note the eye diagram is simply our waveform versus time over one symbol period, plotted again and again (as in infinite hold on a scope) such that we see every possible trajectory from one symbol to the next, which is dependent on the history of all prior symbols over the memory of our channel.

On the right is the constellation diagram for 16QAM, where the large one is what it would look like with a carrier offset (similar to your case), such that it is essentially our 16 points spinning (which doesn't tell us much!), and the other is after removing the spin with $exp^{+j2\pi f_c t}$ or $exp^{-j2\pi f_c t}$ depending on which way the diagram is rotating. In your case with a real signal both spins are present, so multiply by either and low pass filter the resulting higher spin portion. One direction will have an inverted spectrum compared to other, so if you go the wrong way you may have an inverted result, but this is easy to undo by simply swapping I and Q).  The smaller plot on the right is the resulting 16QAM constellation plotted as real vs imaginary of ONLY the correct sampling locations within each symbol were chosen (no timing error as in the QPSK case shown), but we do see evidence of phase error and some amplitude error (hence the spreading of the dots), as well as the possibility of other distortions such as multipath etc causing intersymbol interference; all having the effect of increasing the size of the dots in our constellation (leading to errors!).

Below is another example for QPSK showing the constellation and eye diagram side by side, and also showing the significance of sampling at the correct sample positions (in red) versus the overall trajectory of the waveform.

And finally here is an example eye diagram for 16 QAM:


Understanding FFT size in OFDM modulation

An FFT does not need to be of length 2^N.  However, FFTs of lengths that are a multiple of only tiny primes run faster.  If you don't mind slower performance (or higher CPU utilization), you can use any length.

Note that multiple modulations can be combined, for instance, modulation of each subcarrier, and modulation of a collection of carriers.  Each can be different, thus differently named.  For instance, each subcarrier in an OFDM system can have a different size of QAM constellation.


How to Use the DFT (FFT) to Solve a Least Squares Regularization Problem (Inverse Problem)?

The question really depends on $ f \left( \cdot \right) $.
Yet in order to show how to use FFT we can even use 1D signals.

Let's rewrite the problem:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} \left\| K x - b \right\|_{2}^{2} + \frac{\lambda}{2} \left\| f \left( x \right) \right\|_{2}^{2} $$

The derivative is given by:

$$ g = {K}^{T} \left( K x - b \right) + \lambda f' \left( x \right)^{T} f \left( x \right) $$

Now, $ K $ as a matrix is given by a Circulant Convolution Matrix.
Hence the operation $ K v $ or $ {K}^{T} v $ can be done in the Fourier Domain.
Under some circumstances it might accelerate the operation significantly.


Think about how the output signal $y[n]$ is computed:

$$\begin{align}y[0]&=x[0]h[0]\\
y[1]&=x[1]h[0]+x[0]h[1]\\
y[2]&=x[2]h[0]+x[1]h[1]+x[0]h[2]\\\vdots\end{align}$$

From the first equation above you know that the first coefficients of the FIR filter must equal $h[0]$. The second equation gives you $h[1]$ as the second coefficient, etc. So your intuition was correct: the FIR filter's impulse response is just a truncated version of the original impulse response.

For finding the difference equation, you don't need to compute any zeros or poles, you just need the coefficients:

$$y[n]=x[n]h[0]+x[n-1]h[1]+\ldots +x[n-1000]h[1000]=\sum_{k=0}^{1000}h[k]x[n-k]$$


How is the frequency vector being derived?

linspace(0, 1, fix(numel(FTSignal)/2)+1) $\tag{1}$

"$y =$ linspace(x1,x2,n) generates n points. The spacing between the points is (x2-x1)/(n-1)."

So in the above statement it is generating fix(numel(FTSignal)/2) points between 0 and 1 that are equally spaced

  Frv = linspace(0, 1, fix(numel(FTSignal)/2)+1)*Fn;

This generates Frv with just equation (1) scaled/multiplied by Fn

Edited: the FFT signal can usually be represented by half of the total FFT coefficients, since FFT of a real signal is symmetric around mid point, it is uniquely represented by half the coefficients. Adding 1 is for the DC frequency. Take an example an FFT of length 16, is represented by 8+1 coefficients. The "fix" function is present to handle FFT of odd length signal, so that they can be rounded to an integer value. Since these FFT values have been generated under the normalized scale of 0 to 1, it is then normalized by Fn, which can be $\pi$ or $F_s/2$, where $F_s$ is sampling frequency


What impact have different types of border-handling on high pass or low pass filtering?

HPF or LPF linear operators. It means that an operation at pixel position i, j will only be affected by neighbors within the window size of your filter. Pixels outside the filter will not play a role - in other words, the effects of any filtering operation is localized around its support. Typically, if the image size post filtering is not a concern, we 0 pad an image and crop out the contents near the border with a size equal to half the filter window. To maintain the same size, it is advisable to replicate border values.


You'd generally do the FFT without zero padding, i.e. window size == FFT length.

Having an FFT larger than the window leads to an interpolation, but doesn't add any information to the signal, and that's only useful in certain cases – maybe the paper presents one such case.


Can a software defined radio with a frequency range of 300 MHz – 3.8 GHz detect AM (535 to 1705 KHz) and FM signals (88 to 108 MHz)?

For those given specs, without additional hardware (upconverters) or modifications to the SDR circuit, it's not possible to detect AM or FM (88-108 MHz) range of frequencies. As you might have also seen, there are frequency range extenders (especially on the low frequency range) available as kits.

Also, the sampling rate determines the maximum bandwidth of simultaneous detection (or spectrum analysis if you prefer to use the SDR for that purpose). So, the bigger the sampling rate the larger will be the spectral bandwidth that you can capture instantaneously. This means, more radio channels are observed at once and more carriers are searched without changing the central tuning frequency of the SDR tuner chip. Which makes it faster in spectrum scanning applications.

However, using a large bandwidth (high sampling rate) mode of SDR, might bring the problem of transmitting the data from SDR to PC, especially on a USB 2.0 interface. (this is not a problem on most ubiquitous SDR receivers used at 2.4 MSPS rates) Also more processing and interrupting than necessary (provided a narrow bandwidth radio channel is to be demodulated) is done at those high sampling rates, to decode just a small portion of the analysis bandwidth.


Does delta-sigma ADC also reduce Gaussian noise on input signal to ADC or just quantization noise?

No!  The ADC (delta sigma or not) can not reduce the uncertainty in the input.  It sounds to me that your friend has not made up a real signal flow diagram and then formed equations.  The answer to your second question is affirmative.  In fact it's generically  true that putting as much gain, without introducing extra, in the front end is a good idea.  But... think long and hard about input overloads distorting the signal and possibly overloading the amplifier circuit causing distortion and such; either in the front end amplifier or during subsequent signal processing.  In the "good"  old days overloading due to the input could cause numerous problems; including long term effects.  OTOH: detailed knowledge of the signal and noise spectral (or even time transient) characteristics can reveal ways to optimize your S/N .


Memory efficient filtering with scipy.signal in Python

The way to do this is to break the signal up into chunks and process each one at a time.
I asked in a comment if you could accept one-dimensional filtering, but I guess you can do bidirectional filtering, too; you just have to do it yourself.
For one-directional filtering:

Break up the signal into chunks that can fit in memory  (Is there a convenience function for this?  Maybe np.nditer?)
zi is initially all zeros
For each chunk:

Run it through scipy.signal.sosfilt, with zi=zi
Set zi = zf (from the output of the filter)
Write the chunk of output filtered data to your new array

When you've assembled all the output chunks, it will be as if you had done one continuous sosfilt in one direction.
For bi-directional filtering:

Do the whole thing again, but in the reverse direction.

For this to be exactly equal to sosfiltfilt, you probably need to keep some additional tail past the end of your signal, and then use that as your initial conditions when going backwards?  Look at the code for sosfiltfilt or filtfilt and confirm how they work.
If you are ok with FIR filtering, then that's easier: You can make a zero-phase filter and just do that in one direction (using lfilter instead of sosfilt), with the same process as above, and then throw away the head of the output, so that it looks as if it were done with no offset.


What is the current "state of the art" in future audio sample prediction, given recent developments in "image inpainting" using deep learning?

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


What happens when linear phase is critical?

Causal IIR (infinite impulse response) filters cannot have a linear phase response, so they can only be used if (exactly) linear phase is not required. Unlike IIR filters, FIR (finite impulse response) filters can have a linear phase response (and, consequently, a constant group delay). So in applications where a linear phase is required, we have to choose FIR filters.

A non-linear phase response causes signal dispersion because different frequency components experience a different delay when being filtered. E.g., in many audio applications a linear phase is not required (although in some it might be), whereas in digital communications, a non-linear phase would cause neighboring symbols to interfere with each other (inter-symbol interference, ISI), so a linear phase is important.

Note that it is possible to design IIR filters with an approximately linear phase in the pass band (as shown in this answer). For many practical applications this might be good enough, and it might save a lot of computational complexity and memory compared to a linear phase FIR filter with a comparable magnitude response.

As mentioned in MBaz's comment, more details on applications where a linear phase may or may not be required can be found in the answers to this question .


Let's say you have a signal which is all zeros except for a spike at one point where x(8)=1 (total N=32, for example).  If you perform the DWT on this signal and then calculate the total energy (by taking the square root of the sum of the squares of all the results), you will get a value - call it "E1".

Now, let's take another signal which is still all zeros except at one point, n=12, where x(12)=1 - JUST LIKE THE FIRST SIGNAL.  It is the same signal except it is shifted by 4 points.  If you take the DWT of this signal and calculate its total energy (again by taking the square root of the sum of the squares of the results), you get another value - call it "E2".

You might expect that the total energy of these two would be the same - but they're not!  E1 will NOT equal E2, so therefore the DWT is NOT "shift-invariant" because the energy "varies" whenever you "shift" the incoming signal - even though it's basically the same signal.

So transforms like the dual-tree complex discrete wavelet transform have the property that they are "shift-invariant", which means that you can shift the signal around before calculating the transform and the resulting energy will still be the same.  This is important for all sorts of reasons, but if you think about it, shift-invariance is a really good thing to have because it means that it doesn't matter WHERE in the signal you start calculating - it will still have the same energy coming out.

This page has a pretty good description and demonstrations non-shift-invariant transforms in 1D, 2D, and 3D.

Here's a really simplistic example that (I hope) illustrates the main idea of shift-invariance:

Suppose you are trying to find the wavelet transform of a piece of music you have, but the WAV file you used has 5 seconds of silence before the music starts.  You get the wavelet transform coefficients and graph them nicely on a scalogram (a time/scale graph - which you can think of as a time/frequency graph).  You notice that it indicates 5 seconds of silence (0 values) before the musci starts.

Now you edit the WAV file to cut out the 5 seconds of silence down to 1 second of silence (you are shifting the sample in time by doing this - hence "shift" in shift-invariance).  You take the discrete wavelet transform again and plot it in a scalogram and you expect to see the same thing except without the silence - you expect it to just be shifted over to the left, but...

It won't be the same!  In fact, if you take the total energy of the resultant signal, it won't even be the same, because the DWT by itself is NOT shift-invariant.

But if you use the UDWT (also called "stationary" DWT), then the two images will be the same except the second one will not have all that silence in the beginning.  You could take a screenshot of the first one and crop it in photoshop to get the same image as in the second one.

The dual-tree complex DWT (DTCDWT) is near-shift-invariant, meaning it's very close to being shift-invariant.

The term literally means "does not vary when shifted".

Hope that helps!


what is $1/f$ property of the visual world?

"Essentially, all models are wrong, but some are useful" (attributed to G. Box)

One of the simplest morphological model for 2D images decomposes objects into structure and texture. Both can be used to detect or understand objects. Even though image formation models are still complicated, objects often appear as smooth intensity gradient, filled with decoration, bounded by edges (some being false).
Smooth intensity gradients are often the most energetic, thus low frequencies can be prominent. Edges, though energetic, are sparser and local, thus possibly less present in a spectrum.

So, a simple "low frequencies are prevalent than high frequencies" model is a decreasing spectrum. Now what would be the decreasing simplest model in terms of elementary functions? Possibly $1/f$. And yes, $\exp^{- |f`|}$  could seem better, especially with the singularity at $0$. More involved models involve $1/f^p$ or $\exp^{- p |f|}$ shapes. Indeed, such models are dependent of the transformations used to reshape the data, and Fourier is a very simple reshaping tool. More local tools like wavelets use Generalized Gaussian Distributions: $\exp^{-\alpha x^\beta}$ in the time/space domain.

Simpler can also be seen from a computational point of view, see Frequency Synthesis of Landscapes (and clouds):

  Frequency synthesis is based upon the observation that many "natural"
  forms and signals have a 1/fp frequency spectra, that is, their
  spectra falls off as the inverse of some power of the frequency where
  the power is related to the fractal dimension.

  This leads naturally to a method of generating such fractals:

  Generate a random (white noise) signal.
  Transform this into the frequency domain.
  Scale the resulting spectra by the desired $1 / f^p$ function.
  Inverse transform.

Such simple models can reproduce realistic shapes. From the frequency domain:

one gets:


What is the difference between Linear Interpolation factor and Sampling rate conversion factor?

I'll move my comments to an answer:

With sample rate conversion factor $S$, you want to compute samples at "time" $\frac{i}{S}$ for $i=0,1,2,\dots$. "Time" relating to the index of the previous series. The linear interpolation factor is simply then the fractional part of $\frac{i}{S}$ for each $i$ with the integer part telling which sample index to use (together the next) for the linear interpolation.

Simple example to clarify: upsample a sequence a factor $\frac{10}{9}$. You would like to determine the samples at "time" $0, 0.9, 1.8, 2.7, 3.6, \dots$.  First, start with sample $0$, then, compute the linear interpolation between sample $0$ and $1$ with $\eta = 0.9$. Next sample should be at "time" $1.8$, so $\eta = 0.8$ for samples $1$ and $2$. "Time" $2.7$ gives $\eta = 0.7$ for samples $2$ and $3$ and so on.

Note that you conceptually can compute the "time" $1.8$ sample from sample $0$ and $1$ with $\eta = 1.8$, but using sample $1$ and $2$ with $\eta = 0.8$ is better. So while $\eta$ does not really have to be $0$ and $1$ the results will be better and there is no advantage to not have it in that range.


Your 2-mic array can provide an angle-of-arrival estimate with front-back ambiguity, using cross-correlation to estimate TDOA and then calculating AOA from TDOA and microphone spacing. The front-back ambiguity can be removed with sufficiently directional microphones, or just by constraining the geometry, putting the microphones against walls, for example.
You will also have a multi-path effect. If the space is cluttered, if you lose the direct path, or if the sound source is directional, this multi-path effect will be strong and the errors it produces may be overwhelming.
But even in ideal conditions, the two mics produce a single geo-observable, the TDOA, and thus can calculate only one geometric constraint. In other words, you have only one "known" and thus can solve only one "unknown", not a full 2D position.
If the sound source and the sampling at the microphones can be synchronized, you can also measure propagation delay directly to get distance estimates from each mic - two "knowns". Such absolute distance estimates are also very powerful for resolving the multi-path problem.
If the source and receivers cannot be synchronized, you will need a third microphone to produce a second independent TDOA measurement, a second "known". At this point, you can use the same pseudorange/light-cone cost-minimization calculations that are used in GPS receivers, or just take the two TDOAs and try to intersect the AOA lines.
These are the general principles.
(My comments are based on personal experience with geolocation of radio emitters in several success projects, but the same principles apply to sound.)


What is a speaker behaviour for out of range frequencies?

Imagine an amplifier/speaker combination with no DC blocking, so the speaker cone moves exactly in proportion to the input voltage.
Now imagine a microphone that measures the actual absolute pressure of the atmosphere, again with no DC blocking.
Even with such a speaker and microphone, because of the finite size of the speaker, the air will act as a highpass filter.
It's complicated, and I'd have to do some research myself to be able to put numbers to it, but roughly speaking the actual response of the sound pressure to the speaker position is going to start rolling off at about the point where the wavelength is the same order of magnitude as the speaker's diameter.  There will be much better information out there on the web, but I suspect you can search it as well as I can.
This is why woofers are bigger than tweeters, and its why subwoofers are constructed so (seemingly) oddly.  Subwoofers, in particular, are significantly smaller than a wavelength of the sound they're producing, so the designers need to play all sorts of acoustic tricks to make them work, and they're horribly inefficient, which is why they often need -- or come with -- their own amplifiers.
When the speaker pushes out suddenly in response to the square wave, it tries to generate a permanent rise in pressure -- but for wavelengths longer than the box, the pressure wave more or less goes around the box and into the hole in the back, and very little sound is produced.  What's left is that spike you see.


How to shift input for testing shift invariance in a system?

Time invariance: changing $x$ by $x-x_0$ do not affect anything.

Example 1: "Multiplying the argument by 2" operator
$$H(f(x))=f(2x)$$

Note the operator is very tricky: multiplying the argument by 2.

For example:
$$H(sin(x))=sin(2x), H(x^2+1)=4x^2+1, H(e^x)=e^{2x}$$

The invariance is tested by comparing a shift in the signal vs a shift in the result of the operator:
$$H(f)|_{x-x_0}=H(f|_{x-x_0})$$

For the LHS:
$$H(f)|_{x-x_0}=f(2x)(x-x_0)=f(2(x-x_0))$$
The function is $f(2x)$, which is shifted by $x_0$.

For the RHS
$$H(f|_{x-x_0})=f(2x-x_0)$$
The function is $f(x-x_0)$, which is converted to $f(2x-x_0)$.


Pulse Doppler Radar: how to filter its Range-Doppler matrix?

Two suggestions that may independently, or combined, help to solve your problem.

Consider using a different window/weighting function for your FFT processing.  There are window functions with low uniform sidelobes such as the Chebyshev window. This window function allows for selection of a constant sideline level response.
It sounds like you are using a fixed detection threshold, consider using an adaptive detection threshold so that the threshold can vary with the combination of background noise and sidelobe artifacts.  Mathworks has a pretty good description for a simple case at Constant False Alarm Rate (CFAR) Detection.


How to constrain the object detection problem?

A long time ago I have come across a very good advice: research should be problem-oriented, not approach oriented. That means find a problem you want to solve, and then look for appropriate methods.  Not the other way around. So, I would recommend starting with a specific problem, such as a real-world application.

For example, maybe you want to build a system with a camera looking at your bookshelf, which detects a particular book.  Or maybe you want to build a mobile robot with a camera, which can follow you around.  Or a video player program, which detects whether you are looking at your computer screen, and pauses the video when you look away (some cell phones do that now).  Or a program that can organize your photo collection by detecting and recognizing people, and maybe landmarks, like the Eiffel Tower.  The possibilities are endless.

So the algorithm is as follows:

Find a real-world vision problem that you find interesting.
Make as many restrictive assumptions as you can to scope it down and still keep it interesting.
Start reading papers to find applicable algorithms.
Understand the limitations of existing algorithms, and come up with your own ideas to address them.


Approximate the change in perceived gain from applying transfer function to audio signal

You want to calculate loudness difference between input and output.

Loudness depends on the power spectral density of a signal; knowing only the transfer function of your filter is not enough for calculating PSD.

For example, suppose your filter is a lowpass filter, with attenuation 20 dB in the stop band. If your input signal is completely contained below the cutoff frequency, then no perceived gain (input and output are similar), whereas if your input signal is all above the cutoff frequency, then you will perceive a -20dB gain.

Related question on calculating loudness.


Design of Reed-Solomon Code where $n < q-1$

You have the right idea but are missing some of the fine details.  A cyclic Reed-Solomon code of length $n$, where $n$ is a proper divisor of $q-1$, is a BCH code and the key property is that if a sequence of $2t$ consecutive $\mathbf{n}$-th roots of unity are roots of the generator polynomial, then the minimum distance is at least $2t+1$. Here, $n = 3$ and $\alpha^5$, $\alpha^{10}$ and $\alpha^{15}=1$ are the three cube roots of unity. So, choose two of them to be roots of $g(x)$ (it doesn't matter which ones because you will always get two consecutive roots in this simple case!) and you are done. If you choose $\beta = \alpha^5$ and $\beta^2 = \alpha^{10}$ as the roots of $g(x)$, then $g(x) = x^2+x+1$, and so all the codewords on this $[3,1,d=3]$ Reed-Solomon code are of the form $\gamma +\gamma x + \gamma x^2, \gamma \in$ GF$(2^4)$, or equivalently $[\gamma, \gamma, \gamma]$,   that is, we have a repetition code!

It is also the case here that the cube roots of unity are elements of the subfield GF$(2^2)$ but this is not the case in general. Take $n=5$ which is also a divisor of $q-1=15$ in which case $\alpha^3, \alpha^6, \alpha^9, \alpha^{12}$ and $\alpha^{15}=1$ are the five $5$-th roots of unity which all belong to GF$(2^4)$.


What is the difference between convolution and cross-correlation?

The only difference between cross-correlation and convolution is a time reversal on one of the inputs. Discrete convolution and cross-correlation are defined as follows (for real signals; I neglected the conjugates needed when the signals are complex):

$$
x[n] * h[n] = \sum_{k=0}^{\infty}h[k] x[n-k]
$$

$$
corr(x[n],h[n]) = \sum_{k=0}^{\infty}h[k] x[n+k]
$$

This implies that you can use fast convolution algorithms like overlap-save to implement cross-correlation efficiently; just time reverse one of the input signals first. Autocorrelation is identical to the above, except $h[n] = x[n]$, so you can view it as related to convolution in the same way.

Edit: Since someone else just asked a duplicate question, I've been inspired to add one more piece of information: if you implement correlation in the frequency domain using a fast convolution algorithm like overlap-save, you can avoid the hassle of time-reversing one of the signals first by instead conjugating one of the signals in the frequency domain. It can be shown that conjugation in the frequency domain is equivalent to reversal in the time domain.


Jigsaw puzzle: isolating the pieces (separating "stuck" blobs)

I agree with Abid Rahman K, however eroding only would result in shape shrinkage.  To only get rid of the gap, Abid Rahman suggests a watershed based approach, but for a faster and easier method you might try:

Erode the regions
Apply connected component labeling.
Dilate the regions back.

This way you will have almost the same shape (with some artifacts arising due to morphology) but you will surely end up with single connected components.

Watersheds and distance transform are also powerful and might be good choices if you have a lot of overlap.


Note that the antiderivative of a function is only defined up to a constant. Furthermore, note that if you integrate a periodic function, the result is not necessarily periodic. Let

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0t}\tag{1}$$

Now we integrate $(1)$ with an arbitrary lower integration limit $t_0$. I'll say more about that later.

$$\begin{align}y(t)&=\int_{t_0}^tx(\tau)d\tau\\&=\sum_{k=-\infty}^{\infty}a_k\int_{t_0}^te^{jk\omega_0\tau}d\tau\\&=a_0(t-t_0)+\sum_{k=-\infty\\k\neq 0}^{\infty}\frac{a_k}{jk\omega_0}e^{jk\omega_0t}-\sum_{k=-\infty\\k\neq 0}^{\infty}\frac{a_k}{jk\omega_0}e^{jk\omega_0t_0}\tag{2}\end{align}$$

Clearly, if $a_0\neq 0$, $y(t)$ isn't periodic, so we can't hope to obtain a formula for its Fourier coefficients. So in order to end up with a periodic function $y(t)$, we require that $a_0=0$.

If $a_0=0$, we get from $(2)$

$$y(t)=\sum_{k=-\infty}^{\infty}b_ke^{jk\omega_0t}\tag{3}$$

with

$$b_k=\begin{cases}\displaystyle\frac{a_k}{jk\omega_0},&k\neq 0\\\displaystyle-\sum_{l=-\infty\\l\neq 0}^{\infty}\frac{a_l}{jl\omega_0}e^{jl\omega_0t_0},&k=0\end{cases}\tag{4}$$

So in general the antiderivative $y(t)$ has a non-zero DC component, which depends on the choice of the lower integration limit $t_0$.

The antiderivative

$$\tilde{y}(t)=\sum_{k=-\infty\\k\neq 0}^{\infty}\frac{a_k}{jk\omega_0}e^{jk\omega_0t}\tag{5}$$

is the specific antiderivative of $x(t)$ that has a zero DC Fourier coefficient, and this is the one that is meant by the formula that is usually given in textbooks.

Consequently, the expression

$$y(t)=\int_{-\infty}^tx(\tau)d\tau\tag{6}$$

for the antiderivative of $x(t)$ with Fourier coefficients $b_k=a_k/(jk\omega_0)$, $k\neq 0$, and $b_0=0$,
is at least inaccurate; it is in fact a sloppy way of expressing the specific antiderivative $\tilde{y}(t)$ given by $(5)$, i.e., the one with a zero DC coefficient.


How to name every row of 100*2500 matrix in MATLAB?

You can assign each row to a vector as follows:

row1 = matrix(1,:);

That will store the whole first row of matrix in row1.

To do this "automatically", you can store all the vectors containing each row in a cell array, for example:

for k = 1:size(matrix,1) %For each row
    rows{k} = matrix(k,:); %Store the k-th row in the k-th element of the cell array
end


Practicality of the assumption of i.i.d. for Rayleigh channels

I am a little late but I post my answer anyway so that someone having the same question will find it interesting and discuss.

The discrete baseband multipath channel can be modeled as a FIR, i.e.
$$y[n] = \sum_{l=0}^{L-1} x[n-l] h_l + w[n]$$
where $L$ is the number of channel taps. $L$ depends on the relation between the bandwidth of basis waveform and the delay spread of channel.

The term "Rayleigh fading" channel implies that the channel taps $h_l$ can be modeled as i.i.d circular symmetric Gaussian complex random variables because :

$h_l$ is the sum of a large number of small independent circular symmetric random variables, each random variable is channel of a physical path. This the rich scattering assumption, which is typically vaid for urban environment.
there is not a particular path having gain much significant than others. Otherwise, we have Rician fading.

Let me call this random variables "Rayleigh".

With a sufficient cyclic prefix ("sufficient" means greater than delay spread, thus OFDM receiver captures all delayed versions of OFDM symbol, proof can be found at single-tap OFDM regardless subcarrier spacing), the demodulated data at subcarrier $k$ is
$$z[k] = x[k] \times \sum_{l=0}^{L-1}h_l e^{-j2\pi \frac{l}{N} k} = x[k] \times H[k]$$
where $N$ is DFT size.

The channel taps $h_l$ are i.i.d circular symmetric Gaussian random variables, $H[k]$ are circular symmetric Gaussian random variables, but they are in general not iid.

As pointed out by Maximilian Matthé in the comment, the covariance matrix is   $F \mathrm{diag}(\vec{p}_0) F^H$ where $\vec{p}_0$ is Power Delay Profile zero-padded to size $N$. The frequency bins at $k = u \times N/L, u \in \mathbb{N}$ are independent, if $N/L$ is integer. Other ones are sinc-interpolated thus they are correlated. Note that $N/L \times \Delta f = 1 / L T_s \approx 1/\tau_m$ can be seen as coherence bandwidth.