Yes your understanding is basically correct.
The 1st paragraph (2 lines) expresses the fundamental relation between the DFS and the DFT of a finite-length sequence $x[n]$ while the 2nd paragraph tries to put down the relation between the DFT $X[k]$ of a sequence and the DTFT $X(e^{j\omega})$ of it (assuming it exists).
However this 2nd paragraph shall better be rephrased like this: given a finite length sequence $x[n]$ of length $N$ whose DTFT is:
$$X(e^{j\omega}) = \sum_{n=0}^{N-1} {x[n] e^{-j\omega n}} \tag{1} $$
it will have its $N$-point DFT $X[k]$ defined and computed as the uniform samples of its DTFT $X(e^{j\omega})$ taken at $N$ frequency points $\omega_k = \frac{2 \pi}{N} k$ for $k=0,1,...,N-1$ .
It's important to recognise that some important class of aperiodic signals which are of infinite length, such as $x[n] = a^n u[n]$, will not have a valid $N$-point DFT representation for any (finite) $N$, nor will it have a valid finite DFS representation as it's not a periodic signal. Such a signal requires infinitely many samples for its DFT to exactly represent it.
In most practical cases, for an exponentially decaying sequence, a truncated (windowed) finite length, long enough, portion of it can be used to approximate the infinite length sequence. One can always find a large enough finite $N$ that will reduce the error of the approximation down to an acceptable level. (There's also an exact relation even for finite $N$ as will be shown below.)
For infinite length aperiodic signals, such as the above $x[n]=a^n u[n]$, whose DTFT is $X(e^{j\omega})= \frac{1}{1-a e^{-j\omega}}$ for $|a|<1$, one can think about the following: according to the relation of DFT to DTFT, the DFT $X[k]$ is the uniform samples of the DTFT $X(e^{j\omega})$ of that sequence $x[n]$ :
$$ X[k] = X(e^{j\frac{2\pi k}{N}})= \frac{1}{1-a e^{-j\frac{2\pi}{N}k}} ~,~  k=0,1,...,N-1 \tag{2}$$
And by assuming that these N-samples of DTFT, constitude a proper (valid) DFT representation of some N-point finite length sequence $x_a[n]$, we convert those $N$ DFT samples back into time-domain using an $N$-point inverse DFT. Then what's the relation between the finite length sequence returned by the N-point IDFT (of the N uniform samples of the DTFT $X(e^{j\omega})$ of $x[n]$) and the infinite length sequence $x[n] = a^n u[n]$..?
It can be shown that, the returned sequence, I denote $x_a[n]$, will be a time-aliased version of the infinite length sequqence $x[n]$ defined with the following relation:
$$ x_a[n] = \sum_r { x[n+rN] } \tag{3} $$
Which can be viewed as the periodic extension of the infinite length sequence $x[n]$ by N, however, the tails of $x[n]$ should sum up to finite values (converge) for $x_a[n]$ to exist. This will be the case for $x[n] = a^n u[n]$ with $|a| <1$ .
Note that for the class of complex exponential sequences $x[n]=a^n u[n]$, this aliasing can be corrected, based on an explicit relation between $x_a[n]$ and the parameter $a$ of $x[n]$. The relation between $x[n]=a^n u[n]$ and the aliased sequence $x_a[n]$ indicated by Eq.(3) turns out to be:
$$ x_a[n] = \frac{1}{1-a^N} x[n] \tag{4} $$
For other classes of sequences, such a relation may or may not exist, and an inversion from the aliased sequence $x_a[n]$ into the infinite length sequence $x[n]$, in general is not possible.
From Eq.4, the parameter $a$ of the actual exponential sequence $x[n] = a^n u[n]$ whose DTFT was sampled to obtain the N-point DFT $X[k]$ can be found as :
$$ a^N = \frac{ x_a[0] - 1 }{x_a[0]} \tag{5}$$
Note, however, that it's not a practical attempt.
Another important note about the relation between the sequence $x[n]$, its DTFT $X(w)$, and the DFT $X[k]$ of $x[n]$, is about a practical issue that arise when DFT is computed through an FFT call within a program, such as MATLAB/Octave.
Basically it's about an implicit alignment of the data seqeunce $x[n]$ to begin at $n=0$. DFT assumes that $x[n]$ is a causal sequence and it considers samples of input only for $0 \le n \le N-1$. This will cause a phase error at the computed DFT samples of a non-causal sequence if the DFT input is not properly constructed to yield the expected DFT result.
Consider we want to compute DFT of an even symmetric sequence $h[n]$ with a region of support $-M \le n \le M$ . Let $M = 2$ for the example.
$$h[n]=\delta [n+2] +2 \delta[n+1] + 3\delta[n] +2\delta [n-1] +\delta[n-2] \tag{6} $$
$h[n]$ is not a causal sequence; but a standard FFT function (withing MATLAB/OCTAVE) cannot know this, and treats the first sample of its input to be at $n=0$. This will cause a phase shift on the computation and the resulting discrepancy between the true DFT values of the actual non-causal sequence $h[n]$ and the constructed causal equivalent $h_c[n]$, should be remebered and be explicitly taken care of by the user, if the phase information of the DFT is important (Magnitudes will not be affected).
Lets call the causal extension of the input $h[n]$ as $h_c[n]$. We can define it in two ways as :
$$ h_c[n] = h[n-2]   ~~~,~~~ n=0,1,2,3,4,5 \tag{7}$$
or
$$ h_c[(n-2)_N] = h[n-2] ~~~,~~~ n=0,1,2,4,5 \tag{8} $$
where $(n-2)_N$ indicates a modulus_N operation; i.e. $$(n-2)_N = \text{mod}(n-2,N) \tag{9}$$
Note that Eq.7 treats the sequences $h[n]$ and $h_c[n]$ defined for all $-\infty < n < \infty$, and performs a linear shift on the non-causal sequence to obtain $h_c[n]$. On the other hand, Eq.8 treats them as periodic sequences with a predefined period $N$ and performs a circular shift on $h[n]$ to obtain $h_c[n]$.
$h_c[n]$ according to Eq.7 becomes $h_c[n] = [1,2,3,2,1] $,whereas according to Eq.8 it is $h_c[n]= [3,2,1,1,2]$.
In a more general setting, when the non-causal sequence has its non-zero samples for $a \le n \le b$, $a<0$, and is zero when $n < a$ or $b<n$, the linear shifted causal sequence of length $N \ge b-a+1$ is obtained as:
$$ h_c[n] = h[n+a] ~~~,~~~ 0 \le n < N  \tag{10} $$
and the circularly shifted causal sequence obtained as :
$$ h_c[(n+a)_N] = h[n+a] ~~~, ~~~ 0 \leq n < N \tag{11}$$ Eq.11 can be split into two pieces in the range $0\leq n < N$ as:
$$ h_c[n] = h[n] ~~~,~~~ 0 \le n \le b  \tag{11.a}$$
$$ h_c[n] = h[n-N] ~~~,~~~ N+a \leq n \le N-1 \tag{11.b}$$
Having defined the non-causal $h[n]$ and its N-point causal extensions according to Eqs. 10 & 11, the next step is to compute and relate their DTFTs:
For the linearly shifted sequence $h_c[n]$, it's easy to see that :
$$ H_c(e^{j \omega}) = e^{ j\omega a } H(e^{j \omega}) \tag{12}$$ This follows directly from the DTFT properties of a linear shift.
For the circularly shifted sequence, $h_c[n]$, the relation is more complex at first, but it can be seen that at the DFT sampling frequencies $w_k = \frac{2\pi}{N} k$, the relation simplifies to :
$$ H_c(e^{j \omega_k}) = H(e^{j \omega_k}) \tag{13}$$
This result in Eq.(13) may seem surprising at first but can be derived from splitting the DTFTs into two pieces as follows :
$$
\begin{align}
H_c(e^{j \omega}) &= \sum_{n=0}^{N-1} h_c[n] e^{-j \omega n} \\
&= \sum_{n=0}^{b} h_c[n] e^{-j \omega n} + \sum_{n=N+a}^{N-1} h_c[n] e^{-j \omega n} &\text{h_c[n] =0 for b< n < N+a} \\
&= \sum_{n=0}^{b} h[n] e^{-j \omega n} + \sum_{n=N+a}^{N-1} h[n-N] e^{-j \omega n} &\text{replace hc[] by h[] via Eq.11} \\
&= \sum_{n=0}^{b} h_c[n] e^{-j \omega n} + \sum_{n=a}^{-1} h[n] e^{-j \omega (n+N)} &\text{substitude n-N} \\
&= e^{-j\omega N} \big( \sum_{n=a}^{-1} h[n] e^{-j \omega n} \big) + \sum_{n=0}^{b} h[n] e^{-j \omega n} &\text{reorder the sums} \\
H_c(e^{j\omega_k})&= e^{-j\omega_k N} \big( \sum_{n=a}^{-1} h[n] e^{-j \omega_k n} \big) + \sum_{n=0}^{b} h[n] e^{-j \omega_k n} &\text{w_k =2pi k/N} \\
H_c(e^{j\omega_k})&=  \sum_{n=a}^{-1} h[n] e^{-j \omega_k n} + \sum_{n=0}^{b} h[n] e^{-j \omega_k n} &\text{e^(-j\omega_k N)= 1} \\
H_c(e^{j\omega_k})&=  \sum_{n=a}^{b} h[n] e^{-j \omega_k n} &\text{merge the sums} \\
H_c(e^{j\omega_k})&=  H(e^{j \omega_k}) &\text{RHS: DTFT of h[n] at wk} \\
\end{align}
$$
Hence we see that eventhough the DTFTs of circularly shifted causal sequence $h_c[n]$ and non-causal sequence $h[n]$ are not the same for arbitrary $\omega$, when evaluated at the DFT frequencies $\omega_k = \frac{2\pi}{N}k$ , they become equal :
$$ H_c(e^{j\frac{2\pi}{N}k}) = H(e^{j\frac{2\pi}{N}k}) ~~~,~~~ k=0,1,...,N-1 \tag{14}$$
Now, we can compute N-point DFTs by sampling the corresponding DTFTs at the DFT frequencies. Which indicates the following.
For the linearly shifted causal sequence $h_x[n]$ of Eq.10, we have
$$ H_c[k] = H_c(e^{j\frac{2\pi}{N}k}) = e^{j\frac{2\pi}{N}k a} \cdot H(e^{j\frac{2\pi}{N}k}) = e^{j\frac{2\pi}{N}k a} H[k] \tag{15} $$
and for the cirularly shifted causal sequence $h_c[n]$ according to Eq.11, we simply have
$$ H_c[k] = H_c(e^{j\frac{2\pi}{N}k}) = H(e^{j\frac{2\pi}{N}k}) = H[k] \tag{16} $$
Now in the Eq. 15 &16 above, the quantitiy $H[k]$ on the RHS is the desired DFT of the non-causal sequence $h[n]$, whereas the quantity $H_c[k]$ on LHS is the actually computed DFT by MATLAB/OCTAVE fft() function when feeded with linearly or circularly shifted sequences $h_c[n]$ respectively.
From Eq. 15 it's seen that for the linearly shifted sequence we have a phase error by $e^{j\frac{2\pi}{N}k a}$ and to get the correct (desired) DFT samples of $H[k]$ we have to reverse it; i.e.,
$$ H[k] = e^{-j\frac{2\pi}{N}k a} H_c[k] ~~~, ~~~k=0,1,...,N-1 \tag{17}$$
From Eq.16 we see that the DFT samples $H_c[k]$ of the circularly shifted causal sequence (Eq.11) $h_c[n]$ is exactly identical to the samples of the desired DFT $H[k]$ of $h[n]$. Hence we get the desired DFT result without error, if we feed the fft() function with the samples of the circularly shifted causal sequence.
Below is a MATLAB demonstration of the above discussion based on the given even symmetric sequence $h[n]$. DTFT of such an even symmetric real sequence will be real and even (hence will have zero phase) as given by:
$$H(e^{j\omega}) = \sum_{n=-2}^{2} {H[n] e^{-j\omega n}} = 3 + 4 \cos(\omega ) + 2 \cos(2\omega)$$
Therefore its $N=5$-point DFT $H[k]$ will also be a real, even, and zero-phase sequence as give below :
$$ H[k] =    [9.0000,    2.6180,    0.3820,    0.3820,    2.6180] $$
Lets compute the desired DFT $H[k]$ using MATLAB :
hc = [1 2 3 2 1];   % linearly shifted causal sequence hc[n] acc to Eq 10
fft(hc,5)           % Hc[k] = DFT of h[n] ?

The output is :
  9.0000  ,  -2.1180 - 1.5388i ,  0.1180 + 0.3633i ,  0.1180 - 0.3633i ,  -2.1180 + 1.5388i

This is not the expected result. MATLAB computes not the DFT of the original sequence $h[n]$ but the DFT of the linearly shifted verison $h_c[n]$
To recover $H[k]$ from $H_c[k]$ :
$$ H[k] = e^{ j~ 2 \frac{2\pi}{5} k} ~  H_c[k] ~~~,~~~ k=0,1,2,3,4$$
as the following line of code reveals:
exp( j* 2*pi *[0:4]*2/5 ).*fft([1 2 3 2 1],5)   % output is :
9.0000    2.6180    0.3820    0.3820    2.6180

This final desired result could be obtained directly (without correction) by using the circularly shifted causal sequence $h_c[n] = [3,2,1,1,2]$ as the input to the fft() :
H = fft([3 2 1 1 2],5)
9.0000    2.6180    0.3820    0.3820    2.6180


Preserve specific information in compressed sensing

As a result of the discussion, let me provide an answer here.
Say you want to perform "linear compression" àla Compressed Sensing to a vector, i.e., reduce a length-$N$ vector $\mathbf x$ to a length-$M<N$ vector $\mathbf y$ via a linear transformation $\mathbf \Phi$ such that $\mathbf y = \mathbf \Phi \mathbf x$. Further assume that $\mathbf x$ may be partitioned into $\mathbf x = [\mathbf x_1^T, \mathbf x_2^T, \mathbf x_3^T]^T$, where $\mathbf x_i$ are length $N_i$ vectors and $N = N_1 + N_2 + N_3$. Further, $\mathbf x_2$ is a part of your signal you do not want to compress so that the actual compression is only applied to $\mathbf x_1$ and $\mathbf x_3$.
The most general way of achieving this is by structuring $\mathbf \Phi$ as $$\mathbf \Phi = \begin{bmatrix} \mathbf \Phi_{11} & \mathbf 0 & \mathbf\Phi_{13} \\
\mathbf 0 & \mathbf I_{N_2} &\mathbf 0 \\
\mathbf \Phi_{31} & \mathbf 0 & \mathbf\Phi_{33}
\end{bmatrix} $$ where $\mathbf \Phi_{ij}$ is $M_i \times N_j$ for $i, j \in \{1,3\}$ and $M = M_1 + N_2 + M_3$.
Applying this $\mathbf \Phi$ we can see that $\mathbf y$ will be partioned as well, i.e., $\mathbf y = [\mathbf y_1^T, \mathbf y_2^T, \mathbf y_3^T]^T$ where $\mathbf y_1 = \mathbf \Phi_{11} \mathbf x_1 + \mathbf \Phi_{13}\mathbf x_3$, $\mathbf y_3 = \mathbf \Phi_{31} \mathbf x_1 + \mathbf \Phi_{33}\mathbf x_3$ and $\mathbf y_2 = \mathbf x_2$, as desired.


How we can obtain the time-integration of a sine signal?

No, this improper integral doesn't have a value. Since for the indefinite integral we have

$$\int\sin(\omega_0t)dt=-\frac{1}{\omega_0}\cos(\omega_0t)+C\tag{1}$$

the given definite integral would be

$$\int_{-\infty}^{t}\sin(\omega_0\tau)d\tau=-\frac{1}{\omega_0}\cos(\omega_0t)+\lim_{\tau\rightarrow -\infty}\frac{1}{\omega_0}\cos(\omega_0\tau)\tag{2}$$

but the limit in $(2)$ doesn't exist.


Confusion about collection of bits and symbols

In digital communications, information is transmitted in quantized form. We are constrained to transmit items that belong to a finite, discrete set. The items are called "symbols", and the set is called a "constellation". The symbols may be transmitted by changing a signal's frequency, phase and/or amplitude.

As an example, pulse amplitude modulation (PAM) is a technique where the symbols are pulse amplitudes, and the transmitted signal is a sequence of pulses. Let's say that you are using the 4-PAM constellation $\mathcal{C} = \lbrace -3, -1, 1, 3 \rbrace$, and the pulse shape is $h(t)$. Then, you group your bits 2 by 2 and assign them symbols as follows:

bits | symbol | pulse
-----------------------
 00      -3     -3h(t)
 01      -1     - h(t)
 11       1       h(t)
 10       3      3h(t)

Recall that a PAM signal can be written as $$s(t) = \sum_k c_k h(t-kT_p),$$ where $T_p$ is the pulse interval (the pulse rate is $R_p = 1/T_P$) and $c_k \in \mathcal{C}$ are the symbols. Then, transmission of the bit sequence $00111001$ would be done with the signal $$h(t) = -3h(t) + h(t-T_p) + 3h(t-2T_p) - h(t-3T_p).$$ Note that the information (bits) is conveyed by the "symbols" or pulse amplitudes (hence, PAM = pulse amplitude modulation). The task of the receiver is to recover the transmitted symbols from the received signal.

(Note that there's no unique or "official" way to map bits to symbols in 4-PAM or any other constellation; the mapping above is just an example. The only rule is to use Gray encoding, in order to minimize the bit error rate. And, or course, the mapping must be known at the receiver!)


Of course, the machine could be provided with three, four or more sensors, if the goal of the control unit is perform a balancing of the rotation. Actually three sensors should be the proper selection if the rotation is unconstrained in a biaxial plane, and two sensors, if constrained in a uniaxial direction.


wiener filter (non-causal) calculating covariance matrices

The expectation operator can not be applied in real.In Weiner filter you will require to find the value of R i.e. correlation matrix of input signal which is ideally E[X(n)X'(n)] and also P i.e  cross correlation matrix of input signal and desired signal E[X(n)d(n)]. Actual meaning of this expectation is that the correlation is obtained by using a large number of realisation of the input signal called ensemble .But in reality we have only one realisation and that is the input signal which you are going to apply to the weiner filter to calculate output. hence its exact value cant be find out. and therefore to use weiner filter in reality we require to apply various approximations to calculations of covariance matrix which finally result in various adaptive algorithms .As an example ,LMS algorithm which uses the following approximations for above matrices,
               R=X(n)*X'(n)   and P= X(n)d(n)
Now these values can be calculated using the only available realisation of signal.Although these Algorithms do not provide full convergence to weiner solution but approximates weights close to weiner solution. This error in solution is the result of approximation used. Weiner solution is actually the ideal solution of a adaptive filter hence it can not be attained because the Expectation operation can be applied in theory but practical implementation is impossible.
LMS,NLMS,BLMS,RLS these are some other algorithms that can also be used .each have its own good side and a bad side and have to be applied according to necessity.


Digital Filters deal only with samples, right?

You are almost right: digital filters do deal with samples, but a sample can be any numerical representation of a given signal value at a given instant (so in general, they may accept zeros or ones).

Moreover, a sample is usually represented by a binary word (e.g. 0001), so a digital filter actually deals with 0s and 1s.


The key is to use the "state" argument of the filtering commands. The idea is that filter commands such as lfilter can accept a vector specifying the inital state of the filter, and can also output their state at the end of the filter computation. In pseudocode:

let state = vector of zeros
let filtered_data = 0
for c in chunks
    y, state_new = filter(b, a, c, state)
    append y to filtered_data
    let state = state_new
end

See the docs for more details.


A General Question regarding OFDM symbol

Your question is a little confusing, but I'll give it a shot.

The 64 sample sequence output from the IFFT is a discreet sampled time domain sequence of complex values.  (A signal at complex base band)  The 3.2usec length comes from the sample rate (Fs) of the system.

In this case it sounds like Fs is 20MHz.  That Fs comes from the fact that there are 64 sub-carriers (FFT = 64)  therefore the frequency spacing to maintain orthogonality between sub-carriers is 1/Tu = 1/3.2usec = 312.5KHz.  Tu is the "used time", which in this case is that 3.2usec.  The total bandwidth is then 64*(sub carrier frequency spacing) = 64*312.5kHz = 20MHz.  Since the signal is complex Fs = BW.

Fs = BW = NFFT / Tu = 64 / 3.2usec = 20MHz

In my experience an OFDM symbol (or FFT symbol) refers to the data without the cyclic prefix added with its length equal to the FFT size.  I call the OFDM symbol plus the cyclic prefix and OFDM frame. Other people may have a different definitions.  It looks like this author uses the same definition as me not including the 0.8usec cyclic prefix in what they call an OFDM symbol.  That 0.8usec is also called the guard time and acts like a spacing in between successive OFDM symbols offering multipath protection.

Not sure if that will help or not.


There are many varied impairments in a transceiver that we could call "noise", (see my graphic below), and ultimately I resort to simulations with good noise and impairment models for each effect, that is very dependent on the receiver type and operation.

This question is to broad to give a specific answer. However, let me describe my treatment of colored noise in the receiver in case this adds some insight. In most cases my approach is to first have a thorough understanding of the frequency domain characteristics of the noise process involved, and then to evaluate the equivalent frequency response of the receiver translated to baseband prior to final detection, and integrate the total noise across that bandwidth to get the total noise power in band. The signal bandwidth will have a low frequency limit (typically set by carrier and timing recovery loops, or the duration in time that a packet is observed-- you need to observe a signal for infinite duration to go down to true DC so it will always be a high pass!) and then an upper frequency limit set by the symbol rate if the receiver is optimized (typically a percentage over R/2 where R is the symbol rate and the percentage is based on the pulse shaping factor), or the actual pre-detection bandwidth if not optimized.

Here is an example with phase noise specifically. It is usually clear that the receiver will filter out the higher frequency components beyond the channel bandwidth, but often there is confusion with how to handle a low frequency limit, and knowing that the phase noise keeps increasing as we approach DC, we could get very different answers depending on what low frequency limit we choose!

If there is a carrier tracking loop involved, the low frequency limit is often established by that tracking loop, since it works by tracking out phase rotation (so slow rotations or movements of phase within the loop BW of this tracking loop would be tracked out).  Depending on the loop this could be a first order high pass response as I show in the figure below.

It is the resulting phase noise density after being multiplied by this "filter" of the receiver that I would integrate (and double as half of the band is shown) to assess the phase noise impact on my receiver SNR.

The key point in this example is that with a thorough understanding of the frequency domain characteristics of the noise (the phase noise response, and beyond the scope of this example, how it maps to my receiver channel of interest), and a clear understanding of the receiver both in the low frequency limit and the upper frequency limit of the baseband channel with all operations prior to final detection; one can integrate the resulting colored noise density over the band of interest to establish the total noise in the receiver prior to detection.  I also do time domain histograms to capture and access the effects of different noise distributions, but in most case the combination of multiple independent noise sources makes the Gaussian estimate reasonable.


How to add repeated beeps to sound created with a pause in between each beep in MATLAB

Any pause can be done by using vector/matrix concatenation. Let's say you want to make a signal with pause of length P after it, then do the following:

P = 400;
beep = [sound_array; zeros(P, 2)];

Now if you want to make arbitrary number K of your beeps then use the repmat function as follows:

K = 5;
beeps = repmat(beep, K, 1);

Please read more about matrix concatenation here if you wish.

I believe that with your current code, because you use sound function, signals will be played at the same time. I suggest you to use the pause function to separate playbacks, i.e.:

sound(...)
pause(1)
sound(...)

To create a 1 second pause between playback starts.


Window functions with rippleless spectra

That's an interesting question, I personally don't know if, in theory, that's possible or not (would be interesting to investigate). But as a practical matter perhaps you could design a Chebyshev Type II filter, take the IFFT, and use a sufficiently large number of taps. Chebyshev Type II filters have ripple in the stopband but not the passband (although that is monotonically decreasing). I'm not sure how much the truncation would influence the ripple, or if you can handle that large of a FIR. But it might get you close enough depending on your application.

This page details some of the results.
http://iowahills.com/B2PolynomialFIRFilters.html


The phase response is clearly non-linear, as mentioned in a comment. Note that phase jumps at the zeros of the frequency response occur depending on how the phase is defined, but these phase jumps don't make the phase non-linear. However, in the given example the phase is clearly non-linear at low frequencies.

You could also look at the corresponding impulse response. If the impulse response has a finite length and if the coefficients are symmetrical or anti-symmetrical with respect to the center, then the phase is linear.

Another method is to look at the pole-zero diagram of the transfer function. If the system is causal then all poles must lie at the origin and all zeros must either lie on the unit circle, or they must occur in quadruples $z_0$, $z_0^*$, $1/z_0$, $1/z_0^*$.


How can I transform a time signal into a third octave spectrum?

maybe either the bandpass filters or the peaking EQ filters in the audio EQ cookbook can help.  bandwidth is defined in octaves in that.  the output, $y$, of each 1/3 octave BPF filter should be squared and low-pass filtered to get the mean square, $\overline{y^2}$.  and then

$$dB = 10 \log_{10}\Big(\,\overline{y^2}\,\Big)$$

is applied to that to get dB.

are you planning on doing this to the whole audible spectrum?  31 bands?


Difference between a digital lock-in amplifier and a FFT when extracting phase of a signal?

Is there a difference between the FFT and a lock-in amp?

Yes, two of them:

The FFT assumes the signal at its input to be periodic.

What this means about your FFT integrals is that they are missing a phase variable (let's call it $\phi_r$) which will be random because it depends on the ratio of your signals period (or its frequency), the window of the FFT and the initial phase that the FFT will "catch" the signal on the first frame (and any subsequent phase changes whether desired or not). By the way, Frequency is integrated Phase.
Consequently, your $\theta_{LI}$ will depend on that $\phi_r$ (which is all over the place).
Intuitively now, the FFT assumes the signal to be periodic. That is, it repeats in the same way both to the left and to the right of the observation window. This is "alright" if the window happens to be an integer multiple of the signal's period but this is extremely unlikely because of noise and possible changes in the phase of the signal. It will therefore sound like a whistle with regular "pops" because of the discontinuities and at each "pop" the phase estimation will be disturbed. There are workarounds to this (and this one) and this brings us to the second difference.

There is a concept in the FFT called Spectral Resolution.

Spectral resolution relates the physical frequencies (in Hz) with each of the distinct harmonics that the discrete FT evaluates its integrals at. Therefore, you would have to accurately calculate the length of the window of the FFT with respect to the sampling frequency and the frequency of the input signal so that it lands exactly on one of the bins. And after doing this of course, it would seem a waste to be evaluating the FT for all harmonics when all you are interested in is just one of them.

Having said this, let's try to tackle the rest of the questions:

  So does this mean that I would always get the same result whether I use a digital lock-in amplifier, or a digital FFT on the same data set?

Provided that the above details are taken care of, the result will be almost identical. I say almost because there is another consequence from #1 above that is not exactly obvious. Even if you were to use overlap-add or overlap-save, you would still have no control over the initial phase at which the FFT would "catch" the incoming signal. And once the FFT "starts", it's own local oscillators are not going to adapt to the incoming signal. Therefore, both techniques would provide information about phase but the FFT's estimate would suffer by a systematic error proportional to the initial phase of the incoming signal. The lock-in amp takes care of this with its Phase-Locked-Loop (PLL). For (much) more information, please see this link.

  What is the benefit of using a lock-in amplifier?

Should be becoming clearer by now. A lock-in amp can return amplitude, power and phase information at one specific frequency with the accuracy of its estimate depending on the clarity of the reference signal (the local oscillator).

  Is there no advantage to the lock-in amp in my application?

On the contrary, there is a massive advantage. And this advantage is coming from the use of the PLL. The PLL is an automatic control system that adjusts the frequency (and therefore phase) of a local oscillator to match the frequency and phase of the incoming signal. It therefore can "adjust its frame of reference" and it returns much more accurate phase information.

  What is then an example of an application that the lock-in amplifier is best suited for?

Please see above PDF.

Hope this helps.


Conceptual Question on Least Squares Estimation Method

Gauss did a two-way proof.  In rough terms, he not only showed that least squares gives the best fit to data containing a standard normal distribution of errors, but that if least squares is always giving the best fit, that implies that the distribution is indeed standard normal (also now called Gaussian).

The correlation coefficient, often computed in conjunction with a least square fit, is usually assumed to be parameterizing an underlying Gaussian distribution.  If the distribution is otherwise, it may well infer nonsense.


Compressive sensing: numerical generation of RIP matrices

You can't prove RIP through numerical exploration of all possible cases. If you are interested in numerical analysis I suggest to use Coherence instead,  however Coherence is not as strong condition as RIP. Note that RIP is a tool which is used for Mathematical proofs.
I explain how to check for RIP, in case your are interested.
As I said Coherence is a weaker condition in comparison with RIP, however, it is numerically feasible to check. Coherence of matrix A (which is sensing matrix multiplied in sparsity basis) is defined as below:

where $a_i, a_j$ are columns of matrix A, and $u_A$ is coherence factor. Its computational complexity is of $O(n^2)$ which in many cases is feasible.
For checking coherence, you simply need to compute the given type of correlation between each two different columns of A matrix and among those values chose the max. Based on Welsh Inequality we'll have:

In CS sampling $M<<N$, so the above relation simplifies to $u_A>sqrt(1/M)$, therefore, to satisfy RIP using coherence, number of measurements ($M$) must be of order of $O(K^2)$ where $k$ is number of non-zeros. This number of measurements is so greater than $Klog(N/K)$.

Another conditions instead of RIP, is stRIP which is much much simpler to prove.

Ref: [A.Amini,"Deterministic Compressed Sensing" PhD Thesis, Sharif University of Technology, Iran, 2011]


What is the standard way to handle images of different sizes for classification?

Not sure about your specific application, but the Viola Jones face detection algorithm successively scans the candidate picture for faces in all possible sizes.

I.e., the algorithm search space is both the object (face) location and its size. The beauty of VJ's "cascade architecture" is that it rejects most hypotheses quickly, thus making the search feasible in real-time.


How do I zero pad time samples if I only have frequency samples?

assuming I understand it, your question is interesting.  I'm "shooting from the hip" here.  I wonder if it would be useful to you to interpolate your two freq-domain sequences by a factor of two before the multiplication in the freq domain, and then performing your IFFT of the product?
[-Rick-]


Interpreting the cross correlation results

Yes and no. In principle, you can use the peak of your correlation function. However, it is not at 10000. The correlation function is symmetric around 0, so your peak is actually much closer to zero than you think. This is one of the reasons why xcorr returns two paramters, one for the lags at which the function is calculated. The correct way of plotting the the correlation function would thus be a plot(s,r).
Secondly, you should try zooming in a bit. Since your time signals have 10000 samples, your acf has ~20000, making it hard to see the exact maximum. If you do zoom in you'll notice it's not exactly at zero but two samples to the right.
Why two samples? Because you have 10 samples per period and the shift between sine and cosine is a quarter period, so it's actually 2.5 samples. You notice the problem, with the sampling you've chosen it is a little hard to measure the delay precisely.
Try reducing your fm to a small value, like 10, this makes the correlation function easier to interpret. Here is an example:

In this example, the maximum occurs at sample 247, which is equivalent to 247*ts = 24.7 milliseconds.
In case you wonder about the triangular shape, this has to do with the fact that while at lag 0 (in the middle) we can "compare" the entirety of the two signals whereas for larger lags, we need to truncate both signals to shorter versions. The larger the lag, the fewer samples we have to correlate. For this reason, we typically check the ACF closer to lag zero.


Real-valued ringing when zero-padding odd-length FFT

By zeroing out high frequency bins, you have effectively multiplied the spectrum of the signal with a rectangular function.  Multiplication in frequency is convolution in time and the Fourier pair of a rect is a sinc.  So what you really did was convolve the time domain signal with a sinc with the width of the main lobe of the sinc inversely proportional to length of the rect.  This is why the numerous filter designing techniques like Parks-McClellan design in what's called a "transition region" or "transition" band so that there is not an instantaneous change in the frequency response of the filter.  These filter design techniques are important because the "ideal" filter like you used has such undesirable effects in the time domain.


Is the following system stable and causal?

You are on the right track, but you are confused.
First, regarding the calculations you could simply do
$$y[n]=2\left(\frac{1}{2}\right)^n(u[n]+\delta[n+1])=2\left(\frac{1}{2}\right)^nu[n]+4\delta[n+1]$$
which gives you

$$Y(z)=\frac{2z}{z-\frac{1}{2}}+4z,\,|z|>\frac{1}{2}$$
The impuse resposne is the first difference of the step response:

$$\begin{align}
h[n]&=y[n]-y[n-1]\\
&=2\left(\frac{1}{2}\right)^nu[n]+4\delta[n+1]-4\left(\frac{1}{2}\right)^nu[n]\\
&=4\delta[n+1]-2\left(\frac{1}{2}\right)^nu[n]
\end{align}$$
Obviously, $h[-1]=4$ and the system is not causal.
But here you are mixing up finding the inverse $z$-transform of anti-casual signals with system being non-causal. ROC is $|z|>\frac{1}{2}$ and the system is non-causal full stop.


Easy way to get rid of noise in a hand drawing

You can use opening on your image http://php.net/manual/en/imagick.morphology.php with a kernel big enough to erode the small objects


May I know how to find Time frequency power spectrogram of a signal?

Short-time Fourier transform

Traditional Fourier analysis provides no information about evolution of frequency components over time. However, Short-time Fourier transform (STFT), a method which applies a short-time window to the signal and performs a series of Fourier transforms within this window as it slides across all the data, can overcome this limitation, providing a time–frequency representation of the signal.

In MATLAB, try

doc spectrogram

Continuous Wavelet Transform

Alternatively, continuous wavelet transforms (CWT) provide a useful approach in investigating non-stationary signals, this is often regarded as an “optimal” solution with regard to time and frequency resolution. Reportedly, CWT perform better at detecting abrupt changes in non-stationary signals compared to STFT.

In MATLAB, try

doc cwt

Generally, EEG signals have a low signal to noise ratio (SNR), therefore it is preferred to use ensemble averaging over trials for your time-frequency analysis (TFA), should the study design allow it.


Let me add the following graphic to the great answers already given, with the intention of a specific and clear answer to the question posed. The other answers detail what linear phase is, this details why it is important in one graphic:

When a filter has linear phase, then all the frequencies within that signal will be delayed the same amount in time (as described mathematically in Fat32's answer).
Any signal can be decomposed (via Fourier Series) into separate frequency components. When the signal gets delayed through any channel (such as a filter), as long as all of those frequency components get delayed the same amount, the same signal (signal of interest, within the passband of the channel) will be recreated after the delay.
Consider a square wave, which through the Fourier Series Expansion is shown to be made up of an infinite number of odd harmonic frequencies.
In the graphic above I show the summation of the first three components. If these components are all delayed the same amount, the waveform of interest is intact when these components are summed. However, significant group delay distortion will result if each frequency component gets delayed a different amount in time.
The following may help give additional intuitive insight for those with some RF or analog background.
Consider an ideal lossless broadband delay line (such as approximated by a length of coaxial cable), which can pass wideband signals without distortion.
The transfer function of such a cable is shown in the graphic below, having a magnitude of 1 for all frequencies (given it is lossless) and a phase negatively increasing in direct linear proportion to frequency. The longer the cable, the steeper the slope of the phase, but in all cases "linear phase".
This makes sense; the phase delay of 1 Hz signal passing through a cable with a 1 second delay will be 360°, while a 2 Hz signal with the same delay will be 720°, etc...
Bringing this back to the digital world, $z^{-1}$ is the z-transform of a 1 sample delay (therefore a delay line), with a similar frequency response to what is shown, just in terms of H(z); a constant magnitude = 1 and a phase that goes linearly from $0$ to $-2\pi$ from f = 0 Hz to f = fs (the sampling rate).

The simplest mathematical explanation is that the a phase that is linear with frequency and a constant delay are Fourier Transform pairs. This is the shift property of the Fourier Transform. A constant time delay in time of $\tau$ seconds results in a linear phase in frequency $-\omega \tau$, where $\omega$ is the angular frequency axis in radians/sec:
$$\mathscr{F}\{g(t-\tau)\} = \int_{-\infty}^{\infty}g(t-\tau)e^{j\omega t}dt$$
$$u = t - \tau$$
$$\mathscr{F}\{g(u)\} = \int_{-\infty}^{\infty}g(u)e^{-j\omega (u+\tau)}du$$
$$ = e^{-j\omega \tau}\int_{-\infty}^{\infty}g(u)e^{-j\omega u}du$$
$$ = e^{-j\omega \tau}G(j\omega)$$


How can I spread a modulated signal over walsh code

That's right, you can do spread using the kroon function, but be careful what you spread is not the data, it's the whole possibilities of data which is mod.The data is one of them, which is either, 0.7+0.7i or 0.7-0.7i ...etc.

then you should know how you can do despread in the receiver side too.

Good luck


Why we take Laplace Transform of functions which converged using Fourier Transform

There is a large class of functions for which both the Fourier transform and the Laplace transform exist, and for which one can be obtained from the other by setting $s=j\omega$. (Note that even when both exist, the latter need not be the case). So for this class of functions, obtaining the Laplace transform from the Fourier transform (or vice versa) does not require any additional work.

Example: $$\begin{align}x(t)&=e^{-at}u(t),\quad a>0\\
\text{Fourier transform:}\quad X(j\omega)&=\frac{1}{j\omega +a}\\
\text{Laplace transform:}\quad X(s)&=\frac{1}{s +a},\quad |s|> -a\end{align}$$

What the Laplace transform offers is a description in the complex $s$-plane, such as poles and zeros of transfer functions, from which many system properties can be easily deduced. The Fourier transform only shows the behavior on the frequency axis (i.e., the spectrum).


why difference between this example and matlab result

When you use Matlab's function butter, all frequencies need to be normalized by $f_s/2$, so your second call to the function butter is the correct one, the first one is wrong. And, as usually, Matlab is correct. The difference between the solution in your book and Matlab's solution is that the function butter uses prewarping to guarantee that after the bilinear transformation, the resulting cut-off frequency of the digital filters equals the specified cut-off frequency. For whatever reason, the solution in your book doesn't take the frequency warping of the bilinear transform into account, so the resulting cut-off frequency of the digital filter does not exactly equal the desired $\omega_c$.

The frequency warping caused by the bilinear transform used in your book maps the cut-off frequency $\omega_c$ of the analog filter to the following cut-off frequency of the resulting digital filter:

$$\tilde{\omega}_c=\frac{2}{T}\arctan\left(\frac{\omega_cT}{2}\right)=2\pi\cdot 2906.8\tag{1}$$

instead of $\omega_c=2\pi\cdot 3000$.

You can check the actual cut-off frequencies of the digital filters using the following Matlab commands:

wc = pi/5;                          % normalized cut-off frequency
wc2 = 2*atan((2*pi*3000)/30000/2);  % normalized warped cut-off frequency
[b,a] = butter(2,wc/pi,'low');      % Matlab's solution
b0 = 0.063964*[1,2,1];              % solution from book
a0 = [1,-1.168261,0.424118]
abs(freqz(b,a,[0,wc])).^2           % check |H|^2 at DC and at wc;
                                    % (should be [1,0.5])
abs(freqz(b0,a0,[0,wc])).^2         % [0.99999   0.46637]
abs(freqz(b0,a0,[0,wc2])).^2        % [0.99999   0.49999]

This shows that the filter in the book indeed has cut-off frequency $\tilde{\omega}_c$ given in $(1)$ instead of $\omega_c$.

If you don't understand the frequency warping of the bilinear transform, I'd suggest you study chapter 11 of this excellent and free book by S.J. Orfanidis.


Series connection of digital PT1 filters

You could have mentioned where you got that formula from. But anyway, I'll try to reconstruct it.

A first order discrete-time low pass filter has the following transfer function:

$$H(z)=\frac{1+c}{1+cz^{-1}},\qquad -1<c<0\tag{1}$$

where the constant $c$ determines the cut-off frequency. If you cascade $N$ of these filters, the total transfer function is

$$H_N(z)=\prod_{k=1}^NH(z)=\frac{(1+c)^N}{(1+cz^{-1})^N}\tag{2}$$

Using the Binomial Theorem this can be written as

$$H_N(z)=\frac{(1+c)^N}{\displaystyle\sum_{k=0}^N\dbinom{N}{k}c^kz^{-k}}=\frac{(1+c)^N}{\displaystyle\sum_{k=0}^Na_kz^{-k}}\tag{3}$$

with denominator coefficients

$$a_k=\dbinom{N}{k}c^k\tag{4}$$

Note that the numerator $(1+c)^N$ is just there to make sure that the scaling is correct, i.e., that the transfer function equals $1$ at $z=1$ (DC).

Following MBaz's suggestion, I add the relationship between the constant $c$ in $(1)$ and the filter's $3\,\text{dB}$ cut-off frequency $\omega_c$. For this we have to solve

$$|H(e^{j\omega_c})|^2=\frac12|H(1)|^2=\frac12\tag{5}$$

From $(1)$ we get

$$\frac{(1+c)^2}{1+2c\cdot \cos(\omega_c)+c^2}=\frac12\tag{6}$$

from which we obtain after some elementary algebra

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


Using Matlab/Simulink to build Reader based on Software Define Radio

You probably can. But: the AD9361 is definitely a complex beast, and I honestly think you're setting yourself up for more trouble than your solving when using a model of a device with a couple thousand different settings.

Unless you're planning to build an SDR out of an actual AD9361, I'd recommend not using them in a simulation. It's probably among the most complex models available in Simulink!

You shouldn't be building your own SDR, but rather buy one; as far as I remember, these typically need relatively high bandwidth, and designing the digital side of a PCB that carries multiple Megasamples of complex data per second into your computing unit will not be a beginner's project (aside from being far more expensive to buy and manufacture in single digit quantities than just buying a commercially available SDR).

The question really is whether you need to have a simulation that incorporates the actual SDR hardware – wouldn't for a very long time a simple idealized (IQ) mixer be representative enough? Aren't you starting off at the wrong end by trying to model something that works quite close to ideal very intimately instead of trying to get a system understanding based on idealized models?

For example, your flow graph top row isn't correct, to be honest. You used "QAM modulation" for what should be a quadrature modulator (QAM modulation is a baseband symbol mapping method, not a "modulation" as in "shift of frequency").


Fourier transformed acceleration from Fourier Transform of velocity and position using Omega arithmetic

The problem is that you are analyzing sine and cosine without using a smooth window e.g a Hanning window. Your phase range is from -100 to 100 giving 100 / $\pi$ = 31.83 periods.This means that you have a discontinuity at the end of the FFT time record. So you are not analyzing a sine or cosine as you believe. Using a Hanning window will suppress the discontinuities and give a more meaningful result.


I have a signal s(t) of duration T with limited bandwidth.

No, you don't.

A signal cannot be finite in both time and frequency at the same time. So the whole question about length in both domains comes down to how exactly you define it. Typically you need a criteria like "less than -100dB", "99.99% of the entire energy" or something like this.

These criteria are specific to your application. Maybe -40dB is plenty, maybe you need more than -200dB. Once you have a criteria defined you can apply it to your signals and determine the "lengths", but the result will always be tied to an application specific definition.


To answer your last question, the impulse response having infinite length affects the output length in every possible way. In general, the output sequence of an IIR filter has an infinite length, even for a finite length input sequence. Of course, for stable filters and finite-length input sequences, the magnitude of the output samples decreases, and the output values become virtually zero after a while.

There are, however, specific input sequences that result in cancellations in the output, and, consequently, the corresponding output sequences have only a finite number of non-zero elements. A simple example is the first-order system

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

You can verify yourself that for the specific input sequence $x[n]=\delta[n]-a\delta[n-1]$, the output sequence is finite and simply equals $y[n]=\delta[n]$. However, note that the use of finite precision arithmetic in an actual implementation may still result in non-zero output samples, even in cases when ideally the output samples should cancel.


You are correct. When using uniform sampling, the sampling rate is unrelated to the quantization noise. Only the quantization step and the signal's dynamic range are relevant for calculating the quantization SNR.


Is Entropy Coding in Wavelet equivalent to Huffman Coding in DCT

The goal of entropy coding is to code your data such that the length of the code associated with a symbol is close to $-\log_2 p$, where $p$ is the probability of the symbol (Shannon's limit). This can be achieved through different techniques: Huffman coding is one of them, but arithmetic coding is another option - or if the distribution of the coefficients can be known in advance, a pre-computed code can be specified. The choice of an entropy coding method is motivated by performance/complexity trade-offs, and by the known (or unknown) distribution of the coefficients. If two transforms generate the same distribution of coefficients, the same entropy coding methods can be used for both.

Since the diagram from your question hints that your application is JPEG2000 compression, I suggest you to look into the JPEG2000 specifications to see which entropy coding method they have actually specified.


You are right. If you really want to send the sampled signal $$x_s(t) = {\delta}_{T_s}(t) x_c(t) = ( \sum_k \delta (t- k T_s) ) x_c(t) = \sum_k x_c(k T_s) \delta (t-k T_s)$$ then this weighted impulse train constitudes an infinite amount of (shifted replicas of base period) bandwidth in frequency domain.

Therefore this signal is not transmitted anywhere... Indeed the impulse sampling cannot be realized in practice either. It's just a mathematical model which provides a quantitative description of the conversion process from a continuous time-signal into a discrete-time sequence whose samples have values equal to the continuous-time signal at those sampling impulse instants.

The sampling process is used as the initial stage of conversion from continuous to discrete time domains and it's the sequence of those samples which are processed later on for whatever purpose you have.


hint only, since this is homework and you could get a lot out of this!
The mathematical operation a filter performs between filter impulse response and signal is called __________(1). When we process both filter and signal with an FFT (which is just a fast implementation of __________(2)), then that operation on the originals is equivalent to __________(3) after the (2).
To do (3) on both (2)'ed vectors, these vectors need to be of the same  __________(4), here 256. Often, the signal itself is longer than 256 - (length of filter). This is a problem, because (1) under the (2) is not  __________(5) as the filtering operation normally is, but  __________(6) due to the periodic nature of the (2).
To alleviate this, a segmented approach is chosen.


Tracking a Sine Wave with Kalman Filter - How to Account for Offset (DC Signal)?

Well, in continuous time, a sinusoid with a bias can be seen as the output of the linear system
\begin{align*}
 \begin{bmatrix}\dot x_1\\\dot x_2\\\dot x_3\end{bmatrix} &= \begin{bmatrix}0 & 1 & 0\\-\omega^2&0&0\\0 &0 &0\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}\\
 y &= \begin{bmatrix}1&0&1\end{bmatrix}
\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}
\end{align*}
As the initial condition $(x_1(0),x_2(0),x_3(0))$ vary in $\mathbb R^3$, the system generates all the possible $y(t)$ of the form
$$
y(t) = a\cos(\omega t+\phi) + b
$$
The system is observable for all $\omega\ne 0$,therefore a Kalman Filter is possible. With $x:=col(x_1,x_2,x_3)$ and
$$
A=\begin{bmatrix}0 & 1 & 0\\-\omega^2&0&0\\0 &0 &0\end{bmatrix},\qquad C=\begin{bmatrix}1&0&1\end{bmatrix}
$$
The filter has the form
\begin{align*}
\dot{\hat x} = A\hat x -L(C\hat x-y)
\end{align*}
Where $L$ is chosen according to the Kalman theory, on the basis of the covariance matrices of the noises and the solution of the Riccati Equation. For discrete time systems the procedure is the same.
EDIT:-----
If you only use a model for the oscillator you are actually restricting the $A$ map to the $x_1$-$x_2$ subspace. That's equivalent to take $C=[1\;0\;0]$ in the above system. In this latter case you loose observability (actually you do not have neither detectability) and therefore you cannot observe the bias.


FFT showing evenly spaced signals when tuned to ~56 MHz down to 25.680

These periodic spectral lines look like they might be caused by something like a sparking power supply, DC-DC converter (for LED lighting or USB power, etc.), or high powered digital drivers running at approx. 150 kHz. You might want to walk around the house or neighborhood with a pocket AM radio to see if you can DF this RF noise (determine the direction its coming from).

RTL-SDRs USB dongles based on R820T2 + RTL2832U also have a lower frequency limit somewhere around 52 MHz (where they go deaf or alias or won't lock).  26.68 appears to be near a folding point in the sample rate aliasing, e.g. right around Fs/2, where any spectrum plus its own complex conjugate can partially cancel.


Problem identifying the analytic expression of such determined signal

The simple answer is “Linearity”

The Fourier Transform is linear.  $g(t)$ can be expressed as a sum of signals that have known, as in there are tables of transform pairs, in most books, probably also in your course notes, along with the shift and scale properties. Look at the tables you have and find entries that you can add to get $g(t))$.  You already found rect().  Look for a triangle piece.

A full triangle is the convolution of 2 rects.  A half triangle is a full triangle multiplied by a step.


Mean Square Continuity of Random Process

The whole point of the concept of mean square continuity is to avoid strict assumptions about the continuity of every realization of $X(t)$.

The given proof simply shows that continuity of the autocorrelation function implies that the limit $\lim_{t\rightarrow t_0}E\{[X(t)-X(t_0)]^2\}$ exists and equals zero.

Note that continuity of the mean of $X(t)$ is a consequence of mean square continuity. Since

$$E\{[X(t)-X(t_0)]^2\}\ge E^2\{X(t)-X(t_0)\}\tag{1}$$

it follows from

$$\lim_{t\rightarrow t_0}E\{[X(t)-X(t_0)]^2\}=0$$

that

$$\lim_{t\rightarrow t_0}E\{X(t)-X(t_0)\}=0$$

must hold, i.e., the mean of $X(t)$ is continuous.


Farrow fractional delay filter - Range of delay?

Do you think this link could help you?

http://users.spa.aalto.fi/vpv/publications/vesan_vaitos/ch3_pt1_fir.pdf

I'm not sure I understand 100% of your problem, but IIRC when implementing fractional delays with say a FIR filter of 13 coefficients, you should stick with a group delay between 6 and 7.


Let $y(t) = x(t^2)$, and $a$ be such that:

$$
\begin{align}
  y(at)
   &= x((at)^2) \\
   &= x(a^2t^2) \\
   &= x(bt^2)
\end{align}
$$
In other words, $a^2 = b$.

The Fourier Transform time-scaling property states that given the Fourier Transform pair:
$$
\begin{align}
  y(t) \leftrightarrow Y(w)
\end{align}
$$
and $a>0$, we have the relation:
$$
\begin{align}
  y(at) &\leftrightarrow \frac{1}{a}Y(\frac{w}{a})
\end{align}
$$

Thus,
$$
\begin{align}
  x(bt^2) &\leftrightarrow \frac{1}{a}Y(\frac{w}{a}) \\
          &\leftrightarrow \frac{1}{\sqrt{b}}Y(\frac{w}{\sqrt{b}})
\end{align}
$$

Now, if we define
$$
\begin{align}
Y(w) = G(w^2)
\end{align}
$$
Then,
$$
\begin{align}
Y(w/a) &= G(w^2/a^2) \\
Y(\frac{w}{\sqrt{b}}) &= G(\frac{w^2}{b})
\end{align}
$$

Yielding,
$$
\begin{align}
x(bt^2) &\leftrightarrow \frac{1}{a}Y(\frac{w}{a}) \\
        &\leftrightarrow \frac{1}{\sqrt{b}}G(\frac{w^2}{b}) \\
        &\leftrightarrow \frac{1}{\sqrt{b}}G((\frac{w}{\sqrt{b}})^2)
\end{align}
$$
Or, in prose (1/sqrt(b)) times a function of (w/sqrt(b))^2.


How would one design a (quasi) linear phase adaptive notch filter for a single complex tone?

I might consider doing a complex modulation of the frequency you want to notch down (or up) to 0 Hz.  Then run your favorite deep DC blocker (plus phase adjustment if needed).  Then complex remodulate the result back to the original unshifted spectrum.

DC block: For non real-time, just subtract the entire signal’s average.  Or a long (weighted?) moving average if you can tolerate the latency (Rick Lyons has some notes on this technique).


What is the best way of filtering EEG signal for the specific frequency band?

Perhaps you should try some of the filter design tools that are available online.
You can define the required filter behavior (Cutoff frequency, Ripple ..) and you will get a design solution that fits your needs.
This kind of tools will help you know what is possible to achieve in terms of filtering performance given the constraints that you detailed
Checkout the following ones ( this is not an exhaustive list) :
http://t-filter.engineerjs.com/
https://webench.ti.com/filter-design-tool/filter-response
Best regards,
Mourad


I suggest you do some more reading about the Fourier transform, but the short answer to your question is yes you can calculate the FFT for an arbitrary vector and discard the imaginary part, but it will no longer represent the Fourier transform of that signal. You can then perform the inverse FFT on that purely real spectrum, but it will not give you the original signal back (since it was not the actual spectrum of that signal).

If your starting signal happens to have even symmetry, then the spectrum will be purely real anyway, so you will be able to invert the real part of the FT and get your signal back.


What is wrong with RMS and Maximum delay spread Matlab code?

Not an answer but some guidance how to post this type of questions which is too long for a comment (@mods: feel free to delete if that's inappropriate).

We typically don't debug code. Focus on algorithmic questions. If possible ask about equations, not the implementation of the equations in code.
Define your terms (like "maximum delay spread") using an external reference or preferably with an equation to make sure we are all talking about the same thing.
If you use an equation (which you should), define all symbols used in the equation.
If you post code, make sure it can be run. We don't have your data file so we can't really reproduce your problem.
Do basic debugging first yourself: Run your code with the most simple signal with a known answer. Once that runs ok, add more complicated signals ( but still with a known answer). Only after this is all working, proceed to unknown signals and answers
If you still have a problem, please post why you think this is. What is the answer you expect and why? What is the answer that you are getting and how is different from your expectation? Why does the problem doesn't seem to appear on the basic debugging signals you have used in step 5.


You can unify the definitions for baseband and passband bandwidths by saying that bandwidth is only measured at positive frequencies. So for low pass signals, the lower frequency is zero and the upper frequency is the cut-off frequency. For bandpass signals, you have the difference between upper and lower cut-off frequencies (at positive frequencies).

For complex signals, a reasonable definition of bandwidth should take into account that the negative frequencies are not redundant anymore. So a complex low pass signal extending from $f=-B$ to $f=B$ should have bandwidth $2B$.

Concerning the choice of sampling frequency for bandpass signals, take a look at this answer. It is not true that a sampling frequency equal to twice the bandwidth is always sufficient. Yet you can usually do better than just sampling with twice the upper cut-off frequency.


What is the relationship between a digital audio signal and sound pressure level?

As you say, pressure will be proportional to output level, if the output level is equivalent to voltage. But sound pressure level (SPL) is by convention expressed in the logarithmic decibel (dB) scale, with 0 dB referenced to some standard sound pressure level. The relationship is:

$$\text{SPL} = 20\log_{10}\left(\frac{A}{A_\text{ref}}\right) \text{ dB,}$$

with a reference output level $A_\text{ref}$ that would give an SPL reading of 0 dB. For measures of power such as $A^2$ the formula changes slightly:

$$\text{SPL} = 10\log_{10}\left(\frac{A^2}{A^2_\text{ref}}\right) \text{ dB.}$$

With $A$ and $A_\text{ref}$ non-negative, the equations are equivalent with the square inside the logarithm or outside it as a factor $2$.


Need a similarity metric that describes these two curve as highly similar

While I still don't yet understand all the options in the dtw-python package, it seems as though the relative ordering of my two examples is now as expected using the code snippet

alignment = dtw.dtw(y, yref, keep_internals=True)

This allows me to create the following alignment figure using that package.

Additionally, the relative ordering of my two examples is now as I would hope. The more similar example has a distance metric of 6 whereas the more dissimilar example has a distance metric of ~11. I will have to explore the internals of this package more to see if I could be doing something even better, but this looks promising. Curious what the differences are in implementation that make this ordering so much different from the other two packages used.


Are there algorithms to invert a magnitude spectrum?

Apart from the stochastic literature, there are a set of deterministic algorithms that claim exact(?) recovery of a signal form its partial Fourier description, collectively referred to as:

  signal reconstruction from DTFT phase or magnitude alone

They are generally iterative in nature and exact convergence is hard to achieve in practice. Monson Hayes [et al] published relevant IEEE papers to propose several algorithms that reconstruct a real signal from its DTFT phase function alone, which is based on the phase sufficiency claim. The magnitude is a little more diffcult as I know. Have a look at those papers for your benefit.


Defining an average in time for a density in time

I think it is implied that $\int_{-\infty}^{\infty}|s(t)|^2 dt = 1$ (which is why the author refers to $|s(t)|^2$ as a density).

I agree with your sentiment "I've always thought about an average as dividing by some base", but in this case, the base would be 1, so no division is needed. If you were to plot $|s(t)|^2$ as a function of $t$, and think of this plot as a histogram, the the integral:

$$\int_{-\infty}^{\infty}t|s(t)|^2 dt $$
will tell you the average value, since each value of $t$ is weighted by the density function.


Given a matrix A that consists of zeros and ones only (which you should have after step 4 of your algorithm, you could implement step 5 via the following:

% find column with the maximum number of ones. sum down the rows of the
% matrix and find the column with the largest sum
[ry,xo] = max(sum(A));
% find the corresponding row; do the same thing except sum across the
% matrix's columns
[rx,yo] = max(sum(A,2));
% the values "rx" and "ry" above contain the sums of the row and column
% (respectively) that have the maximum number of ones. divide them by
% 2 according to your algorithm in step 5
rx = rx / 2;
ry = ry / 2;

After the above (which is untested), then xo, yo, rx, and ry should be populated with the values that your algorithm indicated.


Conditions for expressing a system as difference equation

A difference equation does not imply causality.

For example, consider this system, with input $x$ and output $y$:

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

You can solve it (recursively) forwards in time as $y[n] = x[n] + 0.5 x[n-1] - 0.2 y[n-1]$.

Or you can solve it backwards in time as $y[n-1] = 5 x[n] + 2.5 x[n-1] - 5 y[n]$.

These are 2 valid and different results for the same difference equation. These 2 results correspond to the 2 possible regions of convergence of the system response (Z transform). Inner ROC corresponds to the causal solution; outer ROC to the non-causal solution.

Systems of greater order will have more ROCs and more different solutions.

If you impose causality, then you are implicitly selecting the inner ROC.

For the second question, IIR filters can be causal. The example above is an IIR system. The causal impulse response will be $h_{causal}[n] = (-0.2)^n u[n] + 0.5 (-0.2)^{n-1} u[n-1]$. This is causal and IIR.


Increase SNR by using only the first 200 samples

If the SNR is limited by reflections in the room, and the impulse response of the channel from the source to a specific microphone is unknown, then the received signal can be used to determine the impulse response of the channel (together with a known transmitter signal that is spectrally rich in content) and with that equalize the received signal. However 200 samples would likely be far too insufficient for this purpose. In channel equalization we are basically solving for inverse convolution a system of overdetermined equations such that each output (sample) represents a sum of prior samples with unknown weights (which we solve for). These prior samples must extend over the significant duration of the channel impulse response, and then we need several samples to determine a least squared solution: each new sample provides us a new equation, and we would like to have more equations than unknowns as the accuracy of the solution will be improved the more equations we can use -- up to the limits of assumed stationarity in the channel statistics (the channel will change with time, so up to a point further inputs will contribute no further improvement, and then it will actually degrade the result).
200 samples at 48 KSps is only 4.16 ms. I believe typical channel impulse responses for a small room can be in the 100's of ms but I am less familiar with audio engineering to have confidence in that. But to equalize the channel to improve SNR by eliminating reflections/reverb I would recommend determining the impulse response time of the channel and then using at least 5 times more samples than that for the channel estimation with a known spectrally rich signal (white noise such as a maximum length sequence or frequency chirp). The equalization filter itself once the coefficients are determined would be the length of the channel response.
Once all channels are equalized for each microphone, assuming the remaining background noise is white and uncorrelated, further SNR improvement can be gained by coherently summing the microphone outputs: align their delay and sum each signal weighted by the SNR for that microphone for optimal ratio combining. If all microphones have the same SNR, this will provide a $\sqrt{N}$ SNR improvement where $N$ is the number of microphones.
If SNR is not limited by reflections in the room or other coherent sources, then any number of samples from multiple microphones can be added to get the $\sqrt{N}$ SNR gain assuming the coherent signals are aligned and the noise is uncorrelated on each microphone.


Is this a homework problem?

I figured out a way to solve this, but it isn't that simple.

First, you can split your function into three parts:

$$\begin{align}
x_1[n] &= A_1 r_1^n \cos(\omega_1 n) \\
\\
x_2[n] &= A_2 r_2^n \cos(\omega_2 n) \\
\\
x_3[n] &= A_3 \sin(\omega_3 n) \\
\end{align}$$

The $u[n]$ doesn't matter because we are only concerned with $n \ge 0$ and $x_3[0]=0$.

The cosines can be solved like this:

$$\begin{align}
x[n] &= r^n \cos(\omega n) \\
\\
x[n-1] &= r^{n-1} \cos(\omega n-\omega) \\
&= r^{n-1} \big(\cos(\omega n)\cos(\omega) + \sin(\omega n)\sin(\omega) \big) \\
\\
x[n+1] &= r^{n+1} \cos(\omega n+\omega ) \\
&= r^{n+1} \big( \cos(\omega n)\cos(\omega)  - \sin(\omega n)\sin(\omega) \big)\\
\\
r x[n-1] + \tfrac{1}{r} x[n+1] &= 2 r^n \cos(\omega n)\cos(\omega) \\
&= 2 \cos(\omega ) x[n]\\
\\
x[n+1] &= 2r \cos(\omega ) x[n] - r^2 x[n-1]\\
\\
x[i] &= 2r \cos(\omega ) x[i-1] - r^2 x[i-2] \\
\end{align}$$

The sine can be solved like this:

$$\begin{align}
x[n] &= r^n \sin(\omega n) \\
\\
x[n-1] &= r^{n-1} \sin(\omega n-\omega)  \\
&= r^{n-1} \big( \sin(\omega n)\cos(\omega)  - \cos(\omega n)\sin(\omega ) \big) \\
\\
x[n+1] &= r^{n+1} \sin(\omega n+\omega)  \\
&= r^{n+1} \big( \sin(\omega n)\cos(\omega)  + \cos(\omega n)\sin(\omega) \big) \\
\\
r x[n-1] + \tfrac{1}{r} x[n+1] &= 2 r^n \sin(\omega n)\cos(\omega) = 2 \cos(\omega ) x[n] \\
\\
x[n+1] &= 2r \cos(\omega) x[n] - r^2 x[n-1] \\
\\
x[i] &= 2r \cos(\omega) x[i-1] - r^2 x[i-2] \\
\end{align}$$

Notice, whether sine or cosine, the LPC equations are the same, what matters are the initial conditions.  Your parts are multiples of these, but that is also taken care of by proper selection of initial conditions.

Here is some Python code that demonstrates that these equations work:

(Edit:  The code is now in the followup below.)

What remains now is to combine these three parts into one equation.  This is where it gets messy.

What you have is your solution in the form of:

$$\begin{align}
x_1[n] &= a_{11} x_1[n-1] + a_{21} x_1[n-2] \\
x_2[n] &= a_{12} x_2[n-1] + a_{22} x_2[n-2] \\
x_3[n] &= a_{13} x_3[n-1] + a_{23} x_3[n-2] \\
\end{align}$$

You need to convert them to a common set of coefficients ($v_k$'s) in this form:

$$\begin{align}
x_1[n] &= v_1 x_1[n-1] + v_2 x_1[n-2] + v_3 x_1[n-3] \\
&+ v_4 x_1[n-4] + v_5 x1[n-5] + v_6 x_1[n-6] \\
\\
x_2[n] &= v_1 x_2[n-1] + v_2 x_2[n-2] + v_3 x_2[n-3] \\
&+ v_4 x_2[n-4] + v_5 x2[n-5] + v_6 x_2[n-6] \\
\\
x_3[n] &= v_1 x_3[n-1] + v_2 x_3[n-2] + v_3 x_3[n-3] \\
&+ v_4 x_3[n-4] + v_5 x3[n-5] + v_6 x_3[n-6]
\end{align}$$

That way they can be added together:

$$\begin{align}
x[n] &= x_1[n] + x_2[n] + x_3[n] \\
\\
x[n] &= v_1 x[n-1] + v_2 x[n-2] + v_3 x[n-3] \\
&+ v_4 x[n-4] + v_5 x[n-5] + v_6 x[n-6] \\
\end{align}$$

The way I can think of right now is a lot of work, but there may be another slick way to simplify this.  With a bunch of algebra, the following linear algebra problem can be set up:

$$
\left[
\begin{array}{c}
a_{11} \\
a_{21} \\
a_{12} \\
a_{22} \\
a_{13} \\
a_{23} \\
\end{array}
\right]
=
W
\cdot
\left[
\begin{array}{c}
v_1 \\
v_2 \\
v_3 \\
v_4 \\
v_5 \\
v_6
\end{array}
\right]
$$

Where $W$ is a 6×6 array and its elements are functions of the $a$'s.  The $v$'s can then be found my multiplying the left array by $W^{-1}$.

I haven't done the algebra, but there is a pattern there.

Hope this helps, phew.

Ced

Followup:

So there is a much easier way to figure out the $v$ values.  I didn't figure it out myself, I asked on the Math Exchange.  (see https://math.stackexchange.com/questions/2755034).

You simply have to multiply out:

$$ (S^2 - a_{11} S - a_{21})(S^2 - a_{12} S - a_{22})(S^2 - a_{13} S - a_{23}) \\  \qquad \qquad \qquad \qquad = S^6 - v_1 S^5 - v_2 S^4 - v_3 S^3 - v_4 S^2 - v_5 S^1 - v_6 $$

Messy, but it works.  Here it is in Python:

import numpy as np

#==========================================================
def MultiplyTuple( p, q ):

    Np = len( p )
    Nq = len( q )

    r = np.zeros( Np + Nq - 1 )

    for ip in range( Np ):
        for iq in range( Nq ):
            r[ip+iq] += p[ip] * q[iq]

    return r

#==========================================================

x1 = np.zeros( 100 )
x2 = np.zeros( 100 )
x3 = np.zeros( 100 )
x  = np.zeros( 100 )
y  = np.zeros( 100 )

w1 = .1
w2 = .18
w3 = .05

r1 = .9
r2 = .5

A1 = 1.0
A2 = 0.5
A3 = 0.4

a11 = 2.0 * r1 * np.cos( w1 )
a21 = - r1 * r1
a12 = 2.0 * r2 * np.cos( w2 )
a22 = - r2 * r2
a13 = 2.0 * np.cos( w3 )
a23 = - 1

x1[0] = A1
x1[1] = A1 * a11 / 2.0
x2[0] = A2
x2[1] = A2 * a12 / 2.0
x3[0] = 0
x3[1] = A3 * np.sin( w3 )

print 0, x1[0], x2[0], x3[0]
print 1, x1[1], x2[1], x3[1]

for i in range( 2, 20 ):
    x1[i] = a11 * x1[i-1] + a21 * x1[i-2]
    x2[i] = a12 * x2[i-1] + a22 * x2[i-2]
    x3[i] = a13 * x3[i-1] + a23 * x3[i-2]
    print i, x1[i], x2[i], x3[i]

print

for i in range( 0, 20 ):
    y1 = A1 * r1**i * np.cos( w1 * i )
    y2 = A2 * r2**i * np.cos( w2 * i )
    y3 = A3 * np.sin( w3 * i )
    y[i] = y1 + y2 + y3
    print i, y1, y2, y3

p1 = np.array( [ 1, -a11, -a21 ] )
p2 = np.array( [ 1, -a12, -a22 ] )
p3 = np.array( [ 1, -a13, -a23 ] )

p12  = MultiplyTuple(  p1, p2 )
p123 = MultiplyTuple( p12, p3 )

print p123

v1 = -p123[1]
v2 = -p123[2]
v3 = -p123[3]
v4 = -p123[4]
v5 = -p123[5]
v6 = -p123[6]

print

for i in range( 0, 6 ):
    x[i] = x1[i] + x2[i] + x3[i]
    print i, x[i], y[i]

for i in range( 6, 20 ):
    x[i] = v1*x[i-1] + v2*x[i-2] + v3*x[i-3] \
         + v4*x[i-4] + v5*x[i-5] + v6*x[i-6]
    print i, x[i], y[i]


There's an unlimited amount of projects you can do... The extend depends on how much you want to (or you can) turn your home into a lab ?
Audio projects can range from software synthesizers, multimedia players, to aduio compression algorithms, micophone array based advanced imaging, and to all sorts of acoustic magic. You would need alot of electronics and acoustic skills, in addition to DSP, for these projects though.
Hardware projects will also cost extra bucks & mess, compared to pure software projects...
Also, eventhough you've got your masters degree, just still need to polish your DSP skills, and make sure that you don't miss any serious topics in the fundamentals. The best way to do is to code the concepts in a book such as Discrete-Time Signal Processing. Implement (or verify) various algorithms and theorems there. Pay attention to convolutions, LCCDE recursions, filter design, sampling, DFT, and FFT topics. Write your own functions, simulations and libraries. Make theoretical predictions, and try verifying them with simulations... Use Matlab or Octave at this stage.
Have you done some Adaptive Filtering? Why not implement them in realistic applications? Try channel equalisation, echo cancellations, noise reduction, system identification, line enhancement, prediction, compression etc. Try LMS, RLS or their variants.
In the low-level software side, try improving your C / C++ and Assembler programming-debugging skills. Chose platforms such as Intel/AMD or ARM-based CPUs. I bet you will benefit from learning how to program MMX, SSE, AVX SIMD-vector extensions of these CPUs. And get a general glimpse on compilers, switches, libraries, open source development, various tricks and hardware properties and limitations of PCs.
Learn some operating system programming too. Windows,Linux,OS, and lately but more attractively Mobile Android. Learn how to use basic compilers for these platforms and their working mechanisms.
Start with the simplest and cleanest hardware: Arduino. Learn how to code it in its IDE. That's the most user friendly mCU ever. Though, it's a limited one, but you can still do a lot with it. It has an ADC-DAC; you may wish to try sampling some sources and processing them. Best is to go with a simple microphone shield (an attachement board to an Arduino). Record audio, process it, and output some signals based on it. You will learn a lot about USB connection, Serial (RS-232) communication, SPI, I2C communication protocols, motor drive examples, various sensors involving 9DOF acceleromoter-gyro-magnetic sensors, distance measurement sensors, light measuring sensors, temperature, pressure, gas sensors... You will have a lot to do with them... You may also learn practical Kalman filtering for perfecting your measurement results too.
You should support your Arduino with some basic external hardware, such as buffers, filters and actuators. It's good to learn connecting/interfacing them to the mcu, for practical embedded DSP.
Expand your hardware with a Raspberry PI, the most friendly electrically interactable PC platform with ARM based multi-core CPU and a Linux based OS in it.
Add some simple programmable SDR-RTL (software defined radio) kits to your Raspberry Pi to turn it into an RF receiver and data processing machine. You can receive signals from radio stations, Satellites, AirPlanes, wireless phones or any RF source in its range! You can write C-code to implement your analog/digital decoding algorithm. That makes a set of plain advanced DSP projects.
You can also progress with some more industrial MCUs from ST, NXP, Texas, ADI or any similar vendor. They are producing ARM complient cores. HEnce once you learn ARM architecture you will be confident in all of these platforms.
One good thing is, Arduino IDE also supports an STM32 based development card just as it supports ative Arduino UNO cards. So You can transit into ST programming from Arduino IDE and go more professional with ST native development later.
You can further add ADI BlackFin, SHArC, or Texas TMS based DSP processors.
Furthermore, last decade have withnessed great progress in Big-Data processing using GPU computing; so you will benefit form learning some parallel programming using Nvidia CUDA or openCL...
All these will take some time and patience.


What is the difference between the Weiner-Hopf equation and the normal equation?

Weiner-Hopf equation leads to Wiener filter that is optimal filter. For the case of stationarity in some time span it's the only filter minimizing MSE at its output. LMS and RLS algorithms are the adaptive approaches and they converge to Wiener optimal solution (as you can see from their convegence curves). While LMS is the simple steepest descent method of finding such an optimum, RLS solves Least Squares estimation problem at every step recursively.

LS estimation problem can be treated as a practical approach of Wiener filter. Look:

$w = \hat{R}^{-1} \cdot \hat{r}$

is the normal equation in matrix form for the LS problem, $\hat{R}$ and $\hat{r}$ are estimations of autocorrelation matrix and cross-correlation vector respectively. So if these estimations is properly done this solution leads us to the Wiener-Hopf eq., i.e. to Wiener filter. But it's only true if estimations are valid. So in general LS solution is not a Wiener solution, it's an approximation.

RLS approach is derived from LS (as it can be viewed from abbreviation):

$w(k) = S(k) \cdot \hat{r}(k)$,

where $S(k)$ is recursive estimation of $R^{-1}$ and $\hat{r}(k)$ is a recursive estimation of $r$. From the equation above one can see RLS try to solve LS problem at every step.

So the conclusion might be both LMS and RLS try to converge to Wiener solution and Wiener filter theory is in the core of these algorithms. But ways they do this are different.

This is my own understanding, maybe some discrepancies are possible. There is very good reference about adaptive filters: http://books.google.ru/books/about/Adaptive_Filter_Theory.html?id=MdDi_PF7gMsC&redir_esc=y


The term Doppler Shift is actually a bit of a misnomer. The frequencies are not actually shifted but they are scaled (see http://fourier.eng.hmc.edu/e101/lectures/handout3/node2.html for definition of shifting vs. scaling). It's a relative change not an absolute one.

Both time and frequency domains are scaled: when the source is moving towards you, the time axis compresses and the frequency axis expands. If the source is moving a away from you the time axis expands and the frequency axis compresses.

Let's look at a simple example. A car is 100m away from you and moves at a speed of 20m/s towards you. We assume speed of sound to 100 m/s (it's actually 344 m/s but I'm too lazy to do the math that way). Let's now assume that the car honks its horn for 1 second at a frequency of 1 kHz. So the honk consists of 1000 wiggles of air molecules. Let's call the start of the honk time t = 0.

At the beginning of the honk, the car is 100m away so the beginning of the honk arrives at your ear at time t = 1s. It takes one second for the sound to travel 100m. By the end of the honk the car has traveled 20m so it's only 80m away. The end of the honk arrives at time t = 1.8s. The sound only takes 0.8s to travel 80m.

The honk at the listening location is only 0.8s long, so the time has been compressed. It's still 1000 wiggles but they happen in less time. So the frequency is now 1000 wiggles/0.8s = 1250 Hz. The frequency is higher so the frequency axis has expanded.

It's important to note the frequency doesn't go up by a constant shift frequency but by a multiplicative factor. If the honk would have been 2 kHz the resulting Doppler frequency would be 2500 Hz and NOT 2250 Hz.


How to separate Even and odd part (in-phase bit streams I and quadrature-bit streams Q) of the signal?

The OP is mis-understanding how this QPSK signaling scheme operates.
There is a stream of input bits $b_0, b_1, b_2, b_3, \ldots,$ that is coming serially at the rate of one bit every $T$ seconds, so $b_0$ during $[0,T)$, $b_1$ during $[T,2T)$, and more generally, $b_n$ during $[nT, (n+1)T)$.. So, there is a clock signal at $\frac  1T$ Hz at the input.
The first thing that happens is the serial-to-parallel conversion where the serial bit stream gets converted into a stream of dibits (two bits in parallel), one bit on each of two individual wires and each bit lasting for $2T$ seconds (that is, the dibit lasts for $2T$ seconds.  Thus, we have
$$(b_0, b_1) \to \left[\begin{matrix}b_0\\b_1\end{matrix}\right], \quad (b_2, b_3) \to \left[\begin{matrix}b_2\\b_3\end{matrix}\right], \cdots \cdots
(b_{2n}, b_{2n+1}) \to \left[\begin{matrix}b_{2n}\\b_{2n+1}\end{matrix}\right], \cdots$$
So immediately we need to understand that

The output of the serial-to-parallel converter is clocked at with a clock at $\frac{1}{2T}$ Hz, and this clock is synchronized with the input clock of the serial to parallel converter,
There is a delay of $2T$ seconds in that $\left[\begin{matrix}b_0\\b_1\end{matrix}\right]$ is the output of the serial-to-parallel converter during $[2T,4T)$,  $\left[\begin{matrix}b_2\\b_3\end{matrix}\right]$ is the output during $[4T,6T)$, and so on,
The "even bits" have been sorted into appearing on the upper wire and the "odd bits" on the lower of the two wires, just as is envisioned in the OP's diagram. These are the Inphase and Quadrature branches of the transmitter.

The dibit enters into the direct-sequence modulator which says multiplication but this is if the bits ($0$ or $1$) have been converted into real numbers $+1$ and $-1$ respectively, and ditto the PN sequence. The logical operation with $0$s and $1$ is the XOR operation: the bit $b$ lasting $2T$ seconds is replaced by a serial stream of $N$ bits $(b\oplus \text{PN}_0, b\oplus \text{PN}_1, \cdots, b\oplus \text{PN}_{N-1})$ which last for a total of $2T$ seconds.  Thus the PN sequence generators (usually different ones in the I and Q branches) are clocked at $\frac{N}{2T}$ Hz, and this clock is also synchronous with all the other clocks.
Everything is digital up to this point. To get the actual transmitted signal, the $0$s and $1$s are mapped to $\pm 1$, each rectangular pulse (of duration $2T/N$ seconds in the direct-sequence QPSK system) is replaced by the desired baseband pulse shape (e.g. root raised cosine) to get the complex baseband representation $I(t)+jQ(t)$ of the QPSK signal. The actual RF signals in the I and Q branches
of the transmitter are $I(t)\cos(2\pi f_c t)$ and $-Q(t)\sin(2\pi f_c t)$ where $f_c \gg \frac{1}{2NT}$ is the RF carrier frequency. Thus the actual RF QPSK
signal is of the form
$$s(t) = I(t)\cos(2\pi f_c t)-Q(t)\sin(2\pi f_c t) = \text{Re}\big((I(t)+jQ(t))\exp(j2\pi f_ct)\big)$$ and results from modulating a (spread) dibit onto an RF carrier. Dibits consist of two bits; all Americans know that make two bits make a quarter, which explains why they think that the Q in QPSK stands for quaternary instead of quadrature or quadriphase.....


When does a Window function's sidelobe level become out of band attenuation?

Windows with low sidelobe levels have large transition
bandwidth's. Transition bandwidth is inversely proportional to N for a
given window.  Indeed, the ratio of transition width over sidelobe level
in dB for fixed length N is virtually a constant.

We can use the following empirical formula to estimate the required filter length

where,

N is filter length,
A is the out of band attenuation in dB
fs is the sampling frequency


How to plot histogram of difference of two images using MATLAB

In your second line, casting with uint8() will change a negative value to 0. Avoid casting or change to another type that supports negative values.


How to automatically detect blurred document images

Yes, but doing so is difficult. You are basically assuming that when a certain feature is sufficiently present in an image, you can apply a filter and realize an improvement in OCR accuracy. This implies that you can find the right feature, apply the right filter, and evaluate OCR accuracy correctly. Since there are lots of features that might be relevant (edge prominence, high-frequency energy, PSF estimation, unfiltered image OCR accuracy metrics (as @Matt mentions), etc.), lots of filters available (blind deconvolution, among others), and you'll need lots of ground truth images to prove system effectiveness, you are looking at a lot of work and experimentation to create a system that works well.

You also need to make sure you are not working against the OCR software itself, which may do its own image processing to improve its accuracy. In other words, make sure that artifacts introduced by your deblurring do not make OCR accuracy worse, even if the deblurring was done "well" from a deblurring perspective.

All that being said, it should be possible to do this. You might even be able to do this adaptively, by automatically and iteratively adjusting parameters and approaches to optimize the result.


Why the Nyquist frequency is 0.5 of Fs, why not 0.55 or 0.65?, brief explanation

From the definition of sampling, and the definition of periodicity follows that any signal is identical to signal at a frequency higher by an integer amount of sampling frequencies.
Since for real signals, negative frequencies aren't distinguishable from positive ones, it directly follows that for every $0.5f_\text{sample}\le f \le f_\text{sample}$, there is a negative frequency $-0.5f_\text{sample}\ge \tilde f=1- f \ge f_\text{sample}$ that looks the same as a positive frequency.
Hence, only frequencies lower than $|f| < 0.5f_\text{sample}$ can be uniquely represented.


Impulse response of system with $\max\left\{\lvert x[n] \rvert\right\}\ge\max\left\{\lvert y[n]\rvert \right\}$

Not sure what exactly the question is: find some LTI system that meets this constraint (which is trivial) or create a testable condition for any LTI system. (Which is not trivial).
For the latter: Discrete LTI systems are guaranteed to not increase the maximum amplitude if the absolute sum of the impulse response is smaller than $1$, i.e.
$$\sum |h(n)| \leq 1 $$
The proof is a bit tedious, so I'll skip it unless there is need for it.


is it ok to showcase this diagram ?

Sorry, we can't tell you this. This is really a function of your specific application, the physics behind your signal, and how you exactly define "peak".
DSP can offer you methods to suppress "near by" peaks or to de-noise multiple peaks but whether that's the right thing to do or not depends largely on what do want to do with the results and what else you already know about the signal.


Estimation of Two Closely Spaced Frequencies?

This answer is a consideration of the frequency resolution when using the DFT but does not answer the question specifically in terms of what is the best algorithm. (see comments by Amro).

If the two closely spaced signals are similar in magnitude (and sufficiently stronger than any other signals, i.e. high SNR), then the highest frequency resolution (which is what you are interested in) is achieved with a rectangular window. This means simply take an FFT with no windowing, and the signal bandwidth will be 1 bin wide or similarly equal to 1/T where T is the time length of your composite signal. This relationship holds regardless of what algorithm is used; the frequency resolution is related to the time length of the signal with the best resolution of 1/T. If the two closely spaced signals are not similar in magnitude, you will run into dynamic range issues as the rectangular window has the worst sidelobe levels, which can mask lower level signals. In this case we use windowing; specifically we multiply the time series by a window function prior to taking the FFT, which will significantly reduce sidelobe levels, but will have the drawback of increasing the width of the main lobe (and hence decrease frequency resolution).

See this answer for further related info: Specific frequency resolution


No cross-correlation between same signal recorded using two different microphones

This answer is now obsolete.

I still keep it, because future searchers might come here with similar problems.

This is not the signal as it is now, but as it was before, when I answered:

I'll do a hypothesis: the autocorrelation peaks are due to periodicity in your signal. But if that's the case, they should be the same for $R_{11}$ and $R_{22}$.

My guess is that you're comparing two recordings of different sampling rate (for example, 44.1 vs 48 kHz). But xcorr only sees the vector of numbers and doesn't know about the sampling rate!


When does the convolution of $2$ signals equal zero?

Time-domain convolution is frequency-domain multiplication. If at all frequencies at least one of the signals is zero-valued in frequency domain, then the convolution of the two signals will be zero-valued at all frequencies, and at all times. Except for a zero signal, no finite-length signal has a continuous run of frequency domain zeros, so your choices are limited to infinite-length signals such as periodic signals with no coinciding non-zero harmonics and ideal lowpass–highpass filter impulse response pairs.


"plot" command for continous time signal in matlab?uses interpolation?

By default MATLAB uses linear interpolation when creating line plots, which means it simply draws a line from each point to each point, unless there are more points than pixels in which case each point or group of points within each pixel would represent each pixel shown.

Compare the following to see this:

This plots a line from x,y=1,2 to x,y= 5,7

plot([1 5], [2 7])

The above is the same as the following where a line is explicitly declared:

plot([1 5], [2 7], '-')

While in comparison the following will plot a “o” at each sample:

plot([1 5], [2 7], 'o')

Update: To interpolate more smoothly between samples such as polynomial interpolation, you would need to create the interpolated samples and then plot the higher sampled result. The simplest way to do this in MATLAB is to use the resample command where resample(f,M,N) will interpolate the waveform $f$ by $M/N$ samples. Another approach demonstrated below is interpolating using the FFT. Although DSP.SE is not a MATLAB coding site, the following does demonstrate well how to properly do a time domain interpolation using a zero-padded IFFT which is of general signal processing interest. Zero padding in one domain results in interpolation in the other domain, however the technique is to "zero-pad" from the center of the frequency domain waveform with index k from 0 to N-1, such that the high frequency component of the FFT is extended.

N=6;
M=1000;
x= [0:N-1];
y = sin(2*pi*x/N);
fy = fft(y);
padfy = [fy(1:floor(N/2)) zeros(1, M-N) fy(floor(N/2+1:end))];
figure;
plot(N*[1:M]/M, M/N*real(ifft(padfy)));


Is the first derivative operation on a signal a causal system?

If the derivative exists at the given point, then it doesn't matter if you look (infinitesimally) into the future or into the past, you can do both, because both will give the same result:

$$x'(t)=\lim_{h\to 0}\frac{x(t+h)-x(t)}{h}=\lim_{h\to 0}\frac{x(t)-x(t-h)}{h}\tag{1}$$

So a differentiator can be (theoretically) implemented by a causal system.

Clearly, the system is not stable, because a bounded input can cause an unbounded output. E.g., a step at the input causes an unbounded output signal. So if by realizable we understand causal and stable, the ideal differentiator is not realizable. However, good approximations in a certain frequency range can be realized with small errors, in the continuous domain as well as in the discrete domain.


What is a good book reference on the definition of "subharmonic frequency"?

"Highest available frequency" sounds fishy but here's a reference from dynamics:

Cunningham, Walter Jack. Introduction to Nonlinear Analysis. 1st ed. New York, USA: McGraw-Hill, 1958, pp. 191–192.

ISBN 10: 0070149178  ISBN 13: 9780070149175


After my trials, I realize that problem is directly about Warning: Image is too big to fit on screen; displaying at 67% warning. I use imwrite to my matrix and convert it to jpg file. After looking the image using jpg file , there is no problem.


What is the difference between Lumped and Distributed systems?

The elements building a lumped system are thought of being concentrated at singular points in space. The classical example is an electrical circuit with passive elements like resistor, inductance and capacitor. The physical quantities current and voltage are functions of time (only). E. g. the current at a capacitor with capacity $C$ is given by
$$
i(t) = C\frac{\mathrm d v(t)}{\mathrm d t}
$$
Where $C$ is a constant (and so are $R$ and $L$). This leads to ordinary differential equations.

In contrast, the elements in distributed systems are thought of being distributed in space, so that physical quantities depend on both time and space. The classical example is the electrical line where inductance, capacity and resistance are not constant but functions of length $x$. This leads to partial derivatives of $i(t,x)$ and $v(t,x)$ in $t$ and $x$.


Given a complex time-domain signal is it always implied that this is a QAM-demodulated form of a real time domain signal?

No, we generally use complex baseband to describe all kinds of signals, or even passband systems (e.g. frequency-selective or time-variant channels). So, any bandlimited signal has an equivalent complex baseband, and it's not implied this has anything to do with QAM or even digital communication.

QAM is just one example of something that is really easy to write down in complex baseband.

  Are there any other commonly used ways of converting a real time domain signal to a complex time domain signal that don't imply QAM?

I don't even know how QAM comes into this – you transform any real-valued bandpass signal to complex baseband through multiplication with its (complex) carrier and low-pass filtering. QAM really has nothing to do with it.


Feature extraction from spectrum to aid categorisation

First thing to mention is that you must have some kind of Machine Learning algorithm to perform your recognition. You might want to use the following - what's yours?

Neural Networks (easy to implement)
SVM's (should perform best)
k-NN (eventually, but needs lot's of data to store)

Regarding features, it's hard to determine which to use for your application just by looking at spectra - both of categories have peaks at lower and higher part of range. Also there might be some significant information in time domain or a specific method valid for processing this type of signals. You never said what these signals are and what's the application.

The only thing I notice is that class 2 has some harmonic content occurring at frequencies $160 \mathtt{Hz}$, and $200 - 225 \mathtt{Hz}$. I don't think that simple mean and variance will do. Probably you should try different and choose best performing ones. My suggestions of features you've asked for (to start with), are:

Spectral Slope (what's the gradient of the linear regression of your spectrum)
Audio Spectrum Flatness (calculated in frequency banks - tells you if signal is noisy or harmonic)
Audio Spectrum Centroid (kind of 'mean value' you mentioned)
Audio Spectrum Envelope (can be understood as very general spectrum descriptor - envelope)
MFCC (well defined and described - they are describing your spectrum by using cepstral analysis)

For more info about implementations you can refer to this great book (easy to get): MPEG-7 Audio and Beyond.


This is the solution for a linear interpolation between two frequencies as a "weighted average". You would just need to decide what your weighting actually is in terms of magnitude or dB scale for the vertical axis.  X is the interpolated frequency, weighted toward the stronger signal, and Y is the linear interpolation at the new frequency location of the signal strenths a and b.

Here is a simple example demonstrating the results that would be achieved with the formulas given:


Frequency domain to time domain conversion-(varying sample count)

Consider the DTFT of an N-sequence $x_n$:

$$
\begin{align*}
X(e^{j\omega}) &= \sum_{n=0}^{N-1} x_n e^{-j\omega n} & \text{(DTFT)} \\
&= \sum_n \left ( \frac{1}{N} \sum_{k=0}^{N-1} X_k W_N^{-kn} \right ) e^{-j \omega n} & \text{(subs. IDFT)}
\end{align*}
$$

Here $X_k = X(e^{j\frac{2\pi}{N}k}), \, W_N = e^{-j\frac{2\pi}{N}}$, and let $\omega_k = \frac{2\pi}{N} k$, then rearrange:

$$
\begin{align*}
X(e^{j\omega}) = \sum_{k=0}^{N-1} X_k \left ( \frac{1}{N} \sum_{n=0}^{N-1} e^{j(\omega_k - \omega) n} \right )
\end{align*}
$$

Hence the DTFT of an N-sequence is uniquely determined by N frequency samples uniformly spaced between $0$ and $2\pi$. The term in the parentheses is an interpolation function (it's actually the digital sinc).

Check whether your $x_n$ is periodic to 150 samples.

What happens when you zero pad the DFT to length $P > N$ (add $P-N$ more $X_k = 0$ terms), and take the IDFT? You interpolate to a higher sampling rate (since $\omega_k = \frac{2\pi}{P}k$ is less than before).

Note that you're adding points after the first $N/2$ samples, explained here.


The fourier transform of the two dimensional autocorrelation function should do it.


Why FIR uses often direct-form instead of transposed-form?

I don't see how memory or computational burden can possibly be reduced between the two implementations (note this is not implying simplicity, just number of computations required, please read on...) Usually the decision to use transposed-form is for high speed FPGA implementations with a large number of taps as you can eliminate a long adder tree which would otherwise limit the maximum clock frequency to a lower number due to the accumulated effect of multiplier and adder delay, but this can get more complicated for a fixed precision implementation since the bit width must grow as you proceed in the filter. (This consideration is particularly important in IIR implementations where feedback exists but affects FIR implementations with large adder trees, and also a reason to factor longer IIR filters in cascade form). At lower speeds the direct-form is more convenient to implement as you can use a single extended precision accumulator.

I include a summary of the basic forms applicable to FIR/IIR below for benefit of readers that may be less familiar, and see link from Ben in the comments for more advanced structures specific to FPGA implementations:


Magnitude of Power Spectral Density with Different Sampling Frequency

The issue (beside the scaling) here is that your selected frequencies do not align with your FFT grid and that the energy of each complex exponential is smeared out over multiple FFT bins. Your frequency resolution is given by sample rate divide by FFT length. In your case that's 3.906 Hz. Only for frequencies that are an integer multiple of this will all energy fall into a single FFT bin. An easy way to achieve this is to bump the sample rate up to 1024 Hz. Then your frequency resolution is exactly 4 Hz and all your frequencies are an integer multiple of that. Once you fix the scaling as well, it works fine like this:

N=256;
x0=zeros(N,1);
fs=1024; % Fix 1: get to exactly 4Hz frequency resolution
idx = (0:(N-1))'; % Indices
f=20*[1 2 3 4];
A=[1 2 3 4];
K=4;
for ii=1:K
    x0=x0+A(ii)*exp(sqrt(-1)*2*pi*f(ii)*idx/fs);
end
X=fft(x0);

figure;
plot(abs(X)/N); % Fix 2: absolute value only, rescaled


Given a downscaled image, produce Gaussian blur of the original image

I don't think you can do that.
You can probably get a similar result, in that you won't be able to tell by eye that they are different.

The reason I don't think you can do what you have in mind is this:
Imagine that you downscale (downsample, basically throw away pixels) an image, apply a blur, upscale the image.  Suppose now that you apply a blur to the original image.  If those two are the same, then you could deconvolve both images and from both of them recover the original image.

This would be the best compression format ever.  Downscale by 16X, blur.  Small file, send by email.  Upscale by 16X, deblur, original image.  Hell, why stop at 16X?  You could compress all of your images to a single pixel each and be able to recover them.

Unfortunately, things don't work that way.  Downscaling will (must) throw away data, and you can't recover it.

So, the downscaled and blurred image will be different from the blurred original.