Let's find the impulse response for this system:

$$H(z)= \frac{z^{-5}}{1-0.1z^-1} = z^{-5}\frac{1}{1-0.1z^-1}$$

You didn't specify the ROC, so I'll assume that we take it so that the filter is stable. Namely:

$$|z|>\frac{1}{10}$$

I've split the expression of $H(z)$ to notice that the impulse response basically consists of an exponential with a delay. We know that

$$\frac{1}{1-0.1z^-1}\xrightarrow{\mathcal{Z}^{-1}}(0.1)^nu[n]$$

where $u[n]$ is the unitary step. The $z^{-5}$ factor translates to a delay in the time domain. So the impulse response is:

$$h[n]=(0.1)^{n-5} u[n-5]$$

If we want to find the sum of all the terms squared in the sequence:

$$\begin{align}
\sum_{n=-\infty}^{\infty}h^2[n]
&=\sum_{n=-\infty}^{\infty}\left[(0.1)^{n-5} u[n-5]\right]^2\\
&=\sum_{n=5}^{\infty}\left[(0.1)^{n-5}\right]^2\\
&=\sum_{m=0}^{\infty}(0.1)^{2m}\\
&=\sum_{m=0}^{\infty}(0.01)^{m}\\
&=\frac{1}{1-0.01}\\
&=\frac{100}{99}\\
\end{align}$$

Where I've used the formula for the infinite geometric sum.

EDIT:

Note that a smaller NRR implies a worse performance of the filter at reducing noise. In this case, NRR turned out to be almost $1$ (or $0 \ \mathrm{dB}$), so the variance of the noise at the output is practically the same as in the input. As Fat32 pointed out in the comments, this depends on the position of the pole of your filter. The noise is white (as far as I understood from your comments), so it has infinite bandwidth. If the filter's bandwidth is too large, then "a lot" of noise will pass through it, appearing at the output. On the contrary, if it is narrow, then it will "filter out" most of it, leading to a decrease in the variance at the output. In this case, the filter's bandwidth is rather wide due to the pole being near the origin, and then the NRR approaches $1$: not a really good filter to reduce noise.


Why are the bass components of music more audible than other components?

The sound wave is a pressure wave in the air, a mechanical vibration. If you have two rooms separated by a wall or door, there is little or no air passage through which the sound can propagate.

Instead, the wall itself have to absorb and re-emit the mechanical vibration on the other side. Since doors and walls are massive objects, they are resonating on lower frequencies.

The high frequency sound cannot pass through the doors since it is able to make only small particles vibrating (or just surface of the door), hence the wave is not completely absorbed and rather reflected back.

The bass sound wavelenghts are about feet long (the range is quite large depending how you define "bass"), which roughly matches the size of such large objects like walls and doors.

Of course, waves with even larger wavelengths are unable to make objects of this size resonate and hence doors behave like a band-pass filter (as Dilip Sarwate noted).


Adding AWGN noise with a correct noise power to the signal

The kind of AWGN channel you're simulating is sometimes called a "discrete time" channel. That means that you're going to simulate the bit error rate at the output of a decision device, the input to which is produced by a correlator. This is more easily explained by looking into the transmission and reception of a single bit. I'll assume BPSK.

Let's say that you transmit the physical signal $a_0p(t)$, where $a_0$ is either $A$ or $-A$, and $p(t)$ is a pulse of unit energy. To make things easier we'll assume $p(t)$ is symmetrical, that is, $p(t)=p(-t)$ (this is true for all practical systems I know). Note that, in this scenario, the bit energy is $E_b=A^2$.

Assuming that there is no noise, the received signal is $r(t)=a_0p(t)$. The receiver correlates $r(t)$ with the transmitted pulse $p(t)$. The output of the correlator is $$\int_{-\infty}^\infty a_0 p(t)p(t)\,dt=a_0,$$ where I've used the fact that the pulse has unit energy.

Now assume that the input to the correlator is pure Gaussian white noise $n(t)$ with PSD $N_0/2$ and zero mean. Then, the correlator's output is $$n=\int_{-\infty}^\infty n(t)p(t).$$ By the Wiener-Khinchine theorem, $n$ is a zero-mean Gaussian random variable with variance $\sigma_n^2=N_0/2$.

In the case of having signal plus noise at the receiver's input, the output will be $$r=a_0+n$$ because the correlator is linear. Using the value of $r$, a decision can be made on whether $A$ or $-A$ was transmitted.

We define the SNR at the correlator's output as the bit energy divided by the noise variance: $$\text{SNR}=\frac{a_0^2}{\sigma_n^2}=\frac{2E_b}{N_0}=\frac{2A^2}{N_0}.$$

Notes:

When you plot the BER, you'll likely use $E_b/N_0$ in the horizontal axis. The reason is that in quadrature modulations, the total noise power is $N_0$, not $N_0/2$.
Practical systems use a matched filter instead of a correlator. The results are exactly the same.
In a practical system there will be delay, you transmit multiple pulses instead of only one, etc. Those details are not essential to understanding how the SNR is measured in your simulation.


How does viola jones algorithm cope up with a test image size which is very different to training image sizes?

I think you are talking about multiscale object detection. It will rescale the images to different resolution and then use trained classifier over it.


Golay complementary sequences for channel estimation - MATLAB

I would try with this code:

channel =[1 0.2]';
Ga=[+1, +1, +1, -1, +1, +1, -1, +1]';
Gb=[+1, +1, +1, -1, -1, -1, +1, -1]';
y1=conv(channel,Ga);
y2=conv(channel,Gb);
correlation_output1=xcorr(y1,Ga);
correlation_output2=xcorr(y2,Gb);
channel_estimation=correlation_output1+ correlation_output2;
stem(channel_estimation./length(Gb))

Best regards,


Correcting for Doppler shifts at the speed of sound for a 800 MHz signal

frequency shift of about 890 kHz (supposing my math is correct)

$$\begin{align}
\Delta f_\text{Doppler} &= f_0 \frac vc\\
&=800\cdot 10^6\,\text{Hz}\frac{3.4\cdot10^2\frac{\text m}{\text s}}{3\cdot 10^8\frac{\text m}{\text s}}\\
&\approx 800\cdot 10^6\,\text{Hz}\cdot 1.13\cdot 10^{-6}\\
&\approx 906 \,\text{Hz}
\end{align}$$

So, your shift expectation is off by a factor of 1000.

  Supposing we don't have to worry about interference, how should I go about correcting this frequency shift? I was thinking about using an FLL or a PLL, but I'm not sure where whether these could support such important shifts.

You can design the bandwidth of your PLL or FLL to accomodate any shift. It just gets a bit tricky keeping other signals out of the observation.

Even your ~ MHz shift isn't all that much to deal with. Normal observers will simply have a  larger bandwidth than just necessary to observe your communication signal, anyways, and need

to have frequency correction anyways, since no two clocks / frequency generators are perfectly identical, physically (and if someone says "I've built this communication system. For now, I neglected synchronization needs", you know that they omitted the hard part)
When knowing what your communication signals look like, it's often not that hard to find them in the spectrum.
When you observe a communication partner for a while, you can easily write down a Doppler rate, i.e. the rate at which the Doppler shift changes. From that, you can predict future shifts and make your detection even easier (typical Kalman-Filter application)


Circular Convolution Matrix of $ {H}^{H} {H} $

If $ H $ is a matrix form of Circular Convolution then it is a Circulant Matrix.
Being a Circulant Matrix means it can be diagonalized by the Fourier Matrix $ {F} $:

$$ H = {F}^{H} D F $$

Where the matrix $ D $ id a Diagonal Matrix with the Fourier Coefficients of $ \mathcal{F} \left( h \right) $ on it main diagonal.
Also pay attention that we use the Conjugate Transpose $ {\left( \cdot \right)}^{H} $ operator as $ F $ and $ D $ are built by complex numbers (The space is over the complex numbers, hence the Adjoint Operator is the Complex Conjugate).

Now, it is easy to see that:

$$ {H}^{H} = {F}^{H} {D}^{H} F $$

Which is again an operator diagonlaized with the Fourier Matrix hence it is a Circulant Matrix -> Cicrcular Operation.
Basically it is a multiplication in Fourier Domain (DFT) by the conjugate of the DFT of the vector $ h $. By the DFT properties it means it is the flipped version of $ h $.

By the way:

$$ {H}^{H} H = {F}^{H} {D}^{H} F {F}^{H} D F = {F}^{H} {D}^{H} D F $$

Now, pay attention to $ {D}^{H} D $. It is a diagonal matrix multiplied by its conjugate. Namely it is, in the Fourier Domain, the Squared Magnitude.
In time domain it means it can be done by circular convolution with the correlation function of $ h $ (Or just do circular convolution in one direction and then in the other which is basically convolution with the flipped version).


Your $S_{new}=N-F$ requires that you know both $N$ and $F$.
But since $N=F*R$, with $*$ being convolution, and $R$ being the Room Impulse Response (RIR), you're solving the equivalent problem $S_{new} = F*R-F = F*(R-\delta)$, with $\delta$ being the delta dirac impulse or its discrete equivalent.
The hard part here is estimating $R$, which aligns nicely with your initial statement that you want to do echo detection in itself, not signal reconstruction. So, you're building a RIR estimator for itself, not for usage as part in an equalizer. That makes a lot of sense: Maybe a heuristic based on observing the MFCCs as they are, then "cleaning them up" in MFCC domain (simply nullifying things that shouldn't be there), then estimating a time-domain signal that would have these corrected MFCC signals, then using that as estimate $\hat F(t)$ for $F(t)$ and the received time-domain signal $N$ in a classical IR estimation algorithm, giving you a RIR estimate and thus an echo description.
Since convolution with the RIR is an LTI system, trying to do the estimation itself in MFCC domain feels like a bad idea – what is easy in time domain and linear frequency domain suddenly becomes hard: MFCC's are great if your system is not linear, but make working with linear systems harder!
If you insist on integrating MFCCs into your detection process, you could try to use them  as measure of speech quality, and feed that knowledge back into an adaptive time domain filter.
Generally, if I only had MFCCs, I'd try to reconstruct a time-domain signal from that, and then apply the mature theory of impulse response estimation on that – there's plenty of results that achieve optimum performance under some given conditions, and if you have one of these, whatever else you could do can only be worse.


Choice of Gaussian kernel parameters when lowpass filtering before image resampling?

I guess what you want to do is similar to image resize. If so, then multiple build-in functions in MATLAB can be used for this purpose, for example imresize.


Bjorne, I am not exactly sure what you are trying to do but I can assure you that the simplest of DWT code is not more than 10 lines in Matlab. Here is how to do it:

approx = zeros(1, cols);
details = zeros(1, cols);

%for Haar Wavelet, Use filterbanks in Matlab, google it!

hi = ...
lo = ...

for i = 1: Level

   a = conv2(signal, lo);
   d = conv2(signal, hi);

   detail() =  downsample(d);
   approx() =  downsample(a);
   signal = downsample(a);

end

That is pretty much everything you need to compute DWT. As for the values between the parenthesis in detail and approximation, just fill in the indexes.


How to prove a train of sinc pulses in digital communicaton system are orthogonal to each other?

The idea of solving the integral in the frequency domain is good, but you made a mistake rewriting the integral. Note that
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt\tag{1}$$
equals the Fourier transform of $\phi_n(t)\phi_k(t)$ evaluated at $f=0$. As you know, that Fourier transform is given by the convolution of the two individual Fourier transforms of $\phi_n(t)$ and $\phi_k(t)$, respectively:
$$\mathcal{F}\big\{\phi_n(t)\phi_k(t)\big\}=\int_{-\infty}^{\infty}\Phi_n(\xi)\Phi_k(f-\xi)d\xi\tag{2}$$
Evaluating $(2)$ at $f=0$ gives
$$\int_{-\infty}^{\infty}\phi_n(t)\phi_k(t)dt=\int_{-\infty}^{\infty}\Phi_n(\xi)\Phi_k(-\xi)d\xi=\int_{-\infty}^{\infty}\Phi_n(\xi)\Phi_k^*(\xi)d\xi\tag{3}$$
where the last equality in $(3)$ is true because $\phi_k(t)$ is real-valued. Eq. $(3)$ is just Parseval's theorem.
I'm sure you can continue from here and show that the right-hand side of $(3)$ equals zero for $n\neq k$.

Note that the integral you tried to compute equals the inverse Fourier transform of $(\Phi_n\star\Phi_k)(f)$ evaluated at $t=0$, i.e., it equals $\phi_n(0)\phi_k(0)$ which also satisfies
$$\phi_n(0)\phi_k(0)=\delta[n-k]\tag{4}$$


How to resize a Summed Area Table/Integral Image?

Calculation of the integral image is a linear transform on the original image.

Now, assuming you wanting to downscale the integral image means that you want the result to be the same whether you:

downscale the original image, integrate after
integrate the original image, downscale after

and remember that "downscaling" by a factor of $N$ the original image means

applying a low-pass filter to avoid aliasing and
throwing away $N-1$ of every $N$ pixels.

since low-pass filtering is also a linear operation, we can just switch orders of summations and down-scale the result image. Nothing special applies.


The more that sinewaves that have a zero crossing in the same direction at the same point add up in magnitude, the greater the contrast at that point.  e.g. a possible edge.  Phase congruency is another name (actually a superset) for the amount of sinewaves that line up with zero crossings at (nearly) the same point or points.

Whereas a random distribution of phases is more likely to result in something that a human would see as uniform noise rather than an edge.


Implementing Continuous Wavelet Transform

In 1D, some of the standard references are:

Continuous wavelet transform with arbitrary scales and $O({N})$ complexity, A. Muñoz and R. Ertl\'e and M. Unser, Signal Processing, 2002
A fast approximation to the continuous wavelet transform with applications, Berkner, K. and Wells, R. O., Jr., 1997, Proc. Asilomar
Fast Quasi-Continuous Wavelet Algorithms for Analysis and Synthesis of One-Dimensional Signals, Maes, S. H., 1997, SIAM J. Appl. Math.
Comparison of algorithms for the fast computation of the continuous wavelet transform, Vrhel, M. J. and Lee, C. and Unser, M. A., 1996, Proc. SPIE
Fast algorithms for discrete and continuous wavelet transforms, Rioul, O. and Duhamel, P., 1992, IEEE Trans. Inform. Theory

For a GPU implementation:

Effiziente Berechnung der FWT auf Grafikkarten, Richard Maltan, Bachelorarbeit, 2014


For causal FIR linear phase filters the delay you ask for can simply be computed from the group delay and is given by (for an odd length symmetric FIR filter of length $L=2M+1$)
$$\tau = - \frac{d\phi(\omega)}{d\omega} = - \frac{d-M\omega}{d\omega} = M $$

So it will shift every input component by $M$ samples at the output.

NOTE: you can always produce your output sample $y[0]$ as soon as the first input sample $x[0]$ arrives. The delay is related with the information processing and not the sample processing;i.e., your initial $M-1$ output samples will be called as the filter transient response, before the steady output is computed.


How to Extract a Period of a Periodic Pulsed Signal?

Natalia Molinero Mingorance, Welcome to the DSP community.

What you have is basically a shifted periodic signal.
Why? Because what you have is equivalent (Given many samples) of having a periodic signal which you started sampling at $ 0 $ and then we shit it to the real place you started sampling it.

Yea, there is the edge case, since it might not have  its end and beginning to make a complete cycle.

But let's assume our model is perfect, what would you see in Frequency Domain (DFT)?

Well, since shifting in time only changes the phase it means you won't see any difference.

So how would you extract the period of a signal you started sampling at the beginning of a period in the Frequency Domain?
Yep, just looking at the Frequency in the Frequency Domain.

If you post the files of the samples we'll be able to to see this together.


Recover Frame Compression Lossiness from Many Similar Frames

The method to improve resolution by few frames is called Multi Image Super Resolution (As opposed to Single Image Super Resolution).
For compression artifacts you can look for JPEG Deblocking algorithms.


As you pointed out, there are many state-space realizations of one particular transfer function. The reason is that a transfer function only represents the input-output behavior of a system (observable and controllable dynamics) and not the internal states.
That being said, you can directly write state-space realizations from a transfer function with the so-called vertical and horizontal companion forms (equivalently, the observable and controllable canonical forms). The problem is that the physical meaning of the states is not usually guaranteed.
So unless you have a model of your system in terms of differential equations (first principle modeling like Newton's laws or Lagrangian mechanics for mechanical systems), the physical meaning of the states will usually not be accessible. The same holds in the discrete-time domain.


Output of MATLAB's two-dimensional DWT function $\tt dwt2$

The confusion is natural as the notations and practice are misleading, and differ across the different works on 2D wavelets.

First, the order of operations ("row then column" or "column then row") is important. In most cases (to the best of my knowledge), LH means low-pass on rows, then high-pass on columns (so, giving horizontal edges). But in some contexts, people use filtering tools on "columnized signals", so the first filter works on columns, then the image matrix is transposed, the second filter pass is applied, and the result is transposed back.
Second, from an math operator point of view, LH(I) is more often understood as applying the high-pass first, then the low-pass, enhancing the confusion.
Third, one may think about the active derivative: a vertical derivation yields the horizontal features.

In the following image, the rightmost letter seems to denote the first filtering operation:

while the opposite convention is used in second image:

This might stems from a cultural and technical background: does one read images left-right then top-bottom, or  top-bottom then left-right?

The best is to test your implementation, look at the resulting images, and use a consistent notation.


Which image filter can be applied to remove gridded pattern from corrupt jpegs?

You can use a standard inpainting algorithm. These algorithms replace marked pixels in an image with the pixel values that surround these marked pixels. The challenge here is to detect the grid (my tests seem to show that it is not a completely regular grid). So, I came up with this solution:

from PIL import Image
import requests
from io import BytesIO
import cv2

url = "http://i.stack.imgur.com/Ahrnl.jpg"
response = requests.get(url)
img = Image.open(BytesIO(response.content))

plt.imshow(img)
A = np.array(img)
A2 = A.copy()
A_gray = cv2.cvtColor(A, cv2.COLOR_RGB2GRAY)

# Do some rough edge detection to find the grid
sX = cv2.Sobel(A_gray, cv2.CV_64F, 1, 0, ksize=3)
sY = cv2.Sobel(A_gray, cv2.CV_64F, 0, 1, ksize=3)
sX[sX<0] = 0
sY[sY<0] = 0

plt.subplot(221)
plt.imshow(sX)

plt.subplot(222)
plt.imshow(sY)

plt.subplot(223)
# the sum operation projects the edges to the X or Y-axis.
# The 0.2 damps the high peaks a little
eX = (sX**.2).sum(axis=0)
eX = np.roll(eX, -1) # correct for the 1-pixel offset due to Sobel filtering
plt.plot(eX)

plt.subplot(224)
eY = (sY**.2).sum(axis=1)
eY = np.roll(eY, -1)
plt.plot(eY)

mask = np.zeros(A2.shape[:2], dtype=np.uint8)
mask[eY>480,:] = 1
mask[:, eX>390] = 1

A2[mask.astype(bool),:] = 255
plt.figure()
plt.subplot(221)
plt.imshow(A)

plt.subplot(222)
plt.imshow((A2))

restored = cv2.inpaint(A, mask, 1, cv2.INPAINT_NS)
plt.subplot(223)
plt.imshow(restored)

The program output is as follows:

To detect the grid I did a quick-and-dirty solution. It can be improved a lot, but it shows the initial idea. The general flow is:

detect the grid
create a mask that describes which pixels are corrupted by the grid
inpaint the corrupted pixels.

For inpainting I used OpenCV inpaint operation. For detecting the grid, I performed edge detection in X and Y direction using a Sobel filter. Then I add all edge values in the X-direction and Y-direction to find peaks, where the grid lines are. Then, I choose the highest peaks as the coordinates where the grid lines are estimated. It's not working perfect (e.g. strong edges in the image are falsely detected as grid lines), but it shows the idea. It can be improved by e.g. Hough transformation to find lines, kicking out very strong edges etc.

Alternatively, if the grid is really the same for all images, then you can perform the grid detection jointly for all images, which would yield a much better accuracy (just do the technique above, but before choosing the peaks, sum up the results from all pictures). In more detail, you would calculate eX for all images and add all these eX together into a single vector. This vector will have a much clearer peak structure and the thresholding can be done easier.


separating stationary and non stationary parts of univariate signal

All these algorithms are based on the determination of a new basis (i.e. a new set of basis vectors), chosen according to a given constrain (to maximize the variance of each component in PCA for instance).

If you have only one signal, you have only one basis vector. So you cannot use these algorithms because you don't have enough information.

In lower dimensions, you can think of it as the calculation of the arithmetic mean from a set of noisy values corresponding to the same measurement. You need multiple measurement to be sure to get close to the "true" value of the measurement. What you are asking is how to get the "true" value from only one measurement without any repetition. That cannot be done because you don't have enough information even to do anything.

However, depending on you signal and its specific noise, there might be other ways to achieve what you want. Actually, all you have to do is to identify a way to differentiate noise from signal.

Their distributions are obviously different, but these are mere projections, so again, you have too few data to use them.

One working possibility is to filter your data to remove high frequency components (to remove most of the noise) and very low frequency components (to remove pink noise, but usually pink noise is really not a problem). But that would work only if you have an idea of the bandwidth of your signal (the maximum frequency beyond which your signal has no meaning anymore).

What you describe won't work. If I continue my analogy of the average over noisy measurement data, it is like saying that you will split the result of you unique measurement (say, 15 = 12 + 3) and use that result (12 and 3) to calculate the mean (7.5). I hope you understand that this makes no sense (the mean of 15 is at best undefined, or at worst 15). And even if it did, your decomposition would be quite arbitrary (why not use 10+5 or 20-5 or 13264585120-13264585105?). Again, you won't get more information by decomposition of your signal.

Final word:

If you really want us to help you or to give you more details about why this won't work as you expect, please open a new question, provide actual data and explain your setup.


Is the expansion of this expression correct?

Convolution in time domain means product in frequency domain. DTFT of $\delta[n]$ is just $1$, therefore any arbitrary $x[n]$ when convolved with $\delta[n]$ remains same, as in frequency domain you are just multiplying $X(e^{j\omega})$ with $1$.
So, $Y(e^{j\omega}) = 1 \ . \ X(e^{j\omega}) = X(e^{j\omega})$. Hence, $y[n] = x[n]$.
Now, apply same principle for $\delta[n-1]*\delta[n-1]$ to see that the result will be $\delta[n-2]$.


Filterpy Kalman Filter batch processing with multiple measurement sources

I am not use how you can success multiple sensors using batch filter, but, in sequential filtering you need to perform time update before each measurement. You can do it using the timetag information of each sensor measurement. You can modify the batch filter code however, I strongly recommend to write your own function depending on what you need.


Time derivative of signal - effect on noise distribution

Let's call $a[n]$ the samples of the angular measurement. Then, an estimate of the angular velocity (in angle/sec) is given by discrete differences, i.e., $v[n] = (a[n]-a[n-1])/t_0$, where $t_0$ is the sampling time.

Let's assume a signal pluse noise process so that $a[n] = a_0[n] + w[n]$. Inserting, we obtain $v[n] = v_0[n] + \bar{w}[n]$, where $\bar{w}[n] = (w[n]-w[n-1])/t_0$. Certainly, if $w[n]$ is zero mean, so is $\bar{w}[n]$ You can now compute the variance of your effective noise samples $\bar{w}[n]$, namely $$\mathbb{E}\{\bar{w}[n]^2\} = \frac{1}{t_0^2}\left(\mathbb{E}\{\bar{w}[n]^2\} + \mathbb{E}\{\bar{w}[n-1]^2\} -2 \mathbb{E}\{\bar{w}[n]\bar{w}[n-1]\}\right). $$ If your noise samples $w[n]$ are i.i.d. and have  constant variance $\sigma_w^2$, this leads to $\mathbb{E}\{\bar{w}[n]^2\} = \frac{2\sigma_w^2}{t_0^2}$. So there is a doubling of variance, but also some kind of rescaling. This depends of course on the unit of your velocity (whether it's angle/second or angle/sample, in the last case your variance is just $2\sigma_w^2$). NB: If they are correlated, this changes to $2\sigma_w^2(1-\rho)$, where $\rho$ is the correlation coefficient.

All this is true if the velocity is simply estimated from differences of adjacent samples. There are more elaborate schemes available, which would have an effect on the noise variance.


How Does the Non Local Means Filter Work?

Check out the description here.

Basically, it tries to find patterns of pixels all over the image, and average over these self-similarities rather than just using the pixels close to the current pixel.

So to answer your questions:

How is the the value of the non-local mean calculated?

Find all the pixels that are similar to the current pixel (non-local).

Find the average value of these.

How does that value replaces a pixel?

Find the non-local average as above.

Replace the current pixel with this non-local average.

How does that filter sweeps the image?

Well, provided you are creating a new image from the old one so that the mean calculation is being done on the old image and the mean assignment is being done on the new image, it shouldn't really matter.


How do we know that CTFT of the autocorrelation function is the PSD?

$$\begin{align}
R_x(\tau) &=\int_{-\infty}^\infty x(t)x^*(t-\tau)\,\mathrm dt\\
&= \int_{-\infty}^\infty \left[\int_{-\infty}^\infty X(f)e^{j2\pi ft}\,\mathrm df\right]x^*(t-\tau)\,\mathrm dt\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(t-\tau) e^{j2\pi ft}\,\mathrm dt\right]\,\mathrm df\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(t-\tau) e^{j2\pi f(t-\tau)}\,\mathrm dt\right]e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x^*(\lambda) e^{j2\pi f\lambda}\,\mathrm d\lambda\right]e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty X(f) \left[\int_{-\infty}^\infty x(\lambda) e^{-j2\pi f\lambda}\,\mathrm d\lambda\right]^*e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty X(f)X^*(f) e^{j2\pi f\tau}\,\mathrm df\\
&= \int_{-\infty}^\infty |X(f)|^2 e^{j2\pi f\tau}\,\mathrm df
\end{align}$$


Design discrete controller for zero steady state error

You can obtain the steady state error using Finite value theorem (http://www.engr.iupui.edu/~skoskie/ECE595_f05/handouts/fvt_proof.pdf) and then you can compute your controller to make it zero.
or if you have sampling fast enough, you can just discretize your continous version of controller.


How to preprocess signals prior cross-correlation?

What you are describing is called (/ has been solved via) "Dynamic Time Warping".

In DTW, you calculate a series of simple operations between two sequences that if applied on one sequence they would "bring it" directly "over" the other. Of course, there might be more than one ways to do this in which case you are selecting the one with the least cost. This "cost of distortion" is now your similarity.

For more information on DTW, please see this, this and this. More relevant to your application might be this publication. For a practical example, please see this link.

With regards to some of the rest of your questions:

  May cropping the longer signal corrupt the result of the cross-correlation?

Yes, it would...

  What is the best approach for resampling x and y?

...because both cropping and re-sampling assume a linear speed up. This might not be the case across the entire length of your signal.

  As I mentioned, x and y are only looking very similar after filtering. Is filtering signals prior cross-correlation advisable?

It depends a lot on the application. However, filtering two signals $x,y$ with the same impulse response $h$ is expected to increase their cross correlation because after convolution, both $x,y$ time instances would now have the same $h$ imprint on them.

EDIT: Following a brief review of the data provided with the question, it seems that these are some sort of counts (?) or other digital / impulsive data. Therefore, the question "...is filtering signals prior to cross-correlation advisable?" is ill posed. The fact remains that if you filter two signals through the same filter, their cross-correlation will increase. But here, it seems that the signals do not exist in their final form yet. To clarify: We have to agree on which one is the signal. If you were dealing (for example) with a pulse width modulation (PWM) derived sound signal, we would not try to correlate the two PWM pulse trains that produce it but their integrated versions because THAT is what we consider as our signal. Consequently, in this particular case we are dealing with here, we have to understand what is considered as "signal". The only thing we know about it is that it comes out of some process. That's fine but, at the same time, if you have to integrate the data to PRODUCE the signal of interest then the process itself or its context would tell you what sort of integration you have to apply. It might not even be called "low pass filtering". In conclusion: If your signal has to be integrated first, then don't try to compare the "pulse sequences". To get back to the PWM-Audio case, it is like asking "Should I integrate the PWM pulses before cross correlation or would this spoil my signal?". Of course you have to apply the integration because otherwise, you don't have the signal to work with. Decide on the parameters of integration, apply them to your original signal and then, if you still have drifts and compressions and dilations of features in the signal, see if you can use DWT to establish the time delay.

  Is it a must to detrend data prior cross-correlation or might it be enough to remove the mean?

Again, depends on the application and the amount of drift expected in the data but in general, it is a good idea to remove any long term drifts to avoid the cross-correlation being thrown off by the fact that the trends seem to match despite those little "details" on top of them. Where "details", of course, might be your desired signal. This is illustrative but a hypothetical scenario, your data might not even include such a case, it simply shows how detrending might help.

Hope this helps.

EDIT (2):

In the light of more information, I have to say that this does not look like a task for DTW neither cross-correlation the way it is described.

What is attempted here is an estimation of the delay between two waveforms by aligning "landmark" points. NOT the waveforms themselves. Personally, I am struggling to see why the marked points in the signal should align at all. Therefore, what I would suggest is this:

Pick the higher frequency waveform
Mark the peak you are after
Extract $\pm m$ points (or ms) around it. For example, if the peak is at sample 4500, extract x[4500-200:4500+200]. Call this your template $q$.
Subtract the mean of $q$ from $q$
Downsample $q$ so that it matches the sampling rate of the lower frequency signal. Call this $q_d$.

This is now your template. It is shorter than the signal you are looking for and you can "slide" it over the lower rate signal and estimate the cross correlation. At the point that the cross correlation peaks is where this second waveform is, in the lower rate signal. Prior to applying cross correlation, you might want to be subtracting the mean of the signal portion you are about to compare with your template or instead of removing the mean, you might want to be applying detrending. Because the duration of the signal is small, simple linear detrending might work good enough.

This is an attempt to align specific points of interest in the signal. Not the waveform as a whole

Hope this helps.

EDIT (3):
Ideally, a large amount of the information provided in the comments here, should appear in the original question and then perhaps we can delete all of these comments.

If you can control the load of the inverter, maybe you can insert a "sync event". Apply a known load on and off after the inverter has reached its nominal operation to create a "click" with a predictable response. You can use that "click" later to sync the two signals. If you have this capability, then you can average the responses from a number of "clicks" so that you have a good idea of what you are looking for, via cross correlation, as the "sync pulse".

Of course, this covers what happens when the two signals are supposed to have a constant time invariant delay.


A reference request for zero-order hold (ZOH)

I'd personally start with "Signals and Systems" by Oppenheim and Willsky. I have the second edition, and it has a good description of zero-order hold. This book is a must-have if you're interested in signals and systems. There is also a version that focuses on discrete signals. In addition, this book is probably available at most engineering libraries.


Efficient FFT (or SFFT) for adjacent parts of signal with overlap?

One method might be to split your length $N$ overlapping FFTs up, so that you compute the $\frac{N}{2}$-point FFT of each block of samples that come in. Then, in order to get the $N$-point FFT that you want, you can combine the last two $\frac{N}{2}$-point results with appropriate twiddle factors.

This would technically do what you want, but in most cases (especially on contemporary PC-like platforms with good FFT libraries), it's not going to be appreciably (if at all) faster than just doing the overlapping $N$-point FFTs.


Derive Time-Domain Biquad Coefficients from Z-Domain

I'm not entirely sure what you are asking but if your Z-domain representation is
$$H(z) = \frac{b_0+b_1 z^{-1}+b_2 z^{-2}}{a_0+a_1 z^{-1}+a_2 z^{-2}}$$
Then your difference equation is
$$a_0 \cdot y[n] + a_1 \cdot y[n-1] + a_2 \cdot y[n-2] = b_0 \cdot x[n] + b_1 \cdot x[n-1] + b_2 \cdot x[n-2]$$
And solving for $y[n]$ we get
$$y[n] = \frac{1}{a_0}(b_0 \cdot x[n] + b_1 \cdot x[n-1] + b_2 \cdot x[n-2] - a_1 \cdot y[n-1] - a_2 \cdot y[n-2]) $$


First of all, CNNs are typically used with identical sized inputs. So, you resize your input to a constant size, no matter what.

Next, detecting an object and recognizing an object over many samples are different tasks. For object detection, you could be better of with R-CNNs. Its sources are also available and there is a Caffe tutorial for R-CNNs. Of course there are many other works on the same topic. For example, this one.


1) How can we check errors are Gaussian and independent?

Check out this answer for a way to check for independence (whiteness).

This answer gives one way (not necessarily a good way) to check for Gaussianity: calculate the kurtosis and check that it's close to 3.

  2) If the processes are not stationary, could that be a problem?

Only if your signal model (state transition matrix, input matrix, and output matrix, process noise covariance, measurement noise covariance) do not reflect the sort of nonstationarity that your signal has.

  3) If I have the unobservable variables that are latent. How could I get the transition matrix or the covariance matrix for the equations?

How do you get the transition matrix ever? Make it up!

How do you know that there are unobservable states if you don't already have a signal model?


How to determine if an IIR filter is causal or non-causal?

First of all, it's not correct to say "poles should (always) be inside the unit circle for an LTI system to be stable" ; unless it's implied that system is also causal. Otherwise, if the system is noncausal, then its poles should be outside of unit circle for the system being stable.

For IIR systems that are described by LCCDEs causality must be externally imposed on the system (consistent wit the initial conditions), and cannot be derived from the equation alone. For example, the system

$$y[n] + a y[n-1] = b x[n]$$ cannot be concluded to be causal, as the equation can also be written as

$$ y[n] = (-1/a) y[n+1] -(b/a) x[n+1] $$
which is (apparently) indicating a noncausal dependence on the current output on both future input and future output.

Therefore this difference equation can signify both a causal and a noncausal system. Its solution should be derived by assuming a causal system or noncausal system.

This is also understood by the fact that a given transfer function $H(z)$ will have a corresponding LCCDE representation, but will have multiple ROCs. For each ROC the LCCDE should be solved accordingly, yielding a different solution for each assumption. But the LCCDE is the same.

FIR systems do not have poles (except at origin or infinity) and have no issues with stability and they don't posses regions of convergences on the z-plane. Hence their equation will be signifying the causality such as:
$$ y[n] = a x[n+1] + b x[n] + c x[n-1]$$
is a noncausal FIR system and you cannot manipulate the given equation into a causal form.

Or the following system
$$ y[n] = a x[n] + b x[n-1] + c x[n-2]$$
is a causal FIR system and you cannot manipulate the given equation into a noncausal form...

So causality shall be an assumed property for IIR systems.


Wavelet packet analysis, Sampling rate, Decomposition level

"My question is what is the best level of decomposition to see these frequencies?"

You can achieve this through separation of the three components you have into different bands. Since the inner most level is a low passed version of the signal, you should try using at least 3+1 decomposition levels.

With dyadic scaling (the decomposition method typically used with wavelets) outer decomposition level covers the range of frequencies from F/2 to F .. The next level covers F/4 to F/2.. and the next level covers F/8 to F/4....etc.

Assume a signal of two sinusoids (0.4Hz and 1 Hz) with 10Hz sampling rate you get the following plot for its frequency content decomposition into 5 scales dyadic wavelet scales. by changing the sampling ratio, the location of the peaks in wavelets subbands change


Is there an easy way to translate a Fourier transform table from angular frequency $\omega$ to Hertz $f$?

Your confusion comes from the fact that you use $X(\cdot)$ for denoting both functions, the function of $\omega$ and the function of $f$, but they are really two different functions, because

$$X(\omega)=\hat{X}(f)=\int_{-\infty}^{\infty}x(t)e^{-j2\pi ft}dt,\quad\omega=2\pi f\tag{1}$$

That's why in all the correspondences you mentioned you can just replace $\omega$ by $f$ inside the function argument (and note that the two frequency domain functions are actually different).

Note, however, that outside the function you have to replace $\omega$ by $2\pi f$, as for instance in the following correspondence:

$$\frac{d^nx(t)}{dt^n}\Longleftrightarrow (j\omega)^nX(\omega)=(j2\pi f)^n\hat{X}(f)\tag{2}$$

There is one more complication that arises due to the transform as defined by $(1)$ being unitary if $\hat{X}(f)$ is used, but being non-unitary when using $X(\omega)$. This can be seen by the factor $1/(2\pi)$ for the inverse transform when using $X(\omega)$. This implies, for instance, that you need a factor $1/(2\pi)$ in the following correspondence when using $X(\omega)$, but not when you use $\hat{X}(f)$:

$$x(t)y(t)\Longleftarrow \frac{1}{2\pi}(X\star Y)(\omega)=(\hat{X}\star\hat{Y})(f)\tag{3}$$

Also be careful when a Dirac delta impulse is involved, because clearly $\delta(\omega)\neq\delta(f)$. The correct correspondence is

$$\delta(f)=2\pi\delta(\omega)\tag{2}$$

Consequently,

$$\mathcal{F}\{1\}=2\pi\delta(\omega)=\delta(f)\tag{3}$$

and

$$\mathcal{F}\{e^{j\omega_0t}\}=\mathcal{F}\{e^{j2\pi f_0t}\}=2\pi \delta(\omega-\omega_0)=\delta(f-f_0)\tag{4}$$

from which the Fourier transforms of sine and cosine follow easily.


No, if you upsample a signal, its duration (i.e. measured in real-world units like "seconds") does not change.

If your signal $s[n]$ is sampled at a rate $f_s$ Hz, and it is $N$ samples long, then the total duration of the signal is $\frac{N}{f_s}$ seconds.
Say you upsample the signal by a factor $M$, yielding a new signal $s_u[n]$. It will be $MN$ samples long, but it will be sampled at a new rate $Mf_s$. Its duration is therefore $\frac{MN}{Mf_s} = \frac{N}{f_s}$. You can see that the signal duration is the same as $s[n]$.

It's not clear what you mean by a "shaping filter" in your third question.


As you have taken 1024 point FFT, multiply the x-axis with a factor of 2*pi/1024.

For now your peak is at sample number 200.
200*2*pi/1024 approximately equal to 0.4*pi


How can I better extract the hidden images in FEZ's soundtrack?

A few pointers to help you produce more detailed results:
First, you may want to read about a phenomenon called Pareidolia. By tweaking and tuning parameters, you may be simply creating results through coincidence.
There is no such thing as the spectogram.
There is such a thing as a spectogram.
I've looked at some of the results already found that have been posted online, but none actually talk about the specifics of the spectogram used. This makes those results hard to reproduce (outside of running the exact same program with default settings). Why?
Spectograms have many, many parameters. Some explanation below:

Windowing function
First, there is the windowing function, which is technically any function taken from a function space (this stuff gets very complicated). Typically it's taken to be continuous or at least semi-continuous (a finite number of discontinuities).
Theoretically, this kind of space is of the order aleph-3, although practically we're limited by the resolution of the computer, or, the total resolution is N floating-point values big. Typically these have a 53-bit mantissa and are between zero and one, so we're already looking at 53KB of parameters for a 1024 FFT size spectogram.
Window size
This is typically taken to be a multiple of the FFT size, usually W = N. Ideally, at least one full FFT size is inside the window size. Complicated mathematical expressions can identify the resolution of a spectogram, which depends on this parameter.
Window overlap
In order to improve time resolution somewhat, when using non-rectangular windows windows are usually chosen to overlap to reduce spectral leakage losses. High window overlap factors tend to smear out/deform the image. Typically chosen to be 0.5 with most standard windows as that satisfies COLA.
FFT size
Changing the FFT resolution can change the shape and image of the spectogram through rounding effects and also has an effect on the time-resolution capability of the spectogram. You could compare it to the Heisenberg principle: increasing the sample count improves frequency resolution, especially at low frequencies, but reduces time resolution.
Must be a power of two for the traditional FFT algorithm to run quickly enough.
Has an impact on the speed of the visualization. Older programs tend to run with lower FFT sizes to enable real-time visualizations on slower computers. FFTs scale with O(N log N). See here.
Sampling frequency
May not match the recording frequency, in which case the audio is resampled before analyzing. Iff so, then the resampling method used can also impact the result.
Base frequency
Typically not directly configurable (typically the user is asked to select min/max frequency shown in the resulting image), the base frequency of the FFT is equal to fs/N, or the sampling frequency
Filters, transformations, etc.
Post- and pre-processing is typically used to improve the output quality for visualizers. Typically done to remove noise, compression artifacts and the like, they can corrupt hidden images, which are often of very low volume and may be caught by a noise filter.
Image width in px
Finally, the chosen image width and, iff different from the size of the FFT output, the resize filter can make an image fuzzier in quality. Ideally the image width is equal to the FFT output, though, if this is too big, the FFT output is an integer multiple of the image width. Take care to convert the image if needed with an integer ratio as well to reduce artifacting, and to use a lossless format.


Apply Low pass Butterworth filter in Python

You should not be using the analog filter - use a digital filter instead. You want the filter to be defined in Z-domain, not S-domain. Also, you should define the time vector with known sampling frequency to avoid any confusion.

The design of the digital filter requires cut-off frequency to be normalized by fs/2.

Here is a working example:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

fs = 1000  # Sampling frequency
# Generate the time vector properly
t = np.arange(1000) / fs
signala = np.sin(2*np.pi*100*t) # with frequency of 100
plt.plot(t, signala, label='a')

signalb = np.sin(2*np.pi*20*t) # frequency 20
plt.plot(t, signalb, label='b')

signalc = signala + signalb
plt.plot(t, signalc, label='c')

fc = 30  # Cut-off frequency of the filter
w = fc / (fs / 2) # Normalize the frequency
b, a = signal.butter(5, w, 'low')
output = signal.filtfilt(b, a, signalc)
plt.plot(t, output, label='filtered')
plt.legend()
plt.show()


Background subtraction for object tracking

The seminal paper on background modeling is Adaptive background mixture models for real-time tracking by Stauffer and Grimson.

If you have Matlab and the Computer Vision System Toolbox, then you can try using vision.ForegroundDetector, which implements a version of their algorithm.


inverse continuous wavelet transform and [Parm] in cwtft

Why won't cmorwavf work?

http://www.mathworks.co.uk/help/wavelet/ref/cmorwavf.html

Could you then construct a matrix of (complex) Morlet wavelet co-efficients (suitably dilated and translated) and simply invert it for the inverse wavelet transform?


Is a wavelet-based correlation measure worth any additional computational overhead?

First off, you should use whichever tool is appropriate for the job.  Correlation vs coherence vs wavelet-based correlation are all different things, so this question is kind of like asking "Which is better? Screwdrivers or hammers?"  It depends on what you're trying to do, and whether you care about similarity in time, frequency spectra, or both.

Second, I have only a minimal understanding of wavelets, but your assumption that wavelets require more computation might be wrong.  The Fast Fourier Transform takes $O(n\log n)$ operations, while the Fast Wavelet Transform takes $O(n)$.  So the wavelet method might actually require less computation, depending on whether you can use the information you get out of it.

Empirically, producing n outputs from n real inputs, the multi-level wavelet transform in PyWavelets becomes faster than NumPy's FFT when n is greater than about 4096.

However

It's Python, and the two implementations could be very differently efficient.  I don't even know if wavedec() would be considered a Fast Wavelet Transform.  They use the abbreviation DWT in their documentation.  Are Haar DWT and FWT the same thing?
The time varies depending on the wavelet used.  Meyer wavelet takes 6 times as long as Daubechies to produce the same amount of data.
I still don't understand how the FWT tiles the time-frequency plane, or if producing n outputs is sufficient to get the same kind of similarity measurement as an n-point circular cross-correlation using FFTs. (Technically it's a time-scale plane, not time-frequency, but I think they're the same for the complex Morlet wavelet?)  FWT is a "critical sampling" of the plane, and produces the same amount of data as FFT, so it seems fair to compare them.

The main point is that computation time is at least roughly similar for both, so I don't think you should worry about it when deciding which to use.


How does the OFDM system receiver decide the OFDM symbol length?

First of all, the timing, phase and frequency offsets of the channel are removed. After synchronization, the cyclic prefix is removed, to avoid ISI between OFDM symbols. Now because the transmitter and receiver agree upon a common standard before the communication begins, the receiver knows how many samples to consider as a symbol.

Generally, the the number of samples in an OFDM symbol depend upon the bandwidth.


Effect of aliasing when sampling at higher rate

The name of this effect is Spectral Leakage. Remember the relation
$$\frac{k}{N} = \frac{F}{F_S}$$
where $k$ is the bin number, $N$ is the FFT size, $F$ is the continuous frequency in Hz and $F_S$ is the sample rate in Hz. It can be seen that $k$ varies from $-N/2$ to $N/2-1$ (or from $0$ to $N-1$). So there are only $N$ continuous frequencies $F$ for which you will get a single impulse in discrete frequency domain.

So in the first figure, the frequency of the signal contains an exact integer number of periods within the FFT duration and hence the correlation of other DFT sinusoids with the signal is zero. In the second figure, the signal span is not equal to the integer number of periods that violates the orthogonality requirement. Hence, the correlation of different DFT sinusoids appears as leakage around the true frequency.


It's certainly calculating the right thing.  Though instead of

sum(x.*(cos(1000*2*pi*t)-i*sin(2*pi*1000*t)))*2/N;

you might try

sum(x.*exp(-i*2*pi*1000*t))*2/N;

If you need to do something similar, but in-line (not in a batch), you might want to look at the Goertzel algorithm. As the Wikipedia link says:

  .. provides a means for efficient evaluation of individual terms of the Discrete Fourier Transform

And there's nothing particularly "Fast" about this way of doing it, so it's really just calculating one DFT bin.


The total roundoff error for the sum of $N$ numbers is:

$$
S = \sum_{i=0}^{N-1} E_i
$$

The roundoff error for the $i$-th number is represented by the random variable $E_i$. If we assume that the random number generator used by the computer yields numbers $X_i$ taken from a uniform distribution, then the difference between each $X_i$ and the nearest tenth (which is the roundoff error $E_i$) is uniformly distributed on the interval $(-\frac{0.1}{2}, \frac{0.1}{2}) = (-0.05, 0.05)$.

What we're concerned with, though, is the distribution of $S$. Since $S$ is the sum of $N$ independent, identically distributed (iid) random variables, then via the central limit theorem, as $N \to \infty$, $S$ will tend to a Gaussian distribution. If we assume that your case of $N=1000$ is "large enough" for the Gaussian assumption to hold, we can easily estimate the probability that you seek. It's certainly possible to exactly calculate the distribution of $S$, but the Gaussian assumption is likely close enough for most applications with such large $N$.

A Gaussian distribution is characterized by its first two moments, so if we can find those for $S$, then we have all the information we need. These are easy to calculate for a sum of iid random variables. The mean of $S$ is equal to:

$$
\mathbb{E}(S) = \sum_{i=0}^{N-1} \mathbb{E}(E_i) = 0
$$

The variance of $S$ is equal to:

$$
\mathbb{E}\left((S - \mathbb{E}(S))^2\right) = \sum_{i=0}^{N-1} \mathbb{E}\left((E_i - \mathbb{E}(E_i))^2\right)
$$

Recall that the random variables $E_i$ are distributed uniformly. It is well known that the uniform distribution over the interval $(a,b)$ has variance $\frac{1}{12}(b-a)^2$. For this case, that yields a variance $\sigma_{E_i}^2 = \frac{0.01}{12}$. Therefore, the variance of the total roundoff error $S$ is $\sigma_{S}^2 = \frac{0.01N}{12}$.

So in summary, we can approximate $S$'s distribution as Gaussian with mean zero and variance $\sigma_{S}^2 = \frac{0.01N}{12}$. Based on those parameters, you can easily calculate the estimated probability distribution function (pdf), then integrate that result to arrive at whatever probability you seek. The probability that there is a total roundoff error with magnitude greater than one would be:

$$
\begin{align}
P(|S| > 1) &= P(S>1 \lor S < -1) \\
&= 1 - P(-1 < S < 1) \\
&= 1 - \int_{-1}^{1}f_S(s)ds
\end{align}
$$

where $f_S(s)$ is the Gaussian distribution's pdf that we arrived at before.


Entropy when having symbols with same probability

HINT: Compute the entropy $H$ of your source and note that in your case it represents the  number of bits required to represent each symbol. Then from
$\quad H$ [bits/symbol] $\times$ $x$ [symbols/second] = $5$ [bits/second]
compute $x$.


One issue is that, for the first two sampling frequencies, you are not sampling a full wave of the signal.

Below is a plot of the time domain signal (1024 samples) and the resulting FFT.  For the first two, the DC bin (0 index) will capture all the energy.  For the 1kHz example, the first bin (the one after 0 index) will capture it.  Bin 10 is where the peak is for 100 Hz sampling rate.

Perhaps your sampling isn't happening quite the way you expect?

R Code Below

# 43563

N <- 1024

fs <- c(20000, 10000, 1000, 100)
res <- c(0,0,0,0)
idx <- 1

par(mfrow=c(4,2))
for (freq in fs)
{
  res[idx] <- freq / N
  idx <- idx + 1
  t <- seq(0,N-1) / freq
  x <- sin(2*pi*t + 0.892340923)
  plot(x, type='l')
  title(paste('time domain for ', freq))
  x_fft <- abs(fft(x))
  plot(x_fft[1:100], type='l')
  title(paste('100 bins in frequency domain for ', freq))
}


Image processing free open source libraries

For Fourier transforms go get fftw.  I don't know what meshgrids are, so I can't answer that.  For dynamic arrays: that's what std::vector is, it's built into the language.

For a wider variety of image transformations you might look at GEGL and the library underlying Imagemagick.

There's also this question about filtering software: Use Butterworth and Chebychev filters.


Modified Tukey window applied in frequency domain

Regarding the 2:1 frequency ratio recommendation for the transition band (f2:f1 or f4:f3):

The steeper the filter response (i.e. the smaller the frequency ratio between pass band and stop band limits), will necessarily produce a much longer impulse response and produce ringing artifacts in the time domain whenever there is a transient.

These filters, which are applied in the frequency domain, also seem to have linear phase response, which is very desirable for analysing this type of signal.

Butterworth filters have a compact and nice representation as difference equations (for filtering in the time domain), but do not have linear phase response, and as you see from your plots, have not very good transition between pass band and stop band.

(I don't understand how you plotted a symmetric impulse response (time domain) for Butterworth filters; their impulse response is causal and asymmetric in general.)


What criteria should be used for smoothing

You want to look at the innovation sequence; the difference between the filtered value (the estimate or forecast) and the measured value. If you have a perfect filter that removes the noise, the innovation sequence should be white noise (Dirac impulse autocorrelation).

https://en.wikipedia.org/wiki/Innovation_(signal_processing)

Here is some MATLAB code that illustrate these properties for some white sequences.

N = 1e5;
n1 = -1 + 2*rand(1,N); % uniform noise
n2 = n1 + -1 + 2*rand(1,N); % triangular noise
n3 = randn(1,N); % Gaussian noise

[c1,l1] = xcorr(n1); [c2,l2] = xcorr(n2); [c3,l3] = xcorr(n3);
figure(1), plot(l1,c1,l2,c2,l3,c3) # autocorrelation (Dirac delta)
figure(2), clf, hold on
pwelch(n1,kaiser(N/1000,10),N/2000) # power spectral density (constant)
pwelch(n2,kaiser(N/1000,10),N/2000)
pwelch(n3,kaiser(N/1000,10),N/2000)


What is the difference between closed form and opened form in z transform

One possible difference occurs for a pole outside of unit circle. Note that the original system is a truncated IIR filter. As such, it's an FIR filter and have no stability issue for any finite $a$.

However, the closed form representation includes a pole, and implements a pole zero cancellation. Eventhough theoretically (with infinite precision arithmetic) this would not pose a problem, for finite length arithmetic and for large values of $N$, and for $|a| > 1$, you could end up in numerical issues.

Specifcally, the instantaneous dynamic range (the mantissa length of the floating point representation) would be large enough to be able to accurately represent the arithmetic difference indicated in the numerator. Indeed this could also be aproblem with the finite length implementation when $a^N$ is too large to be added to $ax[n-1]$, and that could happen when $|a|>1$. So which one would be more sensitive to this numerical issue? You better test by yourself. Furthermore, the closed-from version could impose an additional problem for those $ a \approx 1$ which would not be the case for the FIR filter.

As an advantage of te closed form, note that the FIR form and IIR form requires the same amount of memory (to hold past values of input up to a delay of $N$ samples), however the closed form (truncated IIR) requires much less computation MACS compared to FIR version for large N.


Finding delay parameters for watermarked signal

Thank you for posting the additional information, it certainly helped.  I'll try to do a prose explanation with less "math" in it, and then point you in what is the right direction.

That's certainly a long winded problem statement, and I believe the solution is actually simpler that your professor may be making it out to be.  From reading the problem statement, a watermark has been made and was supplied to you.  This watermark is just sampled Gaussian noise, such that it lacks coherent frequency content, i.e. the spectral energy of the watermark is not located at any specific frequency (as opposed to say, a constant frequency signal like a sinusoid).

So we've got our original signal, we have our watermark signal, now it's time to actually watermark.  The simplest approach would just be adding the watermark signal to the original signal.  Your professor has gone a step ahead and added 5 copies of the signal.  Each one of these copies has its own time delay, as well as its own amplitude.  So, after we  add in our 5 watermarks to our signal, we have the signal Xm

Since there are 5 copies of the signal at 5 unique points in time in Xm, that means if we did a cross-correlation with our watermark signal and Xm, we would expect to see 5 distinct peaks at different amplitudes. This operation is typically called a matched filter, and you can think of it as using a finding the "locations" (in your professor's words, the delay parameters; think of this as when each watermark "starts" in Wm) of some given template (i.e. the watermark) in your recorded signal.

Your professor has gone ahead and detailed out how you can accomplish cross-correlation in the frequency domain.  This is often the preferred method, since it is typically much more computationally efficient.  Your professor has also detailed common ways to improve cross-correlation detection, specifically the use of windowing functions as well as the use of frequency normalization.

In turns out  those methods he described will be very important for your success.  The watermark signals are pretty low gain (naturally), so we'll need to maximize the signal to noise ratio (SNR) as much as possible.

I think that since this is seemingly an assignment of some sort, I'm not willing to just give you the answer in numeric or code form, but here's a pictorial output of what I got:

You can see that I've marked the 5 points corresponding to the time delayed version of the signal with red X marks.  These  peaks are easily 5x larger than any other point in the output of the cross-correlation, so one can easily  conclude that these are statistically significant, and point to the location of the delayed watermark signals.

Here's some code I used.  Feel free to play around with it and see if you get satisfactory results (note: "data" is your Xm signal, and "wm" is your watermark signal:

% Determine number of samples for the correlation lag vector
nsamps = numel(data);

% define a window function. Taylor windows give great performance
win    = taylorwin(nsamps,7,-10);

% window each signal
xin    = data;
wmw    = win .* wm;

% make the signals complex for phase coherency
xin    = hilbert(data);
wmw    = hilbert(wmw);

% perform cross-correlation in the frequency domain
xf     = ( fft((xin)) .* conj(fft((wmw))) );
% transform out of the frequency domain while normalizing
x      = ifft(xf ./ (abs(xf)+1));

lags   = linspace(-nsamps/2, nsamps/2, nsamps);

figure(1); hold off;
plot(lags,fftshift(abs(x)))
% find the 5 peaks, then plot them
[pks,locs] = findpeaks(fftshift(abs(x)),'NPeaks',5,'MinPeakHeight',0.0075);
hold on;
plot(lags(locs),pks,'rx');
xlabel('Correlation Lags'); ylabel('Scaled Correlation Magnitude (dB)')
grid on;
title('Watermark Cross-Correlation');

So hopefully that gets you going in the right direction.  Like I mentioned before, the watermark signals are pretty low gain, so make sure you do that normalization.  The windowing function is arguably not necessary, but they're pretty standard in practice since  they help your cross-correlation "behave" a bit better.


A signal having a mean-value or DC component of zero is commonly referred to as mean-free or as having no DC component. It does not mean that it cannot be averaged, just that the average comes out as zero.

Might be a little inexact but it is very common.


The term $x(2n)$ is not time invariant: it gives the even samples of $x$. Shift the input by 1 sample and you'll get the odd terms, which are different from the even terms shifted 1 place.

Your reasoning above is not correct. You should compare $T(x(n-k))$ to $T(x(n))|_{n-k}$ (that is, the original output, shifted by $k$ samples).

Edit:
These type of exercises are a common source of confusion because of the following abuse of notation: $x(n)$ is not a signal, but the nth sample of the signal, and having $n$ inside of the expressions is formally wrong. Even Oppenheim's book makes this mistake when working with LTI systems and impulse response.

The correct way for expressing transforms between signals is as follows:

Let $x$ and $y$ be signals; $T\{\cdot\}$ be the transformation (system) under scrutiny; $R_k\{\cdot\}$ be a delay of $k$ samples.

So, in general you would write the relation between input and output signals as $y = T\{ x \}$, where the whole signals are involved. Using $n$ inside these expressions collapses the signal to a single sample and this is wrong. And having the same variable $n$ on both sides of the equation only leads to confusion.

Then, a system is time-invariant if $R_k\{T\{ \cdot \}\} = T\{ R_k\{ \cdot \}\}$ for all integer $k$.

This is why I gave the example above with even-odd samples, in order to avoid the usual confusion with $n$ inside the expressions.

Now, going back to the system $T$ that evaluates to $y[n] = x[2n]$ for each sample $n$:

Let $x_k = R_k\{x\}$ and $y_k = R_k\{y\} = R_k\{T\{x\}\}$.

Evaluating all these signals at sample $n$: $x_k[n] = x[n-k]$; $y_k[n] = x[2n - 2k]$.

Now, $T\{x_k\}$ at sample $n$ is $x_k[2n] = x[2n - k]$ which is different from $y_k$.


Using under-sampled data, what determinations can be made about the Nyquist rate (correct sampling rate)?

This may be a trivial answer, but you can of course use the aliasing information to determine the correct sampling rate, given that you know that the signal is undersampled. In the simplest case, given that the signal appears to be at $k f_s,\ k < 0.5$ one can devise that the frequency of the signal is $(1-k)f_s$. In general it may be at $(N \pm k)f_s$ for some integer $N$.


How to generalise the Fourier transform?

The Fourier Transform is just one of so many different transforms that alter the representation of (usually) a time-series from the time domain, to another domain (usually a frequency domain but other representations exist for other transforms such as time / frequency, time / scale and others).

You can find much more information about transforms in general from this Wikipedia list of articles that lists some popular and often used transforms. (You might want to focus on the Discrete and Integral Transforms at first)

Alternatively, you can check out this recent discussion on how the Wavelet transform, achieves a decomposition similar to that of the Fourier transform.

Finally, when you have the luxury to have acquired simultaneously many different time series from the same phenomenon you can even employ techniques such as Principal Component Analysis (PCA) and Independent Component Analysis (ICA) which go to the point of transforming a signal to a sum of elementary waveforms that are actually extracted from the signal itself (rather than being pre-set as it is done in the Fourier (and related transforms) or Wavelets).


The root locus is a parametric plot. It shows how zeros and poles move in the complex plain depending on the value of a gain $k$. For all points on the plot both the phase condition1 and the amplitude condition2 have to be true.

$$G_0(s) = k \cdot  \frac{\prod^q_{i=1}(s-s_{0i})}{\prod^n_{i=1}(s-s_{i})}$$

In your case, $s_{01} = s_{02} = 2$, $s_1 = 0$ and $s_2 = s_3 = 20$

There are a few rules to construct the root locus by hand, after you added the locations of the zeros and poles to the plot:

symmetry

Everything is symmetric to the real axis.

start and end points

Every line/curve/branch (whatever you want to call it) of the plot starts in a pole and ends in a zero or infinity.

You have $q=2$ as the number of zeros and $n=3$ the number of poles. That means that 1 branch will go to infinity and two others will connect the poles with the zeros.

asymptotes

There are $n-q = 1$ asymptotes that the branches are converging towards for very large gains. The asymptotes have an angle of

$$\phi = \frac{\pi + l\cdot 2\pi}{n-q}$$

and all asymptotes meet in the point

$$s=\frac{\sum^n_{i=1}s_i - \sum^q_{i=1}s_{0i}}{n-q}$$

in your case $s=16$ which is admittedly only of real help if you have more asymptotes.

With $l=0, 1, ..., n- q-1$. For the single asymptote of your system, that means $\phi=\pi$.

From this rule you know that the one pole that goes to infinity will do so by getting close to the asymptote at angle $\pi$ which means it goes to the left of the plot.

real axis

This is one of those rules where you think it might as well be part of a horoscope (at least I did), because it sounds as arbitrary and random as how much the constellation of planets at your birth has an influence on your life.

Points of the real axis can be part of the plot. Those points that are part of the plot have an uneven number of poles and zeros to their right.

Wait, what? To estimate the results for your plot, start from $+\infty$ on the real axis and go to $-\infty$.

The rightmost "thing" in your plot are the double poles at 20. No point of the real axis to the right of those can have any number of poles or zeros to its right, so we conclude that the real axis to the right of 20 is not part of the plot. This means that the 2 branches coming out of the poles at 20 cannot be on the real axis.

Going further to the left, with the 2 poles to the right and then the 2 poles and 2 zeros to the right (4 in total) you know that no point from 0 to $+\infty$ of the real axis is part of the plot.
Note that this also means that the branches ending in the zeros will not arrive there on the real axis, because the real axis around 2 is not part of the plot.

To the left of 0, you now have all your 5 zeros and poles to the right, which means every point to the left of 0 on the real axis is part of the plot.

real axis break in/away points

The branches split or join at points where the following condition is true:

$$
\begin{align}
\sum^q_{i=1}\frac{1}{s-s_{0i}} &= \sum^n_{i=1}\frac{1}{s-s_{i}}\\
2 \cdot \frac{1}{s-2} &= 2 \cdot \frac{1}{s-20} + \frac{1}{s}\\
\frac{2}{s-2} &= \frac{3s - 20}{s(s-20)}\\
2s(s-20) &= (3s - 20)(s-2)\\
2s^2-40s&=3s^2-26s+40\\
0&=s^2+14s+40
\end{align}
$$

that gives $s_1 = -4$ and $s_2 = -10$. At those two points, there's a split or branch.

Now which one's a split and which one's a join?

Remember that this is a parametric plot and the branches grow as the value of $k$ grows. You start with what's essentially a split at 20 (the double pole as starting point), which means the you first have to encounter a join to later split again in order to arrive at 2 (the double zero). We can conclude that the point with a lower value for $k$ is the join and the other is the split.

To calculate the value of $k$ at any given point, the amplitude condition2 can be used:

$$
\begin{align}
k_{-4} &= \frac{|(-4-20)^2(-4)|}{|(-4-2)^2|}\\
&=\frac{|(-24)^2(-4)|}{|(-6)^2|}\\
&=\frac{24^2 4}{6^2}\\
\end{align}
$$

$$
\begin{align}
k_{-10} &= \frac{|(-10-20)^2(-10)|}{|(-10-2)^2|}\\
&= \frac{|(-30)^2(-10)|}{|(-12)^2|}\\
&= \frac{30^2 10}{12^2}
\end{align}
$$

Ok, which one is bigger?

$$
\begin{align}
k_{-4} &\overset{?}{\lessgtr}k_{-10}\\
\frac{24^2 4}{6^2} &\lessgtr \frac{30^2 10}{12^2}\\
\frac{24^2 4}{6^2} \frac{2^2}{2^2}&\lessgtr \frac{30^2 10}{12^2}\\
\frac{24^2 4 2^2}{12^2} &\lessgtr \frac{30^2 10}{12^2}\\
24^2 4~ 2^2&\lessgtr30^2 10\\
3^2 2^{10}&\lessgtr 5^3 3^2 2^3\\
2^7&\lessgtr 5^3\\
128 &> 125\\
k_{-4} &> k_{-10}\\
\end{align}
$$

$k_{-4}$ is the split and $k_{-10}$ is the join.

Of course you could check that more quickly with a calculator, but you said dungeons & dragons style.

1 that is
$$\sum^q_{i=1}\Phi_{0i} - \sum^n_{i=1}\Phi_{i} = (2l+1)\pi$$

with $l = 0, 1, 2, ...$ and $\Phi_x$ being the angle from a horizontal line drawn to the right of the respective pole or zero to the point in question. This is just FYI.

2 that is
$$\frac{|\prod^q_{i=1}(s-s_{0i})|}{|\prod^n_{i=1}(s-s_{i})|} = \frac{1}{|k|}$$


Signal deconvolution from exponential and linear tail

You probably can't. If I interpret your description correctly, this is the unit step response of your system. This looks very much like the unit step response of a high pass filter, that means the system is filtering out exactly the frequency that you want (which is DC or 0 Hz).
You may be better off creating a parametric model of your system response and try to do a weighted least square error fit of the data.


How does the world's simplest sigma-delta DAC work?

This is just a first-order sigma delta DAC with the input as an unsigned binary number. The classical variant with a signed digital number is shown below where the MSB as the sign bit would subtract (Delta) from the input. In this signed variant, consider the input being a scaled digital number ranging from +/- Full Scale, and the feedback would either be +Full Scale or -Full Scale; in this case the accumulator can count up or down with a ramp rate based on the difference between the actual input value and +/-Full Scale with the output toggling when the accumulator output crosses midrange.
In comparison for the OP's case, the input is unsigned and the accumulator grows only in the positive direction with a ramp rate based on the difference between the input and +Full Scale. When the accumulator overflows, it is reset to 0 plus the overflow, which is essentially subtracting Full Scale (which is the Delta in this case).

This has all the drawbacks of any first order Sigma Delta: lower noise supression than a higher order architecture (-9 dB/octave) requiring a higher oversampling ratio to achieve a desired equivalent number of bits, and suffers from pattern noise leading to higher spurious content if the input is not changing.


Spherical equivalent of Nyquist frequency

I'm a little confused about what you are asking and how it relates to the set of spherical basis functions.  Those are the same one that electron orbitals are based on, right?  I have never understood them to be constrained to the surface of a sphere.  Then again, that is not my usual stomping grounds.

Nevertheless, I think I can generalize your question and provide an answer.

$\phi$ is a poor choice for function name in this context.  Yes, I know it is the standard for "potential function" in many contexts, but within spherical coordinates it has a different meaning.  I will use $lon$ and $lat$ for your point location, and retain $\phi$ for the value so there is no confusion.

You have a large collection of points that fit the model:

$$ \phi_n = \phi(lat_n,lon_n) $$

You would like to find $\phi$ as a linear combination of pre-defined basis functions.

$$ \phi(lat,lon) = \sum_{k=0}^{K-1} c_k B_k(lat,lon) $$

Presumably, $B_0(lat,lon) = 1$, but that is not necessary, but I wanted to throw in the "DC" term to help orient folks in relation to a DFT.

Anyway, for each $k$, you have a computable function definition of $B_k(lon,lat)$.

This is now a standard best fit problem.  You want to find the values of the $c_k$s for the best fit of $\phi$ to your data.

Set up the following arrays:

$$
B =
\begin{bmatrix}
B_0(lat_0,lon_0) & B_1(lat_0,lon_0) & B_2(lat_0,lon_0) & ...  & B_{K-1}(lat_0,lon_0) \\
B_0(lat_1,lon_1) & B_1(lat_1,lon_1) & B_2(lat_1,lon_1) & ...  & B_{K-1}(lat_1,lon_1) \\
: & : & : & ... & : \\
B_0(lat_{N-1},lon_{N-1}) & B_1(lat_{N-1},lon_{N-1}) & B_2(lat_{N-1},lon_{N-1}) & ...  & B_{K-1}(lat_{N-1},lon_{N-1}) \\
\end{bmatrix}
$$

$$
C =
\begin{bmatrix}
c_0 \\
c_1 \\
c_2 \\
: \\
c_{K-1} \\
\end{bmatrix}
$$

$$
P =
\begin{bmatrix}
\phi_0 \\
\phi_1 \\
: \\
\phi_{N-1} \\
\end{bmatrix}
$$

$B$ (NxK) and $P$ (Nx1) are known from your data points and basis function evaluations.  $C$ (Kx1) contains the unknowns you are solving for.

Ideally,

$$ BC = P $$

If $N=K$, simply multiply both sides by $B^{-1}$ and you have a solution.

$$ C = B^{-1}P $$

Generally, you want more points than that, but that is the minimum requirement unless you want to be underfitted.

If $N>K$, there is a standard approach to find the best fit $C$:

$$ B^* B C = B^* P $$

$$ C = (B^* B)^{-1} B^* P $$

This is basically an interpolation technique.  Interestingly, on the domain of a sphere surface, interpolation vs extrapolation seems to be meaningless.  You will want your points to be as spread out as possible as to not to leave any large areas of your surface unsampled.  The size of $B^* B$ is KxK, and its inversion is probably your biggest computational load depending on how intensive your basis functions are to evaluate at the surface points.

In a conventional DFT, $N=K$ and $B^{-1}$ is $\frac{1}{N}B^*$ so the inverse is implicitly solved.

The concept of a "Nyquist frequency" doesn't have a direct correlation in this scenario in the general case.  It may in terms of your choice of $B_k$ functions, but it would still be irrelevant because in your situation you have randomly distributed points instead of regularly spaced ones.

Your ultimate answer is then you are going to want to have at least as many sample points as you have basis functions, but are likely going to want significantly more.  The number of points you select will not impact the size of the inverse that needs to be taken.

Hope this helps.

This is in response to Fat32's comment.

Looking at the reference given by the OP, it seems the $B_k$s are defined by equations (18) through (33), though it looks like some others are listed as well.  Of course, the list can be extended.  Since the OP's points are randomly distributed, I don't see any other feasible method than the one I provided.  Note that the first one is indeed a "DC" function.

It does raise the question though if there is a set of definable and spread out locations for which matrix $B^{-1}$ is a multiple of $B^*$ when $N=K$, or $ B^* B $ is a multiple of the identity matrix when $N>K$ .  (Calling Olli) Just because continuous functions are orthogonal does not necessarily mean that discrete sampled versions of them are.  With the trig functions (the DFT) they are, with Legendre polynomials, for instance, they aren't. [Correction:  Being diagonal is sufficient for a trivial inverse.]

Regardless, for a fixed set of sample locations $(B^* B)^{-1} B^*$ can always be pre-calculated.  Then finding the coefficients is merely a simple matrix multiplication, just like with the DFT.  This is not what the OP has.  [Correction:  It is possible to select a set of points with $N \ge K$ so the matrix is not invertible.  With a random selection this is not likely, but the columns of $B$ may end up not being very orhogonal which would lead to poor numerical behavior.  Thanks for the paper reference, Olli.  That looks spot on and the OP should definitely read it.]

The concept of a "Nyquist frequency" is strictly dependent on the spacing of evenly spaced sampling locations on a 1D scale, and completely independent of the number of points in the frame.  We tend to think of it as "the last valid basis function before we get into aliases" (assuming the upper half is considered negative frequencies and are used first).  The OP's situation is similar to unevenly sampled points of a signal, so I really don't see the Nyquist concept as being meaningfully translatable.

I interpreted his question as being more of a the necessary number of points for a good answer and just speculation that it was analogous to the Nyquist.

I do think that it is instructive to understand that the DFT is indeed an example of a best fit solution just like this.


Start by getting rid of the background (Mathematica example):

  map = Import["http://i.stack.imgur.com/HYtbo.jpg"]

  fg = ImageAdjust[map, {1.0, 1.0}]

Then use a graph extraction tool like NEFI. If you want to do it all in Mathematica, that's possible too. You might want to add another preprocessing step to decompose the graph by edge type if you want to treat them separately, like you said.


Cross correlation, when two signals are "uncorrelated"?

Cross correlation is a measure of similarity between two signals, where one signal is allowed to be time-shifted. In this sense, the correlation is not a single number, but a function of the time shift. We say, "these two signal have a certain correlation $R(\Delta)$ for a time shift $\Delta$".

Intuitively, two signals that tend to have the same sign (both positive or both negative) for a given time shift $t$ are similar, and will have large correlation (positive or negative). Signals that are uncorrelated are just as likely to have opposing signs, and then the integral will be small.

As an example, consider two signals that are rectangular pulses of duration $T=1$, but one starting at $t=0$ and the other at $t=10$. Their product is zero (and thus their correlation for time shift 0 is also 0). However, intuitively we could say that they are similar, except for the delay. Thus, the correlation allows us to quantify that intuition: these two pulses are completely uncorrelated for some delays; however, when if the first is time-delayed by 10 seconds, then they become highly correlated.

Note that $R(\Delta)=0$ is a sufficient and necessary condition for the two signals to be "uncorrelated" (meaning they have zero correlation) with time shift $\Delta$.


what is the frequency selective channel and doubly-selective channel

Example. For a channel that can be modeled by LTI system, we send an impulse-like signal and receive channel impulse reponse $h(\tau)$ where $\tau$ is delay.

Frequency-flat channel means $h(\tau) \sim \delta(0)$ or for discrete-time version, $h[n] \sim \delta_0$.

Frequency-selective channel means $h(\tau) \neq \delta(0)$ or for $\tau > 0,h(\tau) \neq 0$.

Now if we do the measure at time $t=t_1$ and $t=t_2$, $h(\tau,t_1) = h(\tau,t_2)$, the channel is called non time-selective (anyone remember the name?).

If we do the measure at time $t=t_1$ and $t=t_2$, $h(\tau,t_1) \neq h(\tau,t_2)$, the channel is time-selective.

Doubly selective is simply both time-selective and frequency-selective.


Let us  stick to traditional functions, considering the almost everywhere equivalence class.  From an actual signal perspective, I am leaving distributions aside.

According to the Paley-Wiener theorem, functions with compact frequency support (corresponding to Bernstein spaces) can be extended
from the real line to the complex plane $\mathbb{C}$ as entire functions
of exponential type; entire functions are holomorphic at all finite points of $\mathbb{C}$ ,  exponential type means that $|f(z)|\le C e^{a|z|}$. They can be called $B$-functions.

So they inherit from all their properties: no bounded support in time, infinite differentiability, Taylor expansion everywhere, controlled decay at infinity...


Why is the carrier suppressed by AM-DSB-SC modulation?

It's a simple trigonometric identity:

$$\cos(\omega_mt)\cos(\omega_ct)=\frac12\left[\cos((\omega_m+\omega_c)t)+\cos((\omega_m-\omega_c))\right]\tag{1}$$

Multiplying two sinusoids of different frequencies results in the sum of two sinusoids with the sum and the difference of the two frequencies.

So this type of AM (DSB-SC) really results in a suppressed carrier, i.e., the modulated signal has no component at the carrier frequency. This is true as long as the message signal has no DC component.

What happens is that the spectrum of the message signal is simply shifted to $\pm \omega_c$, so DC gets mapped to the carrier frequency $\omega_c$, and as long as there is no DC component in the message there is also no carrier component in the modulated signal.

EDIT: Your square wave example is perfect. A square wave is a piecewise constant function, so you could say it is just switched DC, or, paraphrasing your question, it's "made up of sections of otherwise unmodified DC". However, if the duration and amplitude of the positive and the negative part of the square wave are equal, then the corresponding signal's DC component is exactly zero. If you accept that, you must also accept that the signal in your example has no component at the carrier frequency.


Why is the green colour not represented as black in Cr component of YCbCr?

Found the answer with more experimentation: the opposite colour to red in Cr is not green but cyan, RGB (0,255,255), and it does appear black in the Cr channel. This makes sense, as yellow, RGB (255,255,0), is opposite blue.


If your square wave has a mean of zero (this is important!), then a simple accumulator can do the job. Its operation is described by

$$y[n]=x[n]+y[n-1]$$

where $x[n]$ is the input (square wave) and $y[n]$ is the output (triangular wave).

This is a simple Matlab/Octave script showing how it works:

sq = [1,-1,1,-1,1,-1,1,-1]'*ones(1,5);
sq = sq'(:);                 % a few periods of a square wave
tr = filter(1,[1,-1],sq);    % filtered by accumulator
% plot
subplot(2,1,1),stem(sq)
subplot(2,1,2),stem(tr)


Usage of Digital Cosine Transform in lossy compression of images

A rectangular window of data can have a discontinuity between the front and back edge.  An FFT of that window represents that circular edge discontinuity with energy in a bunch of high frequency bins.  Take that segment of signal and mirror it around one edge.  Note that a circular or periodic extension of this no longer has a discontinuity across the edge.  The FFT of that doubled-by-mirroring signal might thus have less information energy in the high bins.  Thus there may be more of these high bins that can be "thrown away" by a compression algorithm without major damage (information loss).

A DCT is an FFT of the doubled-by-mirroring product.  Since the mirrored intermediate vector is symmetric, some efficiency can be gotten by throwing away any computations only required for the sine components since they will be zero anyway.

A DST would be the FFT of a doubled-by-mirroring-the-inverse signal, which would make the discontinuities worse, which would spread even more information/energy into the bins you would really like to throw away, thus making things far worse for simply but lossy compression.


Overlap/Chunk processing possible with Zero-phase?

If you apply a zero-phase IIR (for instance, by filtering both forwards and backwards) to a segment of a signal that has been zero-padded on both sides, and the two zero-padding lengths are equal or greater than the length of the IIR filter's impulse response (until that response drops below your desired noise floor), then you can sum non-overlapping segments, but you have to keep and include the zero-padded tails in the sum (similar to overlap add).  The error will be related to the (infinite length) ends of the IIR tails that have been cut off by finite length zero-padding.


Solving a Linear Mean Square Estimation the Easy Way

So in your case doesn't the relation $x = n+y$ help ?

I mean, assuming your derivation for the mean square estimtor is right, then to compute $E\{xn\}$ you would look for $E\{ (y+n)n\}$ and using properties of $x$ and $n$ you would get

$$E\{xn\} = E\{(y+n)n\} = E\{yn\} + E\{n^2\} = 0.5 + 1 = 1.5 $$


Where is the mother wavelet defined in the Fast Wavelet Transform?

TL;DR: the wavelet appears at the end of the synthesis filter bank, iterated infinitely.
Theoretically founded, practical and fast DSP tools are often derived from continuous theory: think about how the DFT is derived from the continuous Fourier transform, by discretizing both in time (like the Discrete-time Fourier transform) and frequency  (like Fourier series), which are dual variables.
For a function $\psi(t)$ to be a wavelet, very few conditions (admissibility) are required.  However, discretizing the complementary shift and scale parameters $(a,b)$ of $\psi\left(\frac{t-b}{a}\right)$ is more complicated if one want to remain invertible. And even more complicated to obtain a discrete orthogonal basis.
In other words: apart from a handful of well-known cases (Haar, Shannon, Meyer wavelets), if one have a closed-form formula of a nice admissible continuous wavelet, it is VERY unlikely that it can be discretized as a discrete orthogonal wavelet.
The filter bank structure does the job the other way around. If we assume an orthogonal multiresolution framework, we see that an iterated bank of carefully-designed filters can produce a wavelet analysis. Which wavelet? Theoretically, from a two-scale equation, which generally has no closed-form solution.
You can obtain a good approximation of the wavelet shape with the following procedure. Pick a level $L$. In the deepest detail subband (near the approximation), put a one, and zero elsewhere. Then, do the inverse wavelet transform. The highest the level, the better the approximation.Here is an example with a very q&d code.

dwtmode per
nSample = 1024;
data = zeros(nSample,1);
waveletName  = 'db3';
for iLevel = 2:7
[C,L]=wavedec(data,iLevel,waveletName);
C(L(1)+(L(2)+L(3))/2) = 1;
% C((L(1))/2) = 1;
waveletAtLevel = waverec(C,L,waveletName);
subplot(2,3,iLevel-1)
plot(waveletAtLevel(find(waveletAtLevel)));axis tight;grid on ;
xlabel(['At level ',num2str(iLevel)])
end


Get the lag of a signal with respect to the other with cross correlation

The offset is actually 4 samples given the code generated since MATLAB starts indexing at sample 1.

The accumulation from the DC offset (the non-zero average value) creates a larger triangular correlation pattern and together with the random component can sometimes overcome what would otherwise be the strongest correlation at an offset of 4. The solution is to remove any average level prior to doing the correlation.

Given the uniform distribution between 1 and 100, the average will be close to 50 (and random itself so will vary from run to run, but the mean of the random process is 50.5). When both sequences align, the correlation due to the offset alone would on average be $20(\bar x ^2)$, where $\bar x$ is the average value. This result would be close to 50,000 in this case. Since fewer and fewer of the samples are included in the correlation process as the lag increases, this value will linearly decrease to 0 with each shift of the samples.  Due to the random element each sample will vary from this triangular shape. Similarly the correlation of the sequence itself when aligned would on average be $15(\bar x^2)$. By removing the mean we eliminate the DC offset effect and allow the repeated sequence to dominate the correlation result when those sequences are aligned.

This demonstrated in the plots below plotting the 20 correlation runs individually of the OP's code to show the variability, using

[c,lags] = xcorr(a,b);
plot(lags,c)

and

[c,lags] = xcorr(a-mean(a),b-mean(b));
plot(lags,c)

The other option derived from the cross-correlation theorem is to compute a circular cross-correlation using the FFT given as:

$$Corr = \text{ifft}(\text{fft}(a)\text{fft}^*(b))$$

Since it is a circular correlation, the correlation from the DC offset is the same throughout the correlation and not strongest at zero lag as in the linear correlation so need not be removed in this case.

Here is the general MATLAB code for providing the circular correlation results with FFTs given two sequences of equal length with $N$ samples each. Since the inverse FFT will provide a result starting at time =0 and is circular in time (such that the last sample is also the result for the first negative sample), the fftshift command is used to center time = 0 in the middle of the plot similar to the plots above.

# lags is different if odd or even length
if (mod(N,2)); lags = -(N-1)/2: (N-1)/2; else; lags = -N/2:N/2-1; end
c = fftshift(ifft(fft(a).*conj(fft(b))));

As RBJ points out in the comments, this results in time aliasing, as the correlation is effectively done by wrapping around to the beginning rather than correlating to zeros at the ends of the sequence as in the linear correlation that is performed using xcorr.  The sequences can be zero-padded (append 20 zeroes to each sequence) prior to taking the FFT to replicate the linear correlation result, but in this case the mean value would need to be removed for the same reasons given above. Given the correlated samples are random noise I actually prefer the circular correlation personally without zero-padding and removing the mean as it provides a constant noise floor in the result (rather than tapering to zero as you can see in the first two plots above).


Is spectral centroid enough to measure brightness?

'Brightness' is a subjective descriptor of an audio signal. It can't be calculated directly from physical properties of the signal. A panel of subjects can give you their ratings for brightness of one signal relative to another.

That said, the spectral centroid has long been known to be correlated with subjective measures of brightness. So while you are not calculating brightness itself you should be calculating a fairly good indicator.


Deconvolving known impulse response from a sampled noisy signal

This very, very much sounds like the classical problem for which you can derive that the matched filter is the optimal receive filter to maximize SNR.

You say you know the shape of the signal; I interpret this like your system looks like this

TX impulse $p$ -> Pulse shaping filter $g$ -> 500 Hz low pass $h$ -> +Noise $z$ -> Sampling

That implies that the RX signal (given we sample sufficiently fast) is the discrete

$$r[n] = p*g*h + z$$
with $*$ being the convolution operator.

Since you know the order and bandwidth of $h$, I'd argue you actually know what $h$ is precisely. Thus, you can "pre-compute" the combination

$$m=g*h$$

The receive filter that maximizes the SNR is the matched filter – sadly, the matched filter to an IIR is not realizable (when does the time-inverse of something without end start?), you can't directly work with your actual system model.

However, if you can find a FIR filter that approximate your $h$ well enough, you should be able to  bound the error, and still get a reasonably good amplitude/power estimate by calculating the time-inverse conjugate (${}^\dagger$ operator)
$${\hat m}^\dagger =\left(g*\hat h\right)^\dagger$$
 based on the approximative filter $\hat h$. You'd then have
 $$
\begin{align}
\hat r &= p*g*h*\hat m^\dagger &+ z*\hat m^\dagger\\
&= p*m*\hat m^\dagger  &+ z*\hat m^\dagger
\end{align}$$

Since the noise process $z$ should not be correlated to the channel/filter, the second term should become small for sufficient observation.


How to determine brown noise variance from samples?

If the white noise has a Gaussian probability distribution, then it should be possible to find the variance of both noise sources using a power spectral density estimate.

$S_{1}(j\omega) = {\sigma_{1}}^{2}$ and $S_{2}(j\omega) = \frac{{\sigma_{2}}^{2}}{\omega^{2}}$, hence you should be able to minimize
$$
J(\sigma_{1},\sigma_{2}) = \| S_{1}(j\omega) + S_{2}(j\omega) - \hat{S}(j\omega) \|
$$
where $\hat{S}(j\omega)$ is the PSD estimate of your output, in order to find $\hat{\sigma}_{1},\hat{\sigma}_{2}$.

The MATLAB script below shows how this might work in practice.

Ts = 0.1; % sampling time
Nts = 1e6; % number of samples
n = randn(2,Nts); % two Gaussian noise sequences
t = (0:length(n)-1)*Ts; % time vector

sd = 1e3.*rand(1,2); % random standard deviations for noise signals

% filter for the noise
G = sd.*[tf(1,[1 0]), tf(1)];
% zero-order hold
Gd = c2d(G,Ts,'zoh');
% time-response for G(s) by the two Gaussian white noise sources
y = lsim(Gd,n,t);
y = y(1e2:end); % remove initial transient (only want stationary response)

% compute a PSD (of stationary response)
[Pxx,f] = pwelch(y,blackmanharris(length(y)/100),[],[],1/Ts,'psd','onesided');

% scaling using logarithm provides better numerical stability when
% minimizing the function
Jfun = @(p)(log(p(1)*1./(2*pi*f(2:end)).^2 + p(2)) - log(Pxx(2:end)));

p0 = 1e3.*rand(1,2); % random initial values for curve fit
p = lsqnonlin(Jfun,p0); % minimize the cost-function Jfun

sdh = sqrt(p/(2*Ts)); % estimate of standard deviations
est_error = norm(sd - sdh);

disp('Ideal parameters')
disp(sd)
disp('Estimated parameters')
disp(sdh)
disp('Estimation error')
disp(est_error)

fa = 10*log10(p(1)*1./(2*pi*f(2:end)).^2 + p(2));
figure(2), semilogx(f,10*log10(Pxx),f(2:end),fa), grid on
xlabel('Frequency'), ylabel('Power'), legend('PSD','Fit')


What should the time-shift be when implementing a continuous wavelet transform on a computer?

It's indeed wasteful for some signals, but unfortunately not much to do there except downsample afterwards. A fully rectangular representation enables reconstruction and analysis operations that aren't otherwise possible - further, doing 1 shift per sample can be much faster and take less (compute-) memory than e.g. skipping 1, via FFT convolution.
I provide a naive implementation with visuals here, also comparing against Python library implementations. Further, just today I've pushed out a ssqueezepy release, which has the best open-source CWT implementation I know of (I've yet to see if it beats MATLAB's) - can inspect its code.
If you insist on discretizing per-scale, I can only suggest a trial & error approach, as results will vary based on data, wavelet, and wavelet parameters used. You can then plot as a single 2D image by duplicating adjacent time-axis values for higher scales (but that's much effort and guesswork for little or no gain in speed).


Average Power Spectral Density of PAM signals

Question 1:

If we have a PAM sequence like this:

$a[k]=[0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0]$,

then the PAM signal (Waveform) must hold frequency components which cover these fast transitions from 0 to 1 to 0 to 1. The oscillation speed (frequency) of the transitions are basically given by the time $T$ between the PAM symbols in the waveform.

If $G(f)$ has lower bandwidth than $1/T$, these fast transitions are not possible. All the frequency components which could oscillate at that speed are filtered out by $G(f)$.

Here a plot of the above Waveform, its PAM sequence $a[k]$ and the single pulses (when $a[k]=1$) which summed up create the waveform.

And another plot with the above Waveform and its highest frequency component. If the dashed cosine would oscillate any slower, it could not represent the Waveform anymore.

Question 2:

The minimum bandwidth of $G(f)$ is given by $1/T$, but it can be higher and almost always is (see root raised cosine filter and roll-off factor, which let's you control the bandwidth of $G(f)$). BUT, the PSD of a PAM sequence is only defined within $1/T$, and due to the Nyquist–Shannon sampling theorem it is repeated thereafter (aliasing).

Why is this important?

Look at the following plot, it shows the PSD of a PAM sequence (limited by $1/T$), the PSD of a PAM signal and the filter $|G(f)|^2$.

$G(f)$ filters the PAM sequence and the result is the PAM signal, but what happens outside $1/T$? $G(f)$ filters the aliased part of the PAM sequence.

Note: Not sure if the y-axis of the PSD is correct, but we are only interested in the bandwidth.

Question 3:

You are correct!

But there are some tricks (which e.g. are exploited by the (root) raised cosine filter), see my first plot and the single pulses. They obviously overlap in time, but at each symbol a[k] only one pulse is responsibly for all the amplitude of the waveform. All the other pulses related to other symbols are zero. See the Nyquist ISI criterion.


You are correct. In your example, $k$ would refer to the "area underneath" the impulse. Mathematically speaking, the Dirac delta isn't a function in the typical sense of the term. Instead, it is more of a distribution, characterized by the fact that when integrated across any interval that contains $t=0$, the result is unity. That is:

$$
\int_{-\epsilon}^{\epsilon} \delta(t) dt = 1 \ \forall\ \epsilon > 0
$$

The distribution is typically defined as follows:

$$
\delta(t) = \begin{cases} \infty,\ t = 0 \\ 0, \text{ otherwise}\end{cases}
$$

Multiplication by a constant $k$ obviously wouldn't change this definition at all. So instead, the scaling merely changes the area underneath the distribution:

$$
\int_{-\epsilon}^{\epsilon} k\delta(t) dt = k \ \forall\ \epsilon > 0
$$

Extending this concept of multiplication by a constant $k$ to multiplication by a function $f(t)$ yields the following:

$$
\int_{-\epsilon}^{\epsilon} f(t) \delta(t) dt = f(0) \ \forall\ \epsilon > 0
$$

or more generally:

$$
\int_{T-\epsilon}^{T+\epsilon} f(t) \delta(t-T) dt = f(T) \ \forall\ \epsilon > 0
$$

which is known as the sifting property of the Dirac delta and is used extensively in linear systems theory.