DTFT of sine wave using freqz As mentioned in the title, is it possible to use freqz to find the DTFT of a sine wave? I am confused about what the 'a' and 'b' vectors would look like, since there are only impulses in the numerator.

The MATLAB / OCTAVE function freqz(b,a) is typically used for displaying the magnitude and the phase of a Frequency Response function $H(e^{j\omega})$ associated with a discrete-time, linear time-invariant (LTI), system with rational transfer function $H(z)$ (having a corresponding LCCDE representation) given by
$$H(z) = \frac{ \sum\limits_{k=0}^{M} b_k z^{-k} }{\sum\limits_{k=0}^{N} a_k z^{-k} } \tag{1}$$
where the numerator and denominator coefficients $b_k$ and $a_k$ represent such systems.
It's not meant to display spectrum of signals. Of course, an impulse response $h[n]$ associated with an LTI system can also be considered as a signal, and therefore its Frequency Response function $H(e^{j\omega})$ can also be considered as a Fourier transform of a signal ($h[n]$ in this case) too.
However, for such an interpretation to make sense, the signal $h[n]$ should be the impulse response of an LTI system. Stated in other words, to display the spectrum of $h[n] = \cos(\omega_0 n ) $, you should find the coefficients $a_k$ and $b_k$ that constructs a LCCDE whose impulse response solution is $h[n]$.
Consider a first order LCCDE (with initial rest):
$$ y[n] - A y[n-1] = x[n] \tag{2}$$
taking its Z-transform
$$ Y(z)(1 - A z^{-1}) = X(z) \tag{3}$$
and the System Function is
$$ H(z) = \frac{Y(z)}{X(z)} = \frac{1 }{1 - A z^{-1}}  \tag{4}$$
and the impulse response is (assuming causal)
$$ h[n] = A^n u[n] \tag{5} $$
If we let $A = e^{j\omega_0} $ , then the impulse response becomes
$$ h[n] = e^{j \omega_0 n} u[n] \tag{6} $$
From Eq-6, we have a single complex-exponential, by the choice of the coefficient $A =e^{j\omega_0}$. When combined with Euler identity of a causal cosine wave, we can see that
$$ h[n] = \cos(\omega_0 n) u[n]= 0.5( e^{j \omega_0 n} + e^{-j \omega_0 n} ) u[n] \tag{7} $$
Which can also be written with $A$ as
$$ h[n] = \cos(\omega_0 n) = 0.5( A^n + (A^*)^n ) u[n] \tag{8} $$
where $A^*$ is the complex-conjugate of $A$. Using the relation between Eqs.4&5, we can find the Z-transform of $h[n]$ as
$$ H(z) = \frac{0.5}{1- A z^{-1}} + \frac{0.5}{1 - A^* z^{-1} } \\ \tag{9} $$
equivalently
$$ H(z) = \frac{1 - 0.5(A+A^*) z^{-1}}{ 1- (A + A^*) z^{-1} + |A|^2 z^{-2} } \\ \tag{10} $$
we can construct back the differential equation that corresponds to this $H(z)$ as
$$ y[n] - 2\mathcal{Re}\{A\} y[n-1] + |A|^2 y[n-2] = x[n] -\mathcal{Re}\{A\} x[n-1] \\ \tag{11}$$
Plugging $A = e^{j\omega_0 n}$ into Eq-11 yields
$$ y[n] - 2 \cos(\omega_0) y[n-1] + y[n-2] = x[n]-\cos(\omega_0) x[n-1] \tag{12}$$
from Eq-12, you find the corresponding coefficients $a_k$ and $b_k$ as:
$$ b_k = [1, -\cos(\omega_0)] \tag{13}$$
and
$$ a_k = [1, -2 \cos(\omega_0), 1] \tag{14} $$
Now, you can display the Fourier transform of $h[n] = \cos(\omega_0 n)$ by using the freqz(b,a) function, with the coefficeints in Eqs-13,14.
Note that, $h[n]$ is only a causal sine wave valid for $n>0$. It's not the ideal two sided sine wave. Furthermore, I hope you are aware that the ideal sine wave's Fourier transform is a dirac impulse function. And what you can view freqz(b,a) or any other function is actually the spectrum of teh windowed cosine wave of a finite length.
The following MATLAB/ OCTAVE code shows you the implementation of above:
% spectrum of h[n] = cos(w0*n) to be displayed by freqz(b,a) function

N = 128;               % number of samples to evaluate the impulse response

w0 = 3*(2*pi/N);       % frequency of the cosine h[n] = cos(w0*n)
b = [1,-cos(w0)];      % coefficients b[k] of the LCCDE
a = [1,-2*cos(w0),1];  % coefficients a[k] of the LCCDE

h = impz(b,a,N);       % N samples of h[n] in n=[0:N-1] evaluated

figure,stem([0:N-1],h);
title('imp-resp of y[n] - 2 cos(\omega_0) y[n-1] + y[n-2] = (1-cos(\omega_0)) x[n]');

K = 1024;              % number of frequency points to be evaluated.
figure,freqz(b,a,K);
title('Freq-Response of the system: y[n] - 2 cos(\omega_0) y[n-1] + y[n-2] = (1-cos(\omega_0)) x[n]');

The outputs look like :

and

Note that I have chosen a perfect frequency for the cosine wave, so that it displayed almost like an impulse. It will not be the case in general.


Dealing with Random Errors in Sampling Rates? I wonder what generalization I can make about what happens when I introduce some error in my sampler spacing.

For example, If I sample once every second I would want my sampler to take a signal at 1s,2s,3s,4s.

Unfortunately, my real world sampler looks like more like 1.01s, 1.99s, 2.99s, 4.01s.

I want to make a statement about the Fourier domain. I expect my high frequencies to be ambiguous while my lower ones should not to change, too much. Is this correct?

I understand how to do this in terms of a uniform change, for example 1.01, 2.02, 3.03 (which is just a change in sampling rate), but I'm wondering how I can describe this with some knowledge of my error distribution, or for my application even understanding how to describe this knowing the mean of the error would be wonderful.

Can somebody point me in the right direction?

You might read up on "Aperture Jitter". You are correct that higher frequencies are affected more than lower ones. Here is a good reference:

http://www.analog.com/static/imported-files/tutorials/MT-007.pdf


What is a "Hampel Filter" and how does it work? One of the new functions in the MATLAB Signal Processing toolbox for R2015b is the Hampel Filter.  It appears to be used for outlier removal and from the examples, it looks like it might do a better job than a median filter.

Unfortunately, as of yet there appears to be no readily available digestible description of the filter (such as a Wikipedia page).

So, what is a Hampel Filter?  How does it work?  How computationally expensive is it?  When does it work well? When does it not work so well?

The documentation clearly describes its function: "For each sample of x, the function computes the median of a window composed of the sample and its six surrounding samples, three per side. It also estimates the standard deviation of each sample about its window median using the median absolute deviation. If a sample differs from the median by more than three standard deviations, it is replaced with the median."

The $3\sigma$ rule is enforced via the estimator $\sigma \approx 1.4826\ \operatorname{MAD}$.

The most costly operation must be the median of seven elements, possibly achieved by sorting, which can be performed efficiently in a sliding window.

Possibly relevant reference: "Hampel F. R., ”The influence curve and its role in robust estimation,” Journal
of the American Statistical Association, 69, 382–393, 1974."


Is it possible to compute the wavelet transform of a unit step function? Is it possible to find theoretically the continuous or discrete wavelet transform of a unit step function (consntly 0 for t<0 and constantly 1 for t>=0) or of a function starting from zero and approaching 1 in a finite time? I just wonder if this is possible because these signals don't have a finite energy. Thanks. E.

Finite energy functions is a class of choice, but wavelet transforms can exists in more generic cases.

For the continuous case: let us call $H(t)$ the Heaviside step function. Then $|H(t)\psi(t)|$ is dominated by $|\psi(t)|$. If $\psi(t)$ is $L_1$ or integrable (in addition to being  $L_2$), then  $\int_{-\infty}^{\infty}H(t)\psi(\frac{t-a}{b})\mathrm{d}t$ is defined. Moreover, since admissible wavelets have zero-average, the continuous wavelet coefficients can decay at infinity. With wavelet of finite support, they will be zero for all $\psi(\frac{t-a}{b})$ with support in $[0,\;+\infty[$.

For discrete wavelets, this would be about the same. As long as the filters are summable (or $\ell_1$), and this happens  especially for FIR wavelet filters, you can compute coefficients. Special note: if a discrete filter is $\ell_1$, then it is $\ell_2$ as well, so the condition is milder than in the continuous case.


Modulation Index of Conventional AM Signal An AM signal has the form

$$u(t)=100\left[1+\sin(2000 \pi t)+5\cos(4000\pi t)\right]\cos(2\pi f_ct)$$

Where $f_c$ is 800 kHz.

What is the modulation index?

(We can define modulation index as the ratio between the amplitude of message wave and the amplitude of carrier wave)

Modulation index for linear (continuous wave) amplitude modulation may be defined in different ways. But the one I found relevant to your definition is:

$$ \mu = \frac{ | \min \{ m(t) \} | } {A_c} $$

where $A_c$ is the carrier amplitude, $m(t)$ is the message signal and $\mu$ is the modulation index.

For simple demodulation of a conventional AM modulated signal $s(t)$ we look for $\mu < 1$ so that the envelope of $s(t)$ never becomes less than zero.

$$s(t) = [ A_c + m(t) ] \cos(2 \pi f_0 t) $$

Your signal is:
$$ s(t) = [ 100 + 100 \left( ~~\sin(2000 \pi t)+5\cos(4000 \pi t)~~ \right) ] \cos(2 \pi f_c t) $$

Your message signal is:
$$ m(t) = 100 \left( ~~\sin(2000 \pi t)+5\cos(4000 \pi t)~~ \right) $$

whose min value is $-600$ and from which you deduce the modulation index as:

$$ \mu = \frac{ | \min \{ m(t) \} | } {A_c} = \frac{ | -600 | } {100}  = 6 $$

I'm sorry to say but this is severely overmodulated...


What should go to infinity in power's formula N or the repetition KN? I have this question :

I should now calculate power of y[n], where the formula

well in the power's formula N which constant goes to infinity how could this be ?
i see in the question that N=7 which is constant,then should N or KN go to infinity ?

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


Solution verification for this Fourier series problem We have a signal with period $T = 2$

We want to find the continuous time fourier series for this signal.
Since $T = 2$, $\omega = \pi$. All we have to do know is find the frequency domain.
$$x(t) = \sum_{k = -\infty}^{\infty}X[k]e^{-i \omega tk}$$
$$X[k] = \frac{1}{T}\int _{T}x(t)e^{i \pi tk}$$
$$X[k] = \frac{1}{2}\int _{2}x(t)e^{-i \pi tk} = \frac{1}{2} \left( \int _{0}^{1}\frac{3}{2}e^{-i \pi tk}dt  + \int _{1}^{2}\frac{-3}{2}e^{-i \pi tk}  dt \right)$$
$$X[k] = \frac{3}{4} \left( \int _{0}^{1}e^{-i \pi tk}dt  - \int _{1}^{2}e^{-i \pi tk}  dt \right) = \frac{3}{4} \left( \frac{-1}{i \pi k} \right)\left( e^{-i \pi k} - 1  - e^{-2i \pi k} + e^{-i \pi k}
\right)$$
$$X[k] = \frac{3}{4} \left( \frac{-1}{i \pi k} \right)\left( (-1)^k - 2\right) \text{(Because k can only be an integer)}$$
This to me is correct... however the answer is actually
$$X[k] = \frac{e^{-j\frac{k\pi}{2}}3\sin(\frac{k\pi}{2})}{k\pi}$$

You made a small mistake in the final step. The second term in parentheses should be $2((-1)^k-1)$ because you have two terms $e^{-j\pi k}$. So your final result should be
$$X[k]=\frac{3}{2\pi jk}\big(1-(-1)^k\big)\tag{1}$$
With $(-1)^k=e^{-j\pi k}$, Eq. $(1)$ can be rewritten as
$$\begin{align}X[k]&=\frac{3}{k\pi}\frac{1-e^{-jk\pi}}{2j}\\&=\frac{3e^{-jk\pi /2}}{k\pi}\frac{e^{jk\pi /2}-e^{-jk\pi /2}}{2j}\\&=\frac{3e^{-jk\pi /2}\sin(k\pi/2)}{k\pi}\tag{2}\end{align}$$


Link between DFS, DFT, DTFT My understanding of DFT is as follows

For a signal $x[n]$ of finite-length, the DFT is DFS of the periodic extension, $\tilde{x}[n]$, of that signal $x[n]$ and also another way to view DFT is that it’s a sampling of continuous DTFT.

Given that it is possible to reconstruct a original signal from sampled signal, provided the sampling is greater than Nyquist frequency.
We know that the DTFT for sampled signal is a series of replications of the spectrum of the original signal at frequencies spaced by the sampling frequency. Now, since DTFT is continuous and periodic, we can further breakdown DTFT at intervals and still be possible to reconstruct the DTFT and consequently the original signal.  This act of breaking down or sampling the DTFT is called DFT.

Is my interpretation of DFT correct? I would welcome any (true) facts or implications to test my understanding

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


Is the periodogram squared-magnitude DFT or squared-average DFT? How is the periodogram defined? I have been searching the internet a lot, and I find two different definitions:

periodogram is the squared magnitude of the DFT
periodogram is the squared average of the DFT

which one is correct?

The periodogram is simply the squared magnitude of the DFT. Since the periodogram is a rather poor estimate of the power density spectrum of a random process there are methods which use averaging of periodograms to obtain better estimates of the power spectrum. Two such methods are Welch's method and Bartlett's method.


Predicting account balance based on historical data Given the values of a bank account balance over time (see figure below as an example), how can one predict the account balance at a given date in the future ?

Should I just fit one linear regression model to the values of the i'th day of each month, in order to predict the value at the i'th date of future months ?

Given the strong correlation of the downward trend in each cycle, I offer the following that would offer an improved prediction over one linear regression, but not be significantly more complicated (but still would like to see what a precise "best estimator" would be):

Do a linear regression to determine the general upward trend. This would be the Income - Spend = Savings rate.
Do a least mean square linear fit for all the downward trend cycles. This would be the spend rate. (my temptation to do this easily would be detect the upward spike and superimpose the datasets, find the mean downward trend and do a linear regression to that.)
Determine the time of each upward spike to determine the best estimate of the frequency of occurrence of the upward spikes (clock recovery). For this the strongest tone in the FFT should reveal this.  This would be the pay dates.
Determine the mean of the magnitude of the upward spikes.  This would be the income per pay date.

From the information above you can create a future projection of continued upward spikes at the best estimated times along an overall upward trend that will provide a better estimate than one linear regression over the whole dataset.


Is There a Way to Perform a 2D Image Rotation by Matrix Multiplication? In my understanding, unlike shifting and scaling, image rotation actually rotates the coordinates whereby interpolating the pixels at discrete positions. It is like a projection that discrete values are projected to another plane which are along with the pixel value.

Thus, I am wondering whether we could do 2-D image rotation in a matrix multiplication way which directly perform rotation on the pixel matrix of the image?

Any help is appreciated in advance.

Let's say you have an $ M \times N $ image.
If you turn it into a vector in $ {R}^{\left( M N \right)} $ and create a matrix which is $ \left(M  N \right) \times \left( M  N \right) $ by multiplication and rearranging the vector from the multiplication result you can reposition any pixels where ever you want.

Being practical, by padding the original image you can create rotation by any angle you want even with interpolation (By summing few pixels into one position).

It all starting with vectorization of the image.


Computing impulse response for system with feedback terms I'm trying to compute the impulse response for the following system:

$$
y(n) = \frac{1}{2} \left( x(n) + x(n-1) \right) + \frac{1}{2} y(n-1)
$$

and am told to assume that $y(-1)=0$.

My solution: Letting $x(n) = \delta(n)$, we have

$$
h(n) = \frac{1}{2}\left(\delta(n) + \delta(n-1)\right) + \frac{1}{2}h(n-1).
$$

Plugging in a few values for $n$ (and assuming that $h(-1)=0$), I've obtained:

$$
h(0) = \frac{1}{2},
$$

$$
h(1) = \frac{1}{2} + \frac{1}{2} \cdot \frac{1}{2} = \frac{3}{4}
$$

and

$$
h(2) = \frac{1}{2} + \frac{1}{2} \cdot \frac{3}{4} = \frac{3}{8}.
$$

From this, I've observed that

$$
h(n) = \begin{cases}
0, & n<0 \\
\frac{1}{2}, & n=0 \\
\frac{1}{2} + \left(\frac{1}{2}\right)^n, & n > 0
\end{cases}
$$

Does this seem correct? I feel like I made a mistake somewhere and that the expression for $h(n)$ should be in a simpler form.

Note: I computed the impulse response by taking the Fourier Transform and then plugging in values for $n$ and got something different than this. However, I believe that this discrepancy is due to the initial condition $y(-1)=0$.

Your values for $h(0)$, $h(1)$ and $h(2)$ are correct. However, the general formula isn't.
It should be

$$h(n) = \begin{cases}
0, & n<0 \\
\frac{1}{2}, & n=0 \\
\frac{3}{2}\left(\frac{1}{2}\right)^n, & n > 0
\end{cases}$$

A simple way to solve such problems is to use the $\mathcal{Z}$-transform:

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

which gives

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

With the transform pair

$$\frac{1}{1-az^{-1}}\Longleftrightarrow a^nu(n),\quad |a|<1$$

and with the property

$$x(n-1)\Longleftrightarrow z^{-1}X(z)$$

you arrive at

$$h(n)=\frac{1}{2}\left(\frac{1}{2}\right)^nu(n)+\frac{1}{2}\left(\frac{1}{2}\right)^{n-1}u(n-1)$$

which corresponds to above result.


Multiplying the imaginary part of DFT with a linear ramp to get a derivative I am trying to understand the statement in a relatively old publication from 1970s, when Fourier transforms found applications in chemical analysis. The author quotes the derivative theorem citing Bracewell of 1965 (I could not locate this exact statement in later edition of Bracewell).
Edits: As discussed in the comments, the approach as literally suggested does not work universally for example a Gaussian peak will have zero imaginary parts. However the idea is neat to differentiate experimental data while filtering the noise.

I've changed my mind during the answer. Instead of using frequency domain processing to get the derivative signal, I will provide a more direct and clean approach to get the derivative of a signal by using what's called a digital-differentiator.  An LTI filter, which is applied into the samples $x[n]$ of a signal $x(t)$ and computes the samples $y[n]$ of the derivative $y(t) = x(t)'$ of $x(t)$.
Below is the matlab code you can use to compute the samples of the derivative of a bandlimited signal $x(t)$:
N = 512;                   % length of x[n]; samples of signal x(t).

% Assume Gaussian pulses added:
t = linspace(-1,1-2/N,N);  % time is sampled at T = 2/N period.

sigma = 0.001;             % Gaussian pulse shape parameter (adjust accordingly)
x1 = 2*exp(-(t.^2)/(2*sigma) );         % Gaussian pulse one
x2 = 0.5*exp(-((t-0.5).^2)/(3*sigma) ); % Gaussian pulse two
x = x1 + x2;               % total signal

% Generate differentiator impulse response hc[n]:
T = 2/N;                  % sampling period (implicit, defined above)
n = -N/2:N/2;
h = ((-1).^n)./(n*T);
h(N/2+1) = 0;
w = hamming(2*N/2+1)';
hc = w.*h;                % causal,stable impulse response hc[n]

% Filter x[n] to compute samples of derivative signal
y = conv(x,hc);
yd = y(N/2+1:N/2+N);     % extract relevant portion

% display the result:
figure,plot(t,x);
hold on
plot(t,0.1*yd);
legend('x(t)','dx/dt');
grid on
title('x(t) vs dx(t)/dt  (scale adjusted to fit both');

The result is :

For getting the same result using DFT techniques,
the following code can be used :
M = N+N;
H = fft(hc,M);          % generate the RAMP mask to apply on X(w)
Y = fft(x,M).*H;        % obtain ramp multiplication in freq domain
y = real(ifft(Y,M));    % invert back to time domain
y = y(N/2+1:N/2+N);     % extract relevant portions

figure,subplot(2,1,1)
plot(linspace(0,2,M),abs(H)); title('RAMP mask that is applied on X(w) in freq-domain')

subplot(2,1,2)
plot(t,x);
hold on
plot(t,0.1*y);
legend('x(t)','dx/dt');
grid on
title('x(t) vs dx(t)/dt  (amplitude scale adjusted)');

with exactly the same result on the time-domain sequence as shown below:


Can I construct a linear system so that it results in a moving average? I've just learned about linear systems and impulse response functions. I know that in a linear system consisting of $h_1[k]$ and $h_2[k]$, the impulse response of the system is $h_1[k] \ast h_2[k]$.

In an exercise I was given the following system responses:

\begin{align*}
h_1(0) &= a \\
h_1(1) &= b \\
h_2(0) &= c \\
h_2(1) &= d \\
h_1(n) = h_2(n) &= 0 \; \text{otherwise} \\
\end{align*}

The question was whether it would be possible to construct the parameters $a$, $b$, $c$, and $d$ in such a way that the linear system $h_1 \ast h_2$ would result in a moving average of the order 3 being applied to the input.

Well, I started with the convolution $h[k] = h_1[k] \ast h_2[k]$ for $k = 0, 1$:

\begin{equation}
h[k] = [ac, bc + ad, bd]
\end{equation}

I understand that for a moving average of order 3, the required convolution kernel would look like this:

\begin{equation}
[1/3, 1/3, 1/3]
\end{equation}

But that means I have to solve the following equation:

\begin{align*}
ac &= 1/3 \\
bc + ad &= 1/3 \\
bd &= 1/3 \\
\end{align*}

And that's where I'm stuck. From my guts I'd say it's not possible to get to an answer here, because even if $a = 1$ and $b = 1$, $c$ and $d$ would have to be $1/3$ and $c + d$ would have to be $1/3$ as well. Which isn't solvable.

The function you are trying to solve is equivalent to
$\frac{b}{3 \cdot a} + \frac{a}{3 \cdot b}=\frac{1}{3}$

It is also
$a^2-ab+b^2=0$

but this function does not have a real number solution.


Compute phase coherence between two eeg signals I have two EEG signals/channels and I would like to compute phase coherence in MATLAB. I have this paper: https://www.nature.com/npp/journal/v39/n5/full/npp2013330a.html but it does not say what Sij is so I don't know the formula.

Can anyone advise me on how to do this or does anyone have MATLAB code?

Many thanks.

Perhaps the built-in 'mscohere' and 'cpsd' functions may help you.

The mscohere function returns a value between 0 and 1 that measures the correlation between the signals, and the phase delay can be computed using the cpsd function, as per this example from the Mathworks website. It suggests that the relative phase between the correlated components can be estimated with the cross-spectrum phase.

https://uk.mathworks.com/help/signal/examples/measuring-signal-similarities.html#d119e1563

Taken from the Mathworks page above

[Cxy,f] = mscohere(sig1,sig2,[],[],[],Fs);
Pxy     = cpsd(sig1,sig2,[],[],[],Fs);
phase   = -angle(Pxy)/pi*180;


What is the right way to interpolate a 2D grid? Let's say I have a 2D grid of temperature measurements in some area and I want to estimate the temperature at some point between the samples. Or at every point, which would basically amount to upscaling an image. Which interpolation method would give the best temperature estimates?

I'm thinking that the "correct" approach would be something like this: assume a prior distribution of temperature in the area, use the measurements to update it with the Bayesian formula and then calculate either expected value or the most likely value at every point.

What a correct solution constitutes depends on the application. What you have proposed sounds more like an inference or estimation problem (that needs to be fleshed out). It would make sense if you thought your values were noisy.

If that were the case, try looking at Kriging.

https://en.m.wikipedia.org/wiki/Kriging

If you just want to interpolate MATLAB has interp2

Which (not limited to the  2 types mentioned above) is correct depends on what your ultimate purpose is.


IIR filtering vs IFFT filtering by zero binning and making bandpass IIR filters with high sample rates? I am having some issues making thin IIR filters that work at the sample rate of my data and was wondering if there is a way to improve my method? I am also curious about the deficiencies in using the method of filtering of applying and FFT, bin-zeroing and the IFFTing back to get the filtered data.

For example I use scipy.signal.iirdesign with scipy.signal.filtfilt to filter a 300KHz signal with a 5MHz sample frequency with bandwidth of 4KHz and a transition width of 1KHz. If I do this I find I get NaNs in the filtered signal, if I reduce the sample frequency by indexing every $N^{th}$ element of the data (e.g. by Signal[::N]) I get a correctly filtered signal if I use N=4 (i.e. sample frequency of 5/4MHz). However this is not ideal as I am having to reduce the sample rate of my data and can therefore get no higher temporal resolution. I can increase this frequency by widening the bandpass but this is undesirable as other frequencies are often present nearby.

I have also found I can create very thin filters by FFTing, zeroing bins and then IFFTing and this works great for a perfect sine wave signal with some noise. However I am interested in extracting a sine wave with random instantaneous phase jumps, and this appears (in simulated tests I have done where i know the true sine wave a-priori) to result in ripples in the amplitude of the signal with the ripples extending furthur the larger the phase shift between the first half of the signal and the next but this is sometimes a preferable method if I can not construct a filter with scipy that will work with the original sample frequency of the data and requires down-sampling.

The following plots demonstrate what I observe.

The blue is the true signal, the red is the filtered signal.

The first plots are a perfect sine wave with no phase shifts or amplitude variation:

IIR:

IFFT:

The IFFT works much better, but this is not that surprising to me since this is a perfect sine wave and is not particularly realistic to observe.

The next plots are of sine wave with a phase shifts of pi/4 in the center with no amplitude variation:

IIR:

IFFT:

The IIR filter used for this has a bandwidth of 1KHz and a transition width of 0.5KHz, it also has a sample frequency that is 4 times less than the original data. The IFFT filter had a bandwidth of 500KHz and appears to fit what the data is doing much closer, it also benefits from having the same sample frequency as the original data.

It would appear from this that using the IFFT filtering method is vastly superior but I know that it should introduce a sinc function to the time signal and the reason FIR and IIR filters are used to approximate a perfect rectangular window in frequency space is that doing so is a bad idea due to the introduction of this sinc function. What is wrong with IFFT filtering in this way with zero binning and how might I improve my method of bandpass filtering such that I can retain a high sample rate and a reasonable thin frequency window?

An IFFT of a single non-zero bin (or conjugate mirrored pair) produces a single pure sinusoid (with no phase or amplitude changes).  An IFFT of just a few adjacent non-zero FFT bins produces something like a slowly amplitude modulated sinewave (a beat note) possibly very slowly modulated in frequency.  Note that none of those IFFT results resemble a sinewave that includes a sudden instantaneous phase shift.  The information about the location of a possible sudden phase change is located in all the other FFT bins that an ultra-narrow FFT filter throws away.  e.g. you would have to de-zero a vast number of FFT bins to get that change to reappear in the IFFT result.  That is why sharp changes have very wide-band spectrums.

This renders your proposed filter (zeroing all but a few FFT bins) nearly useless for determining the location of any sharp phase shift.


Difference between CTFT and DTFT? I have tried to read different articles but still confused in difference between continuous time Fourier transform and discrete time Fourier transform?

The difference is pretty quickly explained: the CTFT is for continuous-time signals, i.e., for functions $x(t)$ with a continuous variable $t\in\mathbb{R}$, whereas the DTFT is for discrete-time signals, i.e., for sequences $x[n]$ with $n\in\mathbb{Z}$.

That's why the CTFT is defined by an integral and the DTFT is defined by a sum:

$$X(j\Omega)=\int_{-\infty}^{\infty}x(t)e^{-j\Omega t}dt\tag{1}$$

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

Note that in both cases, the frequency variables $\Omega$ and $\omega$ are continuous variables.

There is also the discrete Fourier transform (DFT), which can be applied to finite length (or periodic) sequences:

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\tag{3}$$

It is easily shown that DTFT and DFT are related as:

$$ X\left(e^{j\omega}\right)\Bigg|_{\omega = 2 \pi\frac{k}{N}} = X[k] $$

when $x[n]=0$ for all $n<0$ or $n\ge N$.

Note that the DFT of the length $N$ sequence $x[n]$ is also a length $N$ sequence $X[k]$. The DFT is important because there exist efficient algorithms to compute $(3)$. This class of algorithms is called Fast Fourier Transform (FFT).

The DFT can be used to approximately compute samples of the DTFT and of the CTFT. For finite length sequences, the DFT computes (exact) samples of the DTFT. For longer (or infinitely long) sequences, there is always a truncation error caused by considering only a finite length window of the data. When approximating the CTFT, there are two error sources: the truncation error, and the aliasing error caused by sampling the continuous-time signal. Since a signal cannot be time-limited and band-limited at the same time, at least one of those two errors is non-zero.


What is the Reference in Control Theory? I've started studying Control and I've come across some concepts that I'm finding  a bit difficult to understand.

For example, in the system showed here, there is a signal called "reference" as well as a feedback block. As far as I know, the feedback block has the function of getting the output as close as possible to the reference. But here is what I don't understand: what if the reference is, for example, an impulse, and I want the system to show a step when an impulse is at its input?

I'm getting confused with the difference between reference signal and input signal. Why would someone connect a feedback block in order to get the output to follow the reference? What's the advantage of doing this?

Imagine that you're heating (or cooling) a home with a modern furnace (or air conditioner).

the reference or set point is the temperature that you set your thermostat to be.  the feedback signal is the actual temperature that is measured with some kinda thermometer.

the actual value that you are trying to control, whether it's temperature or the position of a robot arm or the position of a pointer in an asynchronous sample rate converter, that actual value is compared to your reference value.  the difference between those two signals is what drives your "controller".

in the case of heating your home, if the actual temperature is below the reference (what you want your temperature to be), the furnace heat is increased.  if the actual temperature is above the reference, the furnace heat is decreased.


Maximum Likelihood Through a Noisy Channel I have random variables $X_1, X_2, \cdots, X_m $, which can take $n$ values and is distributed iid according to $\Theta=(\theta_1, \theta_2, \cdots, \theta_n)$. That is $X_k$ can take values $\{1,2,\cdots,n\}$ and $P(X_k=i)=\theta_i, \ 1\leq k \leq m, \ 1 \leq i \leq n$.

I observe a random variables $Y_1, Y_2, \cdots, Y_m$, where I have that $P(Y_k=j|X_k=i)=c_{j,i}$.

Let $C$ be the matrix $[C_{j,i}]_{1\leq j,i \leq n}$. I have that $C$ is a full rank matrix.

I want to get an ML estimate of $\Theta$, from observations $Y_1, Y_2, \cdots, Y_m$.

Let $\Psi=(\psi_1,\cdots, \psi_n)$, where ${\psi}_i=P(Y_k=i)$. The ML estimate of  $\Psi$ is clearly the vector of empirical frequencies. That is, we have the ML estimate of $\psi_i$, is $\hat{\psi}_i = \sum_{\ell=1}^{k} I(Y_{\ell}=i)$, where $I$ is an indicator function. Let us call this $\hat{\Psi}$.

I wanted to know if the ML estimate of $\Theta$, $\hat{\Theta}= C^{-1}\hat{\Psi}$. I can prove it if $n=2$, but do not seem to be able to show this for a general $n$.

Another way to put this is that I have a discrete memory-less source, which has a support of $\{ 1,2, \cdots, n\}$. I can observe this though a discrete memoryless channel, with transition probability matrix $C$. I want to estimate the source symbol probabilities. An ML estimate the observed symbol probabilities can be easily derived to be the empirical frequencies. I wanted to know if it can be shown that the ML estimate of the source symbol probabilities is $C^{-1}$ times the ML estimate of the observed symbol probabilities.

I'm inclined to think this is true, but so far I've only gotten a counterexample:

Consider the channel:

\begin{bmatrix}0.5 & 0.5 & 0 & 0\\ 0 & 0.5 & 0.5 & 0 \\ 0 & 0 & 0.5 & 0.5 \\ 0.5 & 0 & 0 & 0.5\end{bmatrix}

This a typewriter channel. This matrix is not invertible (it has zero determinant, and rank 2). So the ML estimate of $\Theta$ cannot be gotten from this channel.

However, for $p \neq 0.5$ I think it could well be true.

EDIT: In fact I think the noisy channel aspect of this may well be a red-herring. The ML estimator of $\Theta$ is as you describe, simply because you have a linear model. You don't need any conditions of $C$ other than it should be invertible (i.e. all rows are linearly independent).


Does the Nyquist Sampling Theorem hold for a triangular wave? Does the Nyquist Sampling theorem hold for a triangular wave produced by a function generator?

For example, a triangular wave with a frequency of 15 Hz sampled at 60 Hz, 180 Hz and 15000 Hz.

The Nyquist sampling theorem holds for bandlimited signals.  Note that a perfect triangle wave is composed of a (convergent) infinite series of harmonics, which requires infinite bandwidth and is thus not bandlimited.

However, the energy remainder after a suffinciently large number of harmonics might be small enough than for a given noise floor that you can treat the portion of the triangle wave that is above that noise floor as bandlimited below some point above that last kept harmonic.


Meaning of the frequencies vector using fft I'm so confused, and it would be nice if someone could clear my head, or at least direct me to suitable references.
Assume that I have two gaussian pulses with the same bandwidth in time, but for 2 different wavelengths. I'm going to handle them separatly- 2 separate graphs of one gaussian in time.
Both graphs are the same because they have the same bandwidth.
Then, I use fft on each graph.
Of course I'm going to have the same graphs after the fft- because I had the same graphs before.
Now, I want to creat the frequencies axis for each signal. In order to do that, I use the exact direction from matlab. df=1/T,  T=N·dt.
So the frequencies axis are also the same!
But how can it be if f=c/lambda?  I'm talking about two different wave length, how is it possible that the both have the same spectrum of frequencies?
Further more, I don't understand what is the different between F and df and why is F=N·df. I do understand that T is the time in which I sample, dt=T/N. but that's it.

You forget that what you did was describing the passband signals (i.e. some waveform on a carrier frequency == on a wavelength) in baseband.
So, of course, after calculating the baseband spectrum, you need to shift it back where it came from - to the carrier frequency. It's as easy as that.
If this isn't clear to you, you might want to revisit what (complex equivalent) baseband and passband are.


Remove background noise from audio file (python or matlab) I am implementing a project for infant cry detection and the audio set contains background noises. So for preprocessing i need to remove the background noise from the audio.I am not able to get a proper output for the code in jupyter notebook. For this code the output file does not contain anything even the baby cry is erased.

import numpy as np
import scipy as sp
from scipy.io.wavfile import read
from scipy.io.wavfile import write
from scipy import signal
import matplotlib.pyplot as plt
get_ipython().magic('matplotlib inline')

(Frequency, array) = read('cry.wav')

len(array)

plt.plot(array)
plt.title('Original Signal Spectrum')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

FourierTransformation = sp.fft(array)

scale = sp.linspace(0, Frequency, len(array))

plt.stem(scale[0:5000], np.abs(FourierTransformation[0:5000]), 'r')
plt.title('Signal spectrum after FFT')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

GuassianNoise = np.random.rand(len(FourierTransformation))

NewSound = GuassianNoise + array

write("New-Sound-Added-With-Guassian-Noise.wav", Frequency, NewSound)

b,a = signal.butter(5, 1000/(Frequency/2), btype='highpass')

filteredSignal = signal.lfilter(b,a,NewSound)
plt.plot(filteredSignal) # plotting the signal.
plt.title('Highpass Filter')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

c,d = signal.butter(5, 380/(Frequency/2), btype='lowpass') # ButterWorth low-filter
newFilteredSignal = signal.lfilter(c,d,filteredSignal) # Applying the filter to the signal
plt.plot(newFilteredSignal) # plotting the signal.
plt.title('Lowpass Filter')
plt.xlabel('Frequency(Hz)')
plt.ylabel('Amplitude')

write("New-Filtered-Sound.wav", Frequency, newFilteredSignal)

Can i know how to debug this code or get any other code.

Since you're adding white noise, the highpass and lowpass filtering will almost  not remove the noise in the frequency band where you want to keep your signal, so you will always have some background noise with this highpass and lowpass filtering strategy.

Not sure if this helps, it depends on the signal-to-noise ratio:
If you can clearly distinguish the noise from the signal in the spectrum (something similar as in the second figure of the Noisy Signal example in Matlab's documentation of the fft), you could set a threshold and make the spectrum with an amplitude below that threshold equal to zero before taking an inverse Fourier transform to get back to the (denoised) time-domain signal.

By the way, GuassianNoise = np.random.rand(len(FourierTransformation)) produces uniformly distributed noise. To get Gaussian noise, you need randn instead of rand.


Why Carrier signal can't have frequency less than Message signal I know that why we use Carrier of high frequency to send message signal over a long distance, when the frequency is low , energy will be obviously low. To increase the energy of the signal we need to increase the frequency. This is achieved by multiplying the message signal with the carrier signal (with high frequency).

suppose,

Carrier signal frequency = 2800KHz
message signal frequency = 3KHz

the two generated sidebands will be ,

2800 + 3 = 2803 KHz
2800 - 3 = 2797 KHz

and eventually, the bandwidth of the signal is,

BW = 2803 - 2797 = 6KHz

this is when the carrier frequency is higher , and we noticed that bandwidth is just twice of the highest frequency of modulating signal.
and the signal will be easily demodulated at the receiver end.

but what will happen if it is reversed ?

i.e

Message signal frequency = 2800KHz
Carrie signal frequency = 3KHz

please explain would happen here ?

suppose, Carrier signal frequency = 2800KHz message signal frequency = 3KHz

Then you will get a signal that looks like this in the frequency plane.

Obviously this is not to scale, but you get the idea.

but what will happen if it is reversed ?

i.e Message signal frequency = 2800KHz Carrie signal frequency = 3KHz

please explain would happen here ?

Then you will get the following.

The positive and negative sides of the signal will almost completely overlap, and thus distort each other to the point that they likely cannot be recovered.  Due to the +/- 3 kHz carrier frequency there will be 6 kHz that is non-overlapping (though they will still be distorted by the other signal's roll-off), but the rest will overlap.

To increase the energy of the signal we need to increase the frequency.

Signal power has nothing to do with increasing the signal frequency.  I think that you are thinking of light photons, which do increase in energy as their frequency goes up.  Communication signals, though they are (at a really, really, low-level) composed of photons, it is a non-issue because adding power just means adding more photons.

Anyway, suffice to say that you don't need to increase the frequency to increase the signal power.  We increase the carrier frequency for the following reasons:

Spectrum Availability.  You send your signal where the FCC (or whatever governing body applies) says you can send it.
Increased bandwidth.  As the answer to your question demonstrated, to send wide-band signals you need higher carrier frequencies.
Practicalities of antennas.  It is difficult (impossible?) to make a good wide-band antenna.  An antenna's bandwidth, though, is proportional to its center frequency, so increasing the carrier frequency makes it much easier to make good antennas with wide passbands.
Channel characteristics.  Different frequencies behave in different ways.  Some frequencies get absorbed by rain and some resist that.  Some frequencies bounce off the ionosphere and so can travel farther than "line of sight".


Linear Prediction Coefficients I want to learn how to compute linear coefficient, but to be honest - I can't.
I know how to use them if i have it and why I want them.

Suppose this toy example:

A signal $x=[-1, -2, 0, 1]$

According to matlab function lpc(X,p) if I want to compute the 2th order coefficients of this signal the result is $a = [1.0, -0.50, 0.50]$ and $ g=1.0$

I Know how to compute the g using $A(k)$ formula and that $a_0=1$ always, but can someone explain me why $a_1=-0.5,  a_2=0.5$ how to compute them without using matlab or any other software.

P.S It there any method which allow me using  $A(k)$ to compute directly the coefficients?

thanks

First you need to compute the signal's autocorrelation:

$$R(k)=\sum_nx(n)x(n+k),\quad k=0,1,\ldots,p$$

In your case with $p=2$ you get

$$\begin{align}R(0)&=(-1)(-1)+(-2)(-2)+0\cdot 0+1\cdot 1=6\\
R(1)&=(-1)(-2)+(-2)\cdot 0+0\cdot 1=2\\
R(2)&=(-1)\cdot 0+(-2)\cdot 1=-2\end{align}$$

Then you have to solve a system of linear equations:

$$\begin{bmatrix}R(0)&R(1)\\
R(1)&R(0)\\
\end{bmatrix}
\begin{bmatrix}a_1\\a_2\end{bmatrix}=
-\begin{bmatrix}R(1)\\R(2)\end{bmatrix}$$

which gives the values $a_1=-0.5$ and $a_2=0.5$. Of course $a_0=1$, as you've already noted.

The matrix that occurs in the system of linear equations is a Toeplitz matrix, i.e. all its diagonals have the same elements. As a result, such systems can be solved by specialized and efficient algorithms, such as Levinson's algorithm.


How does viola jones algorithm cope up with a test image size which is very different to training image sizes? To my understanding Viola jones algorithm is a combination of weak classifiers built on training imageset. Lets say we have training images of 24x24 and algorithm selected 100 haar features and built a classifier. Assume there is a test image which is of size 512x512 and there two persons(just faces) in the image one person(100x100 size in 512x512) closer to the camera and the other(10x10 size in 512x512) at a distance.

How does viola jones detect faces detect both persons having learnt
features      features for faces that fit in 24x24?
Is the algorithm implemented in libraries in such a way that haar
features vary depending on the size of the image

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


Deriving the Langrangian interpolation polynomials in Cook-Toom convolutions I'm working through Blahut's 'Fast Algorithms for Signal Processing'. Trying to develop an intuition for the Cook-Toom algorithm for convolutions as used by Lavin and Gray in their Winograd paper for convnets.

At page 151 (chapter 5.2), Blahut writes that the following is the final result of the Cook-Toom algorithm:

$s(x)=s(\beta_0)L_0(x)+s(\beta_1)L_1(x)+s(\beta_2)L_2(x)$

I'm a bit confused about why one has an $s(x)$ and then has an $s$ instead a few lines later like so:

$s=C[Bg·Ad]$

But that is secondary. I'm guessing that's to do with us no longer treating $s$ as a polynomial after the inverse transform with $C$? But I though convolution was equivalent to polynomial multiplication. Still not entirely clear on that.

My main question is, how are the following Lagrange polynomials derived? For reference, the interpolation points are $\beta_0 = 0, \beta_1 = 1, \beta_2 = -1$:

$L_0 = -x^2 + 1,$

$L_1 = x^2 + x,$

$L_2 = x^2 - x.$

These are stated, but don't follow obviously from the rest of the text. My understanding of Lagrange polynomials is somewhat limited. Any clarification would be very much appreciated!

If you represent a second-order polynomial $s(x)$ with Lagrange polynomials $L_i(x)$ and interpolation points $\beta_i$, $i=0,1,2$, such that

$$s(x)=s(\beta_0)L_0(x)+s(\beta_1)L_1(x)+s(\beta_2)L_2(x)\tag{1}$$

then for the equality in $(1)$ to be satisfied, the polynomials $L_i(x)$ must have zeros at $x=\beta_j$, $j\neq i$, and they must equal $1$ at $x=\beta_i$. This completely determines the polynomials $L_i(x)$.

E.g., the polynomial $L_0(x)$ must have two zeros at $\beta_1=1$ and $\beta_2=-1$, respectively, and it must satisfy $L_0(\beta_0)=L_0(0)=1$. The latter condition just determines the scaling. This results in

$$\begin{align}L_0(x)&=\frac{(x-\beta_1)(x-\beta_2)}{(\beta_0-\beta_1)(\beta_0-\beta_2)}\\&=-(x-1)(x+1)\\&=1-x^2\end{align}$$

The other polynomials can be derived in the same way.

Note that the polynomials $L_1(x)$ and $L_2(x)$ have a scaling factor of $\frac12$, which is later in the derivation absorbed in other constants for reasons of computational efficiency.


Use of complex numbers to calculate Doppler effect? First of all, I'm a beginner in this field, please tell me if any of what I'm asking is too basic / can be found with some google research.

I have a situation with a static TX and a moving RX (v is its velocity vector):

I know the value of $\vec{v}$ and the position of $TX$ and $RX$.

My goal is to simulate the Doppler shift that the $RX$ would experience, with the desired output similar to this:

so I ended up with the following equation:

$$ f' = \frac{c + v\cos(\theta)}{c}f $$

and I ended up with

(This is LOS, so I didn't consider reflectors, so no noise)

Question

This seems to work fine; however in some papers (and I have been told) that to apply the doppler shift to a signal, one must apply a rotation to it, so multiply by $e^{i\theta}$. Is this useful / necessary for Doppler ? In which cases should I use it ? How ?

Bottom Line First

What you suggest will only work if your signal is at baseband and not at a carrier, which is usually the best way to do any simulation work regardless. If your signal is not at baseband first multiply by $e^{-j2\pi f_c t}$ where $f_c$ is the carrier frequency and low pass filter out the higher image (which will now be centered at $-2f_c$) to create the baseband equivalent waveform. With this you can simply multiply by $e^{\pm j2\pi f_d t}$ to introduce a Doppler shift. If you need to have your waveform at a carrier frequency in simulation, then you can shift the Doppler-shifted signal back to the carrier after the Doppler shift operation at baseband described above by multiplying the resulting shifted waveform by $e^{+j2\pi f_c t}$ and then taking the real part if a real signal is desired. Usually however ALL simulation work can be done at baseband with a complex envelope that is equivalent to what is happening at any carrier frequency and much simpler to model.

(Note if the concept of baseband vs passband mentioned above is confusing, see this post that goes into further details on that: modulate complex time domian signal to a carrier frequency)

DETAILS

Simply multiplying any signal by a complex rotation will not work for simulating a Doppler shift in all cases; that applies to when the signal is represented by its complex equivalent at baseband, or if you are only concerned with the complex single-sided equivalent of your signal.

This is because the Doppler will shift a signal either higher in frequency or lower in frequency depending on the relative direction (as you have properly shown), and a real signal is comprised of both positive and negative frequency components, each described by $e^{j2 \pi f t}.$  and multiplying by a single complex rotation will shift EVERYTHING in the same direction.

Multiplying a single complex tone at carrier frequency $f_c$ ( $e^{j2\pi f_c t}$)  by a Doppler frequency $f_d$ ($e^{j2\pi f_d t}$) is simply summing the exponents which results in the expected frequency shift:

$$e^{j2\pi f_c t}e^{\pm j2\pi f_d t} = e^{j2\pi (f_c \pm f_d) t}$$

(And the same would occur is we were dealing with a modulated signal as a band of frequencies around $+f_c$.)

However, look what happens if the signal of interest is a real cosine wave which is comprised of TWO complex conjugate symmetric complex tones:

$\cos(2 \pi f_c t) = \frac{1}{2}e^{j2\pi f_c t}+ \frac{1}{2}e^{-j2\pi f_c t}$

$\cos(2 \pi f_c t)e^{j2\pi f_d t} = \frac{1}{2}e^{j2\pi (f_c+f_d) t}+ \frac{1}{2}e^{-j2\pi (f_c-f_d) t}$

Thus we see that for a cosine wave (or any real signal) which has a positive frequency spectrum and a complex conjugate negative frequency spectrum (same magnitude, opposite phase), all components will shift in the same direction-- NOT what occurs for a Doppler shift.

I would highly recommend doing all such simulations at baseband (remove the carrier frequency and replace with 0 Hz and deal with complex signals expressed as $e^{j\omega t}$, not sines and cosines) in which case you can do such an operation which is significantly simpler. However if it absolutely needed to simulate with a carrier, to shift a real signal at a carrier frequency by a Doppler frequency, you need to first take the Hilbert transform of your signal, which when summed in quadrature with the original signal, effectively removes all the negative frequencies (in practical implementation, this is a quadrature splitter and the resultant complex signal is implemented as two separate real channels), and then with this new complex signal, you can multiply by $e^{\pm j2 \pi f_d t}$ as suggested to shift with the sign chosen dependent on the direction desired. After multiplying you can take the real of the result to re-establish your original real signal with the Doppler shift induced. Note this is exactly how a single-sideband modulator is structured as demonstrated by the slides pasted below:

Note, the operation above is showing how to shift a single-sideband signal at $\omega_{\Delta}$ away from the carrier down to baseband, but the same approach can be used to show how to shift a signal at $\cos(\omega_c)$ to $\cos(\omega_c+\omega_{\Delta})$ where $\omega_{\Delta}$ represents the Doppler shift. Also to note in the last diagram, the Hilbert transform will convert the cosine to a sine, and must be added to the sine with a $j$ term to become the desired exponential rotation:

$e^{j\omega_c t}= \cos(\omega_c t)+ j \sin (\omega_c t)$

In the implementation the required j term shown in front of the sine is implied by carrying that path separately (as is the case with having I and Q paths in a digital modem implementation for example); we need two real numbers to represent a complex number, as is done by having the two signal paths. (either as real and imaginary components such as $I + jQ$ or as magnitude and phase components suchs as $A\angle \phi$, which by the way is equivalent to $Ae^{j\phi}$)

The above approach based on a single-sideband modulator gives excellent insight into how to shift a signal at a carrier frequency in one direction (as happens with Doppler shift), and it would be exactly how a hardware implementation could be done at RF to simulate a Doppler shift. For simulation work my favored approach would be to do all simulation at baseband using the "Bottom Line" approach first described.


SRRC filter -- filter span tradeoffs I am simulating (in Matlab) the transmit portion of a comm system.  I generate random data, then I want to modulate it to BPSK, pulse shape the data, and then potentially shift it up in frequency by a carrier.  My question is regarding the pulse shaping portion.

Atleast in matlab, to generate the root raised cosine filter using the rcosdesign function, you need to know the roll-off factor  ($\beta$) , as well as the samples-per-symbol ($sps$) and the filter span ($span$).  The samples per symbol will be dictated by the sampling frequency $f_s$ as well as the desired length of the generated signal with a certain datarate $R$.  For example, if I have a BPSK signal with a sampling frequency of 8192 Hz, and I want to generate 1 second of data at a datarate of 256 Hz, then

$$
sps = f_s/R
$$

Now we've established 2 of the inputs to creating the filter are set, the rolloff factor $\beta$ and the samples per symbol, $sps$.  Now, my understanding is that the filter span is kind of arbitrary.  I know that a longer filter span creates a longer group delay, but what are the other tradeoffs to having a longer filter span versus shorter?

Lower rolloff means less excess-bandwidth and therefore longer time span.
That is, you will need a longer filter to accurately reproduce the ideal RRC filter response (increases implementation complexity)
Longer filter means more accurate RRC filter and potentially lower ISI at the receiver side (depends on your target SNR)
Longer filter means less noise out of band transmitted (infinite length filter would have zero power outside symbol-rate*(1+rolloff)

In practice, due to the above constraints, RRC filters are usually not designed by just truncating the ideal RRC impulse response. Either you multiply the impulse response with a window function or you directly optimize the filter coefficients for error and rejection.


No cross-correlation between same signal recorded using two different microphones I've recorded the same audio signal using two microphones located at the same approximate location. However, when checking the cross-correlation using MATLAB (using the xcorr function) it is negligible, even though the signals look approximately the same and sound the same. The signals were brought to MATLAB using the audioread function. This scenario was repeated using a different set of microphones, achieving the same results.

Did anybody encounter this kind of problem in the past and has a solution?

EDIT:
First of all, thank you all for the help.
This is a scenario of two sources and two microphones in different locations in space. The signals clearly sound mixed (each signal contains both sources) and the microphones are synchronized and at the same sample rate of 44.1kHz.
The code:

[x1,Fs1]=audioread('channel_1.wav');
[x2,Fs2]=audioread('channel_2.wav');
figure; subplot(2,1,1); plot(x1);subplot(2,1,2);plot(x2);
[R_11,lags]=xcorr(x1,x1);
R_22=xcorr(x2,x2);
R_12=xcorr(x1,x2);
figure;plot(lags,R_11);hold on;plot(lags,R_22);plot(lags,R_12);

The correlations:

The signals:

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!


Can I sample at Nyquist rate if I know that my samples are taken only at the signal's maxima or minima? I know that in general the sampling rate, $f_s$, must be greater than twice the highest frequency of the signal, $f$.
If I sample at the Nyquist rate, it can lead to the following:

However, if the signal is sampled only at maxima and minima, then we get this:

My signal is not pure sinusoidal, it is something like this:

Will I loose any frequency components if I sample it at the Nyquist rate?

when you have a quantized signal (like in your case) you also have a noise floor in the frequency spectrum. This noise floor is proportional also to the sampling rate. If you sample a signal at twice its maximum frequency you can distinguish all the frequency components above the noise floor. In your case you are sampling at $F_p = 2$ because there are 2 samples every second therefore you are assuming your signal has a maximum period of 1.
There is one case when you are interested on sampling at the same frequency of the signal. For example if you know that your signal has an unwanted ripple at the frequency $f_{ripple}$ you can sample at this frequency $F_p = f_{ripple}$ preserving the lower frequencies while the ripple is transformed into a constant error which basically depends on the ripple waveform.


Creating a spectrogram I've been trying to work out the logic for this task, and plan to use the KissFFT source package to perform the fast fourier transform.  Please let me know if this looks about right:

Allocate an FFT structure, ie. kiss_fft_alloc(N,0,NULL,NULL)
Where N is the window size I'm using.  The input buffer will be an array of N elements of type kiss_fft_scalar.  The output buffer will be an array of N/2 + 1 elements of type kiss_fft_cpx.
Decode N (window size) number of PCM samples.
For each PCM sample, average each channel's amplitude (unsigned samples) and scale from 0 to 2 (divide by 65536.0), storing the result into the input buffer.
Perform windowing (ie. Hanning) on the input buffer.
Perform fast fourier transform on the input buffer, storing into the output buffer.  Since I am using real values as input, I can use kiss_fftr().
For the N/2 output values, obtain the squared magnitude of the transformed data and convert the values to the dB scale with the following formula:
10 * log10 (re * re + im * im)
Plot the N/2 values from step 6.
Discard the first half of the input buffer, decoding the next (window size / 2) PCM samples and performing scaling and windowing to the data.  This should effectively slide the input window and avoid having to redo math on processed PCM samples.
Loop to step 5, repeating these steps until all samples are processed.
Free the used memory from kiss_fft_alloc().

It was suggested that I subtract a value from the input window before I perform the FFT, so that the resulting DC value has a magnitude of zero.  Should I be subtracting the mean or the average from the input data?

Also, what are the things I need to consider when I choose a window size?  Besides that it has to be an even number as per KissFFT's instructions, is there a benefit to using a small window size, ie. will it provide for a better graph?  I assume that a large window size reduces the number of FFTs that must be performed, is that the only benefit to using a large window size?

Lastly, when I get to the point that the data is ready to plot, how do I go about plotting it?  When I worked on some waveform graph logic in the past, I just plotted 3 values for each pixel along the $x$-axis (min amplitude, max amplitude, RMS amplitude), but I don't know what I'm supposed to do with spectrogram data.

Thank you in advance for any and all guidance you can provide.

Looks pretty good to me. In step 3, though, you actually want to scale the signal from -1 to 1, otherwise you're adding DC. You mentioned subtracting the mean -- I would not recommend doing this for a spectrogram, since that's effectively filtering out DC, which the spectrogram ought to show if it's there.

Choosing a window size is all about tradeoffs. A larger window will give you sharper frequency  resolution, but blurrier time resolution. A shorter window will give you the opposite: sharper time resolution but blurrier frequency resolution. The appropriate choice of window size will depend on the data you are trying to analyze. Typically it will be a power of 2 just because FFTs tend to like powers of 2. A decent rule of thumb is that your window should be at least roughly twice as long as the period of the lowest frequency you would like to be able to accurately resolve.

You might wonder if it is possible to better deal with this tradeoff, and there are techniques for that: they generally involve computing spectrograms with several different FFT sizes at once, and combining them. There's some good visual information on this webpage: http://www.izotope.com/tech/aes_adapt/

If your window size is too small, two very close frequencies might be indistinguishable from each other since they both end up in the same FFT bin. If your window size is too large, two  close events in time, might be combined, or a sharp transient might turn into a gradual attack. Check out that webpage I posted for some ways to visualize this.

A larger window size doesn't necessarily reduce the number of FFTs. You have chosen to compute a spectrogram using a short-time Fourier transform where there is an overlap of half the FFT size. You could use a higher overlap factor if you wanted. Choosing a window size is much more a matter of making the time/frequency tradeoff than how many FFTs you have to compute. In designing a spectrogram (or any STFT), you can think of choosing your window size, and hop size, the distance between blocks, as independent parameters.

When you plot it, time is typically on the x-axis, frequency is on the y-axis (usually a log scale, Mel scale, etc, rather than a linear scale), and then the magnitudes are represented with color intensity, i.e. very dark colors correspond to small magnitudes and very bright colors correspond to large magnitudes.


When calculating the Nyquist frequency should carrier frequency be included This isn't my field, so please forgive any misnomers on my part.

I'm looking into buying a amateur radio receiver, and I'm thinking about what sampling rate I need. I'm most interested in satellite signals, which transmit on ~400MHz (for example). So that would imply I have a nyquist frequency of ~800MHz, or 800 million samples per second. This seems much higher than any hardware I can find.

My current thinking is that the frequency to be used in the nyquist calculation isn't the carrier frequency. Which frequency should I use?

You are correct: a 400 MHz signal needs to be sampled at least at 800 Msamples/s, and that is either impossible or extremely expensive. However, this is not the entire story:

Under some conditions, you can do bandpass sampling, where the sampling frequency is lower than Nyquist. See more at https://en.wikipedia.org/wiki/Undersampling.
You can do analog bandpass filtering and down-conversion to some intermediate frequency (IF), which can then be sampled at reasonable cost/complexity. See https://en.wikipedia.org/wiki/Intermediate_frequency.

Most amateur digital radio receivers I'm aware of use the IF approach, and some also do bandpass sampling. For example, the RTL2832 can only sample at 2.6 Msamples/s, but it can receive signals at up to 2.2 GHz or so. If you see the device's datasheet, you'll see that it supports several different IFs.


Sampling the Dirac function I would like to ask a theoretical question concerning the Dirac function. The Fourier Transform of the Dirac function is the value 1 (DC) for every frequency. If we consider the Sampling Theorem, we have to find a maximum frequency in the signal $ \ f_{max} $, so that we can sample with $ \ f_s \ge \ 2f_{max}$. But as we can see from its Fourier Transform, the Dirac function contains every frequency, so we cannot find an appropriate $ f_s $ . My question is, from a theoretical point of view, can the Dirac function be sampled?

Edit: Thank you for your helpful answers guys!

I completely agree with Juancho's answer. I would just like to add a few things. I think the main problem is the misunderstanding which becomes evident in the last sentence of the question: "... can the Dirac function be sampled?" The Dirac impulse is NOT an ordinary function having definite values for every $t$, but it is a distribution (even though it is often called 'Dirac function'). One should therefore not try to 'evaluate' (or sample!) it. What is important about the Dirac impulse are its integral properties:

$$\int_{-\infty}^{\infty}\delta(t)dt = 1$$
and
$$\int_{-\infty}^{\infty}\delta(t-t_0)f(t)dt = f(t_0)$$

As Juancho already pointed out, the square of a Dirac impulse $\delta^2(t)$ is not defined. So if you were to sample a Dirac impulse, you would get the undefined result

$$\sum_n\delta(t-nT)\delta(t)=\delta^2(t)$$

Dirac impulses are a convenient tool for analyzing linear time-invariant systems but they should be treated with care because common types of  processing performed on ordinary signals (such as sampling) may lead to undefined and meaningless results when applied to Dirac impulses.


Actual Sample Rate slightly off? When capturing microphone audio, I noticed that my static buffer eventually overflows.
So I timed the samples, and low and behold I'm getting an actual sample rate of 44102-ish Hz (my timer is very accurate).
Is this the actual sample rate recorded by my microphone, or is windows sending me samples twice?
Doesn't this cause all sorts of synchronization problems. Say if you're listening to a live recording, wont a noticeable lag start to accumulate if sample rates aren't exact? How is this problem solved traditionally? Is dropping a sequence of samples, or padding zeros normal?

If you calculate the error, 44102/44100 is only about 45 parts per million. That is well within operating tolerance of many crystal oscillators used in consumer equipment to generate audio sampling rates and USB communication clocks. You can be quite pleased it is not even more off. Another problem is that if you simultaneously use different devices for audio capture or playback, they each have their own local clocks so they are not synchronized in any way and they all run free with their own local clockd, and thus you will have to do some kind of resampling to adapt the slight sampling rate differences.


Help with Time shifting a discrete signal I have the following problem statement which I have no idea how to approach,

I have a IR depth sensor, that produces a distribution for it's depth value. So it kind of looks something like this:

As you can see, it shows that there are two possibilities, one at 10 and one at 20. Taking the expectation of this will give me the real depth.

What i want to do is apply some signal processing to model the motion of this camera in the Z axis. Basically i want it to shift say 5m and spread out to represent the increasing uncertainity. Hopefully to look something like this

I guess for the spreading out part, applying a 1D gaussian convolution would work, but i am not sure about the timeshifting part. Can i convolve the signal below with another signal to actual timeshift it somehow?

yes, you can convolve with $\delta[n - n_o]$ to delay your signal by $n_o$ samples.

I am guessing you know how to represent $\delta[n]$ in MATLAB. And for spreading the signal, you can do by time scaling as $x_{spread}[n] = x[\frac{n}{M}]$. So, basically, the overall operation becomes:$$x_{final}[n] = x[\frac{n}{M} - n_o]$$

In MATLAB you can simulate this by following code:

x = sin(2*pi*(0:0.01:1));
y = sin(2*pi*(0-0.1:0.005:1-0.1));
plot(x);hold on;plot(y);

Here, time-spread is by a factor of $2$ after shifting the signal by $10$ samples.


DFT (FFT) of Non Uniformly Sampled Signal I am trying to perform an FFT of a non-uniformly sample signal. My input comes from Simulink and PLECS which uses variable-time solver.

So far, I have used resample function from Matlab to resample the values to a uniform sample rate and perform the FFT but with this set of values I cannot get it to work no matter what. I guess it's due to the nature of the signal which has a lot of jumps but I am not sure.

In this particular case I'd like to identify the component at 300Hz but for general purpose it can be any component. The code I tried is below if anyone wants to give it a try and help me. For comparison I also have the results from PLECS scope which indicates the 300Hz components

Link to raw data

Link to PLECS FFT data

clear all;
close all;
tab = csvread('data1.csv',1,0);
tout = tab(:,1);
in = tab(:,2);
fSamp = 1e6;
% Resampling
a(1) = (in(end)-in(1)) / (tout(end)-tout(1));
a(2) = in(1);
% detrend the signal
indetrend = in - polyval(a,tout);
[ydetrend,ty] = resample(indetrend,tout,fSamp,'linear');

out = ydetrend + polyval(a,ty);
% FFT
fr_des = 300;
amp = abs(fft(out))/length(out);
phs = unwrap(angle(fft(out)));
Fv = linspace(0, 1, fix(length(out)/2)+1)*fSamp/2;           % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector                                          % Desired Frequency
ampv = amp(Iv);                                         % Trim To Length Of ‘Fv’
phsv = phs(Iv);                                         % Trim To Length Of ‘Fv’
ap = [ampv(:) phsv(:)];                                 % Amplitude & Phase Matrix
ap_des = interp1(Fv(:), ap, fr_des, 'linear');
FFT = ap_des(1)*exp(1i*ap_des(2));

I am trying to perform an FFT of a non-uniformly sample signal. My input comes from Simulink and PLECS which uses variable-time solver.

Assuming that the first column in the data provided is your time variable ($t$) in units of seconds, the variable-time solver introduces a very small jitter.

(All code in Octave)

x = csvread("somefile.csv");
time_variable = x(3:end); %Skipping the first few samples as the solver settles.
plot(diff(time_variable)); %Ideally, this should return a flat plot at the sampling period, but that is not what we have here.

This results in:

Where you are sampling with at least a sampling period of something like $Ts \approx 1.68e-07$ extending to $Ts \approx 7.65e-07$.

The period of a 300Hz sinusoid that you are interested in is something like $0.00333 \ldots$ and the worst case scenario given above introduces an uncertainty of $\approx 0.022 \%$. That is per-cent.

Now, I am not sure which waveforms you are trying to obtain the spectrum of, but the second one is the one with the most discontinuities (the control signal?) and its spectrum looks like:

Fs = 1/7.6587e-07;
NFFT = 512;
FreqVec = Fs.*(0:(NFFT-1))./(NFFT-1);
semilogy(abs(fft(x(3:end,2),NFFT));
xlabel("Frequency (Hz)");ylabel("Amplitude");grid on;

Where as you can see, I have used the Sampling Frequency ($Fs$) of the worst case scenario, again, assuming that the first column contains the time variable in seconds.

Hope this helps.


Circle fitting for a binary image I a new to MATLAB and i am trying to apply circle fitting technique to the following image.

Image

Expected Output

I used the following code with no luck:

[centers, radii] = imfindcircles(W,[30 65],'ObjectPolarity','bright');
figure,imshow(W);% title('Sobel gradient');
viscircles(centers,radii,'Color','b');

Note that the circle size can vary between picture

Any suggestions.

imfindcircles can take in other "Name, Value" pairs as parameters. In your case, a really helpful one is the Sensitivity parameter. I would change your first line to include this (also note that 'bright' is the default value for 'ObjectPolarity', so I took it out). The value `0.95' here can be tuned to get more or fewer circles identified.

[centers, radii] = imfindcircles(W, [30 65], 'Sensitivity', 0.95);

Since you mentioned circles may vary between pictures, it would be beneficial use some other parameters to get better results across the set of images. Check out the MATLAB documentation online for specifics.


simulation of lin-log amplification I struggle with two simple questions.

First. What is the distribution of the signal on the output of the lin-log amplifier if Rayleigh noise is driven to the input of it.

Using an approach with the inverse function it is possible to obtain the expression for the distribution of the output signal. It looks like the First Extreme Value Distribution (i.e. with $\exp(\exp(x))$ pattern).
In other words, if $f_{in}$ is the input of the lin-log amplifier:
$$
f_{in} = \frac{r}{D}e^{\frac{-r^{2}}{2D}}
$$
one can get $f_{out}$ - the  distribution of the output signal:
$$
f_{out} = \frac{e^{-\frac{e^{2x}}{2D}-2x}}{D}
$$

But I found a statement somewhere in engineers forums that a lin-log amplifier works more trickily and "takes a logarithm of an envelope" and "there were a lot of problems with erroneous simulations of lin-log signal transformation". It was just statements without explanations that gave me feed for thoughts.
I found experimental data. It looks more like:
$$
f_{out} = \frac{\log(x)}{K_{1}}e^{\frac{-\log(r)^{2}}{K_{2}}}
$$
I'm afraid that it is very similar to the First Extreme Value Distribution and they are indistinguishable in an experimental environment. For example, the additional logarithm of both of them results in Gaussian-like distributions and the difference between them are again almost completely indistinguishable.

Second. How to simulate the lin-log amplifier (detector) in R or MATLAB.
The problem is that taking the logarithm of a random variable with Rayleigh's distribution results in a distribution with a negative skewness which does absolutely not resemble the above-mentioned alternatives expected. The same distribution with a negative skewness one can get using normally distributed in-phase and quadrature components of the signal.

Is there anybody who understands the problem and could share the idea on how to represent the output of the lin-log detector mathematically and simulate it numerically?

Thanks!

The distribution of the log magnitude output of a log amplifier would be according to the log of the magnitude of the input. Although in implementation the log amp is not simply taking a log of the input signal, the result is indeed the log with reasonable accuracy (typically linear to the log scale within a fraction of a dB). Below shows how this should appear for both Rayleigh and Ricean distributions. Given the OP is not seeing the same skewness as is the case for the Rayleigh at the log magnitude output of the log amp, I would conclude that the input signal (which he does not have access to) was not Rayleigh distributed.

Notice the skewness decreases significantly as the dominant path in the Ricean distribution (A) becomes more prominent, as would be expected.


Complex IIR to Real IIR I have created an IIR design algorithm that generates complex coefficients (there is no symmetry in the poles and zeros). However, the IIRs will be used to filter a real signal. Is there a closed form (analytical) way to convert the IIR to real coefficients. I did think that using the output of the difference equation to solve a system of linear equations might be a known solution?

Another way to put the question: when both complex and real IIRs are acting as a digital filter processing discrete samples, for a given real input they produce the same real output (discarding the imaginary component).

Another way to put the question: Both the complex and real IIRs have the same magnitude response.

This answer shows how to create a real-coefficient infinite-impulse-response (IIR) filter, the output of which equals the real part of the output of a given complex-coefficient IIR filter. Also an example is given that shows that the approach can be useful.

Z-transform $\mathcal{Z}\{h[n]\}$ of a time-domain sequence $h[n]$ has the property that:

$$\mathcal{Z}\{h^*[n]\} = H^*(z^*),\quad\text{where }H(z) = \mathcal{Z}\{h[n]\}.\tag{1}$$

Given an IIR filter transfer function:

$$H(z) = H_0\frac{\prod^M_{k=1}\left(1-q_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)},\tag{2}$$

with real coefficient $H_0$ and $M$ complex zeros $q_k$ and $N$ complex poles $p_k,$ discarding of the imaginary part of its output $h[n]$ to form real output $g[n]$:

$$g[n] = \operatorname{Re}(h[n]) = \frac{1}{2}\left(h[n]+h^*[n]\right),\tag{3}$$

is equivalent to having a filter with transfer function:

$$G(z) = \frac{1}{2}\big(H(z) + H^*(z^*)\big).\tag{4}$$

Equivalently, with $z = e^{i\omega},$ this new filter has a frequency response $\frac{1}{2}\big(F(\omega) + F^*(-\omega)\big),$ where $F(\omega)$ is the frequency response of the full complex filter. Using further properties of complex conjugation:

$$G(z) = \frac{1}{2}H_0\left(\frac{\prod^M_{k=1}\left(1-q_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)} + \frac{\prod^M_{k=1}\left(1-q^*_kz^{-1}\right)}{\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)}\right)\\
= \frac{\frac{1}{2}H_0\left(\prod^M_{k=1}\left(1-q_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right) + \prod^M_{k=1}\left(1-q^*_kz^{-1}\right)\prod^N_{k=1}\left(1-p_kz^{-1}\right)\right)}{\prod^N_{k=1}\left(1-p_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)}\\
= \frac{\mathcal{Z}\bigg\{\operatorname{Re}\Big(\mathcal{Z}^{-1}\left\{H_0\prod^M_{k=1}\left(1-q_kz^{-1}\right)\prod^N_{k=1}\left(1-p^*_kz^{-1}\right)\right\}\Big)\bigg\}}{\prod^N_{k=1}\Big(\left(1-p_kz^{-1}\right)\left(1-p^*_kz^{-1}\right)\Big)}.$$$\tag{5}$

By expanding the products to polynomials of $z^{-1}$, the above gives the procedure to calculate the coefficients of the real filter given the coefficients of the original complex filter, by convolutions of sequences of polynomial coefficients. If the original complex filter with input $x[n]$ and output $y[n]$ is in form:

$$\begin{align}y[n] &= b_0\,x[n] + b_1\,x[n-1] + b_2\,x[n-2] + \ldots\\
         &-a_1\,y[n-1] - a_2\,y[n-1] - \ldots,\end{align}\tag{6}$$

that readily accepts the coefficients $B = [b_0,\,b_1,\,\ldots]$ and $A = [a_0,\,a_1,\,\ldots]$ of the numerator and denominator polynomials $b_0+b_1\,z^{-1}+b_2\,z^{-2}+\ldots$ and $a_0+a_1\,z^{-1}+a_2\,z^{-2}+\ldots,$ respectively, with a hidden coefficient $a_0 = 1$, then the coefficients $B_\text{real}$ and $A_\text{real}$ for the new filter will be given by:

$$\begin{align}B_\text{real} &= \frac{1}{2}\left(B*A^* + B^**A\right) = \operatorname{Re}\left( B*A^*\right)\\
A_\text{real} &= A*A^* = \operatorname{Re}\left(A*A^*\right),\end{align}\tag{7}$$

where $*$ represents convolution (or complex conjugate if superscript). The last $\operatorname{Re}$ is unnecessary, but may be useful in numerical convolution.

Example 1: Minimal example

A minimal example in Octave starts with a complex filter and generates a real filter that matches the real output of the complex filter:

b = [0.5, 0.35 - 0.35*j]; #Coefficients of numerator polynomial
a = [1, - 0.6 - 0.6*j]; #Coefficients of denominator polynomial
length = 1024;
impulse = eye(1, length);
[h_complex, w_complex] = freqz(b, a, length, "whole");
[h_discard_imag, w_discard_imag] = freqz(real(filter(b, a, impulse)), [1], length, "whole");
b_real = real(conv(b, conj(a)))
a_real = real(conv(a, conj(a)))
[h_real, w_real] = freqz(filter(b_real, a_real, impulse), [1], length, "whole");
freqz_plot([w_complex, w_discard_imag, w_real], [h_complex, h_discard_imag, h_real]);

The generated real coefficients of the numerator and denominator polynomials of the new real filter are:

b_real =
   0.50000   0.05000   0.00000

a_real =
   1.00000  -1.20000   0.72000

The frequency response of the new real filter equals that calculated from the real part of the impulse response of the complex filter:

Figure 1. Blue: Frequency response of the complex filter with the imaginary part of its output still present. The frequency response lacks Hermitian symmetry, so no real-coefficient filter can reproduce it. Red: Shared frequency response of the new real-coefficient filter and the complex-coefficient filter with the imaginary part of its output discarded. The shared frequency response has Hermitian symmetry.

Example 2: Sliding discrete Fourier transform

Constructing a real filter from a complex filter can sometimes be useful. One such case is calculation of the real part of a single term of a sliding discrete Fourier transform (sliding DFT) of length $r$ at frequency $\omega = 2\pi m/r,\, m\in N.$ The complex IIR filter is:

$$\begin{align}y[n] &= x[n] - x[n-r]\\
&+e^{-i\omega}y[n-1]\end{align}\tag{8}$$

We get using Eqs. 6 and 7 the coefficients of the numerator and denominator polynomials of the transfer function of the real filter $(A_\text{real},\,B_\text{real})$ from those of the complex filter $(A,\,B).$ Assuming that $r > 1$ (coefficients not given are zero):

$$\begin{array}{rcllll}
k&=&0&1&r&r+1\\
\hline
B[k]&=&1&&-1\\
A[k]&=&1&-e^{-i\omega}\\
B_\text{real}[k]&=&\operatorname{Re}(1)&\operatorname{Re}(-e^{i\omega})&\operatorname{Re}(-1)&\operatorname{Re}(e^{i\omega})\\
&=&1&-\cos(\omega)&-1&\cos(\omega)\\
A_\text{real}[k]&=&1&-e^{i\omega}-e^{-i\omega}&-e^{-i\omega}(-e^{i\omega})\\
&=&1&-2\cos(\omega)&1\end{array}$$

The real filter is thus:

$$\begin{align}y[n] &= x[n] -\cos(\omega)x[n-1] -x[n-r] -\cos(\omega)x[n-r-1]\\
&+ 2\cos(\omega)y[n-1] - y[n-2]\end{align}$$

The recursive part can be recognized as the Goertzel filter.

In Octave:

m = 2;
r = 20;
w = 2*pi*m/r;
b = zeros(1, r+2);
b(1) = 1; b(2) = -cos(w); b(r+1) = -1; b(r+2) = cos(w);
a = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
plot([0:length-1], filter(b, a, impulse), "o");

Figure 2. Impulse response of a real IIR filter that calculates the real part of a single frequency bin of sliding DFT.

The impulse response can be zero phase aligned to the opposite end by multiplying the complex transfer function by $e^{-i\omega}:$

m = 2;
r = 20;
w = 2*pi*m/r;
b = zeros(1, r+2);
b(1) = cos(w); b(2) = -1; b(r+1) = -cos(w); b(r+2) = 1;
a = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
plot([0:length-1], filter(b, a, impulse), "o");

Figure 3. Impulse response of a real IIR filter that calculates the real part of a single frequency bin of sliding DFT, with zero phase aligned to the opposite end than in Fig. 2. This is a circular shift of the impulse response one sampling period to the left, and can sometimes be used to remove one sample of system latency.

In case of $m=0,$ the numerator and denominator polynomials will have an unnecessary common factor, but for other values of $m$ the approach works fine.


Steering Vector for Generic STAP I am interested in understanding the proper terminology for the topic of Space Time Adaptive Processing (STAP) that extends beyond the example in all the text books and papers I could find that reference the Uniform Line Array.  In the case of the ULA, the steering vector is easy for me to understand and is

\begin{equation*}
\mathbf{a}(\phi) = A[e^{i \theta_1}, ...,e^{i \theta_n}]^T
\end{equation*}

where $\theta_i$ is the phase response in the direction $\phi$.  The optimal matched filter is thus:

\begin{equation*}
\mathbf{h}(\phi) = \alpha\mathbf{R}_s^{-1} \mathbf{a}(\phi)
\end{equation*}

where $\alpha$ grabs all of the constant terms and lumps them together

However, my expectation for the more generic case is that the array response on all channels is not equal (more like an array manifold) and is thus:

\begin{equation*}
\mathbf{a}^\prime(\phi) = [A_1e^{i \theta_1}, ...,A_ne^{i \theta_n}]^T
\end{equation*}

and now $A_i$ is the amplitude response in the direction $\phi$.

This leads me to my question.  Is $\mathbf{a}^{\prime}(\phi)$ also a steering vector or is this something different (i.e. elements of the steering vector be other than unity)?  And if the steering vector is still $\mathbf{a}(\phi)$, what is the convention for treating the amplitude response of the array.  I am leaning towards the former.

Yes, $a'(\phi)$ is a steering vector. In fact, the third equation (almost) describes a phased array that is capable of steering lobe patterns at different frequencies.

Each element of the array has a filter behind it which forms a weighted sum of (delayed) samples in Time. The same idea is applied one level above by now summing (delayed) samples of the received signal from the different elements that are in different points in Space.

In STAP, the adaptation in Time maintains the coherence of the probe signal (or carrier) against noise / doppler and other "structural" distortions and the adaptation in Space maintains the coherence of the signal along a specific direction by trying to minimise the covariance between different elements. Think of it as one convolution per element of the array(typical filtering in the Time domain) and one more ACROSS the same time instant (the output) of the element filters (Space).

BUT! With this configuration you can tune the Time filters at some sinusoid and then with the frequency of that sinusoid tune the Space filter to receive from a given bearing. To do the same thing simultaneously, you would need a Fourier Transform behind each element to decompose the received signal into a bank of sinusoids and then apply different space filters ($a'(\phi)$) to those outputs for specific bearings. In this way, you could have the array listening to some low frequency from the left and some high frequency from the right simultaneously. I feel that this is what is meant with "...the more generic case is that the array response on all channels is not equal (more like an array manifold) and is thus..."

Hope this helps.


What does "warp" mean in DSP? I hear the word warp a lot, but cannot quite get a feel for what it really means
I suppose it's literal meaning would be to change in shape, but what does it mean with respect to DSP?
More specifically, I have a DSP chip manual that says a certain parameter adjusts warp for DACs and another one that adjusts warp for ADCs. What does warp mean in this context?
If possible, please explain it without technical jargon.

I've never really heard of "warp" regarding ADCs or DACs, but I imagine it's because the analog-to-digital (or vise versa) mapping function, which is supposed to be linear, a straight line, is curved slightly due to variation in parts.  I don't see how this would happen with sigma-delta ADCs or DACs (unless the summer at the input with the delta feedback is non-linear).  but with an R-2R DAC, if the R's aren't all perfectly the same, the digital-to-analog mapping function will not be a straight line.  For a flash ADC, if all of the resistors in the voltage divider line are not all perfectly equal, that analog-to-digital mapping function will also not be linear and would be "warped".

If, somehow, one can measure or determine the "warpness" of the mapping function, I suppose it can be compensated in the DSP software.


Why do musical instruments have different colour in phone? When I hear a musical instrument in phone, the color of it sounds a bit different than its original 'live' sound. Why?

Typical telephone codecs use a sampling frequency of 8kHz and a bandwidth of 7kHz.  That means that your signal will get cut off at frequencies higher than 3.5kHz.  That's ok for most speech but not sufficient for reproducing the timbre of most music instruments.

You'll likely fare reasonably well with instruments like recorders (which don't have a lot of overtones).  The worst combination is high-pitched and overtone-rich, like violin, oboe, soprano sax.  Also instruments with a full spectrum will sound different, like high hats, snares, gongs.


Detect abrupt changes in signal I'm using MATLAB and I'm trying to find the index of the first signal drop as shown in the figure bellow :

I tried the function findchangepts but it didn't give me the exact index, any solutions ?

A cheap and simple way is to differentiate the signal numerically and threshold the magnitude of differencate to identify the abrupt changes.While deciding the threshold you should take into account the dynamic range of the signal.