Using MATLAB `xcorr()` to Find a Signal Hidden Behind Noise

You're basically doing a bank of hypothesis to find your signal using Matched Filter.

Though you use a slightly different method.

First of all, you should leave the signal in the time domain and calculate the cross correlation or their multiplication at the frequency domain.
Yet, since your signal doesn't have unknown phase (Or delay) multiplication will do (It covers the cross correlation at time zero) -> which is equivalent of what you do at the frequency domain.

Since you're using AWGN the performance depends on the Effective SNR.
The effective SNR is a function of the energy of the signal, variance of the noise and the number of samples.

If you set the signal and the variance of the noise, just use more samples (Of the Signal -> Equivalent of faster sampling rate).


Bilinear Transform (Tustin's Method) applied to the Derivative

It's a matter of perspective if the bilinear transform "fails so miserably" when trying to approximate a derivative. First of all, it's not true that the approximation is quadratic for small $\omega$, it's linear as it should be:

$$D(e^{j\omega})=\frac{2}{T}\frac{1-e^{-j\omega}}{1+e^{-j\omega}}=\frac{2j}{T}\tan\left(\frac{\omega}{2}\right)\approx \frac{j\omega}{T}\quad\textrm{for}\quad |\omega|\ll 1$$

So for small frequencies the discrete-time system behaves as a (discrete-time) differentiator. For larger frequencies we expect the effect of frequency warping, which is exactly what happens here: the infinite continuous-time frequency range is mapped to the interval $\omega\in[0,\pi]$.

In sum, the bilinear transform behaves as expected, but it is also clear that there are better ways to approximate a derivative in discrete-time than applying the bilinear transform to the continuous-time transfer function $H(s)=s$.


Calculation of the salience function (Klapuri 2006)

You are correct that STFT is a function of 2 variables (usually x-axis is time and y-axis is frequency). The author is a bit careless in calling $Y(k)$ the STFT. What he really means is that $Y(k)$ is a fixed column (i.e. fixed time slice) of the STFT. The salience function is more carefully defined by the same author in this paper titled Multipitch analysis of polyphonic music and speech
signals using an auditory mode Eq. 13, where he uses a subscript $t$ to denote a fixed time slice of the STFT.


Why would someone need a Spatial Light Modulator with a maximum phase shift $> 2\pi$

That will be useful at least for pulse shaping, because then it is not just a single wavelength that you are dealing with, and as modulo $2\pi$ is actually implemented as modulo of displacement, it would work only for a single wavelength, creating jumps in the phase for nearby frequencies.


Finding Frequency Differences In Differing Amplitude Signals

Sounds to me like you need to look at the signals in the frequency domain.  Take the Fourier transform of each signal (use an FFT), and look at the spectrum.  I'm not sure if you've got signals with just a single frequency (a tone) or with many frequencies, but in either case if they have the same frequencies in them then the spectra should look the same (except for one being larger than the other).

If you're unfamiliar with frequency domain concepts, I'd suggest you read up on them first.  The wikipedia page on Fourier transforms is pretty large, so I'm not sure how good that would be as an introduction, but there are a large number of books on the topic and many websites too.


Using fft in matlab to retrieve input signal magnitude

I just fixed 3 things in your code and added comments where I did. The reason is basicly just the way the FFT is defined/implemented in Matlab.

fs         = 1024;
fftsize    = 2048;           % Choose as a power of 2 of the sample rate

sim_length = fftsize/fs ;
t          = 0:1/fs:sim_length-1/fs;
x          = 0.5*sin(2*pi*15*t);

data                  = x';
window                = blackman(fftsize);
windowed_data         = window.*data;

fftd_data   = abs( fft(windowed_data, fftsize)/length(x) ); % Normalize the fft

bin_freq = fs/(fftsize) ;

fundamental_bin = ceil( 15/bin_freq );

lower_limit = fundamental_bin-4;
upper_limit = fundamental_bin+4;

2*sum( fftd_data(lower_limit:upper_limit) ) % Multiply by 2, because of the DFT-symmetry

ans =

    0.5000

For a more detailed explanation, check out this thread at Mathworks forum.


Optimal windowing strategy for filtering live signal

If possible, design your filter, and then do what must be done (buffering etc.) to implement that filter properly. When implementing a finite impulse response (FIR) filter such as: (pseudocode:)

out[t] = 0.5*in[t] - 0.6*in[t-1] + 0.7*in[t-2] - 0.8*in[t-3]

where t is an integer time index, you could do it so that you store a few of the past in values in a buffer, like this, stepping t:

t buf[0] buf[1] buf[2]
----------------------
0 0      0      0
1 0      0      in[0]
2 0      in[0]  in[1]
3 in[0]  in[1]  in[2]
4 in[1]  in[2]  in[3]
5 in[2]  in[3]  in[4]
...

That way you always have in[t-1], in[t-2] and in[t-3] available for calculation of out[t]. Perhaps you were worried about the performance cost of moving data in a "rolling" buffer like this. But that extra work can be avoided by using a circular buffer, the length of which is usually rounded up to the nearest 2^n (integer n) so that wrap-around can be done using binary and (2^n-1). The next write position is stored here in an extra state variable i:

t i buf[0] buf[1] buf[2] buf[3]
-------------------------------
0 0 0      0      0      0
1 1 in[0]  0      0      0
2 2 in[0]  in[1]  0      0
3 3 in[0]  in[1]  in[2]  0
4 0 in[0]  in[1]  in[2]  in[3]
5 1 in[4]  in[1]  in[2]  in[3]
6 2 in[4]  in[5]  in[2]  in[3]
7 3 in[4]  in[5]  in[6]  in[3]
8 0 in[4]  in[5]  in[6]  in[7]
...

If you are using single-instruction-multiple-data (SIMD) acceleration, it can be more straightforward to also write the current in[t] into the buffer before calculation of out[t].

If you have an infinite impulse response (IIR) filter (like Butterworth) you can use another circular buffer for the past out values. This would give an IIR implementation structure called direct form I. There are other structures that may or may not have better numerical stability, see http://www.earlevel.com/main/2003/02/28/biquads/.

If you don't care about the phase of the filter's frequency response, a simple two-pole, two-zero notch IIR filter will do the job of eliminating 50 Hz: (pseudocode:)

params:
fs = sampling frequency
f = notch frequency
r = sharpness, 0..1 excluding 1

init:
z1x = cos(2*pi*f/fs)
a0a2 = (1.0 - r)*(1.0 - r)/(2.0*(fabs(z1x) + 1.0)) + r
a1 = -2.0*z1x*a0a2
b1 = 2.0*z1x*r
b2 = -r*r

loop:
out[t] = a0a2*(in[t] + in[t-2]) + a1*in[t-1]
         + b1*out[t-1] + b2*out[t-2]

For such a simple filter you may want to skip implementing an elaborate indexed buffering scheme.


how to extract a radio signal features with python?

Some of the more cutting edge signal classification techniques utilize deep neural networks (DNNs):

https://arxiv.org/pdf/1712.04578.pdf

This particular paper utilizes a convolutional neural network (CNN) and a residual neural network.  The convolutional "front-end" extracts its own nonlinear features prior to connection to the densely connected layers.

DNNs are typically very computationally expensive, so your approach to feature extraction approach might meet your requirements.


"Infinity focus" refers to the Hyperfocal Distance of the camera's lens. This is the distance beyond which all objects will appear acceptably focused.

Therefore, if your objects come closer (to the camera) than its hyperfocal distance they will appear blured and detection of sharp features such as corners will become more challenging.


Which of DFT, DCT and DWT transforms is more robust to noise and geometric transformations in image processing

From the perspective of numerical analysis, the most robust transformations are unitary orthogonal matrices because they have a condition number of one.  Indeed the condition number is the metric that defines robustness.  You can find this in Stuart’s Introduction to Matrix Computations or the more advanced book by Golub and Van Loan.

A related way of looking at perturbation sensitivity are the Holder p norms for Matrices also covered in Stewart


Regularized Least Squares by Laplacian Operator - Optimal Value of the Regularization Factor (Lagrangian Multiplier)

Let me rewrite the problem:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| L x \right\|}_{2}^{2} $$

In order ot say something intelligent about the optimal value of $ \lambda $ one must set an optimization criteria.
Since you don't have one there is no well defined answer.

What might assist you is reformulating the problem in an equivalent form:

$$\begin{align*}
\arg \min_{x} \quad & \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} \\
\text{subject to} \quad & {\left\| L x \right\|}_{2}^{2} \leq \alpha &
\end{align*}$$

Both forms are equivalent in the sense that for any $ \alpha $ of the 2nd form there exists $ \lambda $ such that both solutions are the same.

Now if you can say something about the smoothness level of the image which will let you set $ \alpha $ intelligently than you'll be able to extract $ \lambda $ to match it.

Pay attention that the connection isn't trivial and depends on the data itself among other things (See my answer to Significance of λ in Basis Pursuit).


Derivations here - you can pick any tmin, tmax, and fmin, fmax for some number of samples or sampling rate N.
We adjust the code one line toward the end to rescale phi to end as an integer multiple of $2\pi$ to yield zero phase; this has the effect of nudging fmin and fmax, slightly or greatly depending on all other parameters - see here.
An alternative variant that forces only the end of the chirp to be zero phase will exactly preserve fmin and fmax, by subtracting. Forcing all, zero phase at tmin and tmax without changing fmin and fmax is impossible.
Taking your parameters, with cosine and sine:

Code
import numpy as np
import matplotlib.pyplot as plt

def _lchirp(N, tmin=0, tmax=1, fmin=0, fmax=None):
    fmax = fmax if fmax is not None else N / 2
    t = np.linspace(tmin, tmax, N, endpoint=True)

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

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

def lchirp(N, tmin=0, tmax=1, fmin=0, fmax=None, zero_phase_tmin=True, cos=True):
    phi = _lchirp(N, tmin, tmax, fmin, fmax)
    if zero_phase_tmin:
        phi *= ( (phi[-1] - phi[-1] % (2*np.pi)) / phi[-1] )
    else:
        phi -= (phi[-1] % (2*np.pi))
    fn = np.cos if cos else np.sin
    return fn(phi)

#%%######################################################################
f0 = 7000
f1 = 17000
samplerate = 192000
T = .0013

N = int(samplerate * T)
tmin = 0
tmax = T

t = np.linspace(tmin, tmax, N, endpoint=True)
for zero_phase_min in (True, False):
    for cos in (True, False):
        x = lchirp(N=int(samplerate * T), tmin=tmin, tmax=tmax, fmin=f0, fmax=f1,
                   zero_phase_tmin=zero_phase_min, cos=cos)
        plt.plot(t, x)
        plt.title("cos={}, zero_phase_tmin={}".format(cos, zero_phase_min),
                  weight='bold', fontsize=17, loc='left')
        plt.show()


...how can I extrapolate the pen's possible path across the semi-enclosed area?

There are infinitely many ways that this path could have been generated. Therefore, you will either have to make assumptions or use extraneous information for something that resembles "...the pen's possible path..." (which would not be too far away from the way potrace does it).

In any case, instead of applying edge detection, I would suggest that you apply skeletonisation first to get rid of the stroke's width and then try to work out the "evolution" of a much simpler path. The stroke's width generates what looks like a closed path which is incredibly hard to work out because of ambiguities (do I trace this clockwise or counter clockwise?)

Even when the skeletonised output is given, you still need extraneous information.

For example:

Where did the pen land on the paper for the first time?

You could select the "left-most" pixel for that but cursive writing can result in heavily slanted letters where this assumption is false. Checkout the "y" on the word "yesterday" in this image for example (image taken from the Wikipedia article linked earlier).

"Which process generated this curve?"

This is basically what you are looking for and possibly cursive writing could be expressed as some form of minimisation problem. That is: Find the shortest length curve that links as many of the given points as it is possible. Notice here, "curve", singular, one curve.
An algorithm that does this is Alpha Star. Alpha star operates over a graph data structure and finds a path that is optimal given a heuristic. Usually, this heuristic is an expression of "cost". Say for instance, the "cost" of establishing an edge between two points is proportional to their distance (definitely) but also how close that edge would bring us to the final destination. The algorithm basically calculates this for all transitions and then works out the path with the highest combined heuristics.
In your case, the graph data structure are all the points that result from skeletonisation connected at least with the very simple rule of shortest distance. So, you start at your "starting point", scan its neighbourhood for the nearest ink point, connect it, jump to the nearest and repeat. More than one pixels could be "nearest" in this case, which is fine at this stage, you connect them and then let A* tell you how to resolve a "fork".
In your case, the cost is obvious: The distance between pixels. But the heuristic is possibly where the game is won or lost. You can use the "how close it brings us to the end of the curve" but this means that you have an idea of where the curve ends.... Which brings us to even more assumptions:

"Where did the pen lift from the paper?"

You can possibly identify far enough sequential points and mark them as the "end points". This would give you rudimentary start-end points for tracing the curve but where does this leave the dot above the "i"? Was it added while writing the "i" or at the end of a given "curve"?

For these reasons, I do not think that you would be able to recover the pen's movement without any extraneous data.

It would help even if you could make an assumption on the pen nib that the writing was made with because if it was a flat nib and it was held at an angle, it can give you information about its "orientation" and it is safe to assume that the writer was not rotating the pen in their hands. This gives you an additional criterion to that of the nearest distance. Nearest distance and smooth angle transition.

A more accurate method would be to record cursive writing movements for each letter and then recognise each (visually) and "fit" the movements to the point data with "transition curves" between the letters. This would give you the possible way a given writer could have written a particular script. You can capture this data with lots of different devices (tablets, e-notebooks).

Hope this helps.


I just want to comment on the confusion between $p(t)$ and $h(t)$. What happens is that we transmit a signal $s(t)$ where

$$s(t)=\sum_ks_kg(t-kT_s)$$

where $s_k$ are the data symbols, $g(t)$ is the transmit pulse shape, and $T_s$ is the symbol time.

Transmitting this over a wireless channel with a channel impulse response $h(t)$, the channel acts acts as a filter, and thus the channel output is the convolution between the transmitted signal and the channel, i.e.,

$$y(t) = s(t)\star h(t) = \sum_ks_k\underbrace{\left[g(t-kT_s)\star h(t)\right]}_{p(t-kT_s)}$$

So, the convolution of the transmit pulse shape and the channel is considered in the matched filter design.

By the way, the purpose of the matched filter is not to "remove the noise", but rather to maximize the signal-to-noise ratio (SNR) by maximizing the received signal power at the output of the matched filter.


Filtering Overlapping Frequencies in sound file

In my opinion, the best approach for this kind of spectrally-overlapping noise is to use blind source separation techniques, like independent component analysis (ICA) or time-frequency masking, or sparse source decomposition.  Look up the work of Emmanuel Vincent, Kostas Kokkinakis, or Philipos Loizou.  These methods work best when some statistics of the noise are time-independent, and the SNR is small or negative, as appears to be the case in your example.  Vincent has some MATLAB code available on the web.


Are the samples adhering to Nyquist (properly pre filtered)?
If so, I would just resample to one common rate, prefereable >= max(rate1, rate2). And align them in time. Check out MATLABs resample() function.
If the samples are not adhering to Nyquist, then you need to think about what they represent, how they are sampled and in what respect you need to compare them.


Generating audio clip by superposing two files

We typically model audio to be a pretty linear thing, so yes, just add – you can multiply each input with a constant to weigh them differently.


Doing DFT based interpolation has to keep 3 principles:

Keep the zeros outside the data (Looking form $ -\pi $ to $ \pi $). Namely add zeros to the higher frequencies.
If it is real signal, keep the symmetry (If needed, divide the highest frequency into 2 items).
Preserve the energy (By the division by 2 if needed).

By the way, the interpolation isn't by a Sinc Kernel, it is done with its discrete form - The Dirichlet Kernel.
See The Proper Way to Do Sinc Downsampling (DFT Downsampling) for Uniformly Sampled Discrete Signals with Finite Number of Samples


Use MatConvNet. They have samples for image-net classification, such as this. Pre-trained models already exist. You could also train on your own.

Check cnn_imagenet function in MatConvNet. You will find it useful.
There are tutorials. For example this one.


Intensity change after filtering this image?

For binary images, the filter kernel you have is looking for a $0 \rightarrow 1$ edge in the North-East direction.  Your image raster only has one place where such an edge exists (in the bottom left / South-West corner).

The sign of the output means that the edge is a $1 \rightarrow 0$ edge transition rather than a $0 \rightarrow 1$ transition.


How to derive equation for second order interpolation of soft knee cutoff in a compressor?

If you look at the formula and plot it, you get what you see in the other link on dsp.ee. But if you add the inflexion points and plot the consituent functions, it should become more clear:

The function would have been a linear one, with two slopes, but between the points defined by $T-W/2$ and $T+W/2$, there is a curve uniting the two segments. This means that the conditions become:

a simple, quadratic formula (to avoid too much computation)
the function's value at the first point should be equal to the point, itself, and
the derivatives at those points should be the actual slopes

With these, there are many ways to interpolate, so I can't say exactly how they did, but these must have been their starting points.
For example, considering two slopes of $3/4$ and $1/4$ at $[-1,1]$, the quadratic function $f(x)=ax^2+bx+c$ that would unite those would have to have:
$$\left\{
\begin{aligned}
f'(-1)&=2a(-1)+b&=\dfrac34 \\
f'(1)&=2a(1)+b&=\dfrac14
\end{aligned}\tag{1}
\right.$$
The derivative makes $c$ disappear, but its value at $-1$ is $3/4$, so the solving the equations above results in the function:
$$f(x)=\dfrac34+\dfrac{x}{2}-\dfrac{x^2}{8}\tag{2}$$
and this is how it looks like:

The slopes at $[-1,1]$ are a match, which means the two slopes are now connected by a quadratic formula. The whole graph would be, from left to right: the red slope until $-1$, then the quadratic blue until $1$ and then the green slope; just like in the 1st picture.

I'll add that you can go the other way, just as well: from the derivative, up.
You know that the signal has a slope in the beginning, and one at the end which means that their derivatives give you fixed values (in the example above they were $3/4$ and $1/4$). Uniting the points will be a slope of the form $ax+b$, as seen in (1), starting from $-1$ with the value $3/4$, and ending up at $1$ with $1/4$. Since it's a line, it will pass through the middle point, $m$, and it will have a slope $s$, both given by:
$$\begin{align}
m=\dfrac{\dfrac14+\dfrac34}{2}&=\dfrac12 \\
s=\dfrac{\dfrac14-\dfrac34}{2}&=-\dfrac14 \\
\Rightarrow\quad g(x)=m+sx&=\dfrac12-\dfrac14x\tag{3} \\
\end{align}$$

Integrate this to get the quadratic, to which you add the constant term $3/4$ -- the value of the function at $-1$:
$$\int_x{g(x)\mathrm{d}x}=\color{blue}{\dfrac34}+\dfrac{x}{2}-\dfrac{x^2}{8}\tag{4}$$
(2) and (4) are the same. q.e.d.


What are the parameters of transducer that I can get from an Anechoic chamber?

Mainly frequency response as a function of the 3D direction of incidence. Also known as "polar pattern" or "polar balloon". You can also look at distortion and noise but that's more tricky.


inverse z transform performed on 6th order IIR filter

You can write your transfer function as the sum of three 2-nd order transfer functions (partial fractions), and then the impulse response will be the sum of the impulse response of each.
To help you with the partial fractions you can use this calculator


What is the form of the spectral derivative in the all-positive-frequency notation in DFT?

There is no requirement, it is a matter of choice.

Some conceptual clarity is required.  A sequence of sampled points does not have a derivative.  Derivatives come from continuous functions.

If you treat the inverse DFT definition as a continuous function, the result is a continuous interpolation function that will pass through all your sampled points.  Choosing different ranges (could be 6N to 7N-1 if you wanted) yields different interpolation functions each being differentiable.

This thread should give you some insight:

How to get Fourier coefficients to draw any shape using DFT?

The "simplest" or "most natural" range is using the negative frequencies, and if you have an even number of bins, the Nyquist bin should be split in half, one positive and one negative.


This is achievable with two parallel all pass filters.
The two all pass filters synthesize an odd ordered low pass filter whose pass band extends from -90º to +90º in the z-domain. (I will discuss this below).
$$
G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}
$$
The low pass filter is then rotated by +90º so that its pass band extends from 0º to 180º, which approximates the Hilbert transform. Rotation mathematically is:
$$
H_{Hilbert}(z) =G_{lowpass}(-jz)
$$
As a consequence, one of the all pass filters becomes completely imaginary - this is the Hilbert transformed signal path. By necessity, since the filter is approximating this Hilbert transform, then the imaginary branch must be applying a 90º phase shift relative to the output of $A_0(z)$.
$$
H_{Hillbert}(z)=G_{lowpass}(-jz)=\frac{A_0(z)+jz^{-1}A_1(z)}{2}
$$

Example with Elliptical Filter
The following example was taken from Handbook for Digital Signal Processing pg. 920
You can also use this program that I interpreted to design a new half-band elliptic filter with a higher order.
Starting with a real, half-band elliptical low pass filter $G(z)$ whose frequency response satisfies
$$
1-\delta_1\leq|G(e^{j\omega})|\leq1\>\>\>\>\>for\>\>0\leq\omega\leq\omega_p
$$
$$
\>\>\>|G(e^{j\omega})|\leq\delta_2\>\>\>\>\>for\>\>\omega_s\leq\omega\leq\pi
$$
To obtain half band we choose:
$$
\omega_s+\omega_p=\pi
$$
$$
(1-\delta_1)^{2}+\delta_2^{2}=1
$$
But I just used the examples given seventh order poles or an elliptical low pass half-band filter:
$$
z=0,\>\>\>z=\pm j0.436688,\>\>\>z=\pm j0.743707,\>\>\>z=\pm j0.927758
$$
$$
G_{lowpass}(z) = \frac{A_0(z)+z^{-1}A_1(z)}{2}
$$
Distribute the poles between the two all pass filters like so:

$$
A_0(z)=\frac{z^{-2}+0.190696}{1+0.190696z^{-2}}\cdot\frac{z^{-2}+0.860735}{1+0.860735z^{-2}}
$$
$$
A_1(z)=\frac{z^{-2}+0.553100}{1+0.553100^{-2}}
$$
Rotate by 90º
$$
H_{Hilbert}(z)=G_{lowpass}(-jz)=\frac{A_0(-jz)+jz^{-1}A_1(-jz)}{2}
$$
$$
A_0(-jz)=\frac{0.190696-z^{-2}}{1-0.190696z^{-2}}\cdot\frac{0.860735-z^{-2}}{1-0.860735z^{-2}}
$$
$$
A_1(z)=\frac{0.553100-z^{-2}}{1-0.553100^{-2}}
$$
Phases of $A_0(z)$ and $A_1(z)$:

Magnitude Response of $H_{Hilbert}(z)$:

Phase Difference Between $A_0(z)$ and $A_1(z)$ Branches:


If you can assure that the physical signal being measured by A has no detectable components above 50Hz, then you could take 1 every 10 samples.

However, if you have CPU power available, it is much better to filter the signal (lowpass, with stop-band starting below 50Hz) before downsampling. This way, components above 50Hz won't be aliased down to your working bandwidth. Also, this process will reduce quantization noise from the ADC, which is a good thing.


Processing signals with phantom fundamentals

Methods for estimating a fundamental pitch frequency (missing or not) from an FFT result containing a series of harmonics include finding maxima in a cepstrum, the harmonic product  spectrum method, and looking for clusters in histograms of the distances between frequency pairs of higher order frequency peaks that fit within psycho-acoustic filter bandwidths..  These are all more estimation algorithms that post-process an FFT result, rather than a simple formula.


It is not possible to simply concatenate FFTs to obtain a larger one (that would be a neat trick if it worked). It is possible to combine several smaller point FFT's to create a larger one. To see this, let's start with the original equation,
$$ X[k] = \sum_{n=0}^{N-1} x[n] \exp\left(-\frac{j 2 \pi k n}{N} \right). $$
Now, consider that we have used this formula on two sequences already to obtain two FFTs of length $\frac{N}{2}$,
$$ X_0[k] = \sum_{n=0}^{N/2-1} x[2 n] \exp\left(-\frac{j 2 \pi k n}{N/2} \right) $$
and
$$ X_1[k] = \sum_{n=0}^{N/2-1} x[2 n + 1] \exp\left(-\frac{j 2 \pi k n}{{N/2}} \right). $$
Now, we have two FFTs of length $\frac{N}{2}$, and we want to take the information contained in them to form a larger FFT. Specifically, we would like to form
$$ Y[k] = \sum_{n=0}^{N - 1} x[n] \exp\left(-\frac{j 2 \pi k n}{N} \right). $$
To see how to obtain this, we decompose the summation into two parts, one containing even terms and one containing odd terms,
$$ Y[k] = \sum_{n=0}^{N/2 - 1} x[2 n] \exp\left(-\frac{j 2 \pi k \left( 2 n\right)}{N} \right) + \sum_{n=0}^{N/2 - 1} x[2 n + 1] \exp\left(-\frac{j 2 \pi k \left(2 n + 1\right)}{N} \right). $$
Some simple transformations turn this into,
$$ Y[k] = \sum_{n=0}^{N/2 - 1} x[2 n] \exp\left(-\frac{j 2 \pi k n}{N/2} \right) + \exp\left(-\frac{j 2 \pi k}{N}\right) \sum_{n=0}^{N/2 - 1} x[2 n + 1] \exp\left(-\frac{j 2 \pi k n}{N/2} \right). $$
And this yields our answer. To obtain an FFT that is twice the length of two smaller FFTs, we simply apply a twiddle factor to the later samples and combine the two smaller FFTs via addition. Explicitly, based on the notations above, we can write
$$ Y[k] = X_0[k] + \exp\left(-\frac{j 2 \pi k}{N}\right) X_1[k]$$
where the index $k$ is taken modulo-$\frac{N}{2}$ or else $X_0[k]$ and $X_1[k]$ are considered periodic functions with period $\frac{N}{2}$.

You can arrive at a similar formula through a decimation in frequency approach. I leave the details for that to you or someone else.

Edit: If you are dead set on using $N=5120$ as your output FFT size, then you can use the following variation. First compute 5 FFTs of length 1024. You can compute these by taking every fifth number from the original sequence and applying the FFT from the library. For notational convenience, let's define these as
$$ X'_p[k] = \sum_{n=0}^{N/5-1} x[5n+p] \exp\left(-\frac{j 2 \pi k n}{N/5}\right) \:\: \mbox{for $p = 0, 1, 2, 3, 4$}.$$
Now, combine these components with the equation
$$ X[k] = \sum_{p=0}^{4} \exp\left( -\frac{j 2 \pi p k}{N} \right) X'_p[k].$$
Again, the indices should be considered modulo-$\frac{N}{5}$ in this case.


HINT:

If $z_1[n]$ is the response to input $x_1[n]$ (and $y[n]$ is NOT an input signal), and $z_2[n]$ is the response to $x_2[n]$, then it is straightforward to show that the response to a linear combination $ax_1[n]+bx_2[n]$ is NOT equal to $az_1[n]+bz_2[n]$, which would be a requirement for linearity. And due to symmetry, the same is the case if $y[n]$ is the input (and $x[n]$ is NOT an input signal). The test for time-invariance also fails in both of these cases, because only the input signal is shifted (and not the other nuisance signal).

When both $x[n]$ and $y[n]$ are input signals, the linearity test above still fails, so the system is non-linear (because it is an affine mapping due to the constants in the system equation). It is, however, time-invariant because if both input signals are shifted by some constant, the output is shifted by the same constant (which is not the case if only one of the two signals is an input signal).


How to decide what subcarrier allocation should I use in OFDM?

As is typical, the answer is that it depends. In most cases, you're implementing a system that someone else has designed already, and you just need to match their design. However, if you're in the situation of designing a system of your own, then there are a few high-level considerations:

Data subcarriers are the ones that actually communicate information through the channel. Ultimately, the goal of any communications system is to transfer information with as little energy and time expended as possible, so maximizing the number of data subcarriers used in each symbol helps to reach the second of the two goals.
Null subcarriers don't do anything for information transmission; instead, they are used to shape the spectrum of the transmitted signal. Therefore, we would want to minimize the number of null subcarriers that are in use, as they reduce the information rate of the system. They are often used for two reasons:

Allowing some room for analog anti-aliasing filters. In a practical system, you'll need to leave some portion of the upper end of the frequency band near the Nyquist rate empty to allow for rolloff of pre-ADC and/or post-DAC filtering. Typically, you'll see the upper ~10-25% of the band unused.
Exclusion of the center bin in the band. Common, low-cost direct conversion receivers usually exhibit a high level of interference at band center (due to ADC DC offset after the conversion step). Therefore, it's common in wireless standards that are aimed at low-cost markets to just place a null subcarrier in that slot to avoid the potential interference.

Pilot subcarriers are used in the equalization process. Again, since their contents are known to the receiver a priori (or else they wouldn't be useful as pilots), they don't convey any information through the channel. So, one would want to minimize the number of pilot subcarriers in use. However, you want to have enough of them so that you can effectively estimate the response of the channel across the entire frequency band. This is why you'll often see them spaced across the OFDM symbol. Also in many cases, the pilot positions will rotate over time, allowing the receiver's equalizer to periodically get a view of the entire channel response.

Note that the number of subcarriers doesn't define the channel bandwidth. OFDM transmitters are often implemented using an IFFT structure so you can lay out the subcarriers in the frequency domain first, then transform to the time domain to get a waveform that you can apply to the channel. The sample rate that is used on that waveform, along with the subcarrier allocation and FFT size, is what defines the bandwidth of the transmitted signal. You can increase the FFT size arbitrarily while still maintaining the same signal bandwidth (you will just end up with more lower-rate subcarriers that are spaced more finely in the frequency domain).


Shouldn't the Sampling Theorem imply that there should be no information loss at all after a signal is processed?

For no information to be lost on conversion back to continuous form, the signal would first need to be perfectly band-limited, and you would need an ideal reconstruction filter.  A perfectly band-limited signal is infinite in extent. Since you want this arbitrary signal to be processed by a computer, your computer would need infinite memory.  You would also have to wait an infinite amount of time for your perfect reconstruction filter to settle.

Finite sampling time and filter length problems also completely ignores other potential lossy issues, such as quantization and sampling clock jitter, etc.

One might be able to process an infinite signal that can be represented in certain closed symbolic forms within finite time and resources, but there does not seem to be a method to convert any arbitrary signal into such a form.

Back in the real-world (stuff you can buy or make), one would normally accept an information loss from imperfect band-limiting, finite filter length, jitter (etc.) that is around or below the quantization and numeric noise floors.  This allows the processing to happen using an amount of RAM that one can afford to hopefully finish within one's lifetime.  Thus leading to information loss and imperfect reconstruction.


How to get the background noise of microphone?

For a given TV, if you send a single short pulse (with an arbitrary amplitude of '1') out the speaker. The microphone will receive a slightly distorted pulse a short time later. (Typically you would send a pulse with enough volume that the microphone gets a clear reception of it, then reduce the microphone signal proportionally to simulate an original speaker signal with amplitude of '1'.) That distorted and delayed microphone signal is called the "impulse response" of that particular TV system. That impulse response (including the delay from the instant the signal was sent to the speaker) will need to be stored in your software for that TV.

Then to remove the distorted speaker signal from your microphone signal, you need to take the original electronic speaker signal and multiply it by that impulse response, and then subtract the result from your microphone signal. You do all this in the frequency domain because its far simpler and more efficient than doing it in the time domain. You don't even have to turn the final microphone signal back into time domain since you only wanted the ambient volume, anyway.

Look up "convolve impulse response".

EDIT:
To clarify (per your comment): The impulse response of the speaker-to-microphone path is a measure of the delay and distortion created by that path--and that's what you initially experimentally measure in the lab. But you don't know the impulse response of the noise-source-to-microphone path, so you can't correct for that distortion. All you can do is note the power spectrum of the noise (after subtracting the speaker contribution), and then adjust the speaker volume and maybe spectrum in order to help compensate for that noise. By the way, be sure to do the arithmetic (the multiply and the subtract) in the complex plane, because the FFT data is complex.


Symmetry of real and imaginary parts in FFT

Yes, this is always true if the input to the DFT is real valued. It's called the "conjugate complex symmetry", because
$$
X_{N-n} = {X_n}^*
$$
where $X_n$ is the DFT output and $()^*$ denotes the conjugate. It can be proven by inserting the property into the transformation formula of time domain sequence $x_k$:
$$
X_n = \sum_{k=0}^{N-1}x_ke^{-j\frac{2\pi k n}{N}}\\
X_{N-n} = \sum_{k=0}^{N-1}x_ke^{-j\frac{2\pi k (N-n)}{N}}\\
=\sum_{k=0}^{N-1}x_k e^{-j 2\pi k}e^{j\frac{2\pi k n}{N}}
$$
Using $exp(-j2\pi k) = 1 \:\: \forall \: k$ we find
$$
X_{N-n} = \sum_{k=0}^{N-1}x_k e^{j\frac{2\pi k n}{N}}
$$
Now we exploit that $x_k$ is real for all $k$ to rewrite the above:
$$
X_{N-n} = \sum_{k=0}^{N-1}\left(x_k e^{-j\frac{2\pi k n}{N}}\right)^* = X_n^*
$$
It's now straightforward to show that the same is valid for the IDFT. (real spectrum -> complex conjugate symmetric time domain sequence). Moreover, the opposite statement is also true: if the input to the DFT (IDFT) is complex conjugate symmetric its output is real valued.


Got stack in calculating state-space representation

Model your deterministic input as a $2 \times 1$ matrix :

$$ \begin{bmatrix}
 u  \\
T_a \\
\end{bmatrix}
$$

Then your state space equations will be:

$$ \begin{bmatrix} \dot x_1 \\ \dot x_2 \\ \dot x_3 \end{bmatrix} =  \begin{bmatrix}
 -\frac 13p & \frac 13p & 0 \\
p & -(p+ \frac{c}{10}) & \frac{c}{10} \\
0 & \frac{c}{30} & -(\frac{c}{30}+\frac a3)
\end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} + \begin{bmatrix} 0 & 0 \\ v & 0 \\ 0 & a/3 \end{bmatrix} \begin{bmatrix} u \\ T_a\end{bmatrix} $$


Two bugs in your code.

pol2cart() returns the real and imaginary part, but you only use the real part and discard the imaginary part
You apply ifftshift() on one method but not the other

Try instead

[a,b] = pol2cart(artPhase, artMag);
artInv = ifft(a+1i*b); % Time-domain signal
plot(artInv-invSig);


What does pre-emphasis $\alpha=1$ will do to a signal?

The pre-emphasis filter is a simple high-pass filter with a monotonically increasing magnitude of its frequency response $H(\omega)$. Its value at DC is $H(0)=1-\alpha$, and its value at Nyquist (i.e., the maximum frequency) is $H(\pi)=1+\alpha$. Consequently, for $\alpha=1$ you get a zero at DC, which means that the filter will totally suppress a constant input signal (which is obvious from looking at the difference equation).

Note that a pre-emphasis with $\alpha=1$ cannot be compensated for because the DC component is lost and can't be recovered. In signal processing terms: there is no stable inverse filter.


Trouble with inverse Z-transform and calculating of samples

The poles in the equation are Z=1.5,0.8;hence to be two sided the roc is $0.8<=z<=1.5$
so the 0.8 term is right sided and 1.5 term is left sided let us do partial fractions now $ X(Z)=1+\frac{1.5-0.85z^{-1}} {(1-0.8z^{-1})(1-1.5z^{-1})} $
after doing some maths i got $x(n)=\delta(n)-0.5(0.8)^nu(n)-2(1.5)^nu(-n-1)$
so x(-2)=-0.88888,x(-1)=-1.3333,x(0)=0.5,x(1)=-0.4,x(2)=-0.32.


How to find the period of a noisy signal using MATLAB's $\tt xcorr$?

You are getting that result because you are taking into account other maxima that do not correspond to the points you are actually looking for.

When you call findpeaks() with no other parameter than the input signal, you are taking into account all the points marked with a blue arrow. To determine the period, you just want those inside red circles.

Note that, due to the noise, inside some red circles there are more than just one point considered as a peak. That is undesired as well.

You have to put a restriction to find peaks such that they are at least 160 samples away from each other (in this case, due to your signal and its correlation). You can see that if you replace

[~,locs]=findpeaks(ac);

with

[pks,locs] = findpeaks(ac,'MinPeakDistance',300); %Or any number >160

the period will be the expected one.


How to simulate a signal of a particular frequency according to the number of points in a period of the signal?

I suppose that the $n$ points (hereafter, the chunk) precisely divide the base period (which you don't know yet)  evenly. You are lucky. You can verify that by concatenating  the $n$  discrete samples  (or possibly $n-1$ consecutive ones)
by an integer number $K$ of (periodic) repetitions (copy-paste, several times, hereafter the chain), with $Kn\ge m$ (or $K(n-1)\ge m$. It should be smooth  visually. The FFT of the chunk and the chain should be superposable (with appropriate scaling in time and amplitude).
Then, crop the chain of $Kn$ samples to $m$.
Then, using a non-parametric (discrete Fourier transform, FFT), you can estimate the spectrum of the (underlying) continuous data, and find the minimum Nyquist sampling frequency, which may differ from the original one (the initial sampling).


Optimizing FFT code for 12 bit processing on microcontroller

12-bit is not a native data type on your platform, so the best you can do without loss of precision is 16-bit per sample, in Q15 format.

I don't know the Stellaris line well, but from what I gathered on the TI website, these are either Cortex M3 (LM3xxx) or Cortex M4 (LM4xxx).

The Cortex M4 embeds a floating point unit; so there won't be as much benefit in terms of run-time from converting your code to fixed point as you think. The immediate thing is to switch to 32-bit floats. The next thing is to check your compiler settings and have a quick look at the disassembled generated code to check that it actually generates floating point instructions rather than using emulation. With gcc, you'll probably have to use compilation flags like -mcpu=cortex-m4 -mfloat-abi=hard -mfpu=fpv4-sp-d16 (which are the flags I use with another Cortex-M4 product). There's a good chance you'll get better performance from using a vendor-supplied function rather than your own code. Try arm_cfft_radix4_f32 in the CMSIS DSP library. This should solve your performance problems.

If you use a M3 part instead, and/or are really short on RAM for your samples, you'll have to switch to fixed point (Q15 format). Again, you don't have to write code for that, just grab the CMSIS DSP library provided by TI and use the functions arm_rfft_q15 (for real input) or arm_cfft_radix4_q15 (for complex input, if the input size is a power of 4).


Can FFT convolution be faster than direct convolution for signals of large sizes?

There is probably a bit of a misconception here. In many application the signal is  runnning all the time or is VERY long: modem, audio stream, video etc. In this case you can't really define the "length" of the signal. The relevant metric is here "number of operations per input sample" not the "total  number of operations". If you watch a streaming movie, it's obvious that a two hour movie takes twice the operations of a one hour movie. But that's okay since you have twice the time to process it. The key metric is if your processor can keep up with the stream , so that's why "ops per sample" or "ops per second" are used.
In this case, you don't do an FFT over the entire input signal but you apply "overlap add", i.e. you chop the input signals into blocks that are about the same size as the filter length and then you process one block at a time.
So the effectiveness of the algorithm is only a function of the filter length. And it increases linearly for direct convoltuion but only logarithmically for Overlap Add (which is FFT based). The "break even" point is typically at "many tens of samples".
We can run some numbers assuming that the filter length is $N$ and lets also assume its a power of 2.
An FFT of length $M$ has  $log_2(M) \cdot M/2$ butterflies, each with one complex multiplies and two complex adds. That's a total of 10 real operations. So the price for an FFT is roughly $5 \cdot log_2(M) \cdot M/2$
Overlap add requires zero padding to 2N and that's our FFT size. Hence overlap add costs

Forward FFT: $5 \cdot log_2(2N) \cdot 2N/2$
Spectral Multiply: $3*N$  (N/2 complex multiplies)
Inverse FFT: $5 \cdot log_2(2N) \cdot 2N/2$
Overlap adding: $N$
Total: $N\cdot (10 \cdot log_2(N) + 10 + 3 + 1)$

So overlap add cost per sample is roughly $ 10 \cdot log_2(N) + 14$ as compared to $2N$ for direct convolution. The break even point is a filter length of 32.
DISCLAIMER: I made a lot of simplying assumptions and in a practical implementation the break even point can vary wildly based on processor architecture and smartness of the algorithm. The main cost is in the FFT and for a single channel real system you would deploy a core complex FFT of N/2. I also skipped over bit reverse ordering and other overheads.


How to Remove a Glare / Clipped Brightness from an Image?

Let's assume glare portions are the only saturated areas in the image.
Detection can be performed by thresholding the intensity (code in Mathematica):

saturated = Binarize[ColorConvert[img, "Grayscale"], .9]

Then we need only to replace the portions of the image around the saturation mask (enlarging the mask is done by the morphological function Dilation). Inpainting using texture synthesis (using the function Inpaint) seems to work well in this example, though I can not test it as input to your defect-detection algorithm:

Inpaint[img, Dilation[saturated, DiskMatrix[20]]]


Evaluating the autocorrelation function of a telegraph signal to get exponential decay?

A semi-random telegraph signal is one for which $X(0) = 0$ and $X(t)$ equals $0$ or $1$
according as there have been an even number or an odd number of arrivals in a Poisson Process in the interval $(0,T]$.  A simple physical model would be the least significant
bit (LSB) of a digital counter that is counting the Poisson arrivals.  A random telegraph signal is almost the same except that $X(0)$ is equally likely to have value $0$ or $1$.
(Feel free to translate this to signals having values $\pm 1$ if that is more convenient
or more in accordance with your notation).

With this as background, you are encountering problems because your simulation is
not modeling the random telegraph signal correctly. You should instead simulate a
Poisson process and obtain the random telegraph signal from that.


This reference should help.  If you carry out the matrix-vector multiplication, you will see they are the same.

Also note that it is only possible to express the convolution as a matrix multiplication when the filter impulse response is of finite duration, i.e. a FIR filter.

The equivalence $h_{n,m} = h_{n+q,m+q}$ comes from the "TI" part of LTI, that is, time-invariance. If you delay the input by $q$ samples, the output should be the same, just delayed in time by $q$.


Variable Step Size LMS vs Leaky LMS Adaptive Filter Algorithm

Variable step size LMS is generally used to improve speed of convergence or decrease steady state error.
Leaky adaptation is used to combat problems like potential instability of the filter in a finite-precision implementation. It is closely related to L2 norm regularization technique and results in continuous down scaling of filter coefficients (hence smoothing the extracted filter).
So comparing these two techniques is not that meaningful. In fact you can easily use both techniques together to achieve the above advantages.


You aren't doing anything wrong, but you also aren't thinking carefully about what you should expect to see, which is why you're surprised at the result. For question 1, your conjecture is close, but you actually have things backwards; it's numerical noise that's plaguing your second one, not your first one.

Pictures may help. Here's plots of the magnitude and phase for the first trial:

x = Cos[2.0 \[Pi] 200 Range[0, 1, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

And the second one:

x = Cos[2.0 \[Pi] 200 Range[0, 1 - 1/1000, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

So what's going on here? The second one is easiest to explain. First, the FFT of the second one has zero amplitude everywhere except at the two peaks visible in the magnitude spectrum; this follows since the definition of the FFT using 1,000 datapoints returns frequencies of the form $k/1000$ for $0\leq k\leq 999$, and hence your signal falls exactly on a frequency bin. As a result, at the other 998 points your signal is entirely due to machine noise due to floating point errors, and hence your phase spectrum is nonsense, because it's literally the phase of pseudorandom numbers.

For the first one, however, the definition of the FFT encompasses frequencies of the form $k/1001$ for $0\leq k \leq 1000$, whereas your signal frequency is $200/1000$, which is not of the form $k/1001$. As a result, your signal is broadened by spectral leakage, and will be nonzero almost everywhere. I will not comment on the physical form of the phase plot, but I will say that it admits a closed analytic form.

In general, I think graphs of phase angle alone are a really bad idea for conveying information, precisely for this reason; first, you can't tell if you're looking at the phase of low-amplitude garbage or actual signal, and second, it's not translation-invariant, and it's easy to get utterly bewildering graphs for simple inputs. Much better, if you're still looking for something which conveys phase information, is a graph which simultaneously portrays phase and amplitude information in the same visual manner, such as a plot where phase is encoded as hue and magnitude is encoded as brightness.

ADDENDUM: Here's a couple pictures from Mathematica which illustrate the principle I stated in the preceding paragraph:

hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + \[Pi])/(2 \[Pi]),
Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]},
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
L = 500;
data = Table[Boole[x <= 11 && y <= 11], {x, L}, {y, L}];
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-5, -5}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-4, -4}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {0, 0}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]

All three images are 2D Fourier transforms of the same input signal (an $11\times 11$ square of 1's padded with zeros to a length of $500\times 500$), but the inputs have been cyclically rotated by 5, 4, and 0, and 200 datapoints. The magnitude spectra (encoded by pixel brightness) are identical, but the phase spectra are completely different! Phase encoding is done so that 1 maps to red, $i$ maps to green, $-1$ maps to cyan, and $-i$ maps to purple. This is what I mean when I say that phase spectra are non-shift invariant, and thus not amenable to human visual understanding. For example, with a cyclic shift of 200 datapoints, it's completely impossible to tell what's going on in the phase, since it just looks like static, but the input signal is no more complicated than the other input cases.


Design the coefficients of Wide-band FIR filters

As it stands your question appears to be quite broad, so the answer depends largely on your application and level of knowledge.

If you have no prior experience, you might want to consider picking up a book about digital signal processing, e.g.

Digital Signal Processing by John G. Proakis and Dimitry Manolakis

which may be available in your local library.

For digital filters specifically, there are great resources available online such as this free book:

INTRODUCTION TO DIGITAL FILTERS WITH AUDIO APPLICATIONS by JULIUS O. SMITH III

There is also the Cookbook for audio EQ biquad filters by RBJ which you should also check out, found for example here (although i cannot find the original link right now).


Understanding Asymptotic Equipartition Property

The sentences you've mentioned in your question all have a quite different meaning. Let me try to explain them one by one:

  Clearly, it is not true that all $2^n$ sequences of length $n$ have the same probability.

This refers to the example of a sequence of binary i.i.d. random variables. If $p(1)\neq p(0)$ then it is clear that sequences with different numbers of ones and zeros must have different probabilities, i.e., not all sequences have the same probability.

  One might be able, however, to predict the probability of the sequence that is actually observed. The question is: what is the probability $p(X_1,X_2,...,X_n)$ of the outcomes $X_1,X_2,\ldots, X_n$, where $X_1,X_2,\ldots, X_n$ are i.i.d. according to $p(x)$?

This question is answered by the asymptotic equipartition property (AEP). It turns out that the typical sequences do have approximately the same probability as $n$ gets large, namely $2^{-nH}$.

  The typical set has probability close to 1, all the elements of the typical set have [the] same probability.

That's also what the AEP says. As mentioned before, the typical sequences have approximately the same probability, and the aggregate probability of the typical set approaches $1$ as $n$ gets large.

Finally, the statement

  The concept of a typical set is generally different from that of a high probability set.

is indeed confusing. What is meant is the fact that there can exist individual sequences that are not part of the typical set but which have a higher probability than the typical sequences. Still, the aggregate probability of the typical set approaches $1$ for increasing $n$, and there are too few atypical sequences with high probability to have any significance.


Convolution in DFT is still circular.

Think of the DFT as taking the 1st period (in time and in frequency) of the DFS (discrete Fourier series). In DFS, both the time sequence and the frequency sequence are N-periodic, and the circular convolution applies beautifully.

I personally think all properties in terms of DFS, and then consider the 1st period when speaking DFT.


Which book of Proakis has dealt with continuous time signals and systems?

No, not as far as I know. I can recommend you two books dealing with continuous-time systems, the Fourier transform, and the Laplace transform, among other topics:

Signals and Systems, Oppenheim, Willsky: a classic, used copies can be found quite cheap
The Fourier Integral and Its Applications, Papoulis: less scary for an engineer than the title might make you think; it's actually much broader and very accessible; it deals with linear systems in general, low pass and band pass systems, the Laplace transform, etc.


Finding radius r of the overlappable sphere(s) in 3D image

A simpler solution and much more computationally efficient when compared to Hough Transform is to use the distance transform:

Find the surface of your spheres (i.e. the pixels that have value 1 and have at least one neighboring 0 pixel);
Compute the distance transform with respect to the spheres surface, but constrain the computation only to pixels that are internal to the spheres. The output will be a distance map;
The radius will be exactly the maximum value in your distance map.

Another advantage of this solution when compared to Hough transform is that it provides a much more precise value for the radius.


How would a Linear Gray Level Transform affect an image?

The graph you've shown above describes a pixel-wise operator, basically a function which maps input pixel values $r$ to a new value $s$ in the output image. i.e. each of the functions in the figure you provided describes a transformation $s = T(r)$.

These transformations are typically used to improve image contrast through histogram equalization, basically spreading out the pixel values to make it easier for us to distinguish objects in the image.

To implement something like this, you have to iterate through every pixel, get its grey level $r$, and replace this value with the output $s$ of the function $T(r)$.

For converting color to greyscale, it's the same idea, although this time your function $s = T(r, g, b)$ has three inputs (r = red, g = green, b = blue), since each pixel in a color image will return these three values. Implementation is the same, just iterate through every pixel and replace these three values with just one value $s$ based on the function.

While the identity transform above (given by $s = T(r) = r$) is linear, linear transformations do not necessarily mean identity.

An example of a linear rgb color to greyscale transformation (the one used by MATLAB) is $s = T(r, g, b) = 0.2989 * r + 0.5870 * g + 0.1140 * b$. However, as long as you keep the function in this form (and just change the coefficients to constant values), your transformation will remain linear; this does not guarantee that the transformation will produce good results, though.

EDIT: What I described is only a subset of linear transformations. A linear transformation is any function $T(r)$ where $T(A * r1 + B * r2) = A * T(r1) + B * T(r2)$


Faster way to find the quantization level of a sample

As hinted by Florian, to quantise any $x[n] \in \mathbb{R}$ to a $y[n] \in \mathbb{Z}$ with $K$ levels you can apply the function:

$$y[n] = round\left( \frac{x[n]-min(x)}{range(x)} (K-1)  \right)$$

Where $n \in [0,1,2..|x|]$ and $|.|$ denotes the length of the sequence.

This will map your $x$ from its $min(x)..max(x)$ range to a $[0..K)$ discrete steps range.

For example (in Octave):

Fs = 100; #Sampling frequency (in Hz)
f = 4; # Frequency of a simple sinusoidal signal (in Hz)
T = 1; # Timeframe length (in seconds)
K = 5; # Number of levels to quantise to
t = 0:(1./Fs):(T-(1./Fs));
p = 2.0 * pi * t;
x = sin(f*p);
y = round(((x-min(x))/2.0) * (K-1));
plot(y);
xlabel("Discrete time (sample)");
ylabel("Amplitude");
grid on;

Which produces:

You might also like to see uencode

Hope this helps.


To make it simple, finding user location is done by using four pseudo-range measurements (approximative distance from user/receiver to satellite) to resolve a four unknown equation. As information, these unknown are receiver's (X,Y,Z) coordinates and it's clock offset.

These pseudo-range measurements are done by tracking stage of the receiver. Acquisition results (code phase and Doppler shift) are used to initialize delay and phase tracking loops. I let you the following simple explanation of the receiver behaviour : https://gssc.esa.int/navipedia/index.php/Generic_Receiver_Description

Please let me know if my explanation are not sufficient,

Math


$\rm Q$ number format singed or unsigned?

From the name of the format alone ("Qm.n") you generally can't tell if the numbers are signed or unsigned. There is a convention that would use a prefix "U" to indicate that the numbers are unsigned, but that's not used universally. Often when the numbers are signed, the sign bit would not be included in the number m, so by comparing the register size with the sum m+n you could figure out if there's an extra sign bit or not. But that's also not a 100% indication because sometimes the sign bit is included in the number m, so m+n equals the register size, but the numbers are signed.

So in general you need extra information to be sure if the numbers are signed or not. In some cases, you can get a hint from comparing the register size with the sum m+n.


GPS | Retain Signal Tracking after Noise Bound

Now that the OP has clarified in the comments that this is an IF file for beginners I understand the issue. The point is to not be able to demodulate or recognize the individual chips of the GPS PRN sequence, in fact in typical application the received signal on the surface of the earth can typically be -20 dB SNR (pre-correlation) meaning it is 100x lower than the power of the noise! There is no way in such a signal to recognize any such sequence. However by multiplying with the received signal once translated completely to baseband with an aligned PRN code (the GPS C/A code for a particular satellite of interest) and then summing the result over a complete 1023 chip symbol is a great demonstration of processing gain and spread-spectrum demodulation (whose intent is to "pull signals out of the noise").
I detail this in the slides copied below, starting with the sequence as transmitted from a GPS satellite, representing the phase of the 1575.42 MHz carrier going back and forth between 0 and 180 degrees according to the sign of the C/A code being transmitted. The GPS C/A code, unique to each satellite, repeats every 1023 chips as a pseudo-random sequence as a PRN symbol at a 1.023 MHz rate. After every 20 symbols the PRN can be completely inversed (multiply by -1) according to the sign of the data bit being transmitted from the satellite, meaning the satellite will send 20 symbols to send a data bit "1" or the same 20 symbols all inverted to send a data bit "0" and thus transmits a 50 bps data rate:

Below is live data that I captured (using a Lecroy O-scope as the A/D converter hence shown in the graphic), translated to baseband (all carrier frequency and phase offsets removed) and aligned using acquisition techniques to be exactly in phase with the transmitted code which is also shown in red. The received captured data in blue for all intents and purposes appears to us a additive white Gaussian noise, yet the code sequence is part of that waveform and using correlation we can extract the symbols being transmitted from the satellite (you may think that you can recognize in parts the original PRN, but trust me, you cannot possibly!):

A histogram of the actual captured waveform further confirms the appearance of Gaussian noise (and an FFT would confirm it is "white" over the sampling bandwidth used).

Here are the samples of the received signal on a complex IQ plot after being translated back to baseband and all carrier frequency and phase offsets removed. The red dots shown are all at the chip centers, which when transmitted at the satellite existing only at +1/-1 on the real axis. This plot shows the significance of the added noise to each of those samples, totally swamping out the actual chips that were transmitted.

But after doing a sliding correlation with the same reference C/A code, we are able to clearly make out the boundaries of each of the 1023 chip symbols, and their overall phase associated with the 50 bps data transmitted by the satellite. Each sudden peak below occurs when the received waveform (after having been translated to baseband with all carrier frequency and phase removed) aligns in time with the reference C/A code, which is slid past the received waveform and after every offset a new correlation is computed (cross correlation function of the received sequence with the C/A code). We can see how there are at least 20 symbols in a row before the entire sequence is inverted (to change the data in the satellite between 1 and 0).

Hopefully this makes it all clearer that the intent of demodulation with GPS is nor to be able to extract what would appear to be a noise-free or even low noise C/A code itself, but operate at the symbol level after correlating with a received waveform where the C/A code is unrecognizably buried in the noise.  The extracted C/A code symbols after correlation become the cleaner signal (at a much lower rate) with which to them perform typical receiver processing on.  So the intent of tracking (and acquisition!) is to operate on post-correlation symbols in the latter low noise condition; the actual signals as received in practical application will never have the benefit of the high SNR as presented in the first half of the OP's waveform. Any solution that depends on such high SNR will therefore be flawed in its approach.
For further details of GPS and GPS acquisition in negative SNR condition, where some of these same figures are repeated, please see:
Auto-correlation of C/A code with 30dB noise
Phase transitions in GPS C/A signal
GPS Coarse Acquisition PRN Codes
GPS signal acquisition
GPS Receiver in Frequency Domain
Understanding time delay estimation for GPS?
Isn't Gardner's algorithm and Early-Late gate the same thing?


What is the point spread function and optical transfer function and what uses are they in image processing

PSF (point spread function) is the spatial domain impulse response of the optical imaging system assuming it's modeled an LSI (linear shift-space invariant) system according to its light intensity processing.

In principle, PSF can obtained by inputting an impulse (varry small bright dot - a point) to the system, and measuring the output response where the dot is typically spread into a blurry circle; hence the name point-spread function.

Also for CRT display systems, PSF can be obtained by sending a single pixel dot signal to the display, and observing its resulting brightness spread on the screen.

Once you obtain the impulse response (or PSF) of an LSI system, then you can mathematically compute the output response of the system to any input by convolving the PSF with the input signal. This concept lies at the heart of the signal processing, and that's why PSF is a very important characterisation of optical systems.

OTF is the Fourier transform of the PSF, and therefore carries the same information, but, stated in frequency-domain instead.


How to find derivative of an image along a radial direction with an arbitrary center of rotation?

Essentially, you want to compute the derivative of your image in polar coordinates. Have a look here in particular equation (15):

Your image is given in cartesian coordinates: $I(x,y)$ where $x,y$ are column and row. You want to know the expression

$$\frac{dI(r,\phi)}{d\phi}$$

where $I(r,\phi)$ is the image in polar coordinates. $r,\phi$ relate to $x,y$ by the expressions

$$\begin{align}
x &= r\cos(\phi)\\y &= r\sin(\phi)
\end{align}$$

Now, using the chain rule we get

$$\frac{dI}{d\phi}=\frac{dI}{dx}\frac{dx}{d\phi}+\frac{dI}{dy}\frac{dy}{d\phi}=-\frac{dI}{dx}(r\sin\phi)+\frac{dI}{dy}(r\cos\phi)$$

where $\frac{dI}{dx}$ and $\frac{dI}{dy}$ are the gradients in cartesian coordinates.


Convolutions with changes to the argument

Let me show you how to manipulate the convolution integral.
Let
$$y(t)  = x(t) \star h(t) = \int_{-\infty}^{\infty} x(\tau) h(t-\tau) d\tau \tag{1} $$
Let $a > 1$ that we have a compressed signal $y(at)$ ;
$$y(at)  = \int_{\tau=-\infty}^{\infty} x(\tau) h(at-\tau) d\tau \tag{2} $$
Substitude $\tau' = \tau/a$ so that $\tau = a \tau'$ and $d\tau = a \cdot d\tau'$ into Eq.2 :
\begin{align}
y(at)  &= \int_{\tau'=-\infty}^{\infty} x(a \tau')~ h(at-a\tau') ~a ~d\tau' \tag{3}\\ \\
&= a \int_{\tau'=-\infty}^{\infty} x(a \tau') ~ h(a(t-\tau')) ~d\tau' \tag{4}\\ \\
&= a \int_{\tau=-\infty}^{\infty} x(a \tau) ~ h(a(t-\tau)) ~d\tau  = a \cdot \big( x(at) \star h(at) \big) \tag{5}
\end{align}
Eq.5 shows that the compressed output $y(at)$ is obtained by the convolution of the compressed input $x(at)$ and the compressed impulse response $h(at)$, and (also the result scaled by $a$)....
Note that convolution operator assumes an LTI (linear time-invariant) framework. However, the compressor is not an LTI system and it does not have an impulse response $h_c(t)$, so that $y(at) = y(t) \star h_c(t)$ can not be true.
Assume for a moment that the compressor is an LTI system with impulse response $h_c(t)$. Then we can consider the generation of the compressed output in two (serial) LTI stages:
$$ x(t) \longrightarrow \boxed{ h(t) } \longrightarrow \boxed{ h_c(t)} \longrightarrow y(at) \tag{6}$$
From LTI theory we know that the cascade of two LTI blocks is also LTI with an equivalent impulse response given by $h_e(t) = h(t) \star h_c(t)$ hence Eq.6 can be written as :
$$ y(at) = x(t) \star h_e(t)  = x(t) \star ( h(t) \star h_c(t)) = x(t) \star h(at)\tag{7}$$
Which is contradicting the result we have obtained in Eq.5, hence the assumption that the compressor could have been represented with an impulse reponse was wrong.
Another implication of this is also observed in the time-reversal of the output $y(t)$ :
$$y(-t) = x(-t) \star h(-t) \neq x(-t) \star h(t) \neq x(t) \star h(-t) $$
The reason of the unequal signs is again about reversal operation not being an LTI system.
On the other hand, if you want to apply an LTI operation on the output such as $y_d(t) = y(t-d)$ , then it can be simply be reflected in the convolution as:
$$ y(t-d) = y(t) \star \delta(t-d) = (x(t) \star h(t)) \star \delta(t-d) = x(t) \star h(t-d)$$...


(For c-FastICA) On covariance and pseudocovariance matrix of a complex random vector

You have

$$E[s_ks_k^*]=E[|s_k|^2]=1$$

because the complex random variables $s_k$ have zero mean and unit variance. That means that all elements of the main diagonal of $E[\mathbf{s}\mathbf{s}^H]$ equal $1$. Furthermore, with $s_k=x_k+jy_k$we have

$$\begin{align}E[s_ks^*_l]&=E[x_kx_l+y_ky_l+j(x_ly_k-x_ky_l)]\\&=E[x_kx_l]+E[y_ky_l]+jE[x_ly_k]-jE[x_ky_l]\\&=0,\qquad k\neq l\end{align}$$

because real and imaginary parts are uncorrelated. Consequently, all off-diagonal elements of $E[\mathbf{s}\mathbf{s}^H]$ are zero.

The off-diagonal elements of $E[\mathbf{s}\mathbf{s}^T]$ are zero for the same reason. (There's just a sign difference in the sum of the expectations, but since each of them is zero the result is the same). The main diagonal elements are

$$\begin{align}E[s_k^2]&=E[x_k^2-y_k^2+2jx_ky_k]\\&=E[x_k^2]-E[y_k^2]+2jE[x_ky_k]\\&=0\end{align}$$

because real and imaginary parts are uncorrelated, and they have equal variance, i.e., $E[x_k^2]=E[y_k^2]$.


Generate the Convolution Matrix of 2D Kernel for Convolution Shape of `same`

I cannot test this on my computer because I do not have the convtmx2 function, here is what the MATLAB help says:

http://www.mathworks.com/help/toolbox/images/ref/convmtx2.html

T = convmtx2(H,m,n) returns the convolution matrix T for the matrix H. If X is an m-by-n matrix, then reshape(T*X(:),size(H)+[m n]-1) is the same as conv2(X,H).

This would get the same resulting convolution of conv2(X,H) but then you would still have to pull out the correct piece of the convolution.


Yes this would properly delay the input by 64 samples but note that you can simply do:

$out[n]=in[n-64]$

Which was probably clear to you and you just wanted to show how it related to a 65 tap FIR with all taps = 0 except the last one.

The coefficients of the FIR filter represent the impulse response (strictly in discrete terms we call this the "unit sample response") of the filter. In this case the impulse response is a delayed impulse. The Fourier transform of a delayed (time shifted) impulse has a magnitude response in frequency that is all ones, and a linear phase response where the delay is equal to the phase slope according to:

$\text{delay} = d\phi /d\omega$  where $\phi$ is the phase in radians and $\omega$ is the radian frequency in radians/second.

We see this specifically from the Time Shift property of the Discrete-Time-Fourier-Transform:

$$y[n] = x[n-N]$$
$$Y(\omega)=e^{-j\omega N}X(\omega)$$

You can quickly plot the frequency response using Octave or Matlab as follows:

freqz([zeros(1,64) 1])

With the following result:

With a linear phase filter such as this, the signal is not distorted by the filter, only delayed. For more on that see this post: why is linear phase important


How to find transfer function looking bode plot

locate the 3db point. (The point where the line before and after the curve would meet.) Take note of the magnitude and frequency f.

Find the angular frequency at this point. w = 2 * pi * f

If the graph is going up then it's a zero. H(s) = w + s

If it's going down then it's a pole. H(s) = 1 / ( w + s )

Once you've got that. check the slope. If it's 20dB/decade then you're set. If it is 40, you'll need to square it, 60 and you'll cube it etc.

Finally you'll need to find the gain constant to make it correct at low frequency (or DC)

So you're answers will look like this:

H(s) = Gain * ( w + s )
H(s) = Gain / ( w + s )


FFT Phase interpretation of input signal with non-integer number of cycles in FFT window

There is a circular waveform discontinuity at sample 0 of an FFT input (to sample N-1), if it isn’t exactly integer periodic in aperture.  However, if the waveform is continuous at sample N/2, then the phase can be measured at that point.

You can measure phase at the halfway point by doing an FFTshift, or by flipping the phase of every odd numbered FFT output bin.  If you know, or can estimate the frequency, then if you estimate the phase at sample N/2, you can use the frequency and phase at a known point to compute the phase at any other point in time.  Including at the beginning of the original sample data window.

This answer is for even N Lengths of FFT.  (Dawg may have a better answer for odd lengths.)


Does a .wav file created from a .mp3 file have better quality than the .mp3 file itself?

MP3 is a lossy compression format, so encoding an uncompressed lossless WAV file to MP3 will lose some quality in the encoding process.

Decoding a MP3 file back to WAV file has at most the same quality as the MP3 file has, it is just as good or bad as the contents in the MP3 file was, it is not in any way different to decode it for saving into a file or decoding it for playback. Most importantly, you don't get the lost quality back when decoding, and basically the result actually depends on the quality of the decoder as well.


How to remove the bump from the signal using filters?

A bump like this one is likely to be wide-band, especially with the sharp onset. Plus, the line may be hard to deal with in the Fourier domain.
 Hence, the combination is complicated to remove with a classical linear filter. The problem is very akin to baseline, background or trend removal, answered elsewhere here.

Several options are possible, for instance:

use a non-linear filter, based on a median, or a minimum/maximum statistics,
use morphological operators: a rolling ball, lot of straight segments, etc.
use a knowledge on the data model, like a linear equation: $y=ax+b$, or the fact that the bump is "above",
combine the above in a variational formulation, using appropriate data fidelity and penalty.

In your example, I suspect that a classical linear fit with robust distance (like a least-absolute distortion) could do the job. I will call all the above filters, in the wide sense that you will remplace a value with respect to some sort of combination of the others.

You can also call the following robust regression, LAD fitting. An example at work:

% Standard and Robust fit of a degree 1 polynomial w/ a bump
nSample = 1000;
% Create a similar composite signal
time = linspace(0,5,nSample)';
polyCoef = [0.2  0];
dataLine = polyval(polyCoef,time);
dataParabola = -8*(time-2).*(time-3);
dataParabola(dataParabola < 0) = 0;
data = dataLine+dataParabola;

% Use Matlab curve fitting toolbox
optsRobust = fitoptions('Method','LinearLeastSquares','Robust','LAR');
[fitObject,gof] = fit(time,data,'poly1',optsRobust);
h1=plot(fitObject,time,data);
grid on


Adaptive Gaussian Filter for Image Denoising

It's not clear in your case that Gaussian kernels will help you.

Since you have a noise shape that is closer to salt-and-pepper noise instead of additive white noise, you may want to try other algorithms (that are somehow linked to an $L_1$ best fit instead of a least squares solution):

you can try median filtering, since it is by design more robust to this noise shape,
you can also try to minimize some TV-$L_1$ functional over your image, see Chambolle-Pock's paper on primal dual minimization for a concrete example,
you also can use a bilateral filter. You will have to play with the exponential decay of the pixel similarity function ($\sigma_d$ in Wikipedia page) to get the desired behaviour.


Equation (1) in your question is a pure additive noise channel without any further distortion of the signal. Equation (2) takes a delay into account, but this is usually not relevant. A pure delay does not affect the quality of the communication link in terms of bit error rate. (I'm not talking about the annoying effect of delay in voice communication).

So the model in Equation (2) is often not relevant and can be replaced by Equation (1). A more general and useful model, which includes Eq. (2) as a special case, is

$$y(t)=(x*h)(t)+n(t)$$

where $h(t)$ is the impulse response of the channel, and $*$ denotes convolution. This model represents a linear time-invariant channel with additive noise. Eq. (2) is obtained as a special case with $h(t)=\delta(t-t_0)$.


Determine whether the system is a FIR or IIR by looking the transfer function

The forward $z$-transform is written with negative powers on $z$. So, $z$ transform of your difference equation would be:

$$
Y(z) = \frac{1}{3}\left(zX(z) + X(z) + z^{-1}X(z)\right)
$$

yielding a transfer function of:

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

The easiest way to determine whether a filter is IIR or FIR is to identify its pole locations. For FIR filters, there is a rule for this that is based on the structure of the impulse response:

If the system is causal (i.e. it is zero for all $n < 0$), then it is FIR if all of its poles are located at the origin ($z=0$).
If the system is anticausal (i.e. it is zero for all $n \ge 0$), then it is FIR if all of its poles are located at $z=\infty$.

Your system is neither causal nor anticausal. Instead, it has a causal component ($\frac{1}{3}(x[n] + x[n-1])$) and an anticausal component ($\frac{1}{3}x[n+1]$). So, it can be helpful to look at each component separately:

$\frac{1}{3}(x[n] + x[n-1])$: The $z$-transform of the causal portion of the transfer function is $\frac{1}{3}(1 + z^{-1})$. This expression $\to \infty$ as $z \to 0$, so it has a pole at $z = 0$. This meets the criterion for an FIR response.
$\frac{1}{3}x[n+1]$: The $z$-transform of the anticausal portion of the transfer function is $z$. This expression $\to \infty$ as $z \to \infty$, so it has a pole at $z = \infty$. This meets the criterion for an FIR response.

This would imply that you gave a difference equation for a non-causal FIR filter.


How does causality (i.e. unit step) affect the DTFT of a sine or cosine wave?

Multiplication in the time domain is convolution in the frequency domain. The frequency response of a unit step has 1/frequency shape and it's pretty messy.

What you will see is that the overall frequency response is not a delta impulse anymore but that the lines widens and that the get "skirts". These will bigger and more "wiggly" for the cosine than for the sine as the step function slices right through the maximum of the cosine.


How to Create a Convolution Matrix with a Variable Condition Number (CN)

Let's try to think of it intuitively.
Given the LPF what would make the inverse "Hard"?

The sections we won't be able to recover are where the LPF is 0, since there is no inverse for that we can multiply the result by.
Real world LPF won't reach "Real" zero as usually.
But what would make it hard to recover is where there are big ratios between the biggest and the smallest magnitude.
As close and fast your LPF goes to zero (And has magnitudes which are big) the harder the recovery will be.

Now just build analog LPF which this requirements at different scales, digitize it, create the Convolution Matrix and there you have it...

A practical example is given in my answer to 1D Deconvolution with Gaussian Kernel (MATLAB).


Why do we always characterize a LTI system by its impulse response?

Let's try a simple example.  A discrete-time LTI system has unit pulse response
$$h[0] = 1,\quad h[1] = 1,\quad h[n] = 0 ~~\text{for all}~ n \geq 2.
\tag{1}$$
The response $y$ of this system to a signal $x$ given by
$$x[0] = 1,\quad  x[1] = -1,\quad x[2] = 1,\quad x[n] = 0
~~\text{for all}~ n \geq 3 \tag{2}$$ is readily calculated
via the convolution formula,
and we get
$$y[0] = 1, \quad y[1] = 0, \quad y[2] = 0, \quad y[3] = 1,\quad y[n] = 0
~~\text{for all}~ n \geq  4.\tag{3}$$

============================================================

The unit-step response $g$ of the same LTI system is easily found to be
$$g[0] = 1,\quad g[n] = 2 ~~\text{for all}~ n \geq 1.\tag{4}$$

So, suppose that you are given only the step response $g$ in $(4)$
but not the pulse response in $(1)$.  How will you compute the response
of this LTI system to input $x$?  It is known that the answer is $(3)$,
but you are asked to figure out how to get from $(2)$ and $(4)$ to $(3)$.


What Is the Relation Between Deblurring and Deconvolution in Computer Vision and Image Processing?

Let me present the following Diagram:

So, both Deblurring and Deconvolution are operations within the family of Image Restoration (Which is a subset of Inverse Problem set).

Basically we build the Image Restoration set by different Degradation Models.
The one related to the question are:

Linear Degradation Model
Namely, the degradation is made by a Linear Operator.
Spatially Invariant Model
A Model where the degradation is the same for any place in the image.
We created it as a subset of Linear Model though it doesn't have to be. But for clarity.
Any operator which is both Linear and Spatially Invariant can be defined by a Convolution Operation. Hence it can be reversed by Deconvolution.
Low Pass Operator
A set of images degraded by a Low Pass Filter. Namely by a convolution with a Low Pass Filter.

Now, there is a set of degradation which basically create a blurry image.
Reversing this operation is called Deblurring.
In case the blurring is made by a Low Pass Filter applied by Convolution the Deconvolution of this operation is also a Deblurring process.


How do I calculate correlated color temperature correctly?

The problem is that I don't know whether should I apply white balance process before calculating CCT.

Definitely not. Afterwards, the light source would have a completely different color profile, depending on how you implement white balance, and hence, a different color temperature.

  If it should not, how do I get the conversion matrix that convert the current device-dependent RGB color space to CIE XYZ to calculate the algorithm correctly?

Get the color profile (which might be much more complex than a simple matrix) of your device; usually, that means getting the ICC file for your camera/scanner/sensor, and applying that to convert your image from the device color space to something canonical.

If your camera/scanner/sensor doesn't have such a profile, you will need to calibrate it first.


Discrete time system signal flow graph problem

There is clearly recursion in there and your solution does not reflect it since its purely FIR.

I recommend doing this as a coupled difference equation: Define a new variable $v[n]$ at the point where Peter drew the circle and than derive the transfer function from (1) X->V and then (2) X and V -> Y. Then you can eliminate V in eq (2) by using equation (1).


How does using DCT help in image compression

If you take an image as is (especially a photograph, as opposed to line art, a document scan...), there's very little redundancy in it for Huffman encoding to do its work - because the image contains many small fluctuations and details which are not perceptible to the human eye but nevertheless present in the data. DCT followed by coefficient truncation/quantization eliminates the barely perceptible details and creates redundancy for Huffman compression to work on.


MATLAB flanger...zipper noise, crackling, help with interpolation values?

ok so it turns out setting indices is tricky in MATLAB (who knew :p)

first, a less needlessly complicated interpolator:

function result = dLinTerp(y1,y2,x)

% use weighted sum method of interpolating
result = x*y2 + (1-x)*y1;

and the flanger

% flanger
clear all;
% Creates a single FIR delay with the delay time oscillating from
% Either 0-3 ms or 0-15 ms at 0.1 - 5 Hz
[input,fs] = audioread('pianoriff_s.mp3');
input = input(:,1);
len = length(input);
out = zeros(len,1);

% parameters to vary the effect %
max_time_delay=0.005; % 7ms max delay in seconds
rate=0.3; %rate of flange in Hz

% convert delay in ms to max delay in samples
max_samp_delay= floor(max_time_delay*fs);
delayBuff = zeros(max_samp_delay,1);
amp=0.7;
readIdx = 1;

fetchPrev = 1;

% delaybase = .010;
% new_d_base = delaybase*fs;
% d_base = LIN_INTERP(0,1,0, new_d_base,frac);
%
% dp = (float)(delay_pos - d_base) - (delay_depth * law);
% %     // Get the integer part
% dp_idx = f_round(dp - 0.5f);
% %     // Get the fractional part
% dp_frac = dp - dp_idx;

% for each sample
for i = 1:len,
    delayBuff(readIdx) = input(i);
    sin_ref = 0.97*sin(2*pi*i*(rate/fs));
    sin_next = 0.97*sin(2*pi*(i+1)*(rate/fs));
    cur_sin=abs(sin_ref); %abs of current sin val 0-1
    cur_delay=ceil(cur_sin*max_samp_delay);

    cur_frac=cur_sin*max_samp_delay;

    if (readIdx + cur_frac >= (max_samp_delay))
        fetchIdx = cur_frac - (max_samp_delay - readIdx) + 1;
    else
        fetchIdx = readIdx + cur_frac;
    end

    if (floor(fetchIdx) = max_samp_delay)
        fetchNext = 1;
    else
        fetchNext = ceil(fetchIdx);
    end

    integ = fix(cur_frac);
    frac = abs(cur_frac - integ);
    interp = dLinTerp(delayBuff(fetchPrev), delayBuff(fetchNext),frac);
    out(i) = (amp*input(i)) + (amp*interp);

    if (readIdx+1 >= max_samp_delay)
        readIdx = 1;
    else
        readIdx = readIdx + 1;
    end

end
% write output

soundsc(out,fs);


There is no such proof because it's not always true.  It's a rule of thumb, because I guarantee that you could come up with a situation- an infinite number of situations actually- where higher order splines would do better than cubic splines.  The optimal spline order for any given situation is the exact same order as the system you are trying to model.  If the order is the same and your data points are error free (never the case, of course, except for in theoretical problems), then you should be able to model the system perfectly.

The reason they recommend not to go higher than cubic splines is because overfitting is really, really bad.  Overfitting can greatly magnify errors, while "underfitting" (choosing a spline method with an order lower than the order of the system you are modelling) introduces some low pass filtering that is either not that bad or sometimes even beneficial.


Active Noise cancellation for non periodic signals

Thats wideband and nonstationary interference cancellation. Using LMS adaptive filter and different adaptive filtering structures, you can also eliminate nonstationary and/or wideband interferences too. Not always possible though.


Effect of One Sample Latency in Audio Application

One sample latency is pretty much unavoidable so this is fairly normal. Please note that you need to look at the total system level latency which includes DA and AD conversion. Many audio applications use Sigma Delta converters which have much higher latency than one sample, so your processing latency wouldn't make much of a difference.
The impact of latency depends very much on the specific application. One sample is very low and the only applications that would be impacted by this is anything with a feedback loop (active noise cancellation, active excursion control) and any application that has a parallel analog path (effects loop in an analog mixer for example)
It's not a problem for UI responsiveness, audio/video sync, play/sing/talk along, multi-speaker synchronization, etc.
The most latency-sensitive audio application is active noise cancellation and many dedicated ANC processors run at a significantly higher sample rate, just to reduce latency and increase bandwidth of the feedback loop.


How identify echo patterns in audio samples

The side peaks could be the result of periodic structure in your input signal (say around on the order of 8-10Hz, plus a 2nd harmonic).

If you are creating the sound samples, try playing pseudo-random noise, which has good autocorrelation properties for a long enough sample, notably tight peaks & aperiodic.

2nd John's answer. You can do a back-of-the-envelope estimate of the max echo lag (assume a few reflections) using the speed of sound & room size.


Estimate Variance of Additive White Gaussian Noise (AWGN) Given Multiple Realizations with Different Mean

Basically each pixel is a realization so all you need is to work in the 3rd dimension (Though you can also get better by using the Spatial Data).

Method 1

So the trick here is to use the multiple images to estimate the Mean (True value) of each pixel and then calculate the STD on all samples (numRows * numCols * numRealizations).

Assuming we have single channel image. So we pack all give images into a tensor tI with dimensions: numRows * numCols * numRealizations.

We calculate the mean per pixel (Averaging on the 3rd dimension) and then subtract each image from the calculated average image.
Then we're left with many realizations of the noise (Well, noise and the left over from the estimated mean error).

Method 2

Estimate the STD in the 3rd dimension per pixel. Them average all numRows * numPixels estimations.

Method 3

Just as method 2. But since the property which obeys to linear operatiosn is the Variance calculate the variance along the 3rd dimension, average it over all pixels and then take the sqrt() of the average Variance.

Summary

Here is a graph showing the estimated noise STD as a function of the number of realizations:

As can be seen, for an image with many pixels (Say more than 20,000 or so) Method 3 becomes almost perfect.

The full code is available on my StackExchange Signal Processing Q61273 GitHub Repository.

Remarks

Method 1 - Basically the Maximum Likelihood Estimator for the problem. It is biased as this is the property of the ML for the Variance / STD.
Method 2 - Works nice, really intuitive.
Method 3 - Basically the MLE with the non biased estimation for the variance. As can be seen. It is also the best in real world.

Yep, Method 3 is much better only because it divides the data by N while Method 1 divides by N - 1. But it has significance impact on performance.


As with most things in Engineering the answer is "depends on your requirements". In general you can't do a perfect Hilbert transform unless the signal is strictly periodic (see http://andrewduncan.net/air/ ). Since you can't do perfect you need to define what's good enough for your particular application.

For most typically audio applications, that it's pretty straight forward. The Hilbert transform can be represented as a simple linear time invariant filter. Unfortunately the impulse response infinite and non-causal. For an FIR implementation you need to truncate and delay it. The larger your truncation window and your delay the more precise it will be, however it will also become more expensive computationally and the latency goes up.

If your application requires a really long filter (typically when you need high precision at very low frequencies) than FFT can be useful in an overlap add configuration.

If you just care about relative phase and not absolute phase you can run the signal through a set of differential allpass filters, which is a very efficient way of implementing a band limited 90 degree phase shift.

You can also use advanced filter design methods (least square, Parks McLellan etc.) to optimize any Hilbert approximation for your specific requirements.


Your approach seems awfully complicated. It looks like the frequencies of the "minor cycle" are much higher than the "major cycle" and the ramp so a simple lowpass filter should do the trick. If you don't know the frequencies, a simple FFT analysis can be done.

Subtracting filtered versions of signal is tricky since you need to carefully watch time alignment and phase distortion in the filter.


Second moment ergodicity of gaussian random process

A WSS Gaussian process is also a strictly stationary Gaussian process, but beyond that, Gaussianity or the lack thereof has nothing to do with the matter; the question might as well have asked for a proof for a general discrete-time stationary process that the process is second-moment ergodic if and only if the following condition holds:
$$\lim_{n\to\infty} \frac{1}{2n+1}\sum_{i=-n}^n |R_X[n]|^2 = 0 \tag{1}$$ (Note that I have changed the condition to be more symmetric).   Regardless of all this, I think that this is a question whose answer (the proof of an ergodic theorem) is way too complicated for anyone on dsp.SE to type up for the edification of the OP and any other readers who happen to stray here.

Random processes don't have to be first-moment or second-moment ergodic, and they fail to be ergodic when the autocorrelation function does not decay away fast enough as $n\to\infty$. To see this, note that the sum is one of nonnegative terms, and so if $R_X[n]$ has finite support, there is no problem at all, but if $R_X[n]$ is nonzero for all $n$, then the sum in $(1)$ is increasing with $n$, and if the increase is not slow enough, that $\frac{1}{2n+1}$ will not enough to damp the sum down enough that the limit will be $0$. $R_X[n] = \exp(-|n|)$ or $\frac{1}{1+n^2}$? No problem at all: the process is ergodic.  But consider a process in which all the random variables $X[n]$ comprising the process are the same random variable $Y$. The process is stationary, but all the sample paths are (different)  constant functions of time. If $Y=a$, the corresponding sample path has value $a$ for all $n$. Thus, while $R_X[n] = E[X[m]X[m+n]] = E[Y^2]$ is a constant, it is not necessarily the same constant as $$\frac{1}{2n+1}\sum_{i=-n}^n y[i]y[i+n] = a^2$$ and thus the process is not second-moment ergodic. It's not mean-ergodic either as noted in this answer of mine. That answer also has an example of a process that is not stationary but is nonetheless second-moment ergodic.


how can i get to know the frequency of a speech signal?

A speech signal can contain a lot more than 4kHz bandwidth, however experiments showed that if you removed all but the spectrum below 3kHz (thus allowing the use of cheaper circuitry), it was still reasonably understandable by another human.


Dimension error using findpeaks function in Matlab on Power Spectral Density to calculate SNR of audio data

Firstly, here is how you can fix what you are tryin to do:
noisepsd{k} = psdvectors{k};
noisepsd{k}(signalpeaklocs{k}) = noisepsd{k}(signalpeaklocs{k}) - signalpeaks{k};

Secondly, here is why I think it is not the greatest idea: You are finding the peaks and subtracting them from the spectrum. This is actually equivalent to setting the points where you have local peaks to zero. In effect, you will have a spectrum that is equal to zero at all points where you have local peaks and equal to signal plus noise at all other points. I fail to see why this could be a good estimate of the noise power spectral density.
Thirdly, if you ask me how to do better, well... finding the SNR of a signal is like asking how to get A and B from C = A+B. You cannot! The problem is not unique. If we have a noisy signal, how would we know what is signal and what is noise, how to even properly define "noise" here? In general you cannot, because both could be stochastic. You can only find an approximation and this approximation will only work well if you know something about your signal. For instance, if you had a pure tone in noise, the problem becomes solvable. We can estimate the tone parameters (frequency, amplitude, phase) form the signal, remove it and treat the remaining part as noise. It works in this case because we know how to differentiate signal and noise, the signal should be tone-like. However, in general audio signals need not be tone-like at all.
The only answer I can give is studying your signals closely and trying to come up with some heuristic that describes them well. Maybe your signals are bandlimited, i.e., you have some regions in the frequency range where you expect no signals. Then you could use those to estimate the noise floor. You would need to assume thought that the noise floor is constant (i.e., the noise is white).
All this will only work if you know a thing or two about your signals. If you know nothing it's back to guessing A and B from C = A+B.


Kalman filter using only location information

If you are new, you might use FilterPy and python. It is such a good start point since it has detailed explanation and manuel.


Are both halves of a symmetrical FFT result required to generate the original waveform?

You can perfectly use, $X_h[k]$, one half of the conjugate symmetric DFT $X_f[k]$ which belongs to a real signal $x[n]$ of length $N$.

Of course in doing so, you must symmetrically pad the half length DFT $X_h[k]$ to the full length DFT $X_f[k]$, which is obviously necessary for an unaliased reconstruction of $x[n]$ from its proper length DFT $X_f[k]$ of length N.

Note that during the creation of the mirrored padded full DFT $X_f[k]$ from the symetric half $X_h[k]$, the first sample $X_h[0]$ is not repeated and the remainig samples $X_h[1]$ to $X_h[M-1]$ are mirrored with their phase being negated (i.e. they are conjugated), where $M$ is the length of the half length DFT $X_h[k]$.

The length $M$ of the half sized DFT $X_h[k]$ is given by the following: If the signal's (or its natural DFT $X_f[k]$'s) length $N$ is even then $M= 1 + \frac N 2 = \frac{N+2}{2}$, otherwise, if $N$ is odd $M = 1 + \frac{N-1}{2} = \frac{N+1}{2} $.

Also when $N$ is even, the last sample of the half sized DFT, $X_h[M]$ is not repeated, and when $N$ is odd then the last sample $X_h[M]$ is conjugate mirrored in formation of the full length DFT $X_f[k]$, as examplified below Matlab/Octave code:

x = [1:8];      % N = 8 , a signal of length N = 8 (even)
X = fft(x,8);   % X[k] is the natural DFT of x[n]
Xh = X(1:5);    % M = (8+2)/2 = 5, also matlab indexing starts from k=1
Xf = [ Xh conj(Xh(end-1:-1:2))];    %N:even => Xh[0] and Xh[5] not repeated
xr = ifft(Xf,8); % reconstruct the original signal from Xf[k]
xr = 1     2     3     4     5     6     7     8

as expected.


Synthesize Digital Filter with Arbitrary Frequency-Response

For FIR filter design, common algorithms include Remez exchange or Parks-McClellan :
What is the difference between Remez exchange and Parks-McClellan filter design?

For IIR filters, one possibility is the Differential Evolution algorithm, as described in some IEEE journals.  If you have a realistic desired phase response and can guess the filter order, then Berchin's FDLS method may be suitable.