Why do DFT frequency buckets need to be divided by sample period?

I don't know why you got downvoted.  This is actually a very pertinent question having to do with normalization conventions.

The most conventional DFT formula is this:

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

The $1/N$ normalized form is

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

And the $1/\sqrt{N}$ normalized form is

$$ X[k] = \frac{1}{\sqrt{N}}\sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n k } $$

My very strong preference is to use the $1/N$ normalized version.  This makes the bin values work as you expected them to.  It also makes the expected bin values at whole integer frequencies bin count invariant, and leakage values comparable.

If you do not apply a $1/N$ in the forward DFT, you will need it on the inverse to get your original signal back.  If you put $1/N$ on the DFT, then the inverse needs no rescaling (and the bin values serve as coefficients for continuous Fourier series).

The $1/\sqrt{N}$ is the "proper mathematical" one as it makes the forward and inverse DFT unitary matrix multiplications.  Yeah, nobody does that.

Bottom line, I think the $1/N$ normalization should be the standard convention and I write all my blog articles using that definition.

BTW, you will get a magnitude of $1/2$ for pure whole integer tones on $1/N$ normalized DFTs.  The other half is on the other side of the circle.


Explanation for the quantization matlab function

sigp = (sig_in - sig_nmax)/Delta + 1/2;   % convert into 1/2 to L+1/2 range
qindex = round(sigp);                     % round to 1,2,3..... L levels
qindex = min(qindex,L);                   % Eliminate L+1 as rare possibility
q_out = q_level(qindex);                  % use index vector to generate output
clear
end

sigp is a mid-tread uniform quantizer which is defined using the quantization step size Delta whose details you can find here: https://en.wikipedia.org/wiki/Quantization_(signal_processing)#Mid-riser_and_mid-tread_uniform_quantizers

sigp uses sig_in - sig_nmax instead of sig_in so that the signal can lie between 0 and 2 instead of -1 and 1 since indexing in MATLAB is always greater than 0

1/2 is added to round the first index to 1.

qindex = min(qindex,L);

The highest level can be L so if qindex goes over L, it is truncated to L

q_out = q_level(qindex);

q_out gives the quantization level corresponding to the qindex.


Removal of blocking artifacts in images due to compression

If you talk about JPEG-like blocking artifacts, an interesting appproach consists in the re-applications of the same JPEG algorithm on shifted versions of the image, see Enhancement of JPEG-Compressed Images by
Re-application of JPEG, by Aria Nosratinia, 1999. It has a software implementation. It implement a sort of shift-invariant processing. He did the same with the JPEG-2000 coder, in Post-Processing of JPEG-2000 Images to Remove Compression Artifacts. You can find an overview of these artifacts here.

Improved methods include sparsifying priors like total variation, see for instance Adapted Total Variation for Artifact Free Decompression of JPEG Images. Recent references include Reducing Artifacts in JPEG Decompression Via a Learned Dictionary and Understanding, Optimising, and Extending Data Compression with Anisotropic Diffusion.


Image processing books Gonzalez et al. vs Sonka et al?

I am not sure there are "best" image processing or computer vision books to buy. The topics covered can be very wide, and some can be better on some aspect (morphology, segmentation, denoising, etc.).

Some sites recommend a hanful of such books, and link to their draft versions, as in 8 Books for Getting Started With Computer Vision. Indeed, a lot of authors have a draft on a personal page. This is the case for instance for Computer Vision: Algorithms and Applications by Richard Szeliski.

CVonline: Vision Related Books including Online Books and Book Support Sites also points to classics, like:

D. H. Ballard, C. M. Brown; Computer Vision, Prentice-Hall Inc New Jersey, 1982, ISBN 0-13-165316-4.
A. Blake, M. Isard; Active contours, Springer, London, 1998, ISBN 3540762175.


FFT, padding, IFFT and plot in time domain

user3406207: I'm not sure what you're trying to do because I don't know the meaning of your words "wrong" and "right."  But we can say, if you perform time-domain zero padding on an original unpadded $x1[n]$ input sequence you'll produce some longer-length $x2[n]$ sequence. The $X2[m]$ fft of your $x2[n]$ sequence will be an interpolated version of (more freq-domain samples than) the $X1[m]$ fft of your $x1[n]$ sequence.  And, of course, the inverse fft of $X2[m]$ will be $x2[n]$ and the inverse FFT of $X1[m]$ will be $x1[n]$. I hope that helps. (By the way, you were smart to plot just the real part of your inverse fft.)


DFT of an audio signal stored as a multi-dimensional array

It looks like you are reading a stereo file, so you have a left and a right channel. You can do a DFT on the left or a DFT of the right or the mono or all of the above. It really depends on your specific application and requirements.


Why is anisotropic diffusion useful in image processing?

The anisotropic diffusion algorithm by Perona and Malik is the pioneering work in partial derivatives equations (PDE)-based denoising.

It applies the law of diffusion on pixel intensities to smooth textures in an image. A threshold function is used to prevent diffusion to happen across edges, and therefore it preserves edges in the image. (Unlike for instance gaussian blur filter.) This makes it very interesting if you want to remove noise, but do not want to smooth out the edges of your image, for instance if you want to use these edges to segment the image, without being perturbed by the noise.

Numerous efforts have been made that build on it, improving or extending it.

Now for where it is used, I only have a limited culture. I can cite two

Image analysis in the Life-Science field (where I work): the images that you can get out of a microscope are extremely noisy, and most of the time, it is by construction. Automated image analysis of these data often involve segmentation, for which you sometimes PDE-based algos.
Video games! For instance, try playing Mass Effect (the first one at least).


Using spectral coherence for computing similarity between time series data

One way to check for similarity between time series is to use the cross-correlation:
$$
\rho_{xy}(\tau) = \frac{1}{\sigma_x \sigma_y} E[ (x[n] - \mu_x)(y[n] - \mu_y) ]
$$
where $\mu_{x,y}$ is the mean of $x$ or $y$ and $\sigma_{x,y}$ is the standard deviation of $x$ or $y$.

This gives us a value between -1 (completely anti-correlated) and +1 (completely correlated).

Sometimes, this doesn't work as well as we'd like because there is some filtering operation relating the two time series, so while there is a strong relationship between the two time series, the cross-correlation is not very close to 1.

As a result, the spectral coherence is sometimes used:
$$
C_{xy}(f) = \frac{|G_{xy}(f)|^2}{G_{xx}(f)G_{yy}(f)}
$$
where $G_{km}$ is the cross-spectral density between $k$ and $m$.

This is now a function of frequency, and the closer it is to 1 at a given frequency, the closer the two signals are related.

I've put together a simple example in R which illustrates the difference between cross-correlation and (spectral) coherence.

The top graph shows the two time series.  The middle graph shows the cross-correlation is less than 1.  The bottom graph shows the spectral coherence as a function of frequency.

You can see that the time series are not as similar in the higher frequencies as they are at low frequency.

Sometimes, the integral of these values is used over a limited bandwidth so that a single figure of similarity is arrived at.

R Code Below

  #15466

  x1 <- rnorm(1000)

  h <- c(1,3,2,1)/7
  x2 <- filter(x1,h, circular=TRUE) + rnorm(1000)/10

  par(mfrow=c(3,1))

  plot(x1,type="l",col="blue")
  lines(x2, col="green",lwd=5)
  title("Original and filtered noise signals")

  ccf(x1,x2)

  Ixy <- spectrum(cbind(x1,x2), plot=FALSE,spans=c(5,7))

  plot(Ixy$coh)h)


How does projecting a few IR dots capture an image three dimensionally?

The projected differaction pattern (formed out of Gaussians) has a special assymetric structure, which in turn would help the software resolve the stereo correspondence problem. After correspondence is obtained, the distance can be triangulated, as in a simple stereo setting.

For further info, please check the PrimeSense patent:
http://www.google.com/patents/US20090185274

Other stereo reconstruction techniques might use color coded patterns, blue light patterns, phase shifting projections and so on. Generally, tools projecting fewer dots enjoy muliple pattern projection and encode every pattern differently. This might for example obtained by rotating the projector or programming it to project different light stripes.

Depending on the object to be reconstructed, using multiple cameras might just be enough as the RGB image information might provide you enough data to disambiguate. In this case, you wouldn't even have to project a pattern.

The technique is generally called structured light and for further info you might check:
http://fofi.pagesperso-orange.fr/Downloads/Fofi_EI2004.pdf
http://eia.udg.es/~qsalvi/Tutorial_Coded_Light_Projection_Techniques_archivos/v3_document.html
http://en.wikipedia.org/wiki/Structured_light


Building a spectral envelop of FFT'd audio

Two techniques commonly used for recovering the spectral envelope of an audio signal:

Low-order AR modeling. Estimate your AR coefficients, then plot the magnitude of the frequency response of the estimated all-pole filter.
Cepstral analysis. Take the logarithm of the magnitude of the FFT output. Take the FFT, truncate to keep only the first coefficients, take the inverse FFT, take the exponential. This is equivalent to applying a low-pass filter on the spectrum (with a log-scale magnitude).


Wavelets: How does the admissibility condition imply this $\Psi(0)=0$

It's obvious that Fourier transform of the mother wavelet is zero at the origin, otherwise the integral, whose argument is nonnegative, would be unbounded. Just because both the numerator and the denominator are zero, it does not mean it does not have a limit. So we know that $\Psi(0)=0$, but this is merely necessary, not sufficient for the integral to be bounded.

Consider the function $f: \omega \to |\sin(\omega)|^2/|\omega|$. How does this function's integral from $(-T,T)$ behave as $T \to \infty$. I hope that gives you some ideas.


My best attempt at this so far is kind of inspired by the difference of Gaussians. I'm putting it as an answer, even though I'm not totally happy with it.

Basically I made a blurred version of the image using two different Gaussian kernels. One captures the fine detail (small kernel) and the other captures the broad detail (large kernel). These two were then merged. This new blurred image is then subtracted from the original.

UPDATE: Made the process iterative, added contrast equalisation and changed blur divide to blur subtract.

output_image = get_input_image()

for i in xrange(4):
    # Magic numbers
    sigma_broad = 301
    sigma_fine  = 5
    merge_ratio = 0.5

    # Get blur
    broad_detail = cv2.GaussianBlur(output_image, (sigma_broad, sigma_broad), 0)
    fine_detail  = cv2.GaussianBlur(output_image, (sigma_fine, sigma_fine), 0)
    blurred = merge_ratio*broad_detail + (1-merge_ratio)*fine_detail

    # Subtract blur from original (and correct for average drop in intensity)
    output_image = (output_image + np.mean(blurred)) - blurred

# Contrast equalisation
clahe = cv2.createCLAHE(clipLimit=1.0, tileGridSize=(4,4))
output_image = clahe.apply(output_image.astype(np.uint8))

Result:

The reasons that I'm not totally satisfied with this result:

The information I want is altered. The objects in the ball (e.g. the dark patches) are disproportionally reduced.
There are some artefacts.
There are lots of magic numbers and parameters to fiddle with. Ideally the algorithm should be robust for a variety of balls.

ANOTHER UPDATE:
I had a go at a Laplacian Pyramid filter, which I adapted from here. I got decent results, posting mainly for interest. This strategy allows for targeting frequency bands, but it's very fiddly.


A common practice in your case is to give the mean rank of the correct answer (the ground truth), eventually after querying for only the 5 best matches. An alternative is to give the percentage of tests where the good solution was in the top 5.

Search for related papers in your field, this will give you some hints on what to choose.


I have done very little work using 2-dimensional FFTs.  Most of my work deals with speech, so the short-term Fourier transform applied to a 1-dimensional speech signal is what I understand.

I do not know of a "canned" function, particularly in Matlab, that would realize the time-varying "temporal Fourier domain" characterization that you require.  However, I would suspect you could apply a 2-dimensional FFT to a set of time-slices extracted from the original video sequence that you are working with.  The time resolution of when the slices are extracted would depend on the nature of the events captured by the video.  The magnitude and phase response of each 2-dimensional FFT computation could then be determined and examined to see how the Fourier domain response changes across time.

From what I understand, viewing images in the Fourier domain is not readily intuitive.  I believe the Fourier domain is used to realize efficient image filtering.

I hope this helps.


Continuous time double exponential filtering in state space form?

I think I figured it out, dumb mistake in writing the system in continuous time.

$\dot{l} = r_1(u - l)$

$\dot{s} = r_2(\dot{l} - s)$

$\ \ = r_2(r_1(u - l)  + r_1 s - s)$

$\ \ = r_2(r_1(u - l)  + (r_1 - 1)s)$

$\ \ = r_2 r_1(u - l)  + r_2(r_1 - 1)s$

So $$
A = \begin{bmatrix}
-r_1 & 0\\
-r_1r_2 & -r_1r_2\\
\end{bmatrix}
$$

which is stable for all positive $r_1, r_2$!


It's great that you are interested in signal processing at that early stage of your educational path.

The best path to get there is to read some introduction books on the topic. There are a lot of good and free online resources to get you started. [Note to the esteemed editor: good introduction books may a really good topic for a "sticky"].
I sometimes use

http://www.dsprelated.com/dspbooks
http://www.dspguide.com/pdfbook

One of the most important mathematical concepts you will need to get your arms around are “complex” numbers”. It’s clearly a misnomer since it’s really not that complicated and clearly makes nearly all engineering math a lot simpler. Another great free resource for all things math related is http://www.khanacademy.org and in this case specifically
http://www.khanacademy.org/video/complex-numbers--part-1?topic=core-algebra

Back to your first question: There are actually four different flavors of the Fourier Transform: Fourier Series (most likely to show up in high school), Fourier Transform, Discrete Fourier Transform and Discrete Fourier Series. All of them of them use a combination of both sine and cosine (or a complex exponential, which is essentially the same thing). You will need both.

Let’s say you calculate the sine and cosine Fourier coefficients of an input sine wave. (Under certain conditions) you’ll find that all Fourier coefficients will be zero except for one cosine and one sine coefficient. However, depending on the phase of the input sine wave you these two numbers will move around. You may get [0.707 0.707], or [1 0], or [0 -1], or [-0.866 0.5] etc. You will see that the sum of the squares of those two numbers will always be 1, but the actual values depend on the phase of the input sine wave.

If you want to deep dive, try this:
http://www.dsprelated.com/dspbooks/mdft/


How to derive the $\log$ of fraction of two Gaussian distribution

Your final condition is correct, but it doesn't show explicitly for which values of $x$ it is satisfied. This can be worked out as follows. From what you have it follows that $p_1(x)\ge p_2(x)$ is equivalent to

$$0\ge x^2(\sigma_2^2-\sigma_1^2)-2x(\mu_1\sigma_2^2-\mu_2\sigma_1^2)+\mu_1^2\sigma_2^2-\mu_2^2\sigma_1^2-2\sigma_1^2\sigma_2^2\log\left(\frac{\sigma_2}{\sigma_1}\right)\tag{1}$$

Let's first consider the case $\sigma_1=\sigma_2$. From $(1)$ it follows that in this case the condition on $x$ becomes

$$x\ge \frac{\mu_1^2-\mu_2^2}{2(\mu_1-\mu_2)}
= \frac{\mu_1+\mu_2}{2}\tag{2}$$

where we assume $\mu_1\neq\mu_2$ because otherwise $p_1(x)=p_2(x)$, which makes no sense.

In the case $\sigma_1\neq\sigma_2$, we have to consider $(1)$ as a quadratic inequality in $x$, which is satisfied with equality for the following values of $x$:

$$\begin{align}x_1&=\frac{a}{2}-\sqrt{\frac{a^2}{4}-b}\\
x_2&=\frac{a}{2}+\sqrt{\frac{a^2}{4}-b}\end{align}\tag{3}$$

with

$$\begin{align}a&=\frac{2(\mu_1\sigma_2^2-\mu_2\sigma_1^2)}{\sigma_2^2-\sigma_1^2}\\
b&=\frac{\mu_1^2\sigma_2^2-\mu_2^2\sigma_1^2-2\sigma_1^2\sigma_2^2\log\left(\frac{\sigma_2}{\sigma_1}\right)}{\sigma_2^2-\sigma_1^2}\end{align}\tag{4}$$

Depending on the sign of $\sigma_2-\sigma_1$ we have two conditions equivalent to $p_1(x)\ge p_2(x)$:

$$\begin{align}\sigma_2>\sigma_1:\quad & x_1\le x\le x_2\\
\sigma_2<\sigma_1:\quad & x\le x_1\text{ or }x\ge x_2\end{align}\tag{5}$$

The inequalities in $(5)$ state explicit conditions on the variable $x$ that can be easily evaluated.

EDIT: As suggested by Dilip Sarwate in a comment it can be instructive to distinguish another special case: $\sigma_1\neq\sigma_2$ and $\mu_1=\mu_2=\mu$. In that case the limits $x_1$ and $x_2$ of the decision regions given by $(3)$ lie symmetrical to the common mean $\mu$, and $(3)$ simplifies to

$$x_{1,2}=\mu\pm\sqrt{\frac{2\sigma_1^2\sigma_2^2\log\left(\frac{\sigma_2}{\sigma_1}\right)}{\sigma_2^2-\sigma_1^2}}$$

Of course, the conditions $(5)$ remain valid.


If you sample your signal with $f_s=40\,\text{kHz}$ then the band between 20kHz and 30kHz will be folded back to the band between 10kHz and 20kHz. So aliasing will affect the upper half of your spectrum after sampling, i.e. the band between 10 and 20kHz. The lower band from DC to 10kHz will not be corrupted by aliasing.


Design one of 4 type of FIR filter by window function?

The topic your referring to goes by various names such as "FIR Filter Design Using Windows", "Windowing Method of FIR Filter Design", "Window Designed FIR Filters", etc. I'll bet your class textbook covers this topic (most DSP books do). Meaningful answers to your questions require more space than is available here. I suggest you have a look at the following web page (and its following three pages) for a decent introductory discussion of that topic:

http://www.bores.com/courses/intro/filters/4_fir.htm

Then have a look at Section 2.2 of the following web page:

http://www.mikroe.com/chapters/view/72/chapter-2-fir-filters/


What Is the Best Color Space to Work with for Hand Detection: RGB, HSV, or YCrCb?

Intuitively, HSV is the place to easily define Skin Color Hues.

Yet there is a broad work on that and even articles about the optimal Color Space for Skin Detection.

Yet, you should have a look at OPTIMUM COLOR SPACES FOR SKIN DETECTION (Alternative at IEEE - Optimum color spaces for skin detection).

According to them there is no difference in the performance as long as you create the Optimal Detector.
This actually makes sense, since all Color Spaces in the article are One to One functions it means the real thing here is to develop the correct classifier.


How Does Gaussian Blur Affect Image Variance

In Gaussian blur the value of each output pixel is calculated as a weighted sum of all input pixels:
$$\text{out}(x, y) = \sum^\infty_{j = -\infty} \sum^\infty_{i = -\infty} \frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\text{in}(x+i,y+j).$$
We want to calculate the variance $\text{Var}[\text{out}(x, y)] = \sigma^2_f$ based on the variances $\text{Var}[\text{in}(x + i, y + i)] = \sigma^2_0$. We require that $\text{in}(x + i, y + i)$ are independent, which you should mention in your question. Based on the identity (variance of linear combination):
$$Var\left[\sum_i c_i X_i\right] = \sum_i c^2_i Var[X_i],$$
where $c_i$ are constants and $X_i$ are independent random variables,
$$\sigma^2_f = \sigma^2_0\sum^\infty_{j = -\infty} \sum^\infty_{i = -\infty} \left(\frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\right)^2.$$
For a large $\sigma_G$, the squared Gaussian is rather smooth and the sum can be approximated by an integral as:
$$\sigma^2_f \approx \sigma^2_0\int^\infty_{-\infty} \int^\infty_{-\infty} \left(\frac{1}{{2\pi \sigma_G^2}} e^{-\frac{i^2 + j^2}{2 \sigma_G^2}}\right)^2\,\text{d}i\,\text{d}j
= \frac{\sigma^2_0}{4\pi \sigma_G^2}\\
\implies \sigma_f \approx \frac{\sigma_0}{2\sigma_G\sqrt{\pi}},$$
which is the same result as in Wikipedia. You can calculate the integral using Wolfram Alpha.

Note that the noise needs not be Gaussian for the above to apply.


A discrete-time signal represents a continuous-time signal that is band-limited to at most half the sampling frequency. Sampling a continuous-time signal that contains higher frequencies than half the sampling frequency will result in aliasing: Samples of those high frequencies are identical to and cannot be discerned from samples of other frequencies below half the sampling frequency. To avoid contamination of the signal by the aliased frequencies, you would first band-limit (lowpass filter) the continuous-time signal and only then sample it.


Inverse Fourier transform Of a triangular impulse

It's nice to see alternative paths, like this one: First of all, properties:

Frequency shifting: $x(t)e^{j2\pi f_0 t} \longleftrightarrow X(f-f_0)$
Duality: if $x(t) \longleftrightarrow X(f)$ then $X(t) \longleftrightarrow x(-f)$

You have correctly expressed the depicted transform: $$X(f) = \frac{1}{2}\mathrm{tri}\Big(\frac{f+f_0}{B}\Big) - \frac{1}{2}\mathrm{tri}\Big(\frac{f-f_0}{B}\Big)$$

As you said, $\mathrm{tri}\Big(\frac{t}{B}\Big) \longleftrightarrow B\mathrm{sinc}^2(fB)$. From the duality property, you have $$B\mathrm{sinc}^2(tB) \longleftrightarrow \mathrm{tri}\Big(\frac{-f}{B}\Big)$$ but since $\mathrm{tri(\cdot )}$ is an even function, you can write $$B\mathrm{sinc}^2(tB) \longleftrightarrow \mathrm{tri}\Big(\frac{f}{B}\Big)$$

Now you have your FT-pair you need. If you apply the frqeuency shifting property on $\mathrm{tri}\Big(\frac{f\pm f_0}{B}\Big)$, you can easily get... what? :-) You can continue from this point.


Determining the final value of the output of a discrete system

Your final equation is incorrect: it should be

$${\lim_{k\to\infty}f\left(kt\right)=\lim_{z\to1}\left(\left(1-z^{-1}\right)H\left(\text{z}\right)\right)} $$

Where $H(z)= F(z)G(z)$ which is the $\mathcal Z$-transform of the output of the system that has a transfer function of $F(z)$ with $G(z)$ applied to the input.

From that you already have $F(z)$ from the first line, multiply it by $G(z)$ to get the $\mathcal Z$-transform of the function that would be at the output of the system, multiply that product by $(1-z^{-1})$, and then take the limit of the resulting equation as $z$ goes to $1$.


Phase rotation without Fourier Transform

You can use the Hilbert transform, then multiply the real part with the sine of the angle you wish to transform, and the imaginary part wth the cosine. A quick code in Octave:
t=[0:0.01:4];
s=hilbert(1-2*mod(t,1));
a=real(s);
b=imag(s);
phi=[0:0.5:6];
test=a.*sin(phi)+b.*cos(phi);
plot(a+4,"",test');

I know there are Hilbert transforms made with IIRs, but their phase is not exactly linear, so a FIR is a better, but costlier choice. Its order will be a function of the lowest frequency of interest, which will count as the transition width. If your bandwidth is large, the order will be large.

A FIR is a filter, which means it acts in time domain (but, implicitly, the frequency domain will be affected). To give an example, the signal you're showing of 1 Hz, sampled at 32 Hz, filtered with a FIR with a Kaiser window with As=60 and 1 Hz transition width:

The output isn't sampled because I am brute-adding the two inputs, but that's how the result would look like. The green number is the order. It's a symmetrical FIR so you'll need that many delays and half as many multipliers. You could try a polyphase approach and that will reduce the number of delays, but not the multipliers.


What is the DTW result (applied to the MFCCs of an audio file) proportional to?

DTW is not proportional to anything, including frame length. It's not a good idea to compare DTW values between two different pairs like you are doing.

The issue is that DTW is not just the number but also an alignment. If alignment does not match, you can not relate numbers. You only can relate numbers if you align to the same thing.

Normalization is often suggested, but mathematically it doesn't have much sense.

In classical tasks like speech recognition you never compare the DTW values between pairs, instead, you compare a sample to several templates and find the best one. The fact that you try to match to same sample allows you to compare DTW values. In keyword spotting using DTW the threshold must depend on every template you are matching, you can not use single threshold for all the templates. Even more, for the two different templates of the same length DTW thresholds must be different.


Any better way to detect Signal with ADC Overflow problem?

How do you transform overflowed signal back to normal? Just by adding or subtracting multiplies of 256. So with every new sample you just determine whether overflow occured - you need to define some threshold and compare it with first derivative of your signal, e.g. x[i+1]-x[i] > 200. If this happened, you know that you need to subtract 256 from x[i+1] to get continuous signal. Opposite case x[i+1]-x[1] < -200 means that you need to add 256. So the pseudocode would be:

base = 0
thresh = 200
for each i in x
  d = x[i] - x[i-1]
  if d > thresh then base = base - 256
  if d < -thresh then base = base + 256
  y[i] = x[i] + base
next


What does it mean that the DFT equals the Complex Fourier Coefficients for even frequencies?

okay, i am using these definitions for the continuous Fourier Transform

$$ \mathscr{F}\Big\{x(t)\Big\} \triangleq X(f) \triangleq \int_{-\infty}^{+\infty} x(t) \, e^{-i 2 \pi f t} \, \mathrm{d}t  $$

and inverse

$$ \mathscr{F}^{-1}\Big\{X(f)\Big\} \triangleq x(t) = \int_{-\infty}^{+\infty} X(f) \, e^{i 2 \pi f t} \, \mathrm{d}f $$

$f$ is "ordinary" frequency, not angular frequency (which for the continuous-time signals, the notation we use is $\Omega \triangleq 2 \pi f$).  this is a preferred representation of the Fourier Transform for electrical engineers because it is the easiest manner to apply the Duality property (lotsa similarity between the FT and the inverse).

in this answer is shown how sampling in one domain (e.g. the "time" domain) causes periodic extension (repetition and overlapping and adding) in the reciprocal domain (e.g. the "frequency" domain).  this is pointed out without regard to bandlimiting, so if the input is not sufficiently bandlimited, there will be aliasing.  given the assumption of a real input to start with (which means the Fourier transform or two-sided spectrum has symmetry about $f=0$) frequency components above Nyquist (which is $\tfrac12$ the sample rate) are folded back to locations below Nyquist and added to whatever frequency components existed there.

now, your continuous-time input

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} c_k \ e^{i 2 \pi k t} $$

has a spectrum (which is not necessarily bandlimited, and is actually not bandlimited since the $c_k$ coefficients are not zero above some finite $k$) of

$$ X(f) = \sum\limits_{k=-\infty}^{\infty} c_k \ \delta(f - k) $$

first defining the Dirac comb:

$$ \operatorname{III}_T(t) \triangleq \sum_{n=-\infty}^{\infty} \delta(t-nT) $$

now, the ideally sampled signal (still represented in the continuous-time domain) is

$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \sum_{n=-\infty}^{\infty} \delta(t-nT) \\
 &= T \sum_{n=-\infty}^{\infty} x(t) \ \delta(t-nT) \\
 &= T \sum_{n=-\infty}^{\infty} x(nT) \ \delta(t-nT) \\
 &= T \sum_{n=-\infty}^{\infty} x[n] \ \delta(t-nT) \\
\end{align}$$

and also this

$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \cdot \sum_{m=-\infty}^{\infty} \frac1T e^{i 2 \pi m \frac{t}{T}} \\
 &= \sum_{m=-\infty}^{\infty} x(t) \ e^{i 2 \pi m \frac{t}{T}} \\
\end{align}$$

now from the latter, we know that the spectrum of the ideally sampling signal is this:

$$ \mathscr{F}\Big\{x_\text{s}(t)\Big\} = X_\text{s}(f) = \sum_{m=-\infty}^{\infty} X\left(f - \tfrac{m}{T} \right)
 $$

which is periodic in the frequency domain with period $\frac{1}{T}$ where your sampling period $T=\tfrac{1}{N}$.  that means your sample rate is $\tfrac{1}{T}=N$ and your Nyquist frequency is $\tfrac{N}{2}$.  that is the "folding frequency" or "foldover frequency".

Okay, winding up, you plug the expression of $X(f)$, the non-bandlimited Fourier series, into the equation for $X_\text{s}(f)$, which is the spectrum of the ideally-sampled signal.  And you plug in $N$ for $\tfrac1T$.

$$ X(f) = \sum\limits_{k=-\infty}^{\infty} c_k \ \delta(f - k) $$

$$\begin{align}
 X_\text{s}(f) &= \sum_{m=-\infty}^{\infty} X\left(f - \tfrac{m}{T} \right) \\
  &= \sum_{m=-\infty}^{\infty} X(f - N m) \\
  &= \sum_{m=-\infty}^{\infty} \sum\limits_{k=-\infty}^{\infty} c_k \ \delta\big((f - N m) - k\big) \\
  &= \sum\limits_{k=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} c_k \ \delta\big(f - (k+Nm) \big) \\
  &= \sum\limits_{k=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} c_{k-Nm} \ \delta(f - k) \\
  &= \sum\limits_{k=-\infty}^{\infty} \tilde{c}_k \ \delta(f - k) \\
\end{align}$$

where

$$ \tilde{c}_k \triangleq \sum_{m=-\infty}^{\infty} c_{k-Nm} $$

so the coefficients of the Fourier series of the sampled periodic signal $x_\text{s}(t)$ are different from those of $x(t)$.  even in the baseband where the DTFT (and following, the DFT) sees the spectrum.


How do I apply a function window to a signal?

How can I apply a window function like Hamming or Lanczos to a signal, using its coefficients?

Just multiply, point-by-point.


Relationship between Nyquist and Shannon channel capacity

Well, I think both are unrelated unless you have some constraints.

The Nyquist channel capacity says you can transmit $2B\log_2(M)$ bits per second for given channel bandwidth $B$. It does not say if these bits would be received reliably at the receiver in the presence of noise or channel distortion. I am assuming here that bandwidth limit is imposed by channel. Your transmitter can afford to have higher bandwidth but channel does not.

The Shannon Channel capacity theorem says, in order to receive the transmitted bits reliably (that is Probability of Error $\text{P}_e \rightarrow 0$) in the presence of additive noise, the maximum rate at which you can transmit is $B\log_2(1+\text{SNR})$ bits per second.

So if your requirement is reliable communication and the channel has a bandwidth limit of $B$, then $C_n \le C_s$.


It seems you mixed up CIC filters and recursive moving-average

A recursive moving-average uses a comb-filter and an accumulator.

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

The numerator corresponds to a comb filter (i.e a delay line of N elements and 1 adder).
The denominator corresponds to an accumulator. The 1/N factor is to get a moving-average, otherwise you'd have a moving-accumulator. You only need 2 additions, a delay line, an accumulator and a divider. With a non-recursive moving-average filter, you'd need N additions and you'd need a delay line anyway. A recursive moving-average is much
less computionally expansive than a non-recursive moving-average filter.

Now the CIC filter combines the recursive moving-average filter with decimation and interpolation.

For a CIC decimation filter, you insert the rate transition block between the accumulator and the comb filter. Say you have a 16-sample moving average-filter (N=16) and you want to decimate by 8 (R = 8). You'd need a delay line of 16 elements. However, by inserting the decimation between the accumulator and the comb filter, you can reduce the number of elements in the delay line. You'd only need N/R elements, i.e. 2. You'd still need to subtract 2 numbers, however you would do it at sampling rate of fs/8.

I suggest you consult this site

https://www.researchgate.net/publication/228648728_Understanding_cascaded_integrator-comb_filters

Update: see this nice write-up by Richard Lyons-
https://www.dsprelated.com/showarticle/1337.php


ROC of this LTI system given $x[n]$ and $y[n]$

Let's look at the pole locations and the corresponding regions of convergence (ROCs) of the $\mathcal{Z}$-transforms of the input and output signals:

$X(z)$: poles at $z=\frac12$ and $z=2$ with ROC $\frac12<|z|<2$ (two-sided).
$Y(z)$: poles at $z=\frac12$ and $z=\frac34$ with ROC $|z|>\frac34$ (causal).

Now we know that the ROC of the multiplication of two $\mathcal{Z}$-transforms equals the intersection of the ROCs of the individual $\mathcal{Z}$-transforms, unless there is pole-zero cancellation, in which case the resulting ROC can be bigger than the intersection.

$H(z)$ has a zero at $z=2$ and a pole at $z=\frac34$. Since the zero cancels the pole of $X(z)$, the restriction $|z|<2$ for the ROC of $Y(z)$ is removed. The two theoretically possible ROCs of $H(z)$ are $|z|<\frac34$, corresponding to an anti-causal and unstable system, and $|z|>\frac34$, corresponding to a causal and stable system. In the first case, the ROC of $Y(z)$ would need to be $\frac12<|z|<\frac34$, which is not the case, and which would correspond to a two-sided unstable signal. $^1$ Consequently, the ROC of $H(z)$ must be $|z|>\frac34$, which means that the system described by $H(z)$ is causal and stable.

1. I use the term unstable signal for a signal that, if interpreted as the impulse response of an LTI system, would correspond to an unstable system (i.e., the ROC does not contain the unit circle).


Phase vocoder, pitch shift in MATLAB and C

alright, conceptually, for one frame (which is centered at $n=0$ and has non-zero width of $N-1$ samples):

$$\begin{align}
 X[k] &= \mathcal{DFT}\Big\{ x[n] w(n) \Big\} \\
\\
  &= \sum\limits_{n=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} \big( x[n] w(n) \big) \, e^{-j \frac{2 \pi}{N} nk} \qquad |k| \le \tfrac{N}2 \\
\end{align}$$

here, our frame hop is $\tfrac{N}2$ samples, which is half of the window width.

$w(t)$ is a continuous-time complementary window with non-zero width of $N$ samples and centered at $t=0$, like a Hann window:

$$ w(t) = \begin{cases} \tfrac12 + \tfrac12 \cos\left(\tfrac{2\pi}{N}t \right) \qquad & |t| < \tfrac{N}2 \\
 0 \qquad & |t| \ge \tfrac{N}2 \\
\end{cases} $$

the output spectrum is a scaled copy of the input spectrum.  this is what shifts the pitch.

$$ Y[k] = \begin{cases}
 X(\alpha^{-1} k) \qquad & |k| \le \tfrac{N \alpha}{2} \\
 0 \qquad & |k| > \tfrac{N \alpha}{2} \\
\end{cases} $$

where $\alpha$ is the output-to-input frequency ratio.  $X(f)$ is a continuous-frequency function that is interpolated from the discrete-frequency spectrum $X[k]$ in such a way that equality exists when the argument is an integer.

$$ X(f) \bigg|_{f=k} = X[k] \qquad k \in \mathbb{Z} $$

the discrete-time-domain output for this frame (before adjusting for stretching or scrunching the window) is:

$$\begin{align}
 y[n] &= \mathcal{iDFT}\Big\{ Y[k]  e^{j \phi[k]} \Big\} \\
\\
  &= \tfrac{1}{N} \sum\limits_{k=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} \big( Y[k]  e^{j \phi[k]} \big) \, e^{j \frac{2 \pi}{N} nk}  \qquad |n| \le \tfrac{N}2 \\
\end{align}$$

and the frame of output (centered at time $n=0$) is

$$ y[n] \frac{w(n)}{w(\alpha \, n)}  $$

the $\frac{w(n)}{w(\alpha \, n)}$ factor is undoing the stretched or scrunched input window and re-applying the original window at the same scale (since we're not time-scaling and the frame hop length remains unchanged from input to output).

$\phi[k]$ is a phase adjustment (with odd symmetry in $k$) that is constant for each sinusoidal component (a la Miller Puckette) rather than a changing phase adjustment for each FFT bin (a la Portnoff).  that phase adjustment for each frequency component is what is required to make each sinusoid component continuous because of the frequency shift.  that's how the phase vocoder does glitch-free pitch shifting, even if the input is not periodic or a single harmonic note.

this is why this phase adjustment is needed.  consider a single sinusoid with normalized angular frequency $\omega_0$.

$$\begin{align}
 x[n] &= \cos(\omega_0 n + \theta_0) \\
  &= \tfrac12\big(e^{j (\omega_0 n + \theta_0)} + e^{-j (\omega_0 n + \theta_0)}  \big) \\
  &= \tfrac12\big(e^{j\theta_0}e^{j \omega_0 n} + e^{-j\theta_0}e^{-j \omega_0 n} \big) \\
\end{align}$$

now let's consider only the positive frequency complex component.

$$ \hat{x}[n] = e^{j\theta_0}e^{j \omega_0 n} $$

this has an instantaneous angle of $\theta_0$ at time $n=0$.  now the instantaneous value of this sinusoid in the center of the previous frame is

$$ \hat{x}[-\tfrac{N}{2}] = e^{j\theta_0}e^{-j \omega_0 \frac{N}{2}} $$

and the instantaneous angle is $\theta_0 - \tfrac{N}{2} \omega_0$.  and it shouldn't surprise us that halfway between the two frame centers (which is exactly in the middle of the crossfade from the previous frame to the current frame) the instantaneous angle $\theta_0 - \tfrac{N}{4} \omega_0$

now if the output spectrum is mapped from the interpolated input spectrum as

$$ Y[k] = X(\alpha^{-1} k) $$

then this component frequency at the input, $\omega_0$ gets mapped to $\alpha \omega_0$.  and the current frame output (before the phase adjustment) is

$$\begin{align}
 y[n] &= \cos(\alpha\omega_0 n + \theta_0) \\
  &= \tfrac12\big(e^{j (\alpha\omega_0 n + \theta_0)} + e^{-j (\alpha\omega_0 n + \theta_0)}  \big) \\
  &= \tfrac12\big(e^{j\theta_0}e^{j \alpha\omega_0 n} + e^{-j\theta_0}e^{-j \alpha\omega_0 n} \big) \\
\end{align}$$

and the positive frequency component is

$$ \hat{y}[n] = e^{j\theta_0}e^{j \alpha \omega_0 n} $$

and the value at the left edge of the current frame where it is joined to the right edge of the previous frame

$$ \hat{y}[-\tfrac{N}{4}] = e^{j\theta_0}e^{-j \alpha \omega_0 \frac{N}{4}} $$

now the instantaneous angle of the output sinusoid at the center of the previous frame is the same instantaneous angle of the input sinuosoid at the center of the previous frame which is $\theta_0 - \tfrac{N}{2} \omega_0$.  this makes the output sinuosoid (the positive frequency component) of the previous frame, when "justified" or expressed in terms of the time index of the current frame

$$ \hat{y}[n] = e^{j(\theta_0 - \frac{N}{2}\omega_0)}e^{j \alpha \omega_0 (n + \frac{N}{2})} $$

(note that when $n=-\tfrac{N}{2}$, which is the center of the previous frame, the instantaneous angle is $\theta_0 - \tfrac{N}{2} \omega_0$.)  now at the output sample midway between the previous frame and the current frame (which is $n=-\tfrac{N}{4}$), the sinusoid of the previous frame is at angle

$$\theta_0 - \tfrac{N}{2}\omega_0 + \alpha \omega_0 (-\tfrac{N}{4} + \tfrac{N}{2})
 = \theta_0 + (\alpha - 2)\omega_0\tfrac{N}{4}$$

but, at that same sample, the angle of the sinusoid of the current frame is

$$ \theta_0 - \alpha \omega_0 \tfrac{N}{4} $$

note that when there is no pitch shift and $\alpha = 1$, then the two angles are the same and the splice is seamless.  but when $\alpha \ne 1$, the phase of one of the frames must be adjusted to make the splice seamless.  the previous frame is already a "done deal" and will not be modified, so it's the phase of the sinusoid of the current from that must be adjusted to make the splice seamless.  that phase is modified for all DFT bins associated with this particular sinusoid as

$$\begin{align}
 \theta_0 - \alpha \omega_0 \tfrac{N}{4} + \phi[k] &= \theta_0 + (\alpha - 2)\omega_0\tfrac{N}{4} \\
 -\alpha \omega_0 \tfrac{N}{4} + \phi[k] &= (\alpha - 2)\omega_0\tfrac{N}{4} \\
  \phi[k] &= (\alpha - 1)\omega_0\tfrac{N}{2} \\
\end{align}$$

so the phase adjustment for the current frame necessary to align it with the previous frame is the normalized angular frequency of the sinusoid, $\omega_0$, times the frame hop displacement $\tfrac{N}{2}$ times the factor $(\alpha - 1)$.  you should apply same phase adjustment to all DFT bins of $Y[k]$ that are associated with this particular sinusoidal component at frequency $\omega_0$.  the values of $k$ will be around $\tfrac{\omega_0}{2 \pi} N$.

so Miller Puckette says

$$\phi[k] = (\alpha - 1) \omega_0 \tfrac{N}{2} \qquad \text{for } k \approx \tfrac{\omega_0}{2 \pi} N$$

and Portnoff would say that the adjustment would be

$$\phi[k] = (\alpha - 1) \tfrac{2 \pi k}{N} \tfrac{N}{2} = (\alpha - 1) \pi k$$

but, in my opinion, Puckette is correct and Portnoff is wrong.


Upsampling signal using convolution-based interpolation filters

Starting from the equation with the convolutions and assuming unit distance between sampling grid points, for a one-dimensional signal you'll have just:

$$I(\mathbf{x}) = I((x_1))= I_s((p_1))*h(x_1),$$

where $\mathbf{x}$ is the vector of continuous coordinates of which there is only one: $x_1 \in \mathbb{R}.$ Likewise there is just a single discrete coordinate $p_1 \in \mathbb{Z}.$ For simplicity, you could drop the vector notation and subscripts altogether and have just:

$$I(x) = I_s(p)*h(x).$$

Without having fully read the paper, I don't think the authors dare to venture into defining what is meant by convolution between discrete-time $I_s(p)$ and continuous-time $h(x)$ which results in continuous-time $I(x)$. Such a convolution can be defined as:

$$I(x) = \sum_{p=-\infty}^{\infty}I_s(p)h(x-p)$$

For example, if you want to obtain $I(x)$ within the inter-sample interval $0 \le x < 1$ and you use the linear interpolation kernel, then the above simplifies as:

$$\begin{eqnarray}I(x) &=& \sum_{p=-\infty}^{\infty}I_s(p)h_{Lin}(x-p)\\
&=& \sum_{p=-\infty}^{\infty}I_s(p)\times\cases{
1-|x-p|, &if$\quad 0 \leqslant |x-p| < 1,$\\
0, &if$\quad 1 \leqslant |x-p|$
}\\
&=& \sum_{p=0}^{1}I_s(p)(1-|x-p|)\\
&=& I_s(0)(1-x) + I_s(1)(1-(1-x))\\
&=& I_s(0)(1-x) + I_s(1)x.\end{eqnarray}$$

Or likewise for any inter-sample interval $x_0 \le x < x_0 + 1:$

$$I(x) = I_s(x_0)(1-(x - x_0)) + I_s(x_0 + 1)(x - x_0).$$

Another way to understand such convolutions is that $I_s(p)$ is the product of the continuous-time signal multiplied by a Dirac comb, and to use a continuous-time convolution integral.

As you can see from the above, the interpolation filter is not a regular discrete-time filter, but one that converts a discrete-time signal into a continuous-time signal. The convolution filter definitions cited in the question are directly applicable, without further "sizing", when the sampling interval is 1.


Amount of exactly approximated samples with Padé approximation

The problem and its solution is consistent with the Monson Hayes' book that describes Padé approximation.

Your mistake is that you've taken $p=2$ but $q=3$ which is wrong. Correct is $q=2$ and there are 3 unknowns in $b[k]$ , that the problem defines a second order system , i.e, $a[k] = [1, a_1, a_2]$ and $b[k] = [b_0, b_1, b_2]$ , according to the book.

Since there are $p+q+1 = 2+2+1 = 5$ (that +1 is for $b_0$) unknown coefficients in the system transfer function $H(z)$, then from its impulse response $h[n],$ only the first 5 samples will match to the data $x[n]$ for $n = 0,1,2,3,4$.

And from the given solution as you can see there are 5 equations for 5 unknowns, of last 2 is used to solve for the unknowns $a_1,a_2$ and the first 3 used to find $b_0,b_1,b_2$. The data set used is the $x = [x_0, x_1,x_2,x_3,x_4] = [2,1,0,-1,0]$.

Note that for the first two equations $x(-1)$ and $x(-2)$ are implictly referenced, which are assumed to be zero for the causal modeling.


Your professor is using sloppy notation. When one writes the sample value $x[\cdot]$, it is generally accepted that the argument is an integer (which $nT$ might not be) and we relate this to the continuous-time signal $x(t)$ via $x[n] = x(nT)$. So, what your professor should have said is that $y[n]$ has value $T\cdot x[n] = T\cdot x(nT)$, and not that $y[n] = T\cdot x[nT]$ (even if $T$ happens to be an integer).  An outline of a proof that $y[n] = T\cdot x[n] = T\cdot x(nT)$ is given below.

If $x(t)$ is a continuous-time signal bandlimited to $\left(-\frac{1}{2T}, \frac{1}{2T}\right)$, that is, its Fourier transform $X(f)$ has the property that $X(f) = 0$ for all $|f| \geq \frac{1}{2T}$,then it can be reconstructed from its samples spaced $T$ seconds apart. So, let's consider the signal $x(t)\operatorname{sinc}\left(\frac{t}{T}\right)$ where
$$\operatorname{sinc}(t) = \begin{cases}\frac{\sin(\pi t)}{\pi t},& t \neq 0,\\1, & t = 0.\end{cases}$$
The Fourier transform of $\operatorname{sinc}(t)$ is $\operatorname{rect}(f)$ and so the Fourier transform of $\operatorname{sinc}\left(\frac{t}{T}\right)$ is $T\operatorname{rect}(fT)$, and the Fourier transform of $x(t)\operatorname{sinc}\left(\frac{t}{T}\right)$ is the convolution  $X(f)\circledast T\operatorname{rect}(fT)$ of the Fourier transforms $X(f)$ and $T\operatorname{rect}(fT)$ of $x(t)$ and $\operatorname{sinc}\left(\frac{t}{T}\right)$ respectively.  Now, in principle, the convolution of a function and a rectangular pulse is easy to compute but I refuse to write out the answer because it is not needed: what we need to find is the value of $\int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t}{T}\right) \mathrm dt$ which, if you think about it a bit, is actually the value of the Fourier transform of $x(t)\operatorname{sinc}\left(\frac{t}{T}\right)$ evaluated at $f=0$. That is,
\begin{align}
\int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t}{T}\right) \mathrm dt &= X(f)\circledast T\operatorname{rect}(fT)\bigg |_{f=0}\\
&= T\int_{-\infty}^\infty X(\lambda)\operatorname{rect}((0-\lambda) T) \, \mathrm d\lambda\\
&= T\int_{-\infty}^\infty X(\lambda)\operatorname{rect}(\lambda T) \, \mathrm d\lambda & \scriptstyle\text{rect is an even function}\\
&= T\int_{-\frac{1}{2T}}^{\frac{1}{2T}}X(\lambda) \, \mathrm d\lambda &\scriptstyle\text{rect has finite support}\\
&= T\int_{-\infty}^\infty X(\lambda) \, \mathrm d\lambda &\scriptstyle{X(\lambda)=0 ~\text{for}~|\lambda| > \frac{1}{2T}}\\
&= T
\cdot x(0).
\end{align}
Thus, $y[0] \overset{\text{def}}{=} \int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t}{T}\right) \mathrm dt$ equals $T\cdot x(0)$. I leave it as an exercise for the reader to work out the Fourier transform of $\operatorname{sinc}\left(\frac{t-nT}{T}\right)$ and show that $y[n] \overset{\text{def}}{=} \int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t-nT}{T}\right) \mathrm dt$ works out to be $T\cdot x(nT)$ and so if we define the samples of $x(t)$ as $x[n] = x(nT)$, then we have that $y[n] = T\cdot x[n]$ for all $n$.


Nonnegative or positive band-limited interpolation

If we require the impulse response of the interpolation filter to be non-negative then we can obviously guarantee that the interpolated signal is also non-negative. However, it can be shown that the frequency response of a low pass filter with non-negative impulse response is always far from ideal. More specifically, it can be shown that if $\epsilon$ is the maximum negative pass band deviation of the squared magnitude $|H(\omega)|^2$ of the filter's transfer function from its value at $\omega=0$

$$|H(\omega)|^2\ge |H(0)|^2-\epsilon,\quad |\omega|<\omega_c \tag{1}$$

where $\omega_c$ is the cut-off frequency, then the following inequality holds:

$$|H(\omega)|^2\ge |H(0)|^2-4\epsilon,\quad |\omega|<2\omega_c\tag{2}$$

This inequality imposes a severe restriction on the minimum stop band attenuation if the ripple size $\epsilon$ is chosen sufficiently small (as should be the case for a good low pass filter).

As an example, if $\epsilon$ is chosen such that the pass band ripple is $1$ dB, the attenuation at $\omega=2\omega_c$ cannot be better than $7.5$ dB.

Of course, if $\epsilon\ge |H(0)|^2/4$ then $(2)$ represents no restriction, but in that case the maximum pass band deviation of the magnitude of the transfer function is greater than usually required for a good low pass filter.

Since the inequality $(2)$ is independent of the filter order, an interpolation filter with arbitrary accuracy cannot be realized by a filter with non-negative impulse response.

The proof of $(2)$ can be found in Attenuation Limits for Filters with Monotonic Step Response by A. Papoulis, IRE Transactions on Circuit Theory ( Volume: 9, Issue: 1, March 1962).


What is the right way to extract particular frequency phase information from FFT?

To estimate spectral component phase angle relative to shaft angle using an FFT of length N, first gather a sequence of N samples such that the shaft angle is known with accuracy at sample number N/2, then window with a Nuttal window and do an FFT shift before the FFT.  The  atan2() of the FFT bin result closest to the desired spectral component will produce an estimate of the phase relative to the center of the original sample window.  Adjust that phase by the known shaft angle at that N/2 sample point to get a spectral phase relative to shaft angle.

If you know the spectral frequency or use a frequency peak estimator, then the phase can also be interpolated between bins for increases phase estimation accuracy.  Interpolate phase by using a Sinc interpolation seperately on the real (even) and imaginary (odd) FFT result components before taking the arctangent.


You can probably use Matlab's filter command. If you call it like this:

[Y,Zf] = filter(B,A,X,Zi);

the variable Zf contains the state of the filter after filtering, and Zi contains the filter's state before filtering. You could proceed roughly like this (in pseudocode):

Acquire several samples and store them in vector Samples
[Y,Z] = filter(B,1,Samples); % assuming FIR filter stored in B
while 1
    acquire new sample and store it in Sample
    [Y,Z] = filter(B,1,Sample,Z);  % Y contains the updated filter output
end

To design the filter, I'd start by looking at the fir1 command.


Suppose you're imaging a square with sides of length $d$, and you're a height $h \gg d$ above the square. You're tilted in one dimension by $\theta$. Then the center of the square appears displaced by $-h \sin \theta$ and its apparent width along that dimension is $d \cos \theta$, not $d$. So, what object would be imaged onto a single pixel of size $d_p$? It would need to be a rectangle, displaced from the image center by $h \sin \theta$. Its width along the tilted dimension would have to be $d / cos \theta$, but along the untilted dimension it should still be $d$.

So after correction, the image consists of rectangular pixels of width $d / \cos \theta$. There's no transform on the raw image data that will implement this directly. However, it's easy to write a function that will return the value of the image at a given lat/long.


confusion about poles and zeros of Lead compensator?

You're right, the given pole and zero are wrong. They should be

$$s_0=-\frac{1}{\tau}$$

and

$$s_{\infty}=-\frac{1}{\beta\tau}$$

because for $s=s_0$ the numerator becomes zero, and for $s=s_{\infty}$ the denominator becomes zero.


Out of memory problem with filter design tool in Matlab

Your sampling frequency is currently set to 5000000000 MHz, which is 5000 terahertz!  Notice that the units on the frequency specification should all be in MHz. Change the Fs value to be 5000 (MHz) and you will be much better off.

However, even once you fix this problem, you still may run into trouble.  You are specifying transition bands that have a width of less than 0.01% of your sampling rate.  This is not going to be practical, especially as an FIR filter.  I recommend that you consider a good low-pass filter followed by lots of down sampling.  Then you will be able to realistically extract the band of interest.


To add here are two diagrams showing common transceiver architectures: (1) a super-heterodyne where the down-conversion is done first to an IF frequency and then to baseband and (2) a zero-IF receiver where RF is translated directly to baseband. Note in both architectures it is arbitrary (technical / technology choice) where the ADC/DAC boundary is as either approaches could be done with baseband, IF or direct RF sampling to the extent technology allows. Ultimately, at baseband we desire a complex IQ signal to support most modern modulations which require asymmetrical spectrums about the RF carrier for maximum spectral efficiency (and hence a complex baseband signal with I and Q components). We can separate into I and Q in the analog using quadrature local oscillators and then use a dual ADC/DAC, or we can use a single ADC/DAC and then separate into I and Q digitally using a quadrature NCO.

The transmit and receive architectures need not match (Can use a ZIF transmitter and Super-het receiver for example).

For high efficiency with a digital IF/RF approach, quadrature sampling can be implemented at a $f_s/4$ IF frequency with a single ADC datapath interleaving the I and Q channels as as follows:

Channel I: 1 0 -1 0 1 0 -1 0

Channel Q: 0 1 0 -1 0 1 0 -1

This concept is further illustrated in the figure below. Sampling at Fs/4:

A subsequent phase rotator would be required to remove carrier offsets, or the ADC clock itself can be in the carrier tracking loop.

Placing the digital IF at $F_s/4$ or $N F_s \pm F_s/4$ also has the advantage of simplifying ADC anti-alias filter design as the images will all be equidistant.


I have a PRN generator that I have validated with live captured signals that is available on the Mathworks Exchange site at this address and equally runs in Octave (Update: I also pasted the core of this in a code block below):

https://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator

The two tap coder is as given in the diagram in that spec and copied in the figure below:

What they mean by “two tap coder” is the selection and addition of two outputs of the G2 generator in the diagram and that this is just one way to generate the sequences (a very convenient way). An alternate approach which I have done is to implement the generator with one LFSR. In this case you simply convolve the coefficients for the two generator polynomials. This results in a length 20 shift register and non-maximum length sequence (since it clearly can be factored), and will generate 1023 different sequences each 1023 long, depending on how you seed it. (So if you change the state of the 20 elements with any 20 consecutive chips from a PRN of interest, it will continue to run the code for that SV).    So when they mention it will only generate a limited selection of the Gold Codes, they mean it won’t generate all 1023 possible codes (because with the selection of two taps there are only $9!=42$ combinations -if I did that right, and if so then I highly suspect that the selection of which of the 1023 different codes to use was based on the assumption of using this 2-tap generator). The two-tap generator shown won’t generate all possible Gold Codes for these polynomials used but it will indeed generate all the GPS C/A codes.

This is because the Gold Code (which is the GPS C/A code) is formed by adding two maximum length sequences in GF(2) (basically x-or'ing the output). Doing that is identical to what could be done with one generator where the polynomials are multiplied together, and you multiply polynomials by convolving their coefficients.

The advantage of doing that with just one long LFSR would be the creation of a flexible code generator capable of creating any maximum length sequence up to $2^{20}-1$ including all the GPS C/A codes and many other compound codes.

The interesting point about the Gold Code is that if you delay one of the generators relative to the other by even 1 chip, the result will be a completely new 1023 length sequence. So since there are 1023 possible delays, there are 1023 total possible sequences that can be generated, each 1023 long.

Further an interesting point about maximum length sequences in general, is if you GF(2) add the sequence to a delayed version of itself, it will create the exact same sequence just at a completely different delay! Thus you see what is happening in the Gold Code Generator with the "two-tap coder" which is used to select the satellite; in the figure I gave we select taps 2 and 6 resulting in SV1 (which should match the first part of the table that you didn't post). The table you have shows the equivalent delay of the G2 coder for each two tap selection that is made.

Since a login is required for the MathExchange site, here is a copy of my code with the front-end error-trapping removed (for brevity):

function g=cacode(sv,fs)
% function G=CACODE(SV,FS)
% Generates 1023 length C/A Codes for GPS PRNs 1-37
%
%
% g: nx1023 matrix- with each PRN in each row with symbols 1 and 0
% sv: a row or column vector of the SV's to be generated
%    valid entries are 1 to 37
% fs: optional number of samples per chip (defaults to 1), fractional samples allowed, must be 1 or greater.
%
% For multiple samples per chip, function is a zero order hold.
%
%
% For example to generate the C/A codes for PRN 6 and PRN 12 use:
% g=cacode([6 12]),
% and to generate the C/A codes for PRN 6 and PRN 12 at 5 MHz use
% g=cacode([6 12],5/1.023)
%
%
% Dan Boschen 12-30-2007
% boschen@loglin.com

% table of C/A Code Tap Selection (sets delay for G2 generator)
tap=[2 6;
    3 7;
    4 8;
    5 9;
    1 9;
    2 10;
    1 8;
    2 9;
    3 10;
    2 3;
    3 4;
    5 6;
    6 7;
    7 8;
    8 9;
    9 10;
    1 4;
    2 5;
    3 6;
    4 7;
    5 8;
    6 9;
    1 3;
    4 6;
    5 7;
    6 8;
    7 9;
    8 10;
    1 6;
    2 7;
    3 8;
    4 9
    5 10
    4 10
    1 7
    2 8
    4 10];

% G1 LFSR: x^10+x^3+1
s=[0 0 1 0 0 0 0 0 0 1];
n=length(s);
g1=ones(1,n);   %initialization vector for G1
L=2^n-1;

% G2 LFSR: x^10+x^9+x^8+x^6+x^3+x^2+1
t=[0 1 1 0 0 1 0 1 1 1];
q=ones(1,n);    %initialization vector for G2

% generate C/A Code sequences:
tap_sel=tap(sv,:);
for inc=1:L
    g2(:,inc)=mod(sum(q(tap_sel),2),2);
    g(:,inc)=mod(g1(n)+g2(:,inc),2);
   g1=[mod(sum(g1.*s),2) g1(1:n-1)];
   q=[mod(sum(q.*t),2) q(1:n-1)];
end

%upsample to desired rate
if fs~=1
    %fractional upsampling with zero order hold
    index=0;
    for cnt = 1/fs:1/fs:L
        index=index+1;
        if ceil(cnt) > L   %traps a floating point error in index
            gfs(:,index)=g(:,L);
        else
            gfs(:,index)=g(:,ceil(cnt));
        end
    end
    g=gfs;
end


Is ideal LPF a LTI system? And is it time invariant?

If your system is defined by its time-invariant impulse response $h(t)$, then it is by definition a linear time-invariant system (LTI). The output of the system is the convolution of the input signal with its impulse response. This holds, since convolution is translation-invariant (https://en.wikipedia.org/wiki/Convolution#Properties).

A time-invariant impulse response can be recognized, if it only depends on one time-variable, i.e. $h(t)$. For time-varying multipath channels, the impulse response can become time-variant, and is in this case expressed as $h(t,\tau)$ and the output of the system is given by

$$y(t)=\int x(t-\tau)h(t,\tau)d\tau.$$

However, these systems are more advanced and not often studied in fundamental courses.


How to apply dictionary learning to several images?

Dictionary learning is usually done on more than one image. So yes, I would normalize and concatenate multiple images together. What I usually do is create 8x8 patches from images and vectorize them so that they represent one column (data-item) from the training set (X in your case).


Diagonal vs non diagonal processing for audio denoising

I find the diagonal term a bit confusing. Sometimes, for wavelets, people separate scalar and vector or block thresholding. [EDIT] Now I believe I understand it. If you put all your samples or coefficients in a long 1D vector, then multiplying it by a diagonal matrix yield a scalar operation. With non-diagonal, you allow cross-terms.

For structured signals,  or images, or transformed signals (in a wavelet or a time-frequency domain), neighboring samples or coefficients often possess some local correlations, that you want to exploit in denoising. So assume you have an observation $r(k)$ at location $k$ (could be multidimensional), and look for an estimator $\hat{s}(k)$ of some underlying signal $s(k)$, under a noisy model  $r(k) = s(k)+n(k)$. In the original domain, you could for instance apply a scalar mask $f(k)$ on the observation: $\hat{s}(k) = f(k).r(k)$, especially if you know that at some location you have signal, at other noise. Windowing is a typical example. More generally, you can apply a function $f$: $\hat{s}(k) = f(r(k))$. If you assume local correlations, you can borrow information for close samples. For instance, average filtering can be $\hat{s}(k) = \frac{1}{3}(r(k-1)+r(k)+r(k+1))$.

In other words, $f$ now apply over a domain, gathered in a vector $v(k)= [r(k-1),r(k),r(k+1)]$. It can be called Reference Observation Vector (ROV). Then  $\hat{s}(k) = f(v(k))$. So average filters are sort of non-diagonal processing.

Traditionnaly with wavelets or time-frequency, the transform compacts or sparsifies noisy data, and data features tends to become relatively higher than  the noise, compared to the original domain. So the previous idea of a scalar mask works better. You can choose $f$ to keep the value of $r(k)$ if $|r(k)|$ is above a threshold $t$, and set it to zero if below. This $f$ is discontinuous.  let us write it $f(x) = \mathrm{sign}(x) (|x|-t)_+$
This is called scalar wavelet shrinkage or thresholding, scalar because you apply it on a single wavelet coefficient, illustrated in the center plot.

This idea was already used in the Fourier domain for noise reduction in speech. You can have smoother transistion, if you attenuate a bit the higher values, which give you soft thresholding, and you have many more shrinkage functions (garrote, firm) for scalar thresholding.

The next idea is to combine all of the above. You can modify one coefficient with respect to a certain norm of a surrounding $\|v(k)\|$, as:

$\hat{s}(k) = f(\|v(k)\|)r(k).$.

In wavelets, the ROV $v(k)$ can be more complicated. It can gather neighbors of the coefficients, but also coefficients from other scales.

The rationale is that, if a coefficient is small, it might be useful to preserve it anyway. If its ROV has high norm (think of energy), the coefficient might small by mistake, and is not noise.

The nice think is that the Stein Unbiased Risk Estimator theory allows you to estimate something on the signal $s(k)$ that you do not know by using only the observation $r(k)$, using some technical conditions, which allows you to find, automatically, the best threshold $t$, or shape paramaters of the function $f$. Such methods yielded almost the state-of-the-art on wavelet denoising. I am not familiar with audio. However, for multispectral images, here is an example:

a satellite image of the port of Tunis:

its noisy version, with heavy Gaussian noise:

the denoised version:

While not perfect, it was able to recover details, almost invisible from the noisy image, with very few parameters to set (only the shape of the ROV). This answer is a short version of what you can read in  A Nonlinear Stein-Based Estimator for Multichannel Image Denoising, 2008.

Why can  block denoising be better that scalar?  Because scalar is a specific case of block (non-diagonal) processing, if you take the ROV to be restricted to the coefficient $r(k)$ itself. If you choose the time-frequency or time-scale method well, and the other parameters (block shape and size), the block processing is generally (much) better.


This is the core subject of linear modulation theory that's considered in analog communication systems. It's analysis is based on the modulation property of the Fourier trabsform.

Given a bandlimited baseband signal $x(t)$ such that its CTFT $X(\omega)$ is zero outside of the band  $-W < \omega < W$, its multiplication with the signal $\cos(\omega_0t)$ has the following spectrum:

$$ x(t) \cos(\omega_0 t) \longleftrightarrow \frac{1}{2\pi} X(\omega) \star \mathcal{F} \{ cos(\omega_0 t) \} $$ where $\star$ stands for convolution.

The Fourier transform of $\cos(\omega_0 t)$ is
$$  \cos(\omega_0 t) \longleftrightarrow \pi \delta(w-\omega_0) + \pi \delta(w+\omega_0)$$

Using the convolution property $ f(x) \star \delta(x-a) = f(x-a)$ i.e., signal shifts to the impulse position, we deduce the spectrum of the modulated signal as:

$$ x(t) \cos(\omega_0 t) \longleftrightarrow \frac{1}{2\pi}  \left( \pi X(\omega-\omega_0) + \pi X(\omega+\omega_0) \right) = 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0)$$

Now if you filter (convolve) this modulated signal with a (n ideal) BPF with a center frequency of $\omega_c = \omega_0$ and a bandwidth of at least $W_{BW} = 2 W$ then you would get the modulated signal spectrum alone. This can be shown by utilizing the filtering property of Fourier transform which says
$$ h(t) \star x(t) \longleftrightarrow H(\omega) X(\omega) $$

Now the defintion of the (ideal) BPF is such that its frequency response is:
$$
H(\omega) =
\begin{cases}
1 &, |\omega-\omega_0| < W \\
0 &, o.w. \\
\end{cases}
$$

Then the multiplication of $H(\omega)$ with $0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0)$ would yield:

$$ H(\omega) \left[ 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0) \right] = 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0) $$

Which means:
$$ \mathcal{F} \{ \text{BPF} \{ x(t) \cos(\omega_0 t) \} \}= 0.5 X(\omega-\omega_0) + 0.5 X(\omega+\omega_0)$$

In this example the whole spectrum is the modulated spectrum so you don't actually filter out anything, but in other situations the BPF would indeed perform some useful operations.


Transfer function: from poles and zeros to polynomial coefficients

In short: $(0.937)^2=0.878$

A bit more detail:

$$\begin{align}(z-re^{j\theta})(z-re^{-j\theta})&=z^2-r\left(e^{j\theta}+e^{-j\theta}\right)z+r^2\\&=z^2-2r\cos(\theta)z+r^2\end{align}$$


I find it surprising that it didn't matter where you started.  Normally incorrect symbol alignment will completely change your results.  The usual method for determining the correct alignment (assuming you don't have other information like a Frame Alignment Word to help you figure it out) is to try both possible phases and see which produces a lower path metric.

You can, by the way, do the same when the code is punctured, you just have more phases to try.


This is not in matlab but in R.  It seems to do what you need.  The plot below mimics your original in colors (red is the original; circles are peaks and troughs; orange are the interpolating points; and blue is the interpolating signal).

R Code Below

#31438

#install.packages("quantmod")
require(signal)
require(quantmod)
T <- 100
filt <- butter(5,0.1)
red_line <- filter(filt, rnorm(T, 0, 1))
peaks <- findPeaks(red_line)
valleys <- findValleys(red_line)

plot(as.numeric(red_line), col="red", type="l")
points(peaks + 1, red_line[peaks], col="red")
points(valleys + 1, red_line[valleys], col="red", lwd = 3)

peak_idx <- 1
valley_idx <- 1
counter <- 1

values <- array(0, c(2,T))

while ( (peak_idx <= length(peaks)) && (valley_idx <= length(valleys)))
{
  mid_idx <- floor((peaks[peak_idx] + valleys[valley_idx]) / 2 )
  mid_value <- (red_line[peaks[peak_idx]] + red_line[valleys[valley_idx]])/2
  points(mid_idx, mid_value, col="orange", lwd = 5)

  values[1,counter] <- mid_idx
  values[2,counter] <- mid_value
  counter <- counter + 1

  if (valleys[valley_idx] > peaks[peak_idx])
  {
    peak_idx <- peak_idx + 1
  }
  else
  {
    valley_idx <- valley_idx + 1
  }
}

values[1,counter] <- T
values[2,counter] <- red_line[T]

interpolated <- spline(values[1,], values[2,], n=T)
lines(interpolated, col="blue")


What is the Effect of Multiplying a Function by the Unit Impulse Function in the Frequency Domain?

I think there is a slight typo in robert bristow-johnson's answer. Should be

\begin{align}
x(t) &= \frac{1}{2 \pi} \int\limits_{-\infty}^{\infty} X(\omega) e^{j \omega t} \, \mathrm{d}\omega\\
\\
 &= \frac{1}{2 \pi} \int\limits_{-\infty}^{\infty} H(\omega_0) \delta(\omega - \omega_0) e^{j \omega t} \, \mathrm{d}\omega\\
\\
 &= \tfrac{1}{2 \pi} H(\omega_0) e^{j \omega_0 t}
\end{align}

What this actually means is that during the inverse transform, $\delta(\omega - \omega_0)$ selects point frequency $\omega_0$ from the frequency response $H(\omega_0)$, and the value of $H(\omega_0)$ at that point determines the magnitude and phase of the time domain sinusoid $e^{j \omega_0 t}$.


I think the Wikipedia definition is limited, and anybody using the terms keying and line coding interchangably is wrong.

A line code is an encoding (or mapping) from bits to pulse amplitudes (voltages) whose main purposes are spectral shaping and receiver synchronization.

These are some of the things you can acomplish with line coding:

Remove DC. Since long transmission lines do not transfer DC very well, and since DC uses power but transmits no information, many systems use line codes that remove it. The simplest example is AMI (alternate mark inversion). In AMI, a bit 0 is mapped to 0 V, and a bit 1 is mapped to +1 V or -1 V in an alternating fashion.
Spread the signal energy across the entire bandwidth evenly. This is achieved with a scrambler.
Avoid long runs of consecutive '1's or '0's. This helps the receiver synchronize to the received signal; sync algorithmis typically rely on signal transitions. One example is the 8b10b code. This line code substitutes runs of  more than 5 symbols by sequences with transitions.

Note that the last two are different ways of achieving similar results.

It is true DC removal is more commonly used in baseband wireline systems. However, line codes, in one shape or another, are used in all communication systems.

Contrast line encoding with other types of codes: source coding compresses the information before transmission, and error control codes add extra information in order to correct and/or detect channel errors. A typical communications system will include all three types of codes.

Regarding keying: keying is the process of taking the output of the line coder and mapping it to a carrier. In this sense, keying does imply a passband channel.


Digital filter performance with noncontinuous signal?

It is not surprising that the filter output resembles the filter's step response at discontinuities of the input signal. It's like applying a (modulated) step at the input. Apparently the cut-off frequency of the high-pass filter is higher than the frequency of the sinusoid, so the output goes to zero after each discontinuity. However, at each point of a discontinuity, the input signal contains frequencies above the cut-off frequency of the high-pass filter and these frequencies are passed to the output.


Using the second half only (chopping) is the correct way. Note that setting the first half to zero is actually the same thing, up to a constant delay. The 100ms delay you observe might be due to the use of the "same" flag for the numpy convolution. It seems to me that the right truncation for this application would be to use numpy.convolve(signal, window, 'full')[:len(signal)] - without the kind of centering that 'same' does.

As you have noticed, the peaks in the smoothed envelope function do not exactly coincide with the peaks in the original envelope and there is no magic trick to make them coincide (the delay will depend on the sharpness of the attack) - so the resulting curve is not ideal for very accurate beat detection and chopping.


System Stability: Can we derive stability of a discrete system (Frequency domain, Z-transform) by applying analogous methods?

If you use the bilinear transform to map from s-plane poles and zeros to z-plane poles and zeros, the left-hand half-plane of the s-plane will be mapped to the inside of the unit circle on the z-plane. This preserves stability because on s-plane the requirement for stability is that poles are on the left-hand half-plane and on the z-plane it is that they are inside the unit circle.

I did not examine too carefully your s and z-domain formulas, but indefinite integration is marginally stable with a pole at the origin of the s-plane, or on the unit circle of the z-plane. Implemented as a nonrecursive filter:

$$\text{out}[i] = \text{in}[i-N + 1] + \text{in}[i-N + 2] + \dots + \text{in}[i-1] + \text{in}[i],$$

integration over $N$ samples is stable. Implemented as a recursive filter:

$$\text{out}[i] = \text{out}[i-1] + \text{in}[i] - \text{in}[i-N],$$

it is stable in fixed-point arithmetic, but can drift if implemented in floating point arithmetic due to accumulation of numerical error. Even then it can be stabilized by moving the pole slightly inside the unit circle, or by periodically resetting the filter state by calculating it using the nonrecursive filter. But sample and hold is about taking an input sample and keeping the output at that value for some time. That can be programmed directly and gives no stability problems.


understanding image error metrics, mean squared error

First of all, the mean squared error should not have a square root. As the name implies, it is the mean of the squared difference. What you are referring to is the root mean squared error. Second, the this is a pixel-wise comparison between to images. Therefore the number of pixels is the dimension, and the sum is performed on it. There is even a comment on that page about color images.


Consider a WSS (wide sense stationary) discrete-time random process (RP) with i.i.d. samples $X[n,s]$ whose auto-correlation sequence (ACF function) is $\phi_{XX}[m]$. We call a random process as white-noise if its ACF is the following:

$$ \phi_{XX}[m] = \begin{cases}
{\sigma_X^2 ~~~~~, \text{ for } m=0 \\ 0 ~~~~~~~~, \text{ otherwise} }\end{cases} $$

This also means that the RP is uncorrelated. If a RP does not have this property then it's a colored noise. And it means its samples are correlated.

When a RP is colored, it may be possible to find an LTI filter which will try to remove the correlation between the samples of the noise and hence make it whitened. Such a filter is called either as a decorrelator or equivalently as whitening filter.

It's best described in the context of AR-p (Auto-Regressive) random processes generation through all pole filtering, where the whitening filter will be the inverse of that all-pole AR filter which generates the colored noise from an input white noise.


What should go to infinity in power's formula N or the repetition KN?

For periodic signals you can use
$$P_y=\frac{1}{N}\sum_{n=0}^{N-1}|y[n]|^2\tag{1}$$
where $N$ is the period of the signal. I'm sure you can take it from here.
Also take a look at this question and its answer where the equivalence of $(1)$ and the formula in your question is discussed.


Range and motion in received signals by radar

You will need more than 2 received pulses to get the Doppler between your sources. What is your ADC sampling rate and PRF?


Turn your guitar to MIDI and Guitar Synthesizer

There is a sensor for each string which individually records the string vibration. This signal is converted to digital, analyzed for estimating its fundamental frequency and its amplitude (envelope detector). Whenever a sharp increase in amplitude is detected, a MIDI note on message is sent, at the detected frequency. Whenever the amplitude falls below a threshold, a MIDI note off message is sent.


I ran your code and got the below plot as well as an almost zero vector of:

[-1.19209290e-07, -2.42083020e-09, -8.94069672e-08,  2.23517418e-07, -6.33299351e-08]

This means numpy successfully curve-fit a line to your data.

An order of 6 works a lot better, though:

[-2.79396772e-09 -2.04281037e-14 -1.27897692e-13  4.33431069e-12 -2.27373675e-11]

numpy won't give you an infinite number of solutions.  It will give you an optimally curve-fit solution to your data for the given polynomial order.


In principle you can do that. You should, however, note that the two noise signals (both of which you called $n(x,y)$) are of course not identical. The two models are

$$g(x,y)=h(x,y)*f(x,y)+n(x,y)\\
g(x,y)=h(x,y)*\left[ f(x,y)+z(x,y)\right]$$

where the two noise components $n(x,y)$ and $z(x,y)$ have different properties. So you have to be careful with methods that assume certain noise properties. E.g. if it is assumed that $n(x,y)$ is white, then $z(x,y)$ is not white anymore and this has to be taken into account in the noise removal algorithm. The basic question you have to ask yourself is what you gain by changing the model.


Online materials on Support Vector Machine

Video Lectures : http://nptel.iitm.ac.in/courses/117108048/

Notes : http://www.nptel.iitm.ac.in/courses/106108057/ (Read from Lesson 25. These notes helped me very much)


The problem with your example is that $\infty\cdot 0$ isn't necessarily equal to zero. The only way to judge what is happening in the limit $K\to\infty$  is to divide the original equation by $K$:
$$\frac{D(s)}{K}+N(s)=0\tag{1}$$
Now it is obvious that for $K\to\infty$ the actual value of $D(s)$ is irrelevant, as long as it is finite. Consequently, the only necessary condition for $(1)$ to be true as $K$ becomes large is that $N(s)=0$.
This problem is a bit similar to the problem of determining the limit
$$\lim_{x\to\infty}\frac{x}{x+c}\tag{2}$$
The value of the limit $(2)$ is independent of the choice of the constant $c$, which becomes negligible compared to $x$. The same is true in the first equation of the question: the actual value of $D(s)$ becomes irrelevant compared to $KN(s)$.


Is there any algorithm to decompose a wave into a set of given frequencies?

You could use the Goertzel algorithm to compute the discrete Fourier transformation for just the frequencies you want to detect.

The Goertzel algorithm is used to detect single bins from the DFT.  You implement a bunch of them, each set to one of your desired frequencies, and calculate them for each audio sample.

Anything that doesn't fit your selected bins is ignored.  If you implemented it for the frequencies given in your example (5Hz, 14Hz, 37Hz, 42Hz, and 59Hz) but your audio was a sine wave at 100Hz, then your software could not detect that 100Hz signal.

The Goertzel algirithm can be "tuned" to make the bins narrower.  You would probably need to do that to separate the frequencies given in your example.  They are very close together, so you would need to use very large $N$ to calculate the coefficients for each bin.  This has the side effect of making the algorithm respond slowly - quickly changing frequencies would not be detected properly.


What is the amplitude of the impulse resonse

In case of discrete signals, the amplitude of unit impulse is 1.

It does not mean that your impulse response should have unit amplitude. The impulse response can have any real valued amplitude in it.

Note: The area under unit impulse in continuous time domain is 1 not the amplitude. The amplitude of a unit impulse approaches infinity as the area under remains constant (=1).


Theoretically it's possible to do this, although it will often not be practical.

Let's consider this in polynomial space. For a filter of order N you have 2*N+1 independent variables (N for the denominator and N+1 for the numerator). Let's look an arbitrary point $z_{k}$ in the z-plane and let's say the value of the transfer function at this point is H($z_{k}$). The relationship between the transfer function and all the filter coefficients can be written as equation that's linear in all filter coefficients as follows:
$$\sum_{n=0}^{2*N}b_{n}\cdot z_{k}^{-n} - H(z_{k})\cdot \sum_{n=1}^{2*N}a_{n}\cdot z_{k}^{-n}=H(z_{k})$$
So if you pick M different frequencies $z_{k}$ you'll end up with a set of M complex linear equations or 2*M real equations. Since your number of unknowns is odd (2*N+1) you probably always want to pick one frequency where z is real, i.e z = 1 or $\omega$ = 0.

If M is larger than N than the system of equations is linearly dependent. You can find the filter order by starting at N=1 and increase N until the equation system becomes linearly dependent. The largest N at which the system is linearly independent is the actual filter order. For this approach, it doesn't even matter what frequencies you pick. As long as they are different, any set of frequencies will work.

However, this a numerically very tricky problem. Polynomial representation for larger filter orders is numerically very fragile and the smallest amount of noise or uncertainty leads to very large numerical errors. For example, if you determine the values of the sampled transfer function through measurement, the required measurement accuracy will be prohibitive unless it's very benign low order filter.


What is the difference between i.i.d noise and white noise?

Independence is 'stronger' than whiteness. I believe that independence between the random variables implies whiteness but whiteness does not imply independence. Whiteness means that the random variables are uncorrelated but not necessarily independent.

The use of i.i.d noise is seen very often when formulating probabilistic models because it makes inference much easier. For instance if two random variables $X_1$ and $X_2$ are i.i.d it means that the joint pdf $p(X_1,X_2)$ factors into the product of the individual pdfs $p(X_1)p(X_2)$. If the two random variables are uncorrelated this factorization is not valid.

In ML estimation typically the log of the product is considered $\log(p(X_1)p(X_2)) = \sum \log p(X)$ because then differentiation with respect to the parameter of interest is much more straight forward.

I'm not sure where to find i.i.d noise in the real world but I believe that the assumption about i.i.d observation noise is made more of convenience than because it is realistic.


Using FFT to convert to frequency domain, then IFFT back to time domain C++

A couple of things:

Your Data[] array seems to start at on offset of "1". That's unusual for C++, where "0" is the first index. Matlab starts at 1. This may be part of the calling convention for fft_basic, but it's a bit strange
Doing the inverse Fourier transform of the Power Spectrum (magnitude squared of the forward transform) gives you the auto correlation. That has nothing to do with phase compensation at all. Wrong algorithm.
You cannot calculate the phase as atan(imag/real) since you'll get half of the quadrants wrong. The correct functions is atan2(imag, real) which takes into account the relative signs of real and imaginary part.
Just adding to the choir: solve it in Matlab first and then port to C++


How @Pichenettes said, we need more information to help properly

I'll take a guess:

If your Frame Rate are 44100, you will have 1587600 samples for 36 second of audio (36*44100), for the calculation example:

Pxx, freqs, bins, im=pylab.specgram(sound_info, Fs=frame_rate)

pylab.specgram seems use FFT Size = 256 and Overlap = 128 (50%) its give me the follow dimensional array:

floor(1587600 /  128) = 12403

(128, 12403)

You need set properly FFT Size and Overlap in your souce code.


Why do we need 4 points for homography but 7/8 points for fundamental matrix calculation?

It is because in the case of fundamental matrix, each correspondence point relates to only one constraint(i.e it maps a point to a line in other image). Hence 8 correspondence points are required.

But in the case of homography, each correspondence solves two constraints. Hence only 4 correspondence points are sufficient.


Why a Convolution Matrix Is Called a Kernel?

In general, a discrete (1D) transform can be defined as:

$$ y[k] = \sum_n \phi[k,n] x[n] $$

where the input function $x[n]$ is mapped into the output function $y[k]$, by the kernel, $\phi[k,n]$, of the transform. For example, the kernel for DFT (discrete Fourier transform) is  $\phi[k,n] = e^{-j \frac{2\pi}{N}kn}$.

In essence the kernel acts on the input and produces the output by modifying it, like in the sum above. Therefore the kernel defines the character of the mapping.

Programmatically, a kernel is the core portion of a function which accepts a sequence as its input and returns a computed output whan called.

The input-output relationship for a discrete-space 2D LTI system (typically an image filter) is expressed as a convolution sum :

$$ y[n_1,n_2] = \sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} h[n_1-k_1, n_2-k_2]x[k_1,k_2] $$

The convolution operation is similar to the transform considered above with its kernel being the impulse response $h[n_1, n_2]$ of the LSI system.

For a discrete-time 2D LSI system whose impulse reponse $h[n_1, n_2]$ has a finite domain of support; i.e, finite size, then programmatically the system can be represented as a matrix as well. And the convolution operation than is performed by calling a function which contains the impulse response matrix as its kernel and producing the output pixels per call. And here is the most typical meaning of a (filter) kernel in image processing.

Furthermore, a matrix called mask operates on a block of pixels (sample wise) on the image being processed. They are more of a programming concept than of signal processing, but also find applications such as windowing or frequency domain filtering. Masks find extensive usage for image processing effects. Their operating matrix can also be called as a kernel.


How to show Image FFT with origin as DC?

The quadrants of the FFT should be shifted Q1->Q3, Q2->Q4 ...
For a [RxC] matrix X this may be accomplished by shifting as:

tmpX = shift_row(    X, floor(R/2) )
Y    = shift_col( tmpX, floor(C/2) )

where shift_row and shift_col is your shift operator in each dimension and floor is the round towards zero operator.

For the inverse FFT shift you should shift in the same way but with a value -ceil(R/2) and -ceil(C/2), i.e rounded towards infinity


In the general case you have

$$H(z)=\frac{P(z)}{Q(z)}$$

where $P(z)$ and $Q(z)$ are polynomials in $z$. If - as is the case in your example - $Q(z)$ just has one single term, $H(z)$ is definitely FIR, because you can simply divide each term of $P(z)$ by that respective power of $z$, and the number of terms of $H(z)$ equals the number of terms of $P(z)$.

If $Q(z)$ is a polynomial with finite zeros away from $z=0$, then $H(z)$ is generally IIR, except if all zeros of $Q(z)$ are cancelled by zeros of $P(z)$. One example of such a case is a recursive moving average (where $N$ is the window length):

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

The zero of the denominator is cancelled by one of the zeros of the numerator, so the transfer function $(1)$ can equialenty be written as

$$H(z)=\frac{1}{N}(1+z^{-1}+\ldots+z^{-N+1})\tag{2}$$

which is obviously FIR.


Your understanding is correct. If you sample at rate $f_s$, then with real samples only, you can unambiguously represent frequency content in the region $[0, \frac{f_s}{2})$ (although the caveat that allows bandpass sampling still applies). No additional information can be held in the other half of the spectrum when the samples are real, because real signals exhibit conjugate symmetry in the frequency domain; if your signal is real and you know its spectrum from $0$ to $\frac{f_s}{2}$, then you can trivially conclude what the other half of its spectrum is.

There is no such restriction for complex signals, so a complex signal sampled at rate $f_s$ can unambiguously contain content from $-\frac{f_s}{2}$ to $\frac{f_s}{2}$ (for a total bandwidth of $f_s$). As you noted, however, there's not an inherent efficiency improvement to be made here, as each complex sample contains two components (real and imaginary), so while you require half as many samples, each requires twice the amount of data storage, which cancels out any immediate benefit. Complex signals are often used in signal processing, however, where you have problems that map well to that structure (such as in quadrature communications systems).


Why not use zero-signal in digital transmission?

All implementation aspects aside, the constellation you propose performs worse than QPSK in an additve white gaussian noise (AWGN) channel.

I claim this based on simulations that I have run with Matlab calculating the symbol error rate (SER) as a function of signal-to-noise ratio (SNR). Here is the result:

As you can see, for a given SNR, the proposed constellation has an higher SER than a QPSK constellation. Additionally, I expect the proposed scheme to perform even worse in terms of bit error rate (BER) as -- in contrast to QPSK -- Gray labelling is not possible. Some details concerning the simulation:

1 million bits per SNR value
maximum likelihood decision
additive, Gaussian distributed noise
memoryless channel (no inter-symbol interference)

It might well be that I made a mistake writing the simulation script, so I uploaded it in case you would like to check my implementation.

This might not be the mathematically rigorous explanation that you have been hoping for. To answer your question theoretically, one would have to calculate the SER analytically. For that, one should first find the decision thresholds, which is the main problem (for me, at least). It's far from obvious, what the optimal decision thresholds would be. They are not straight lines as with QAM constellations but probably a kind of zero-centered disc with three "rays" emerging from its center. Calculation of the SER should probably involve a transformation of the problem to polar coordinates.

With that kind of irregular constellation it also doesn't help much to compare the minimum distance, i.e. the length of a straight line between two given constellation points. It's true, that for a given mean signal energy, the outer points of the proposed contellation have greater distance to the inner point than any two points of the QPSK constellation. However, the inner point of the proposed contellation is surrounded by three other points which  decreases its decision area and consequently increases the symbol error probability for that point. It's not obvious what is better. According to the simulation, QPSK is better.

Update

The areas of the maximum likelihood decision can also be simulated. In the following image, each color corresponds to a distinct, decided symbol. The input to the decision device is a very noisy (0 dB SNR) signal with the proposed constellation.

We can clearly see how the inner point is "squeezed" between the outer points. In my opinion, that is the reason why the proposed constellation performs worse than QPSK. Now that I see this plot, I find the SER calculation not so difficult anymore. With a little time it could be done.


Tikhonov Regularization Alternate Formulation

Just expand the term to be minimized, take the derivative with respect to $x$ and set it equal to zero:

$$\begin{align}||Ax-b||_2^2+\lambda^2||x||_2^2&=(Ax-b)^T(Ax-b)+\lambda^2x^Tx\\&=x^TA^TAx-2b^TAx+b^Tb+\lambda^2x^Tx\tag{1}\end{align}$$

Taking the derivative w.r.t. $x$ and setting it equal to zero gives

$$2A^TAx-2A^Tb+2\lambda^2x=0\tag{2}$$

from which the given equation follows.


Grayscale images are typically stored as 8 bits per pixel in files. Perhaps the network was trained to work with such data. $2^8 = 256,$ so 8 bits can represent any number from 0 to 255. It may well be that the author was converting from range $[0, 1]$ to $[0, 255]$. Note that $0.8681 - (-0.1319) = 1.$ Multiplying by 255 to stay within the range would have been more appropriate, but multiplying by 256 is appealing because it does not introduce any rounding error in finite precision floating point arithmetic.

Subtracting the mean does not change the story as it can be equivalently done before or after the multiplication.


Phase of signals in real and complex DFT

And all of these waves have phase equal to 0.

No. The phase will be zero for $X[0]$ (DC -zero frequency component) but the phase is given by $arctan(Im[i]/Rx[i])$, which usually doesn't result in zero.

Yet if we use the forward DFT formula
$$X[r] = {\sum}_{k=0}^{N-1} x[k]e^{-j(2 \pi k r / M)} $$
to compute ($M$ is usually considered to be equal to $N$), the output
is a sequence consisting of $N$ complex numbers, for a input sample
with size $N$.

N output complex numbers requires just as many input complex numbers (which means you need N real numbers and N imaginary numbers (all zeros) as input, which together make up the complex input).
If, the input is ,however, purely real (N real numbers), and you don't provide an imaginary part, your real fft algorithm will simply treat your input as complex (for instance, the even numbers will be treated as real and odd as imaginary) but now you'll get only $N/2+1$ complex numbers at the output.
I

s it the case that the frequency components have amplitude $A =
> \sqrt{a^2 + b^2}$, phase $\varphi = \arccos{(\frac{a}{A})}$ and they
are shifted sine waves?

Yes, but the the phase fomula is actually $\varphi = \arctan{(\frac{b}{a})}$
EDIT: As mentioned in a comment below, one generally needs to use $atan2$.


What is the intuition of "averaging is a low pass filter"?

Consider a sequence of numbers

x = [1, 100, 1, 100, 1, 100]

If we send this sequence through an averaging transformation

y(i) = 0.6*x(i) + 0.4*x(i-1)

we get

y = [1, 60.4, 40.6, 60.4, 40.6, 60.4]

Original signal x had a lot of sudden changes 1 -> 100, 100 -> 1 where as signal y has smaller changes 60 -> 40, 40 -> 60

When we say that a signal has high frequency components we mean that the values change rapidly with time. So x had rapid changes in amplitude, while y does not have that much of rapid changes in values. This is the intuition behind why averaging is equivalent to low-pass filtering (disallowing high frequencies).

--

Another way to look at it

x = [1, 100, 1, 100, 1, 100]

could be written as sum of two sequences, x1 and x2

Low(zero?) frequency component

x1 = [50, 50, 50, 50, 50, 50]

High frequency component

x2 = [-49, 50, -49, 50, -49, 50]

After averaging transformation, we get

y = [1, 60.4, 40.6, 60.4, 40.6, 60.4]

Which has low frequency component of

y1 = [50, 50, 50, 50, 50, 50]

and high frequency component

y2 = [-49, 10.4, -9.4, 10.4, -9.4, 10.4]

As you can see that y2 has much smaller absolute amplitudes than x2, meaning that we have diminished the high frequency components in x by doing averaging.


Sampling multiple sensors at different rates

My suggestion is to use a single Arduino for all the sensors, create timestamps for the data and transmit it via I2C to whatever sink you want. Your sampling rates are so low, that a single uC should do the job.
Then do a generous oversampling to make data from different sensors easy to synchronize.


There is a trade-off between modulation order, required energy and bit-error rate (BER). As you increase the modulation order, you need to increase the average energy per bit to keep the same BER.

As an easy example, consider BPSK vs 4-PAM. Say you have BPSK with distance between constellation points $d_{min}=1$; this works out to an average bit energy of 0.25 joules per bit.

Now, you want to double the bit rate by going to 4-PAM, but you want to keep roughly the same BER, which means having the same minimum distance $d_{min}$ between constellation neighbors. Now, your average symbol energy is $$\overline{E_s}=\frac{2\cdot1.5^2+2\cdot0.5^2}{4}=1.25$$ joules, or an average bit energy $\overline{E_b}=0.625$.

As you increase the modulation order, you'll see the energy required increase exponentially. Exactly the same thing happens when using QAM.


Have a look at the Python code below. It calculates the auto-correlation of the data, after the DC part of the signal was removed. The autocorrelation shows nice peaks with much less noise than the original signal. Then, we detect the local maxima and measure the distance between them. This distance is then the estimated period of your signal:

import csv
import requests
import numpy as np
from io import BytesIO
from scipy import signal

url = 'http://pastebin.com/raw/5ghKbBRm'
response = requests.get(url)
A = np.genfromtxt(BytesIO(response.content)) # read the CSV file

# begin signal processing

A = A - np.mean(A)  # remove the DC part of the signal
plt.subplot(121)
plt.plot(A)

plt.subplot(122)
C = np.correlate(A, A, 'same')  # calculate autocorrelation to smooth out the noise
plt.plot(C)

plt.subplot(122)
P = signal.argrelextrema(C, np.greater)[0]  # caluclate local maxima (this step can be improved)
P = P[C[P]>15000] # filter maxima, where the autocorrelation values is too low
plt.stem(P, 30000*np.ones_like(P))  # plot where the detected peaks are

print ("Peak distances: ", np.diff(P))  # calculate distance between peaks and the mean of the distances
print ("Estimated period: ", np.mean(np.diff(P)))

Program output:

Peak distances:  [117 118 117 117 117 117 118 117]
Estimated period:  117.25

There is room for improvement. For example, the detection of local maxima is not very robust and the threshold was set manually, but this should enable you to find out the rest on your own.


Frequency Domain Equivalent to Zero-Phase Filtering (`filtfilt`)

In essence filtfilt does the the following
$$y[n] = x[n] \ast h[n] \ast h[-n]$$
It's convolution with the impulse response and then the time reversed impulse response (in any order since the convolution is commutative).
So in the frequency domain it's simply
$$Y(z) = X(z) \cdot H(z) \cdot H'(z) = X(z) \cdot |H(z)|^2 $$
which obviously is zero-phase since $|H(z)|$ is real.
For an FIR filter that's easy enough: If your filter length is $N$, the cascade length is $2*N-1$ and your FFT length needs to be at least $4N-3$. Take your impulse response, zero pad it to the FFT size, take the magnitude squared and then mutliply with $e^{-j2\pi n/N} $ to get it circulated to the proper time location. Then proceed as with any other overlap add.
Alternatively you can just calculate the aggregate impules respones as
$$ h_{agg}[n] = h[n] \ast h[-n] $$
IIR filters can be done using the same approach but it requires turning the IIR into an FIR by suitable truncation of the impulse response. filtfilt() does that implicitly as well but uses the length of the signal as the truncation length.


Is the Hilbert transform computable "losslessly"?

The loss you are seeing is due to the use of a finite length rectangular window.  If your signal is long enough that you can ignore the Hilbert filter transients at the both ends of the window (e.g. if it is infinite in length or so long that you don't care about the ends below some noise floor), then the process won't have the particular visible losses that your example points out.