What does a star shaped autocorrelogram mean?

I'll go ahead and post a practical answer for you, and we can both see if someone can post a theoretical explanation.

The shape you see comes from the offset in your data.  If you were measuring voltage, I'd say you have an AC signal on a DC offset.

I don't know what you are measuring, but it spends most of its time at around 100.  It has short intervals where it changes drastically, and makes small changes almost constantly.

You have two ways to get rid of the triangle.

Take the average of all data points and subtract that average from all the data points.
Use a high pass filter to remove the offset.

If this were an audio signal, you could pick a lower cutoff frequency and be happy - there's some point where the audio frequencies are too low to be interesting.

Since you don't mention the source of the data (could be the sum of loans and borrows of elephants in an elephant lending library for all I know,) I'd go with the first option.

If it is some physical quantity, then you may want to consider the high pass filter, after all.

If the process generating the signal has some natural low frequency, then you could use that as your cutoff.

Or, you set the cutoff by the longest time period you want to consider.

If your data set represents 1 second, then autocorrelation for more than 1second is pointless - you could set the cutoff for one hertz.   At any rate, you would use 1/(time period) as the cutoff frequency.

If you have all the data in blocks at hand, you can use filtfilt to do your filtering.  That preserves the phase and the absolute timing of your signal.


The best performance you will be able to squeeze out of a PC is with an audio interface that supports ASIO.

To be able to get near 10ms of total delay (input and output) you will also need a very fast machine. Not only in terms of CPU but more importantly in terms of memory performance.

On top of this delay, you are going to have to add the delay in processing your signal and generating the output signal.

1ms (total delay (?) ) on the hardware you are describing is unrealistic.

The way ASIO works will handle issues #1,#2. I am still not really sure what you mean by #3. If you want to use the PC to drive a device with some signal and collect the response of the device to that signal then that is relatively easily achievable with an ASIO enabled audio interface too.

Hope this helps.


Signal Processing and Emotion Recognition

For the latest improvements read the papers from the latest speech conferences, for example here are two tracks on Emotions in latest ICASSP:

Emotion recognition I

Emotion recognition II

You can also check reviews:

Speech emotion recognition: Features and classification models

Emotion recognition from speech: a review

Survey on speech emotion recognition: Features, classification schemes, and databases

To find the relevant research use Google Scholar.


Identify large solid color areas of image

Perhaps you could classify the difficulty by assuming that higher spatial frequency implies higher difficulty.  Make sure all the images are the same size, do a 2D FFT, and look at the spatial bandwidth in x and y dimensions, would be my suggestion.


Find the relation between two noisy periodic signals of equal length

All right, we can cope with this one.

So first, we should like to have some knowledge about the relation of  $i$ and $\phi$,
$$
  \phi\colon\ \mathbb{Z} \to S^1,\qquad i\ \mapsto\ \phi_0+\frac1{\nu_{\!s}\cdot\omega_0}i
  =: \varphi + \lambda i\quad {}_{^{(\operatorname{mod}2\pi)}}
$$
Where we can essentially forget about the $\operatorname{mod}2\pi$, just view it as $S^1=\mathbb{R}/\!{\sim}$.
There are multiple ways to determine the most likely values of $\varphi$ and $\lambda$; what I'd do is a least-square fit in those variables and $A_M$, perhaps also an offset $A_O$ such that
$$
  \zeta := \sum_i \bigl(A_{\!M}\cdot\sin\,\phi(i) + A_O - A_i\bigr)^2
$$
is minimized. A fitting algorithm like this can be somewhat slowish if you want it to work reliable, I don't know if that's a problem.

Once we have this relation $\phi(i)$, we can just read out the $B_i$ for which $i$ is close to $\tfrac\pi2$ resp. close to $-\tfrac\pi2$. We can simply average over all next-neighbours to these real values, this way we're quite immune from nonlinear response, but rather more subject to noise than with using some kind of interpolation. You want to do it this way if you have a very large amount of data, spanning many periods of $\phi$. If you don't, better use some interpolation between the $B_i$ with $\phi(i)\approx\pm\tfrac\pi2$. I'd try a Friedrichs function as convolution kernel.


Is $y[k] = y[k-1] + x[k]$ an integrator?

The system

$$y[n]=y[n-1]+x[n]\tag{1}$$

is an ideal accumulator, i.e., it computes the cumulative sum of the input samples:

$$y[n]=\sum_{k=-\infty}^nx[k]\tag{2}$$

It is in a way analogous to a continuous-time integrator, but this doesn't mean that you will necessarily obtain an ideal integrator by transforming the discrete-time system to a continuous-time system. There are several methods to do the conversion, and I'm not sure which one you used when calling d2c.

In any case, the properties of the transformed continuous-time system not only reflect the properties of the original discrete-time system, but also the properties of the transformation. So the properties of a discrete-time system should be investigated directly in the discrete-time domain, not by transforming it to a continuous-time system.


Sequence Length for a Linear Feedback Shift Register

The duration of the linear feedback shift register (LFSR) in number of "chips" is $2^N-1$ where N is the number of states in the shift register (and order of the generator polynomial), and chips refers to each unique output that is generated in the sequence. This is assuming that you are using a maximum length sequence in the implementation of your shift register, which means the generator polynomial is both primitive and irreducible in GF(2).

Primitive: Repeated multiplying by x (representing a shift of one position in the shift register) will generate every possible state (except for the all 0 state). If the generator is not primitive, possible states (combinations of the register values) will be skipped and the sequence will be less than $2^N-1$ in length before repeating.

Irreducible: Cannot be factored, see example:

Reducible: $1+x+x^2+x^3 = (1+x^2)(1+x)$

Irreducible: $1+x+x^3$

Note the above polynomials are in GF(2) so the only elements are 0,1. 1+1=0, 0-1=1, etc

All primitive polynomials are irreducible but the converse is not true!

Primitive and irreducible polynomials are well tabulated. For a handy table of primitive and irreducible polynomials in GF(2), see Peterson's Table of Irreducible Polynomials

The graphic below shows the two common implementation approaches for the LFSR from the generator polynomial, in this case N=4, using generator polynomial $1+x^3+x^4$ which would generate $2^4-1=15$ chips before repeating:


Negative frequency in the Fourier Transform

I think you may be confused about the use of negative frequencies and the sign of the integration coefficient.

Let's take a look at the inverse Fourier Transform

$$ x(t) = \int_{-\infty}^{+\infty}F(\omega)\cdot e^{ 2 \pi i  \omega} d \omega $$

That basically means that you can construct any time domain signal as the from a set of complex exponential. So you can break down a complicated signal into a number of constituent parts that are hopefully easier to work with.

That's typically called an orthogonal base. Complex exponentials are not the only choice, you could use a orthogonal polynomials, wavelets, etc. Complex exponentials are popular since they are easy to work with, specifically when you are dealing with linear time invariant systems

Note two things

The integration goes from $-\infty$ to $+\infty$, i.e you need both negative and positive frequencies to generate any possible time domain signal
The integral variable is positive, but that does only mean than you first integrate the negative frequencies and then the positive ones.

Of course you need to know what the weights for the complex exponentials, so you need to invert the formula above. You get

$$F(\omega) = \int_{-\infty}^{+\infty}x(t)\cdot e^{ -2 \pi i  \omega} d t$$

That's just how the math works out. The integration variable is negative, but that doesn't mean anything. The formula needs to produce results for both positive and negative values of \omega since we need both positive and negative frequencies to represent a time domain signal.


VHDL fixed point implementation of IIR filter

Filter output has to be divided by the first term of denominator.
In this case, division by 16384 can be effectively done by logical shift 14 bits to the right.


Exponential weighted moving average time constant

If I understood you correctly, you want to compute the value of $\alpha$ that results in a specified 3dB cut-off frequency for an exponentially weighted moving average filter. If you square your last equation, you get

$$\frac{\alpha^2}{1-2(1-\alpha)\cos(\omega_c)+(1-\alpha)^2}=\frac12\tag{1}$$

which can be rearranged into the following quadratic equation:

$$\alpha^2+2\alpha(1-\cos(\omega_c))-2(1-\cos(\omega_c))=0\tag{2}$$

with the positive solution

$$\alpha = \cos(\omega_c)-1+\sqrt{\cos^2(\omega_c)-4\cos(\omega_c)+3}\tag{3} $$

So, e.g., for a desired cut-off frequency $\omega_c=0.2\pi$, you obtain from $(3)$ a value of $\alpha=0.455886780102867$. The figure below shows the magnitude of the frequency response of the resulting exponentially weighted moving average filter, from which you can see that the desired cut-off frequency is achieved.

EDIT: The formula for $\alpha$ in your question should actually be

$$\alpha=1-e^{-T/\tau},\qquad \tau=1/\Omega_c\tag{4}$$

Note that unlike $\omega_c$ in Eq. $(3)$, $\Omega_c$ in Eq. $(4)$ is not normalized by the sampling frequency. So we have $-T/\tau=-\Omega_cT=-\omega_c$.

Eq. $(4)$ is an approximation, and it comes from applying the impulse invariant transformation to the continuous-time transfer function

$$H(s)=\frac{1}{1+s\tau}\tag{5}$$

which has a 3dB cut-off frequency $\omega_c=1/\tau$. Applying the impulse invariant transformation to $(5)$ gives

$$H(z)=\frac{T}{\tau}\frac{1}{1-e^{-T/\tau}z^{-1}}\tag{6}$$

Comparing the denominator of $(6)$ with the denominator of the discrete-time transfer function of an EWMA filter

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

results in the given formula.
Note, however, that this is only an approximation. Especially for cut-off frequencies close to Nyquist, the error of formula $(4)$ becomes relatively large.


Determine reflections from received signal

This is a classical deconvolution problem in that your received signal consists of multiple delayed and scaled copies of your transmitted signal, no different than what could be modeled with an FIR filter (which is the sum of delayed and scaled copies of the input signal), and your goal is to reconstruct the original transmitted signal from the received signal.

The received signal is the convolution of transmitted signal with the effective channel given by the "filter" described above which in this case would be the multiple different reflections from the surface of the wall. So we have $r[n] = t[n]\ast h[n]$ (received signal $r[n]$ is the convolution of the transmitted signal $t[n]$ and the channel $h[n]$), with two knowns $r[n]$ and $t[n]$ and we wish to solve for $h[n]$. This is a classical deconvolution problem typically solved using the Wiener-Hopf equations resulting in a least squared solution when there are sufficient samples for an overdetermined equation. This is explained in detail along with MATLAB code at this post here which was specific to an audio signal but would equally apply to the OP's case of a broadband microwave signal from 3.2 GHz to 4.8 GHz as he explained in comments:

Compensating Loudspeaker frequency response in an audio signal

As described above this would provide for a solution to the actual channel modeled as an FIR filter, but the process can equally be used to solve for an equalizer $c[n]$ which would optimally in a least squared sense return the received signal to the closest copy of what was transmitted simply by swapping $t[n]$ and $r[n]$ since in the receiver, the copy of the transmitted signal $t[n]$ is $r[n]\ast c[n]$.

Some precautions when using "sounding" signals to determine a channel:

The transmitted signal would ideally be white over the frequency band of interest. Certainly if there is any frequency portion with no signal content, no solution can be determined for that portion of frequency (since there would be no input/output relationship). If the signal is not white, the least squared solution would be weighted according to the respective spectral density over the entire bandwidth of interest.

If the delay spread is long compared to the inverse of the bandwidth (so becomes a greater problem for wideband signals), deep nulls can occur in the received spectrum due to the cancellation effects of the various delay paths. This is known as "Frequency Selective Fading" and makes equalization and deconvolution challenging as noise enhancement would occur at the locations where deep spectral nulls exist. The received spectrum should be reviewed to confirm that such deep nulls do not exist. For communication systems this is typical motivation for use of non-linear equalizers such as decision-feedback equalizers which would perform better for equalization under such conditions. (see non linear equalizer vs linear equalizer)

In this related post I showed how the channel estimation formulas can be used to determine the group delay of a channel:

How determine the delay in my signal practically


Yes you are violating! But not the Nyquist rate as a matter of fact :))

First of all Nyquist (Shannon) sampling theorem is based on the analog signal's bandwidth, and states that you must sample the signal at a rate twice the bandwidth of it. Now, what's the bandwidth of your signal? You must compute this one first.

Given a DAC (digital to analog converter) modeled as D/C block; i.e., an ideal signal reconstructor with an ideal interpolation (image rejection) filter with $h_r(t)$ whose cutoff frequency is at a half of the DAC sampling rate, then the output of the DAC given at a sampling rate of $1$G samples per second (1 GHz) will be bandlimited to $500$ MHz.

Now given a continuous-time (analog) bandlimited baseband signal $x(t)$ with a bandwidth of $500$ MHz, the required Nyquist sampling rate will be twice the bandwidth of the signal which is $1$ GHz, or stated as $1$ G samples per second.

Hence in principle (ideally) there's nothing wrong with (Nyquist is not violated by) using a DAC-ADC communication pair with the same sampling frequency. Of course due to several nonidealities (please refer to @MarcusMüller 's answer for those) you will find it quite impractical to achieve. Anyway, Nyquist is ok ;-)

But, then, you may be confusing the sampling theorem with channel capacity (information rate) for AWGN channel. So given your DAC-ADC AWGN channel system with a bandwidth of $B = F_s/2$, and given an SNR of the line in between; your maximum rate of information transfer is given by Shannon-hartley law which states that $C = B \times \log_2( 1 + SNR)$ b/s...

Assuming an SNR of $90$ dB (for about $16$ bits of an ADC,DAC with LSB noise) this means in principle you should be able to reach a data rate of $C = 500 \times 10^6 \times 15$ bits per second. This is about $7.5$ Gbps.


Most natural interpolation to simulate increasing distance?

One typically interpolates going up, 16 $\times$16 to 32 $\times$ 32.  Going down is low pass filtering (anti aliasing) and decimation.  As you have acknowledged, this isn’t the same as the effect of distance.  As long you understand that your original image was captured with a camera with a specific MTF, and sensor resolution, and are going to neglect scattering (which would require something like MODTRAN),  your original image has been acquired with some angular resolution, $\theta_0$, your image has a linear spatial resolution $\Delta_0$ such that $$
\Delta_0 = r_0 \theta_0
$$ where $r_0$ is the distance that the image was acquired at.  At some greater distance $r_g$, $\theta_0$ can be considered the same so,
$$ \Delta_g= \theta_0 r_g. $$
Given some coordinate system associated with the original system, you could use something like a Gaussian blur with a standard deviation related to $\Delta_g$, and resample. This doesn’t really help with a large depth of field.


change the central frequency of a wavelet

MATLAB documentation indicates that discretization of the continuous and infinite-support Meyer wavelet is done simply by regular sampling. I doubt that there are any unwanted consequences if you sample it over a longer interval than the nominal -8 to 8 effective support, to have a different central frequency.

The standard values for the first two parameters to meyer(LB,UB,N) are LB = -8, UB = 8, defining the interval to be covered by sampling. Multiply those parameters by the same factor by which you want to multiply the central frequency. If factor < 1, you should double LB and UB to ensure that you sample the whole effective support of the Meyer wavelet and scaling functions. Then you will find your wanted frequency on the next scale level.


Logarithmic vs Nth root transformation of Image

Mathematicaly, inverse logarithm (or exponential) should indeed grow faster than powers, and logarithm is slower than roots... at infinity. Here, their range is limited to finite interval, so with proper scaling, you may observe such behaviors. I will come back later with numerical examples.


Why use Wavelet Denosing instead of Band pass Filtering

If you use frequency selective filters for denoising you implicitly assume that the desired signal and the noise occupy different frequency bands. You only remove noise in frequency regions where there is no desired signal, whereas in the passbands of the filter the noise is left unchanged. For some applications the assumption that noise and signal occupy different frequency regions is justified, for some it isn't. Wavelet based denoising methods compose a noisy signal into different scales and attempt to remove the noise while preserving the signal, regardless of the frequency content.

In practice you have to look very carefully at the characteristics of the noise. It is always advisable to apply a frequency selective filter to suppress all out-of-band noise. For narrow band noise sources a notch filter is often sufficient. In cases where you can get a reference signal of the noise (i.e. an independent recording of the actual noise source) the best strategy is using an adaptive noise cancellation algorithm. Only if these more basic methods do not lead to a satisfactory result, and if signal and noise largely overlap in frequency and time, you should look into more sophisticated methods such as wavelet denoising, or blind source separation, which is discussed in the context of your application in this paper.


Bin sizes for non-uniform discrete Fourier transforms

The notion of a "frequency bin" is somewhat misleading. In general a sine wave of a single frequency will show up in ALL frequency "bins" of the Time discrete Fourier Transform. It will generate the highest value where $f_k$ is closest to the signal frequency and but the values at any other Fourier frequency are not zero, just smaller.

The only exception is if the signal frequency is an exact integer multiple of $f_k$ and not windowing is applied.


As you're probably starting to realise this is a very big field with a number of different methods which you can try.

If you want to try to use colour, you can try transforming from RGB to HSV space, which attempts to decouple the colour from the "lightness", this may improve your lighting inconsistencies.

The next problem is how to track your object. If you know the object's colour distribution (histogram) then you can use backprojection to assign a likelihood to each pixel being that colour, then perhaps you could simply track the centre of mass of the back projection.

Alternatively, you could look into template matching, where a template is scanned across the image and areas which match this template closely are scored more highly than non-matches, giving you likely areas for where the object is.

If you have any more specific information about what you are trying to acomplish, it will be easier to help you: do you have an example input image, with targeted highlighted?


What is an Intuitive Explanation of the Phase of a Signal

Compare $y= sin(\omega t)$ to $y=sin(\omega t-\phi)$ as shown in the plot below. $\phi$ is the additional phase term in radians where $\omega$ represents the frequency in radians per second. Thus the phase term shifts a sinusoid along the horizontal axis. So at a given frequency this will result in a time delay for that frequency, although to be noted that a fixed time delay will cause a phase shift that is proportional to frequency as a linear variation (The phase will increase linearly as the frequency increases). Thus our interest in "linear-phase" filters as this will have a constant "group delay" meaning the group of all frequencies in a waveform will be delayed the same amount in time, and thus not cause destructive interference between the different frequencies in the group. (Review Fourier Series and decomposition of waveforms into individual frequencies and the resulting reconstruction for more insight on that).

More pertinent (and useful) is complex signal representation with $Ae^{j \phi}$ representing the magnitude A and phase $\phi$ of a signal. This can be expanded using Euler's identity which may help give further insight: $$Ae^{j\phi}=Acos(\phi)+jAsin(\phi) = I + jQ$$

From which we see that in order to completely describe a complex signal we require two real signals, this could be magnitude and phase, A, $\phi$ as in $Ae^{j\phi}$ or the real and imagninary components I, Q as in I+jQ.

Refer to this post which may help give further insight into the importance and utility of complex representation for signals: Frequency shifting of a quadrature mixed signal


Is Smith's book still usable and up-to-date?

Do you think that there are a lot of new things have been introduced in the last 10 years? :-)

I like this book. Book is quite accurate and quite pragmatic.  I think this book is a great introduction to DSP. But feel free  to use another books, especially if you want to known about the modern hardware equipment.


How to create a control loop using a 100Msps digital signal as input

No Microprocessor is able to deal with such extreme rates as 100MS/s. That is clearly the domain of custom ICs and FPGAs, bringing the data to a very capable CPU, not something small and embedded.

Regarding USB:

I don't know your sampling bit width, but 40MB/s is not very much if you need to transport 100MS/s – you get 12bit per sample, in the ideal, zero-overhead case (which doesn't exist).

As a matter of fact, I know a series of devices very well that can deliver 100MS/s and more to a PC – these work through expensive interface such as a proprietary external PCIexpress interface, or better even through 10Gbit Ethernet (which you'd normally only find in data center / ISP equipment).

With a single 100MS/s stream you can easily flatten a single CPU core of a pretty capable x86_64 computer – really, that MUCH MUCH more data than what "microcontrollers" handle.

This might really be a case of insufficiently discussed system design – Nyquist and my understanding of the term "Ultrasound" say that you'd never need more than twice the highest signal frequency as sampling rate. I wouldn't expect ultrasound sonar to have frequencies higher than maybe 4 MHz, but I might be wrong.


Sampling of a continuous function: Kronecker's or Dirac's delta?

Modeling of the sampling process via multiplication of a continuous-time signal by a train of Dirac impulses is the most common interpretation in my experience. If you dig into it deeply enough, you'll find some disagreement about the mathematical precision of this approach*, but I wouldn't worry about it; it's just a convenient model for the process. There are no impulse generators inside your cell phone's ADC generating periodic bolts of lightning that multiply their analog inputs.

As you noted, you can't calculate the continuous-time Fourier transform of the Kronecker delta function, as its domain isn't continuous (it is limited to the integers). The Dirac delta function, in contrast, has a simple Fourier transform, and the effect of multiplying a signal by a train of Dirac impulses is easy to show due to its sifting property.

*: As an example, if you're going to be mathematically precise, you would say that the Dirac delta isn't a function at all, but a distribution instead. But at an engineering level, these issues are really just semantics.

Edit: I'll address the comment below. You gave your mental model of the sampling process as:

$$
f_s(t) = \sum_{k=1}^N \int_{t_k-\epsilon_k}^{t_k+\epsilon_k} f(t) \delta(t-t_k) dt.
$$

The problem with this interpretation is that the typical ideal sampling model does not have that integration built in. Instead, it's a pure multiplication of the input signal by a Dirac impulse train. If you look more closely at the equation you showed for $f_s(t)$, you'll see that the right hand side actually has no independent variable; $t$ is the dummy variable of integration. For any $\epsilon_k > 0$ in the above, according to the Dirac impulse's sifting property, you would get:

$$
f_s(t)=\sum_{k=1}^N f(t_k),
$$

which isn't correct. Instead, the model for the sampled signal is:

$$
f_s(t)=\sum_{k=-\infty}^{\infty} f(t) \delta(t-kT)
$$

Which is very similar to the above, except generalizing for an infinitely long impulse train along the time axis and assuming that the data is uniformly sampled at time instants $t_k = kT$. The Fourier transform of the resulting signal is:

$$
\begin{align}
F_s(\omega) &= \int_{-\infty}^{\infty}f_s(t) e^{-j \omega t} dt \\
&= \int_{-\infty}^{\infty} \sum_{k=-\infty}^{\infty} f(t) \delta(t-kT) e^{-j \omega t} dt \\
&= \sum_{k=-\infty}^{\infty} \int_{-\infty}^{\infty} f(t) \delta(t-kT) e^{-j \omega t} dt \\
&= \sum_{k=-\infty}^{\infty} f(kT) e^{-j \omega kT}
\end{align}
$$

If we define the discrete, sampled version of the signal $f(t)$ to be $x[n] = f(nT)$, then you are left with:

$$
F_s(\omega) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n}
$$

which is exactly the definition of the discrete-time Fourier transform.


Relation between poles and cut off frequency of a filter

The Butterworth filter is in fact a continuous-time filter, and what you are seeing is very probably a bilinear transformation applied to this filter in order to get a discrete-time filter.

Software such as Octave or Matlab do this when you ask for a digital Butterworth filter.

The original Butterworth filter, in continuous time, $H(s)$ has all poles equally distributed around a semicircle in the left of the $s$ plane. The radius of this circle is the cutoff frequency (in $rad/s$).

The relationship between frequencies in discrete-time ($\theta$) and continuous-time ($\omega$) given by the bilinear transformation is: $\omega = \frac{2}{T} \tan (\theta/2)$ (where $T$ is an arbitrary sampling period used during design).

You can imagine the left $s$ semiplane as mapped into the unit circle in $z$, or the imaginary axis mapped to the unit circumference.


Producing images with subpixel shift in region of interest

I tried to distill your question down to a 1-d problem. Does this "moving sine" do ~what you want to do?

t = 0:50;
f = (1/4)/5;
x = sin(2*pi*f*t);
nsteps = 50;
figure(1),
clf
ylim([-1 1])
xlim([0 nsteps])
stepvect = linspace(0,3,nsteps);
stepvect = [stepvect stepvect((end-1):-1:2)];
for step = 1:length(stepvect)
    y = interp1(t, x, t+stepvect(step)*sin(pi*t/50), 'linear');
    plot(t, y)
    grid on
    drawnow;
end


How much deviation is acceptable in FFT calculations (using computer algorithms)?

I wouldn't like to give a trivial answer as "it depends on the application" but indeed it does depend on the application!

Considering the fact that FFT computation will not be the only source of error in many applications it is reasonable to allow some acceptable deviations but generally being bounded by such as the required precision of the overall system.

On the other hand very strict precision bounds may be defined on such as IEEE floating point hardware (or soft simulation) units so that a prescribed set of computations will not generate an error larger than very small epsilon values.

I don't know if such a set of normalized limits exist specificaly for FFT or like algorithms. If you find out that there are mismatches from ideal DFT values and the corresponding fast algorithms you must make your own decision as to whether you should allow this error in favor of some other gains.


The AWGN variance is:

sig= 1/sqrt(2*SNR);

Then, you define a random sequence:

n1 = randn(NN,1);

So the noise signal is:

noise1=sig*n1;

And the PSD is:

N0=2*(sig.^2);

Have a look to this example for more info.


What do different Cepstral or MFCC coefficients represent intuitively?

When you take the log of some measurement value, you are simply changing the nature of the scale of what you are measuring.  Where the dividing line between positive and negative values on the new scale is depends only where the "1" was defined on the original scale.  There is no qualitative distinction between positive and negative per se, so what you are really asking is why are slowly changing components (lower coefficients) larger than other ones.
If a sinusoidal is enlarged proportionally, then the amplitude is inversely proportional to the frequency.  If your signal indicates this (regardless of context) it means similarity is being preserved across levels.
So this doesn't answer the title question, but hopefully helps with the enclosed ones.


How can I reconstruct a Time series using it AR coefficients in MATLAB?

If you drive your AR model with white noise, you can generate a time series that has the same power spectral density as the original time series.  But it won't be the same time series.


Complex numbers in $\tt ifft$ of an MMSE amplitude estimator

Note that the MMSE estimator gives you a real-valued gain function with which the noisy signal amplitude is modified. Equation (7) is actually

$$\hat{A}_k=G_kR_k\tag{1}$$

where the gain $G_k$ is a function of $v_k$ and $\gamma_k$. You want to reconstruct the original signal as

$$\hat{X}_k=\hat{A}_ke^{j\phi_k}\tag{2}$$

where $\phi_k$ is the phase of the noisy signal $Y_k=R_ke^{j\phi_k}$. Combining (1) and (2) gives

$$\hat{X}_k=G_kR_ke^{j\phi_k}=G_kY_k\tag{3}$$

So you can simply multiply the complex spectral component of the noisy signal with the gain $G_k$ to get the estimation of the spectral component of the input signal.

In Matlab this is easily done. You just have to extend the gain vector to satisfy the symmetry condition for real-valued signals:

X = Y.*[gain; gain(N-1:-1:2];

where I assume that the vector gain has elements $G_k$, and that N is its length. The relation between the FFT length and N is N = Nfft/2+1.


What happens, when an unlimited bandpass-signal gets mixed?

what happens, if the RF-signal from the antenna isn't band- or lowpass-filtered?

If the system is analog, lots of uncontrolled harmonics. If the system is digital, then yes, you will get aliasing.

Aliasing itself is not a form of mixing. But sampling is.

When an analog signal $x(t)$ is sampled at distinct times $\frac{1}{Fs}$ seconds apart, it is effectively mixed with a train of impulses at frequency $Fs$. Sample-and-hold turns these spikes into "steps". This is what is causing the "mirroring" of the spectrum.

In fact, by tuning your $Fs$ carefully, you can "shift" a signal into a baseband by undersampling it.

  I know for example from decimation (lowpass-filtering and downsampling) that there would be aliasing, which isn't a form of mixing but nonetheless.

Probably it is clear by now that it is. Take an analog signal $x(t)$, multiply it with a train of impulses $u_s(t)$ at some frequency $Fs$, now multiply the result by another train of impulses $u_d(t)$ that is equivalent to "taking out" every other pulse. If you don't first low pass $x(t)$ for the spectrum you are trying to "squeeze it in" post mixing by $u_s(t) \cdot u_d(t)$, you will get aliasing. It's the same concept applied twice.

Hope this helps.


How to calculate the gain in a bivariate fft in R?

I am/was shaky on the correct terminology. I have two signals, and I want to compare them to see how much the 2nd matches the first (coh), how much the 2nd lags behind the 1st (phase)

Normalized cross-correlations are good at that.  I don't know what the relevant command would be in Statistica, but in Matlab you would use the command xcorr.

  and the extent to which the 2nd is under/overshooting the peak amplitudes of the 1st.

To answer that in a useful way it would help a great deal to see the signals that you are talking about.  Are you just interested in the peaks, or the signal in general?  If you are talking about the signal in general then you could time-align the signals using the information from your cross-correlation, and then divide one of the signals by the other.  If you are just interested in the peaks you could do the same and then just extract the results at the peak locations.


RMS calculation in frequency domain after windowing

For a reasonably stationary time domain signal you can assume

rms(signal*windows) = rms(signal)*rms(window)

and can estimate the RMS of the un-windowed time signal this way. For a hanning window that means calculate the RMS using Perceval and then multiply with sqrt(8/3). It's not exact but will be fairly close if the signal is stationary with respect to the window length.
If it isn't stationary the interpretation of RMS isn't particularly well defined and it would helpful to learn more about your application and what exactly you are trying to do.


You have pairs of features-target's which is the first step.

I'm no expert in texture classification, there are probably specialized and/or well-known good methods.

To make it easy for a beginner, I would take a look at the examples for the Python library scikit-learn:
http://scikit-learn.org/stable/auto_examples/index.html#general-examples
and see what looks interesting, for example a few different random forest examples http://scikit-learn.org/stable/auto_examples/ensemble/plot_forest_iris.html

What you need to run these is a matrix with the filter responses features (the same shape as iris.data in the example (number of examples on the first axis and number of features on the second)) and a binary vector of targets (the same shape as iris.target (number of examples on the first axis))
basically you are replacing the "iris" dataset with your own dataset.

to load text data

import numpy as np
features = np.loadtxt('data.csv')

(as a side note: numpy works much like the builtin matrix stuff in matlab)

Try out a few different classifiers with a few different parameters and see what happens with the results http://scikit-learn.org/stable/modules/generated/sklearn.metrics.classification_report.html.

Make sure you split the data into train and test data http://scikit-learn.org/stable/modules/generated/sklearn.cross_validation.train_test_split.html or use cross-validation http://scikit-learn.org/stable/modules/cross_validation.html to avoid overfitting, again there is examples of use under the scikit-learn homepage.


Discrete Time Processing of Continuous Time Signals LTI Concept

Assuming a system like this:

x(t) -> C/D -> DSP -> D/C -> y(t)

where $x(t)$ and $y(t)$ are sufficiently bandlimited and the C/D and D/C blocks are ideal, then the system behaves as if it were LTI. Note that the system is not truly LTI, since it is not LTI for all inputs.

In particular, if the process of discretizing $x(t)$ produces aliases, then the output $y(t)$ can potentially contain frequencies not present in the input. Since an LTI system cannot create new frequencies, then in conclusion the system is not LTI.


Geometric explanation of a methodology in the article about Image Denoising

The explanation below the section you quote seems to be clear to me:

It's forming a vector field made up of the largest gradient direction, smallest gradient direction, and the normal to them both.


FIR filter design with frequency sampling method - setting proper phase response for IDFT

Your procedure is of course correct. It is in the details where things usually go wrong. One important thing is how to extend the desired frequency response beyond Nyquist taking the required symmetry into account. In order for the filter coefficients to be real-valued, the desired frequency response must satisfy

$$D[k]=D^*[N-k],\quad k=0,1,\ldots,N-1\tag{1}$$

where $*$ denotes complex conjugation, and $N$ is the FFT length (= filter length). From (1) it follows that for even $N$, the elements of the complex desired frequency response are

$$D_0,D_1,\ldots,D_{N/2},D^*_{N/2-1},\ldots,D^*_1$$

For odd $N$ you get

$$D_0,D_1,\ldots,D_{(N-1)/2},D^*_{(N-1)/2},\ldots,D^*_1$$
Note that in order for (1) to be satisfied, $D_0$ and, for even $N$, $D_{N/2}$ must be real-valued.

This little Matlab/Octave code works for even and odd filter lengths.

% frequency sampling design of linear phase FIR filter

N = 64;                         % FFT length = filter length
np = floor(N/2) + 1;            % number of independent frequency points
n = 0:np-1;
w = n*2*pi/N;                   % frequency vector
M = sin(n*pi/(np-1));           % some desired magnitude response
D = M.*exp(-1i*(N-1)/2*w);      % desired complex frequency response (linear phase)
D = [D,conj(D(N-np+1:-1:2))];   % append redundant points for IFFT
h = ifft(D);                    % compute impulse response
max(abs(imag(h)))               % should be very close to 0
h = real(h);                    % remove numerical inaccuracies

% check result
[H,w2] = freqz(h,1,4*N);
plot(w/2/pi,abs(D(1:np)),'.',w2/2/pi,abs(H))


Why the calculated value of power is not match?

Discrete signals are by definition periodic and you need to calculate the power in both domains over one period with proper normalization.

In practice this means, you need to divide the spectrum by the FFT length. See https://en.wikipedia.org/wiki/Parseval%27s_theorem


Combining different Likelihoods - Particle Filter

Let $\mathbf z = [z_1, z_2]^T$ be the observations and $\mathbf x = [x_1, x_2]^T$ be the hidden states, and the observation model is $z_i = x_i + n_i$ for $i=1,2$. Here $n_i$'s are independent (but not identical, since the noise has different distributions).

  So I would like to know under which restrictions the first step is allowed (I for example know that it its correct for a multivariate normal distribution)

By definition, the density function of a random vector is given by the joint density of its components. So $f(\mathbf z | \mathbf x) = f(z_1, z_2 | \mathbf x)$ is always true.

  and also how I could enhance the second approximation

The second approximation is actually an equality under the observation model with independent components. Since $z_2 = x_2 + n_2$, conditioned on $x_2$, we have $z_2 \perp (z_1, x_1)$. This implies that $f(z_2 | \mathbf x, z_1) = f(z_2 | x_1, x_2, z_1) = f(z_2|x_2)$. And similarly, $f(z_1|\mathbf x) = f(z_1 | x_1)$.


MVDR Beamformer, minimize but do not change pattern in main beam

Use a white noise gain constraint on the filter coefficients, which is a quadratic inequality constraint. When you have sources in the main beam that are not exactly on the unity gain constraint, the filter coefficients get large.  You can control the magnitude by diagonal loading the covariance matrix. There is a full derivation in "Detection, Estimation, and Modulation Theory, Volume 4, Array Processing" by Harry VanTrees.  It solves the problem but I'm not aware of a nice online algorithm. There is a paper by Cox and Zeskind that has an online algorithm but it isn't optimal no matter what Harry Cox thinks.

Another approach people use is to use a direction finding algorithm prior (before) to set the unity gain constraints (actually LCMV ) and then beam form. I prefer the white noise gain constraint.

Bob Zeskind suggested using an unmodified MVDR but forming about 3 times as many beams as you would in a conventional beamformer in a conversation I had with him. It's not elegant but does seem to have uses.


Comparing MFCC Features ,What do they represent?

Like already said, an eye inspection in this case is not possible since the data is just to high-dimensional. Thus what you want to do is to inspect the difference magnitudes of the frequencies per time frame.

I would suggest that you have a look at the cosine distance in this case. It's not the only way but a good possibility to compare MFCCs without distortion due to the different overall loudness of the word. The cosine distance is 1 - cosine similarity and is frequently used in text mining. In this area it might be important to be independent of the document length (document A with 100 words, 10 * "Hello" and 90 * "Bye" should be similar to document B with 200 words, 20 * "Hello" and 180 "Bye"). Another way of explaining this distance is perform a normalization step (project each data point on the unit sphere) and then measure the euclidian distance.

The cosine distance looks at the angle between two vectors (in your case you have one vector per time/frequency-bin) and therefore gets rid of the overall magnitude. This would by one way of comparing your MFCCs.


Generating a time domain signal using the fundamental frequencies

No way. You cannot recover a time domain signal from its fundamental frequency. Although, it is not clear what does a plural means here, "the fundamental frequencies and amplitudes", so maybe I am not quite understanding you.
Also, the power spectral density (PSD) alone, whether in dB or else, is not sufficient to recover a time domain signal: you need an energy spectral density for this deed. Citing a Wikipedia article:

The above definition of energy spectral density is suitable for
transients (pulse-like signals) whose energy is concentrated around
one time window; then the Fourier transforms of the signals generally
exist.

The Fourier transform of the signal is what you need in order to recover a time domain signal, because the PSD lacks the phase information. Moreover, a periodogram, if it is what you hint by mentioning "fundamental frequency", does embrace averaging, either in time or frequency, which makes the periodogram still less suitable for a time domain signal recovery.
Although it seems that your question is about PSD, you may be interested in retrieving phase information from energy spectral density. There is a number of techniques, the classical one is use of the Kramers Kronig relations for analytical signal. See also a review of three techniques (including KK) used in optical measurements, in Deep learning as phase retrieval tool for CARS spectra. Section 2.2 of this article presents the formula for deducing the phase from the squared modulus (formula 15), with stipulations made about a noise and discrete nature of measurements.
A word about implementation: in order to write a C++ implementation, you need a thorough understanding of signal processing algorithms. Searching github.com for open source implementations of phase retrieval techniques, I've found a Python implementation of phase retrieval for Raman spectroscopy, with Jupyter notebook, explaining the implementation details, see CRIKit2: Hyperspectral imaging toolkit. You may browse this presentation, at least out of curiosity.


The theory behind sweep-sine measurements of LTI systems requires a signal with constantly changing the frequency. You cannot simply playback few tones - the whole frequency range is necessary.

So that if you want to identify your system with the impulse response $h[n]$, you feed the sweep sine signal $s[n]$ into it and record the output. Obviously output will be given by the convolution:

$$y[n]=h[n]\star s[n]$$

In order to obtain the impulse response of the system, it is required to convolve the output $y[n]$ with the inverse filter $f[n]$. Generally, the inverse filter is a sweep signal, inverted in time and with scaled amplitude. So in the end:

$$h[n]=y[n]\star f[n]$$

There are two major types of sweep signals: linear and logarithmic (exponential) frequency change. The latter one is more commonly used, since its spectrum is similar to a pink noise. Personally I always use the logarithmic sweep sine:

$$s(t) = \sin\left[ \dfrac{2\pi f_1 T}{\ln\left( \dfrac{2\pi f_2}{2\pi f_1} \right)}\left(e^{\dfrac{t}{T}\ln\left(\dfrac{2 \pi f_2}{2\pi f_1} \right)} -1\right) \right]$$

where:

$f_1$,$f_2$ - initial and final frequency of a sweep

$T$ - sweep duration in seconds

There is no strict rule on what should be the length of a sweep. For example in acoustics, the rule of thumb is that it should be at least as long as predicted reverberation time.

For more literature on the topic, please refer to the following sources:

  Meng Q., et al. - Impulse Response Measurement With Sine Sweeps and
  Amplitude Modulation
  Schemes

  Farina A. - Advancements in impulse response measurements by sine
  sweeps

  Farina A. - Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique


Decimation with IIR filter altering amplitude of signal

First, 0.3dB is generally considered small potatoes in signal processing.

Second, the documentation for scipy.signal.decimate states that it uses an 8th-order Chebychev filter by default (without specifying the ripple, dangit!).  Chebychev filters have some ripple in the passband; if the bulk of your signal's energy is falling into the troughs of the ripple, then you should expect some attenuation.  In particular, it appears that your filter has a lot of DC content, and even-order Chebychev filters tend to have a gain of less than unity at DC.

So I suspect that it's doing exactly as expected.


How do I determine whether I have pattern-free noise?

If I were you, I would simply take a large number of no-signal, where you aren't measuring any real signal. Run it through whatever identification algorithms you have, and see if you see anything. If you do, then you need to worry about correlations.

I think what you might be missing is that correlation does not always mean a false detection, especially if you have a robust algorithm to this kind of noise.


How is channel frequency response/channel state information related to free space path loss?

In the OFDM system, such as Wi-Fi, there are infinite paths exist, assume it is in free space, and only one path exist,

Exactly. So if you use OFDM in a scenario where FSPL describes your transmission sufficiently, then you're doing something wrong, because you don't need OFDM, as you have a flat channel.
We can still describe large scale fading in multipath channels using models similar to your FSPL formula, but neither of your equations hold! Then, we're describing distance-dependent stochastic properties of the random variables that make up your CSI.
You need to cross over from your deterministic "this channel has a defined attenuation" model to a stochastic "and here's how we describe the properties of random variable, whose realizations are channels".


Convolving complex exponential with box function (discrete)

You can use the convolution theorem to convert the convolution to the product of Fourier transforms of the given signals.
$$h[n]*x[n]=\mathcal{F}^{-1}\{X(\omega)H(\omega)\}$$

$$H(\omega)=\mathcal{F}\{h[n]\}=\delta(\omega-\frac{\pi}{2})$$
$$X(\omega)=\mathcal{F}\{x[n]\}=e^{-j\frac{3 \omega}{2}}X_r(\omega)$$ where $$X_r(\omega)=\left\{\begin{matrix}
\frac{\sin(2\omega)}{\sin(\omega/2)} &\omega\neq 0 \\
4 & \omega=0
\end{matrix}\right.$$

However since $X(\omega)\delta(\omega-\frac{\pi}{2})=X(\frac{\pi}{2})\delta(\omega-\frac{\pi}{2})$ and $X(\frac{\pi}{2})=e^{-j\frac{3\pi}{4}}\frac{\sin(2\pi/2)}{\sin(\pi/4)}=0$, the inverse Fourier transform and therefore the covolution is zero.

You can also assume $x[n]$ is the impulse response of an LTI system and $h[n]$ is an input signal. Since the complex exponential is eigenfunction of the system, the output whould be an scaled version of the input i.e. $Ah[n]$ where $A=X(\frac{\pi}{2})$ and $X$ is the FT of $x[n]$.


Coefficients of the Fourier Series expansion of physical signals having finite power

A finite power periodic signal will have the following property:

$$ P_x = \frac{1}{T_0} \int_{0}^{T_0} |x(t)|^2 dt  < \infty \tag{1} $$

where $P_x$ is the power averaged over a period of the signal. Then from Parseval's theorem, the total power can also be shown to be (for a real signal I'm assuming):

$$P_x = |a_0|^2+2 \sum_{n=1}^{\infty} |a_n|^2 < \infty \tag{2}$$

where the $a_n$'s are the continuous-time Fourier series coefficients for $x(t)$.

From the theory of power series (calculus), it's known that in order for (2) to converge (sum being less than infinity) a necessary condition is

$$ \lim_{n \to \infty} |a_n|^2 = 0 \tag{3}$$

(3) implies therefore that, for finite power periodic signal's CTFS coefficients $a_n$ also goes to zero as $n$ goes to infinity.


How to deal with low fundamental when using AMDF for pitch extraction?

This is what we call in the pitch-detection biz, the "octave problem".
First of all, I would change the AMDF to ASDF.  And I would not reduce the window size as the lag increases.  (Also, I am changing notation to what I consider to be more conventional.  "$x[n]$" is a discrete-time signal.)
The Average Squared Difference Function (ASDF) of $x[n]$ in the neighborhood of sample $x[n_0]$ is:
$$ Q_x[k, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 $$
$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $.
Now, expand the square and consider what the summations look like as $N \to \infty$ (not that $N$ is going to infinity, but to give you an idea if $N$ is large).  The ASDF is directly related to the autocorrelation.  It is essentially the autocorrelation turned upside down.  These steps I will leave to you. take a look at this answer.
So now consider this finite-length "autocorrelation" (in the neighborhood of sample $x[n_0]$) defined from the ASDF:
$$ R_x[k,n_0] = R_x[0,n_0] - \tfrac12 Q_x[k, n_0]  $$
where
$$ R_x[0, n_0] \triangleq \frac{1}{N} \sum\limits_{n=0}^{N-1} \Big(x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor]\Big)^2 $$
This value $R_x[0,n_0]$ is a measure of the mean power of the signal $x[n]$ in the neighborhood of $n \approx n_0$. Since $Q_x[0,n_0]=0$ and $Q_x[k,n_0] \ge 0$ for all lags $k$, that means that $ R_x[k,n_0] \le R_x[0,n_0] $ for all lags $k$.
Another useful way to look at this autocorrelation taking place in the neighborhood centered at sample $x[n_0]$ is to normalize $R_x[k, n_0]$ with $R_x[0, n_0]$:
$$ r_x[k,n_0] \triangleq \frac{R_x[k,n_0]}{R_x[0,n_0]} = 1 - \frac{Q_x[k,n_0]}{2 R_x[0,n_0]}  $$
This normalized autocorrelation has $r_x[0,n_0]=1$ and $r_x[k,n_0] \le 1$ for all other $k$.
Suppose for a minute that $x[n]$ is periodic with period $P$ (and $P$ happens to be an integer), then
$$ x[n+P] = x[n] \quad \forall n $$
and  $Q_x[mP, n_0] = 0$ and $R_x[mP, n_0] = R_x[0, n_0] \ge R_x[k, n_0]$ for any integer number of periods ($m$ is an integer).  So you get a peak at $k=0$ and at $k$ equal to any other multiple of $P$ if $x[n]$ is periodic.  If $x[n]$ is not perfectly periodic, what we might expect is the biggest peak at $k=0$, another peak (but slightly smaller) at $k=P$ (the period we are looking for) and progressively smaller peaks for larger multiples of $P$.
We can then expect that the value of the normalized autocorrelation, $r_x[k,n_0]$ evaluated at a lag of $k=P$ or other multiples of $P$ should be pretty close to 1.  That value $r_x[P,n_0]$ can be thought of as a measure of the degree of periodicity (sometimes called the pitch confidence) of the estimated period $P$ for the quasiperiodic $x[n]$ in the neighborhood of $n \approx n_0$.  If $r_x[P,n_0]=1$, we can say that $x[n]$ is perfectly periodic with period $P$.  If the best $r_x[k,n_0]$ you can get (with $k$ that's not close to $k=0$) is very small, then $x[n]$ shows no periodicity and your pitch confidence is low.
So the octave problem comes about because of a couple of reasons.  First of all, $P$ is not necessarily an integer.  That is an interpolation problem, not a big deal.
The second reason and more difficult problem is that of subharmonics.  Consider that you're listening to a nice periodic tone at exactly A-440 Hz and it sounds like an A that is 9 semitones above middle C.  Now suppose someone adds to that tone a very tiny-amplitude (like down 60 dB) A-220?  What will it sound like and mathematically what is the "true" period?

Choosing the "right" peak for the period.
Let's say you run your note through a DC-blocking filter, so that the mean of $x[n]$ is zero.  It turns out that causes the mean of the autocorrelation $R_x[k, n_0]$ for every $n_0$ to also be zero (or close to it if $N$ is large).  That means $R_x[k, n_0]$ must sum (over $k$) to be about zero which means there is as much area above zero as below.
Okay, so $R_x[0, n_0]$ represents the power of $x[n]$ in the vicinity around $n=n_0$ and must be non-negative.  $R_x[k, n_0]$ never exceeds $R_x[0, n_0]$ but can get as large as it when $x[n]$ is periodic.  $R_x[P, n_0] = R_x[0, n_0]$ if $x[n+P]=x[n]$.  So if $x[n]$ is periodic with period $P$ and you have a bunch of peaks spaced apart by $P$ and you have an idea for how high those peaks should be.  And if the DC component of $R_x[k, n_0]$ is zero, that means in-between the peaks, it must have negative values.
If $x[n]$ was "quasi-periodic", one cycle of $x[n]$ will look a lot like an adjacent cycle, but not so much like a cycle of $x[n]$ farther down the signal in time.  That means the first peak $R_x[P, n_0]$ will be higher than the second at $R_x[2P, n_0]$ or the third $R_x[3P, n_0]$.  One could use the rule to always pick the highest peak and expect the highest peak to always be the first one.  But, because of inaudible subharmonics, sometimes that is not the case.  sometimes the second or possibly the third peak is oh-so-slightly higher.  Also, because the period $P$ is likely not an integer number of samples but $k$ in $R_x[k, n_0]$ is always an integer, so the true peak will likely be in-between integer values of $k$.  Even if you were to interpolate where the smooth peak is (which I recommend and quadratic interpolation is good enough), and how high it really is between integer $k$, your interpolation alg could make a peak slightly higher or slightly lower than it really is.  So choosing the absolutely highest peak can result in spuriously picking the second over the first peak (or vise versa) when you really wanted the other.
So somehow you have to handicap the peaks at increasing $k$ so that the first peak has a slight advantage over the second, and the second over the fourth (the next octave down), etc.  How do you do that?
You do that by multiplying $R_x[k, n_0]$ with a decreasing function of $k$ so that the peak at $k=2P$ is reduced by some factor, relative to an identical peak at $k=P$.  Turns out that the power function (not the exponential) does that.  so compute
$$ k^{-\alpha} \ R_x[k, n_0] $$
So, if $x[n]$ were perfectly periodic with period $P$, and ignoring interpolation issues for non-integer $P$, then
$$ R_x[2P, n_0] = R_x[P, n_0] $$
but
$$\begin{align}
 (2P)^{-\alpha} R_x[2P, n_0] & = \\
 (2P)^{-\alpha} R_x[P, n_0] & < P^{-\alpha} R_x[P, n_0] \\
 \end{align}$$
The factor by which the peak for a pitch of one octave lower is reduced is the ratio
$$ \frac{(2P)^{-\alpha} R_x[2P, n_0]}{P^{-\alpha} R_x[P, n_0]} = \frac{(2P)^{-\alpha}}{P^{-\alpha}} = 2^{-\alpha} $$
So if you want to give your first peak a 1% boost over the second peak, which means you will not choose the pitch to be the sub-harmonic pitch, unless the sub-harmonic pitch autocorrelation is at least 1% more than the first peak, you would solve for $\alpha$ from
$$ 2^{-\alpha} = 0.99 $$
That is the consistent way to weight or de-emphasize or handicap the peak corresponding to the subharmonic pitch one octave below.
It still leaves you with a thresholding issue.  You have to choose $\alpha$ well. But this is a consistent way emphasize the first peak over the second, which is an octave lower, but not so much that if the note really is an octave lower, but the energy in all of the even harmonics was strong, compared to the odd harmonics, this will still leave a possibility for the second peak being chosen.


I found the answer in Akos, D.M., Stockmaster, M., Tsui, J.B., & Caschera, J. (1999). Direct bandpass sampling of multiple distinct RF signals. IEEE Trans. Communications, 47, 983-988.

The extra condition between the 2 bands is: $|F_{IF1} - F_{IF2}| >= \frac{BW_{I1} + BW_{I2}}{2}$ where $F_{IF}$ is the intermediate frequency and $BW$ the bandwidth.

There are then 4 equations that describe all these constraints (eq $(1)$ to $(4)$) in the paper. Unfortunately, the first constraint contains a non-linear combination which cannot be solved analytically (by something like mixed-integer programming).

This is confirm in the paper: "An iterative procedure, based on these equations, was developed to determine sampling frequencies that would meet the requirement, the results of which are displayed in Fig. 5.". I conclude there is no better solution than iteratively going through the frequencies under consideration.


The amplitude of an image represents the intensity of the different frequencies in the image. Therefore, it holds the geometrical structure of features in the image (i.e. changes in the spatial domain).
The phase on the other hand, represents the locations of these features (which helps our human eye to better comprehend the image).
Here's a visualization I made in Matlab:

As you can see at the bottom, when we combine amplitude and phase from different images and perform inverse FFT, the resulting image looks much more like the original from which we took the phase. This means that phase holds more information about the image than amplitude.


With reference to $N_o$ this usually is the symbol for the power spectral density of thermal noise, where $N_o = kT$, where k is Boltzmann's Constant and T is the temperature in Kelvin. The thermal noise signal is a complex, white Gaussian distributed noise, with half of the power in the real component and half the power in the imaginary component (but power in itself is not a complex quantity). Therefore if the total power density is $N_o$, the power density of the real component would be $No/2$ and the power density of the imaginary component would also be $N_o/2$.

Considering the complex signal we can represent the spectral density as "Single-Sided" or "Double-Sided". With a Single-Sided power spectral density, the frequency axis starts at 0 and extends to infinity, making use that the the negative frequency components have the same density as the positive. We do not need to show this redundant information, but we do need to double the quantity that is shown to account for the power that extends into the negative frequencies. If we include the negative frequency axis, as done in a Double-Sided power spectral density plot, then the density is given as $N_o/2$.

In the end in either case the power density is $N_o$ Watts/Hz. In the Single Sided case the power in 1 Hz of bandwidth is completely shown on the positive axis, while in the Double Side case, half the power is in the positive axis and half is in the negative axis.

Thus in this we have made two considerations in how the power may be distributed; either in the positive and negative frequency axis if we present Single or Double-sided spectrums, and in the real and imaginary axis if we consider how the power in the full complex signal ($N_o$) is distributed between real and imaginary components. For example, if we were only to consider the real component, the power spectral density would be $N_o/2$ and then if we viewed that as a Double Sided PSD, the PSD would be $N_o/4$.

Single and Double Sided power spectral densities are frequency encountered when working with phase noise, which may provide a further intuitive way of seeing this: The double sided PSD for phase noise is given as $\scr{L}_\phi(f)$ and represents what we would read directly off of a spectrum analyzer, where we can observe both the upper and lower sidebands of the spectrum (although we typically only report one side on the plot which causes confusion in if it should be called single or double sided; one side is shown but it is a double sided power spectral density quantity--- we just don't need to show the other side since it is identical). The single sided power spectral density $S_\phi(f)$ represents the power from both sidebands as given on a plot with a positive frequency axis only (since the negative sideband is equivalent to the positive, we can put it all on the positive frequency axis simply by doubling in power what we read off of the spectrum analyzer).  So $S_\phi(f)$ represents all the power and is used in related calculations where we integrate over a band of frequency offsets, while $\scr{L}_\phi(f)$ is convenient in that we can read it directly off of test equipment (spectrum analyzers).


Difference between Gabor filtering and Discrete Wavelet Transform

Per se, a  Gabor filter  in image processing is one linear filter at a certain scale and 2D frequency used for orientation filtering and texture analysis.
It would be easier to compare Gabor representations and discrete wavelet transforms. Both are related to a linear decomposition of possibility multidimensional data at different scales with somehow oscillating functions. Main differences are:

Discrete wavelets: critical or non-redundant scheme (stable, invertible), discretize some families of wavelets, more scale-based  though wavelets range from weakly oscillating (Haar) to oscillation (Shannon), often applied separably across dimensions so weakly directional, a lot of statistical and approximation properties (moments, regularity). Their redundant counterparts: discretized continuous wavelets, shift-invariant or stationary wavelets.

Gabor transforms: (highly) redundant decomposition, mostly one shape: a modulated Gaussian, more frequency-based though computed at a couple of scale (Gaussian spread), inherently non-separable or directional. Their non-redundant counterparts: modulated or orthogonal lapped transforms, Malvar wavelets.


If you open the determinants along the middle column, we get for the first determinant:

$ -(bf-ec) + (af-cd)  -(ae-bd) = -3$ $\hspace{1cm}$ $ [1]$

And for the second determinant, we get

$-(bf-ec) + 2(af-cd)  -3(ae-bd) = -2 $ $\hspace{1cm}$ $ [2]$

The required determinant is

$  2(af-cd)  -4(ae-bd) $ $\hspace{1cm}$

Which is just $2*([2] - [1])$, which boils down to $2$.


Matlab results.

The filters coefficients have a big dynamic variation. So you have those results. For second filter, coefficients b are about zeros.

The first filter

b = [-1.40374673978609e+008 1.40374675978609e+008];
a = [1 1.000000235752411];
freqz(b, a)

The second filter

b = [2.22044577169459e-016 2.22044632680604e-016];
a = [1 -1];
freqz(b, a)


How can I lowpass filter by only reducing peak data?

There is indeed a 2D version for your attempt #2 - it is similar in theory, but it cannot be decomposed into two 1D operations. Please read about "2D grayscale morphological filtering". It is faster than curve fitting.

http://opencv.willowgarage.com/documentation/image_filtering.html
http://en.wikipedia.org/wiki/Dilation_(morphology)

Median filtering might also be useful if you are trying to remove speckles. A more advanced form of median filtering is "ordinal filtering".

In all cases, requirement #1 can be met trivially by taking the pixelwise minimum between the output and the input. It is an important quality criteria, but it will not limit the choice of algorithms.

Gaussian filtering (and a number of other useful filters) can be decomposed (first from 2D to 1D operations, then via Fourier transform), but there are many other useful image processing techniques that are not decomposable, which make them slow but does not diminish their usefulness.


Shape of Kaiser windowing function FFT response

The Kaiser window approximates the discrete prolate spheroidal sequence (DPSS) window which I think tends to a Gaussian function in an appropriate limit of large window and large $\alpha$, because the DPSS window concentrates energy both in time and frequency domain, as does the Gaussian function. The the zero-phase Kaiser window of length $N$ with shape parameter $\alpha$ is:

$$w_0(n) = \frac{I_0\left(\pi\alpha \sqrt{1-(\frac{2 n}{N-1})^2}\right)}{I_0(\pi\alpha)},$$

where $I_0$ is the zeroth-order modified Bessel function of the first kind. Substituting $n=\frac{(N-1)t}{2\sqrt{\alpha}}$ appears to normalize time so that the shape of the window does not change at large $\alpha$:

$$= \frac{I_0\left(\pi\sqrt{\alpha^2-\alpha t^2}\right)}{I_0(\pi\alpha)}$$

I was unable to calculate the limit $\alpha \rightarrow \infty$ symbolically but numerical testing indicates that at large $\alpha$ the window tends to a peak-normalized Gaussian function:

$$\approx e^{-\displaystyle\frac{\pi t^2}{2}}.$$

The maximum difference between the window and the Gaussian appears to be in the order of $\frac{1}{\alpha}$. The Fourier transform of a Gaussian function is a Gaussian function.


What is the principle behind my voice echoing perpetually when 2 separate devices are used for playback and recording on video calls?

On each individual device, the speaker output can get subtracted from the microphone before it gets sent to other locations. This prevents others from hearing themselves through your microphone. When using two devices in within audible range of each other, the devices cannot subtract the speaker audio from the microphone audio because the information path does not exist. The delay time is the determined by the network delays between the two devices and if the device with the speaker audio is very loud, you can get a gain of more than 1 to the device with the microphone. If the device with the microphone is using automatic gain control and you are not the loudest thing in the room, then you will eventually only hear the loudest thing in the room. In your case, the echo must have been very loud.


Distinguish between sinusoidal oscillations and rythmic pulses

You are going to want to take your DFT (Discrete Fourier Transform) on a whole number of repeat cycles.

Here is a good way to identify boundaries for your repeat pattern:  Use the difference method I describe in my blog article "Exponential Smoothing with a Wrinkle".  Smooth your signal heavily and locate your boundaries by finding positive to negative zero crossings.  Select a duration of at least three cycles, but keep the number limited so you know your cycles have a similar shape within your duration.

Take a DFT of your original signal on your duration.  A FFT is a fast version of this.  For real valued signals, only the first half of the return values are of interest.  The second half are a mirrored complex conjugate of the first and don't add any new information.

For a pure tone, only one bin will have a value in it.  The rest will be zero or near zero.  For a more complicated waveform, the bin values at multiples of the number of cycles you included will give the data to calculate Fourier coefficients for your waveform.

By a "not pure sinusoidal oscillations", I assume you mean a pure tone with a little harmonic distortion.  In this case, the magnitudes of the bins of the harmonic tones will be smaller compared to a sharp "rhythmic pulse" like the ones in your graph.

If you are using MATLAB, be aware that it is one based, and DFTs are zero based so you will have to subtract one from the MATLAB index to get the DFT bin number.

Hope this helps.

Ced

Followup

If your rhythmic signals are like the one you showed, you can simply count zero crossings of your signal within the bounded intervals defined above.


Finding the maximum phase and the maximum frequency deviation of an angle-modulated signal

I think this is treating the $10^8\pi\cdot t$ as the carrier frequency and the $5 sin(2\pi\cdot 10^3 t)$ as the modulation.  So the maximum phase difference of this signal to the carrier without modulation would be $\pm5$ radians.
To compute the change in frequency you need to find the compute $\frac{d}{dt}$ of the argument of the cos (divided by $2\pi$) and subtract the carrier frequency.
$$\frac{d}{dt} 10^8\pi\cdot t + 5 sin(2\pi\cdot 10^3 t) = 10^8\pi + 2\cdot 5\cdot 10^3\pi cos(2\pi\cdot 10^3 t)  $$
Divide by $2\pi$ to get Hz instead of radians.
$$Frequency = 5\cdot 10^7  + 5\cdot 10^3 cos(2\pi \cdot 10^3t)$$
Subtract off the carrier of 50MHz and you have a maximum frequency excursion of $\pm 5$ k Hz


Transformation of random variables vs shift of functions

I think you've actually answered your own question by observing that the pdf of the random variable (RV) $Y=X-a$ is given by $f_Y(y)=f_X(y+a)$, which corresponds to a left shift for $a>0$. This is also immediately clear from the fact that by subtracting $a$ from $X$, the mean of the transformed RV is also shifted to $m_Y=m_X-a$.
Shifting a general function $f(t)$ is just the same, and you actually did exactly that with the pdf of the original RV to obtain the pdf of the transformed RV. For a shifted function $f(t-t_0)$, the origin (i.e., the value $f(0)$) is shifted to the value $t=t_0$. So for $t_0>0$ this corresponds to a shift to the right, and for $t_0<0$ this corresponds to a shift to the left.


MATLAB phase of 2D rectangular pulse's Fourier transform

The issue is that the phase value found is between $-\pi$ and $\pi$ (or $0$ and $2\pi$) but it needs to be "unwrapped" to be continuous.

In 1D, the unwrap function will help. Your mileage may vary when applying it to a 2D Fourier transform.

Also, there is a bug in your code that forces the R variable to be real valued:

R=re+imm;
phase=angle(R);

is not going to get you the right phase at all.


Looks like your data is virtually free of noise. That, combined with a very high sampling frequency would mean that at the jumps the data is exactly at the threshold between two quantized values. Set up nodes at the middle points of the vertical jumps and construct splines that connect the nodes. The easiest is to just draw straight lines between successive nodes, which gives a piece-wise constant differential. I wonder if that is good enough, or if you already tried that. If you need the velocity real-time, this approach is problematic because occasionally you might have to wait for a new node for some time.

You can further low-pass filter the interpolated data. If you use a filter with an impulse response that is nowhere negative, such as a Gaussian function, then there will be no overshoot.

With linear interpolation, everywhere between successive nodes, the speed will be simply the position difference between the nodes divided by the time difference between the nodes. You can run the smoothing filter on that piece-wise constant speed data and the result will be the same as if you'd run it on the linear position ramps and then differentiated (associative property of convolution, as both differentiation and filtering are convolution).


What you are asking about is typically called "pitch detection". There are 100s of articles written on the topic. A good starter would be this one: https://ccrma.stanford.edu/~pdelac/154/m154paper.htm
If you need fairly high precision (for example a tuner application), I recommend a phase locked loop (PLL) or delay locked loop. There is a fundamental trade off between precision and tracking speed and with a PLL you can adjust this easily.


How can we know the geometry of a blob detected by OpenCV

Instead of connected components you may also use blob detector mentioned here

It provides some information on the shape and geometry of blob

For geometric moments you will have to compute contours.

I am not sure what do you mean by connected components being a black box.


How to find the fundamental frequency of a discrete signal using partial autocorrelation?

Very simple: your sampling rate is 512Hz and your sample time is the inverse of this: about 1.953 ms. Multiply with 506 and you get the time length of the fundamental: 0.988s.

So you have one beat every 0.988s which is 1.0119Hz or 60.711 beats per minute.


Z-domain transfer function to difference equation

$$\begin{align*}\dfrac{Y(z)}{X(z)} &= \dfrac{1+z^{-1}}{2(1-z^{-1})}\\
\\
2(1-z^{-1})Y(z)&=(1+z^{-1})X(z)\\
\\
Y(z) -Y(z)z^{-1}&= \frac{1}{2}X(z) +\frac{1}{2}X(z)z^{-1}\\
\\
y[n]-y[n-1]&=\frac{1}{2}x[n] + \frac{1}{2}x[n-1]\\
\\
y[n]&=y[n-1]+\frac{1}{2}x[n] + \frac{1}{2}x[n-1]\\
\end{align*}$$

You'll need at least one register for the feedback of the output (the y[n-1] term).

You might as well call that register the accumulator, as this is essentially an integrator, integrating up the average of the two most recent input values.  Integrators are not stable, so you may want to clip the output at some point instead of letting it roll over.


Why do we integrate the modulating function in FM?

You need to integrate the modulating signal because frequency is the time derivative of phase. Therefore, the typical relationship from introductory calculus holds:

$$
\phi_i(t) = \int_{-\infty}^{t} \frac{d\phi_i(\tau)}{d\tau} d\tau = \int_{-\infty}^{t} 2\pi f_i(\tau) d\tau
$$

For causal signals, the lower limit on the integral changes to zero:

$$
\phi_i(t) = \int_{0}^{t} 2\pi f_i(\tau) d\tau
$$

$$
\phi_i(t) = \int_{0}^{t} 2\pi\left(f_c + k_fm(\tau)\right) d\tau
$$

$$
\phi_i(t) = 2\pi f_ct +\int_{0}^{t} 2\pi k_f m(\tau) d\tau
$$

which is the answer you were trying to get to.


Modal analysis and frequency response function - interpretation needed

This may be a better question for https://physics.stackexchange.com/ but I'll give it a shot.

I assume you show the FFT of the time waveform at the sensor. In order to get the actual frequency response you would need to divide by the FFT of the hammer response. Since the hammer is reasonably smooth, we can eyeball this. It also looks like your measurement is "good": sufficient signal to noise ratio and nothing unexpected in there.
Force is a vector: it's unclear whether you are only measuring the normal force or a full 3 component vector and calculate the magnitude.
At a broad brush, the frequency response can be interpreted as "reasonably flat" with a bunch of local wiggles at specific frequencies
The flat part of the frequency response (FR) simply means the piece itself is moving (or compressing) uniformly: i.e. it's moving as one without the sensor surface deforming.
The wiggles are modal behavior. Here the piece (or surface) has different movement at different locations.
The shape of the modes is determined by the geometry (and material) of the piece. That's why the modes tend to be very narrow band.
Not all modes will be visible in a measurement like this. Modes have "nodes" and "peaks" as a function of position. In order to see mode both excitation and sensor need to be far enough away from a node. If the hammer is close to the node, the mode will not be excited and if the sensor is close to the node, the sensor won't pick it up (despite it being there and active).
Modal movement implies different areas of "positive" and "negative" movements.
The modal movement simply adds to the "move as one piece" movement. If both go in the same direction you see a peak and if both go in opposite directions you see a dip.

  Does it mean that forces measured at these frequencies are damped?

Not really: what it means is that movement due to modal behavior and movement of the entire piece as a single chunk are in opposite directions at this particular point and frequency. So the sum is smaller than each individual movement.


Subgradient Method for K-Means Like Problem

You should use something like Robust Convex Clustering:

Assume that we are given $ n $ data points and each data point
is described by a $ d $ dimensional feature vector. We collectively
represent the data by a data matrix $ X \in \mathbb{R}^{d \times n} $.
The convex clustering method clusters data points into groups via solving
a convex optimization problem:

$$ \arg \min_{P} \frac{1}{2} \left\| X - P \right\|_{F}^{2} + \alpha \sum_{i < j} {w}_{i, j} \left\| {P}_{i} - {P}_{j} \right\|_{p} $$

Where $ \left\| \cdot \right\|_{p} $ is the $ {L}_{p} $ norm, $ P \in \mathbb{R}^{d \times n} $ is the matrix consisting cluster centers and assignments, $ {P}_{i} $ is the $ i $ -th column of the matrix $ P $ and represents the centroid of cluster that the $ i $ -th data point is assigned, $ \alpha $ is a non-negative regularization parameter that controls the number of clusters. $ {w}_{i, j} $ is the weight between $ i $ -th and $ j $ -th data point specified by the user. Clusters are implicitly given by $ P $ matrix and $ {P}_{i} = {P}_{j} $ means the two data points are in the same cluster. When $ \alpha $ is zero, each data point forms a cluster since $ P $ equals to $ X $ when Eq. (1) reaches the
optimum. As $ \alpha $ increases, $ \left\| {P}_{i} − {P}_{j} \right\|_{p} $ becomes smaller for all pairs so that the whole function can reach the minimal value. Hence more data points will be grouped into a cluster.

Another great reference is Splitting Methods for Convex Clustering.

As the formulation you chose is not convex and of top of that I don't think anyone has ever calculates its Sub Gradient.

Update

I encountered another simple yet effective of K-Means - Revisiting K-Means: New Algorithms via Bayesian Nonparametrics.
It formulates the problem in as, intuitively, the Lagrange form of the K-Means with continuous parameter $ \lambda $ instead of the discrete parameter $ k $.


What factors go into choosing DCT type 1 over type 2 or type 3 etc?

As you described correctly, DCT is closely related to DFT. The two main motivations behind DCT are it's being a real transform and its compaction of signal energy. The first makes it computationally simpler, while the second makes it a prime tool for transform based lossy data compression applications. As you understood correctly.

Then what remains is where the types come from and why is it the DCT-II as the most common choice for data compression? The answer is as follows:

Given a signal $x[n]$ of length $N$, The different DCTs are defined by forming even-symmetric extensions of shifted copies $x[n]$ in a number of possible ways, each different combination yielding a different type of DCT. Also the odd-symmetric extensions are called as the DST (Discrete Sine Transform).

The DCT-II , by virtue of its type of symmetry, has the most energy compaction among other DCT types and therefore is naturally the choice.

The DCT-II signal formation (based on A.Opp DTSP-2e) is as follows:

Consider a signal $x[n]$ of length N, whose DCT-II we want to compute. First extend $x[n]$ by padding $N$ zeros and create the extended signal $x_e[n]$

Then the DCT Type-II signal $\tilde{x}_2[n]$ is a periodic signal with period $2N$ and related to $x_e[n]$ as:
$$\tilde{x}_2[n] = x_e[((n))_{2N}] + x_e[((-n-1))_{2N}]$$
whose base period $x_2[n]$ extends from $n=0$ to $n=2N-1$. And where the $((\cdot))_{N}$ operation here is the modulus operator, providing the periodicity of $\tilde{x}_2[n]$.

Eventually, it is shown that the $N$-Point  unnormalised  DCT-II, $X_{c2}[k]$, of the signal $x[n]$ is related to the $2N$-Point DFT, $X_2[k]$, of the signal $x_2[n]$ as:
$$X_{c2}[k] = e^{-j\pi k/{2N}} X_2[k]$$
for $k=0$ to $k=N-1$

The following Matlab script creates the Type-II DCT signal $x_2[n]$ from $x[n]$ and finds DCT of x[n] from DFT of x2[n] for comparison

N=8;                            % Period for x[n]
n = [0:2*N-1];                  % A set of sequence indices of at least length 2N
x = [1:N];                      % A ramping test sequence for x[n] of period N
xe = [x zeros(1,N)];            % Extended sequence from x[n], by padding N zeros
x2 = xe(1 + mod(n,2*N)) + xe(1 + mod(-n-1,2*N));    % DCT-II type periodic signal
DFTx2 = fft(x2,2*N);            % 2N-Point DFT of x2[n]
Dft2Dct = real( exp(-j*pi*[0:N-1]./(2*N)).*DFTx2(1:N) );   % Create N-DCT of x[n] from 2N-DFT of x2[n]
DCTx = sqrt(2*N).*dct(x,N);     % N-POint DCT of x[n], for test purposes
DCTx(1)= sqrt(2)*DCTx(1);       % UNnormalise the Matlab's default normalised DCT-II


Gaussian Pyramid - How is Subsampling Rate Related to Sigma?

I've spent a bit of time searching for an answer to this as well.

I suspect since the Gaussian kernel has infinite support there will always be a little bit of aliasing, as there will always be frequencies passed that are greater than 0.25 of the sampling rate. That is, in the upper half of the spectrum.

So one would choose $\sigma$ to keep the amount of energy in the aliased bands to an acceptable level. Using $\sigma = 0.25 f_s$ would give 68%, $2\sigma = 0.25 f_s$ would give 95%, and so on.

From wikipedia:


What is the effect of a sawtooth phase mask in the frequency domain?

let

$$ \phi_0 \triangleq \phi - 2 \pi \left \lfloor \frac{\phi}{2\pi} \right \rfloor $$

where $\lfloor \cdot \rfloor$ is the floor() function, the largest integer not exceeding the argument.

it shouldn't be too hard to show that

$$ e^{j \phi_0} = e^{j \phi} $$

for all real $\phi$.


Laplace Transform of Cosine, Poles and Mapping to Frequency Domain

You're comparing the transforms of two different functions. You consider the Fourier transform of the function $x_1(t)=\cos(\omega_0 t)$, but you took the Laplace transform of the function $x_2(t)=\cos(\omega_0t)u(t)$, where $u(t)$ is the unit step function:

$$X_1(j\omega)=\int_{-\infty}^{\infty}x_1(t)e^{-j\omega t}dt\\
X_2(s)=\int_{0}^{\infty}x_2(t)e^{-st}dt=\int_{0}^{\infty}x_1(t)e^{-st}dt$$

Note the difference in the lower integration limits.

The (bilateral) Laplace transform of the function $x_1(t)=\cos(\omega_0t)$, $-\infty<t<\infty$, does not exist, whereas the Fourier transform of the function $x_2(t)=\cos(\omega_0t)u(t)$ does exist:

$$\mathcal{F}\{\cos(\omega_0t)u(t)\}=\frac{\pi}{2}[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]+\frac{j\omega}{\omega_0^2-\omega^2}\tag{1}$$

However, it cannot be obtained by simply replacing $s$ by $j\omega$ in the expression for $X_2(s)$ because $X_2(s)$ has poles on the imaginary axis, and replacing $s$ by $j\omega$ only gives the correct expression for the Fourier transform if the imaginary axis is inside the region of convergence.

Note that in general for causal functions there are three cases concerning the relationship between the Laplace transform and the Fourier transform:

If the region of convergence (ROC) contains the $j\omega$-axis (i.e., all poles are in the left half-plane), the Fourier transform is simply obtained by the substitution $s=j\omega$.
If there are poles on the $j\omega$-axis, but no poles in the right half-plane, then the Fourier transform contains Dirac delta impulses (as in the example above) plus a term that can be obtained from the Laplace transform by setting $s=j\omega$. Note that in $(1)$ the right-most term is simply $X_2(j\omega)$.
If there are poles in the right half-plane, the Fourier transform of the corresponding causal signal does not exist.


Golay complementary sequences in 802.11ad

For part (a), I discovered that the IEEE 802.11 protocol doesn't recommend adding auto-correlations of $Gu_{512}$ and $Gv_{512}$. Instead the correct method is to autocorrelate the first 256 bits of each of these sequences with the next 256 bits. So, $Gu_{512}$ and $Gv_{512}$ are not complementary Golay sequences themselves. But, essentially, both $Gu_{512}$ and $Gv_{512}$ consist of two 256-bit Golay complementary pairs.

~ryan


How is that after random numbers with their negatives alternating passing through 5, 4, 3 order moving average filter, become quasi sinusoidals?

There isn't a mystery here. Your input signal isn't random; you chose 5 numbers at random, perhaps, but then you repeated the sequence of 5 numbers several times. This is a periodic signal.

In continuous time, a periodic signal has a Fourier series representation; its spectrum contains discrete lines (i.e. sinusoidal tones).

In a finite-length discrete-time signal, such as the one you proposed, you see a similar phenomenon. The spectrum of a signal that is periodic within its duration will have a line structure to it.

Your moving averages are simple lowpass filters. As you apply more and more averages, you're filtering out the higher harmonics, leaving a combination of relatively fewer sinusoids. Thus, output of the cascade of moving average filters has a readily discernible structure as a sum of a relatively small number of sinusoids, as you described.


Can I use more than one Window on an FFT?

That makes little sense; the FFT (which is just an implementation of the DFT) is a linear operation; summing/weightedly averaging the same signal windowed with different windows is mathematically identical to just windowing with the summed/weightedly averaged window to begin with. It's literally the same.
So, what you really want is to design / choose a window that's appropriate for your use case instead!


Phase difference calculation for non divisible sample rate to frequency

Re-reference your FFT phase measurement reference point to the center of your data window by doing an FFTShift (N/2 vector rotate) before the FFT.  That moves the portion of data where there is a circular discontinuity (due to non-integer periodicity) away from the phase reference point of a typical FFT.


A dB scale is simply a logarithmic scale.  Bels are base 10, and decibels are base $ 10^{1/10} $.  Because $log(a/b) = log(a) - log(b)$, a decibel scale converts ratios to intervals.  3 dB corresponds roughly to a factor of 2.  So, what you are really asking when you say looking for a dB value that is 3 dB down is equivalent to asking for where the raw value is half of the original raw value.

Added:  The units of the raw values don't matter since you are taking a ratio.  Simply take $ dB = 10 \cdot \log( raw ) $ to convert to a dB scale.

One thing you have to be careful of when taking the log of PSD values is that the PSD may have a zero value which makes the log undefined, but can be thought of as negative infinity.

Hope this helps,

Ced


Inverse Fourier Transform of a squared sinc function

For $T_o>0$, this integral can be tackled with a Schwinger parameterization and Fubini's Theorem.
(The equivalent of the intermediate integral you specifically asked about is at the 7th equal sign below, with the trig functions written in terms of complex exponentials.)
$$\begin{align*}h(t) &= 2AT_o\int_{-\infty}^{\infty} \frac{\sin^2(2\pi fT_o)}{(2\pi fT_o)^2}e^{2\pi ift} \; df\\
\\
&= 2AT_o\int_{-\infty}^{\infty} \frac{\sin^2(2\pi fT_o)}{(2\pi fT_o)^2}\cos(2\pi ft) \; df\\
\\
&= 4AT_o\int_{0
}^{\infty} \frac{\sin^2(2\pi fT_o)}{(2\pi fT_o)^2}\cos(2\pi ft) \; df\\
\\
&= 4AT_o\int_{0}^{\infty} \sin^2(2\pi fT_o)\cos(2\pi  ft)\int_0^{\infty}ue^{-u2\pi fT_o} \; du\; df \\
\\
&= 4AT_o\int_{0}^{\infty} \int_{0}^{\infty} u\sin^2(2\pi fT_o)\cos(2\pi  ft)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{i2\pi fT_o}-e^{-i2\pi fT_o}\right)^2\left(e^{i2\pi  ft}+e^{-i2\pi  ft}\right)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{i2\pi f 2T_o}-2+e^{-i2\pi f 2T_o}\right)\left(e^{i2\pi  ft}+e^{-i2\pi  ft}\right)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{i2\pi f (2T_o+t)}+e^{i2\pi f (2T_o-t)}-2\left[e^{i2\pi  ft}+e^{-i2\pi  ft}\right]+e^{-i2\pi f (2T_o-t)}+e^{-i2\pi f (2T_o+t)}\right)e^{-u2\pi fT_o} \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty} \int_{0}^{\infty} u\left(e^{-f(2\pi uT_o - i2\pi  [2T_o+t])}+e^{-f(2\pi uT_o - i2\pi  [2T_o-t])}-2\left[e^{-f(2\pi uT_o -i2\pi  t)}+e^{-f(2\pi u T_o+i2\pi  t)}\right]+e^{-f(2\pi uT_o + i2\pi  [2T_o+t])}+e^{-f(2\pi uT_o + i2\pi  [2T_o-t])}\right) \; df\; du \\
\\
&= -\dfrac{AT_o}{2}\int_{0}^{\infty}  u\left(\dfrac{1}{2\pi uT_o - i2\pi  [2T_o+t]}+\dfrac{1}{2\pi uT_o - i2\pi  [2T_o-t]}-\dfrac{2}{2\pi uT_o -i2\pi  t}-\dfrac{2}{2\pi u T_o+i2\pi  t}+\dfrac{1}{2\pi uT_o + i2\pi  [2T_o+t]}+\dfrac{1}{2\pi uT_o + i2\pi  [2T_o-t]}\right) \; du \\
\\
&= -\dfrac{A}{4\pi}\int_{0}^{\infty}  \dfrac{u}{u - i\left[2+\frac{t}{T_o}\right]}+\dfrac{u}{u - i\left[2-\frac{t}{T_o}\right]}-\dfrac{2u}{u -i\frac{t}{T_o}}-\dfrac{2u}{u+i\frac{t}{T_o}}+\dfrac{u}{u+ i\left[2+\frac{t}{T_o}\right]}+\dfrac{u}{u + i\left[2-\frac{t}{T_o}\right]} \; du \\
&= -\dfrac{A}{2\pi}\int_{0}^{\infty}  \dfrac{u^2}{u^2 +\left[2+\frac{t}{T_o}\right]^2}+\dfrac{u^2}{u^2 + \left[2-\frac{t}{T_o}\right]^2}-\dfrac{2u^2}{u^2 +\left[\frac{t}{T_o}\right]^2}\; du \\
\\
&= -\dfrac{A}{2\pi}\int_{0}^{\infty}  \dfrac{u^2}{u^2 +\left[\frac{t}{T_o}+2\right]^2}+\dfrac{u^2}{u^2 + \left[\frac{t}{T_o}-2\right]^2}-\dfrac{2u^2}{u^2 +\left[\frac{t}{T_o}\right]^2}\; du \\
\\
&= -\dfrac{A}{2\pi} \left(u - \left[\frac{t}{T_o}+2\right] \tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}+2}\right]+u - \left[\frac{t}{T_o}-2\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}-2}\right]-2u +2\left[\frac{t}{T_o}\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}}\right]\right)\bigg{|}_0^\infty\\
\\
&= \dfrac{A}{2\pi} \left(\left[\frac{t}{T_o}+2\right] \tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}+2}\right]+\left[\frac{t}{T_o}-2\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}-2}\right]- 2\left[\frac{t}{T_o}\right]\tan^{-1}\left[\dfrac{u}{\frac{t}{T_o}}\right]\right)\bigg{|}_0^\infty\\
\\
&= \begin{cases} \dfrac{A}{2\pi} \left(-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]+2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right) & \dfrac{t}{T_o} \le -2 \\
\\
\dfrac{A}{2\pi} \left(\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]+2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right) &-2 \le \dfrac{t}{T_o}\le 0  \\
\\
\dfrac{A}{2\pi} \left(\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]-\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]-2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right)&0 \le\dfrac{t}{T_o}\le 2 \\
\\
\dfrac{A}{2\pi} \left(\dfrac{\pi}{2}\left[\dfrac{t}{T_o}+2\right]+\dfrac{\pi}{2}\left[\dfrac{t}{T_o}-2\right]-2\dfrac{\pi}{2}\left[\dfrac{t}{T_o}\right]\right)&\dfrac{t}{T_o}\ge 2\\
\end{cases} \\
\\
h(t) &= \begin{cases} 0 & \dfrac{t}{T_o} \le -2 \\
\\
A\left(1+\dfrac{1}{2}\cdot\dfrac{t}{T_o}\right) &-2 \le \dfrac{t}{T_o}\le 0  \\
\\
A\left(1-\dfrac{1}{2}\cdot\dfrac{t}{T_o}\right)&0 \le\dfrac{t}{T_o}\le 2 \\
\\
0 &\dfrac{t}{T_o}\ge 2\\
\end{cases}\\
\\
\end{align*}$$


With a real cosine a phase shift is equivalent to a time delay.  If you delay the signal by $\phi_0$ then you will get the result you want.

The number of samples that you will need to delay by will likely not be an integer number.  You can delay by a non-integer number of samples via a fractional delay filter.


You're dealing with digital signals. Power hence doesn't have a physical unit, just as the amplitude of the signal doesn't have one – where it originated, i.e. the ADC, it might have had one, but usually, that unit just doesn't matter in DSP and the fact that the digital numbers are proportional with some constant factor to a physical signal is totally sufficient.

If you later need to actually map a digital power to a physical power, you'd just calibrate (ie. measure with a known source of amplitude or power) and calculate that constant factor.

TL;DR: This is DSP. We don't care about physical units, so (mentioning the) resistance is futile.


How do you alter the phase relationship between I and Q in a quadrature signal?

As you already know, the quadrature error, also called quadrature skew, describes how far off the actual angle between I and Q is from the ideal 90 degrees.
It is one component of modulation error ratio.

There seem to be lots of people talking about how to measure quadrature skew,
but few people talking about how to compensate.
The two major approaches to compensation seem to be:

During initial factory tuning, some human watches an o'scope measuring the skew while manually tweaking the analog parts, or
During normal operation, the receiver internally estimates the quadrature skew and compensates digitally.

Often hardware designers deliberately build "off-tuned recievers".
With such hardware, when the transmitter repeatedly transmits "the same" constellation point, the raw I/Q point sampled by the reciever -- rather than staying close to a single fixed location in the I/Q diagram as with a perfectly tuned direct-conversion receiver -- instead rotates in a circle.

Such a receiver somehow estimates (typically in a FPGA or in software) how far the constellation has been rotated by that analog hardware, and corrects for it by rotating the I/Q sample in the opposite direction.
It is possible to simultaneously estimate the phase skew between I and Q and correct for that at the same time.
(In effect,
if you know the skew is S and you measure raw x,q coordinate,
you can deskew and find the true x,y in perfectly orthogonal coordinates
with y = qc - xt,
where c is some approximation of 1/cos(S) and t is some approximation of tan(S); often c = 1 and t=S (in radians) is adequate).

Further reading:

"Blind Equalization of Receiver In-Phase/Quadrature Skew in the Presence of Nyquist Filtering"
a patent search for "quadrature skew" compensation turns up a few.
"Effects of physical layer impairments on OFDM systems"
"Compensation for In-Phase/Quadrature Imbalance in Coherent-Receiver Front End for Optical Quadrature Amplitude Modulation"
"Correcting Imperfections in IQ Modulators to Improve RF Signal Fidelity"
"Estimation and Compensation of Impairments caused by Quadrature Modulators and Demodulators"


Fourier Transform Computer Vision Textbook

Thats funny,  I was exactly trying to fill the same gap.
I found the several sites which helped me

this is one of them :
http://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm

Moreover , I recommend on creating a black image and playing with it.

e.g.

black = false(500,500);
black(10:40:end,:) = true;
black = imrotate(black , 45 , 'crop');
ffim = fftshift ( abs ( fft2 ( black ) ) );
lffim = log ( ffim ) ;
fs(lffim,[]);


Is interpolation (interp1) better than FIR filtering when rational integers are close to 1?

You are looking at the wrong metric of "correctness".  Nearest neighbor is introducing significant discontinuities that are showing up as massive quantities of noise in the result.

The problem is that you should be comparing to the result you would have gotten if you had sampled at 26.25MHz in the first place.

Let's try it (sample a 12Hz sine wave at 26 Hz for 2 seconds):

numResamples = 52;
resampleRate = 26;
t1 = (0:1:numResamples - 1) * 1/resampleRate;
s1correct = sin(t1 * 24 * pi);
plot(t1,s1correct);
plot(t1,abs(fft(s1correct)));

Ok.  Now let's see what happened if we had sampled a 12Hz signal at 25Hz, and then interpolate up to 26Hz using matlab's linear interpolation.  (Actually, I'm using Octave because I'm on a machine without matlab today)

numSamples = 50;
sampleRate = 25;
t0 = (0:1:numSamples - 1) * 1/sampleRate;
s0 = sin(t0 * 24 * pi);
s1linear = interp1(t0, s0, t1, 'linear', 'extrap');
plot(t1,s1linear);
plot(t1,abs(fft(s1linear)));

Pretty good.  Now let's look at nearest neighbor interpolation.

s1nearest = interp1(t0, s0, t1, 'nearest', 'extrap');
plot(t1,s1nearest);
plot(t1,abs(fft(s1nearest)));

It's completely different, and filled with noise at all frequencies!!!  (The "beat pattern" in the nearest neighbor result is the beat pattern from sampling at 25Hz instead of 26Hz.)  The noise is because nearest neighbor is compressing the 12Hz signal to 12.48Hz, and then inserting an extra sample every 26th sample to get us back to 24 peaks.  That extra sample is the discontinuity that is creating all that extra noise.

I think you'll find that cubic interpolation is even better than linear.


Minimum possible sampling frequency for continuous time signal

This is one of those problems which ask for finding the sampling rate that avoids aliasing in the band of interest but don't care aliasing for the rest of the spectrum.

In this problem the real, baseband analog signal has a bandwith of $B = 1800$ Hz, which would require a minimum Nyquist sampling rate of $F_s = 2 B = 3600$ Hz when we are to avoid aliasing on the complete band of the signal; $-1800< f < 1800$ Hz which is depicted as the condition on figure-a below: (without amplitude scaling)

On the other hand the condition on the figure-b tells that if you want to preserve the signal only in the band of $-1000< f < 1000$ , then you can relax the aliasing condition to lower the sampling rate so that aliasing is prevented only in the band of ineterest; $-1000< f < 1000$ . That sampling rate you look is found by requiring on the shifted spectrum that, the leftmost cutoff frequency of right shifted replica is above $1000$ Hz of the unshifted replica, which is dispalyed on the figure-b.

$$ F_s - 1800 > 1000 \implies F_s > 2800$$

Hence the minimum sampling rate that will avoid aliasing in the band $-1000< f < 1000$ but otherwise will create aliasing in the range $ 1000< |f| < 1800 $ is $$ F_s = 2800 $$

Note that what you considered would be true if the signal was lowpass filtered to $1000$ Hz before being sampled. But in the problem no such antialiasing bandlimiting filtering is mentioned.


How to plot the periodic digital spectrum?

As a simple example, consider the discrete-time signal

$$x[n]=a^nu[n],\qquad |a|<1\tag{1}$$

where $u[n]$ is the unit step function. The discrete-time Fourier transform (DTFT) of $x[n]$ is given by

$$X(e^{j\omega})=\frac{1}{1-ae^{-j\omega t}}\tag{2}$$

Note that $\omega$ is a normalized angular frequency:

$$\omega=2\pi f/f_s\tag{3}$$

where $f_s$ is the sampling frequency.

Also note from $(2)$ that $X(e^{j\omega})$ is a $2\pi$-periodic function in $\omega$, i.e., it is periodic with period $f_s$ in the frequency variable $f$. This periodicity of the spectrum is a general property of discrete-time signals.

Now you can plot the spectrum of $x[n]$ (i.e., $X(e^{j\omega})$) as a function of the frequency $f$, and depending on the chosen range of $f$ and the chosen sampling frequency, you will plot a certain number of periods. This is shown in the following Octave/Matlab script:

a = 0.9;
fs = 1;
fmin = 0;
fmax = 2.5;
f = linspace(fmin,fmax,1000);
X = 1 ./ (1 - a * exp(-1i*2*pi*f/fs) );

plot(f,abs(X))
xlabel('f')
title('|X(f)|')
grid on