Frequency Translation in MATLAB I have to design a generic baseband filter. I have a C++ code which gives me filter coefficients based on the input parameters such as no. of taps, no. of bands, band edges, amplitudes and weights(code uses remez exchange algorithm).
I import the filter coefficients in MATLAB, do fft and plot the graph of it. I have to shift the graph such that its center is the center of baseband. I know that if I multiply filter coefficients with an exponential function and then do fft the resulting frequency domain output will be shifted in frequency $$x(t) e^{j 2\pi f_0} \implies X(f-f_0)$$. But I dont know how to do it practically. Do I just have to multiply my filter coefficient vector with an exponential like exp^(j*2*pi*f0) and then take fft and plot? I tried doing that but the results of both the graphs are exactly the same.

You are partly right. However It's not the filter coefficinets $a$,$b$, but the filter impulse response $h[n]$ which you have to multiply with the exponential term $e^{j w_0 n}$ to shift its frequency response to a desired center of frequency $w_0$. Note that for FIR filters the coefficients are such that $a=1$ and $b[n]=h[n]$; i.e., the impulse response equals the coefficients $b[k]$, but that's not so with the IIR filters.

Assume that you have an FIR filter with coefficients $a=1$ and $b[k]$ so that its impulse response is: $$h[n] = b[n]$$ And let the frequency response of the filter be $H(e^{j \omega})$, then the frequency response of the new filter $h_+[n] = e^{j w_0 n} h[n]$ will be
$$ h_+[n] = e^{j w_0 n} h[n] \implies H_+(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) $$

and similarly the frequency response of the new filter $h_-[n] = e^{-j w_0 n} h[n]$ will be
$$ h_-[n] = e^{-j \omega_0 n} h[n] \implies H_-(e^{j \omega}) = H_(e^{j( \omega + \omega_0)}) $$

Note if you want a real impulse response filter, say $h_0[n]$, then you should add those two sub shifted filters:
$$ h_0[n] = h_+[n] + h_-[n] \implies H_0(e^{j \omega}) = H_+(e^{j \omega}) + H_-(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) + H(e^{j (\omega + \omega_0)}) $$

You can achieve the same result by the following multiplication then;
$$ h_0[n] = 2\cos(\omega_0 n) h[n] \implies H_0(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) + H(e^{j (\omega + \omega_0)}) $$

MATLAB implmenetation would be as follows:

h = fir1(64, 0.1 );     % prototype lowpass filter
L = length(h);

w0 = 0.5*pi;            % bandpass filter new frequency
n = [0:L-1];            % discrete-time index vector

h0 = 2*cos(w0*n).*h;    % new bandpass filter

% Display prototype lowpass FIR filter
figure,subplot(2,1,1)
stem(n,h);title('Prototype Lowpass filter h[n]');
subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('DTFT magnitude |H(e^{j \omega})| of lowpass filter h[n]');

% Display new bandpass FIR filter
figure,subplot(2,1,1)
stem(n,h0);title('Prototype Lowpass filter h_0[n]');
subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h0,1024))));
title('DTFT magnitude |H_0(e^{j \omega})| of bandpass filter h_0[n]');


Determining the features of the DSP I need to sample a signal of 40MHz with a sampling frequency of 120MHz at least with a 8/12 bits ADC and send them through the USB serial port of a PC. That means the serial communication bandwith should be at least 960/1440 megabits per second.

What kind of DSP should I use to satisfy this bit-rate, and what signal-processing-methods should I learn to optimize this bit-stream to use the same data simultaneously with a mainstream PC. (The bitrate seemed a bit huge for me.)

I have no experience with DSPs before. I was seeking for answers for a while but not even my teachers, or their assistants could help me for a little long period of time.

I just need a beacon to find the answers of the questions, thanks in advance.

Edit: I wrote 40Hz instead 40MHz, I am really sorry.

This doesn't sound like you need a DSP at all.

You'll need something that does USB3 towards your computer. Typical solution to that is the fx3 chip from cypress. I'm not aware of DSPs that integrate a USB 3 phy or even just the logic to speak to one - not a typical interconnect in their world.

Then you need something that orchestrates the ADC, keeps timing correct, and hands over the data to the usb3 chip, but doesn't do much of actual signal processing.

The whole "glueing one high speed digital thing to another" aspect screams for an FPGA. So, what I'd recommend: pick an FX3 development board and design an FPGA/ADC add-on board to that.

May I also interest you in other interconnects? Maybe gigabit Ethernet is sufficient for your use case; in that case, you could basically use the Ettus N200 as reference design (schematics are public). Or buy one - it comes with a dual-100MS/s ADC, and with the BasicRX daughterboard, you can directly interface with that ADC in baseband.

Full disclosure: I'm affiliated with Ettus, but to be honest, I'd like to see you succeed in building your DAQ device :)


Is interpolation of an audio signal to increase frequency resolution possible? I apologize if some of what I ask is not entirely correct, I'm new to this field, but extremely interested.

I have an Audio signal of sample rate 44.1 kHz that I want to segment into 30 frames, and get the DFT of each frame to find the magnitude of certain frequencies for that frame. However, this means that I have a frequency resolution of 30 Hz in each bin, which isn't narrow enough.

Is it possible to interpolate the data to attain more data samples? As far as I'm aware, doubling the number of points would give a 88.2 kHz sampling rate, but still give a frequency resolution of 30. Would it be possible to treat the interpolated data as still having a sample rate of 44.1 kHz?

I do not understand why you have 30Hz resolution so I will focus only to the principle of the question "does interpolation increase resolution?".

Short answer is no, no new data, no new information.

A longer answer needs the spectrum visualization below, with time domain, continuous frequency domain and discrete frequency domain from left to right.

The interpolation technique here is to preserve the information of spectrum. DFT works on discrete frequency domain which is the part from $-f_s$ to $f_s$ of the continuous frequency version.

First, look at the continuous frequency domain, if you upsample your signal correctly, it is equivalent to changing the sampling frequency. You wish to double the number of data, but it is just removing one-half the spectrum replicas of sampling process.

Now, look at the discrete frequency one. This version is normalized from $-f_s$ to $f_s$ of the continuous frequency counterpart. $f_s$ is doubled, the spectrum is then shrinked by a factor 1/2. If we call $0 < \alpha < 1$ the proportion of non-zero frequencies before upsampling, this proportion is $\alpha/2$ after upsampling. Before upsampling, DFT gives you $N$ bins for $\alpha$ then $N\alpha$ bins for the spectrum; after upsampling it is $2N$ for $\alpha/2$ then always $2N \times \alpha/2 = N\alpha$ bins for the same spectrum. No, your resolution does not change at all.

To have smoother resolution, the only way is to add more data. In your example, instead of dividing to 30 frames, divide your audio file to 15 frames.


What is the difference between object detection and target detection? Are they the same thing? I was reading some articles regarding object detection and target detection and have a hard time understanding the difference between object detection and target detection. i want to know whether they are different from each other? and if so, how?

"Target" is a military term; historically, that's why it's the word used for "identifiable thing from a radar system".
So, I wouldn't use "target detection" outside of radar systems, whereas the word "object detection" is generally used in image processing (for anything that isn't a person). I might be a bit alone with this, but if I look at an image and someone asks me to perform "target detection" on that, I'd assume they want to shoot at that.
I think "target" in radar has a more "physical meaning", as in, it's actually a single contiguous physical object which emits/reflects radio-frequency emissions that should be recognized / tracked as such.
"Object" can be anything, from a round ball, to scattered remnants of a white blood cell; it's slightly less "physical", more "appearance-based".  (Of course, reflecting light is still reflecting EM waves; the difference lies in the number of sensor elements in a camera vs radar array.)


DAPSK bit assignment for the amplitude difference I am studying DAPSK(Differential amplitude phase shift keying) at the moment. I have a 64-DAPSK. We are using 4 concentric rings and 6 bits per symbol. The first two bits contain the information of amplitude difference and the rest contain the information of the phase difference. I am having trouble understanding how the table on the bottom left of this example is filled out:

[]

Each circle is identified by its amplitude, either 1, $a$, $a^2$, or $a^3$. The table shows the "amplitude transition factor"; that is, you need to multiply the previous amplitude by this factor to find the next amplitude. In other words, it is not the amplitude difference that is encoded, but the amplitude ratio.

Let's follow the example above the table. You start with $d_{n-1}$ with amplitude $a^3$. The next symbol, $d_n$, has amplitude 1, so the factor is $1/a^3$, which corresponds to bits $01$ (indicated also as the first two bits under "info bits" in the figure, top right).

Now, $d_n$ has amplitude 1 and the next symbol, $d_{n+1}$, has amplitude $a^2$. The factor is $a^2$ and the corresponding bit sequence is $11$.

Notice how the table is designed to allow transitions from any of the four amplitudes to any other amplitude.


Periodic signals in Continuous and discrete time Is there any signal which is periodic in Continuous Time but not in Discrete Time?
I have this doubt prevailing in me for a long time.
Are all CT periodic signals periodic in DT?
If so, how is sin wave a periodic signal?
It is periodic only at intervals of pi radians, which is an irrational number and not an integer.
But for a DT signal to be periodic shouldn't the signal be periodic at discrete intervals of time?

A periodic continuous-time signal satisfies $x(t)=x(t+T_0)$ for all $t$. The period $T_0$ doesn't need to be a rational number. A periodic discrete-time signal satisfies $x[n]=x[n+N]$ for all integers $n$. The period $N$ is an integer.
If you sample a periodic continuous-time signal, you don't necessarily get a periodic sequence. E.g., sampling the periodic signal $x(t)=\sin(\omega_0t)$ results in the sequence $x_d[n]=x(nT)=\sin(\omega_0Tn)$, which is only periodic if $\omega_0T=2\pi k/l$ with integers $k$ and $l$.
In general, sampling a $T_0$-periodic continuous-time function with sampling interval $T$ results in a periodic sequence only if
$$kT_0=lT,\qquad k,l\in\mathbb{Z}$$
is satisfied.


Frequency Domain Interpolation: Convolution with Sinc Function I am reading a paper, and I came across the following definition of sinc interpolation.

Warning. I don't have a strong background in signal processing. Also, I have no clue what that bar on $\bar{F}$ means. I don't know why that would be the conjugate. Would it actually be the conjugate?
Context. In this part, for a given kernel $f_k$ of size $k$, we want to expand $F_k = DFT\{f_k\}$ to  $F_N = DFT\{f_N\}$ have size $N$, without going to the spatial domain, padding and back to the spectral domain.
Let's discuss on a 1D scenario for notation simplification.
$$
G(u) = F(u) \ast
\left [
e^{\frac{-j 2 \pi u}{N}(\frac{k-1}{2})}
\frac{sin(\frac{\pi u k}{N})}{sin(\frac{\pi u}{N})}
\right ]
$$

Q1. I am assuming that $G(u)$ and $F(u)$ are discrete, having the respective sizes of $N$ and $k$. Is this assumption right?
Q2. Would that second term actually be the $sinc$ kernel, i.e. would the expression below be true? If not, what is their between that term and the sinc function?
$$
sinc(u) = \frac{sin(u)}{u} = \left [
e^{\frac{-j 2 \pi u}{N}(\frac{k-1}{2})}
\frac{sin(\frac{\pi u k}{N})}{sin(\frac{\pi u}{N})}
\right ]
$$

Another way I know to avoid using the DFT to compute $F_N$, is by representing $f_N$ in terms of $f_k$ as
$$f_n = \sum_{i=0}^{k-1} \delta(n-i) f_i$$
Thus, having that
$$DFT\{\delta[n]\} = 1$$
and
$$ DFT\{x[n-n_o]\} = e^{-j \frac{2 \pi}{N} k n_o} DFT\{x[n]\}$$
Hence
$$F_N(k) = \sum_{n_o=0}^{k-1} e^{-j \frac{2 \pi}{N} k n_o } f_{n_o}$$
I am not sure if that the convolution on the frequency domain $G(u)$,
and the impulse based composition $F_N(k)$ are similar approaches of doing the same thing. I did notice the similarity between
$$
e^{\frac{-j 2 \pi u}{N}(\frac{k-1}{2})}
\sim
e^{\frac{-j 2 \pi k}{N}n_o}
$$
But, if this is the case, then,
$$
\frac{sin(\frac{\pi u k}{N})}{sin(\frac{\pi u}{N})} \sim F\{\delta(k-n_o)\}
$$

Q3. Are these the same approaches? Is there any relation between these expressions?
$$\frac{k-1}{2} \sim n_o$$

The "discrete sinc kernel" is also known as a Dirichlet Kernel.  It looks like a you got a 2D (u and v) going there.  I'm not that familiar with that, but we just did a pretty good number on the 1D case.

This post explains the inverse DFT can be interpreted as "epicycle drawing":

How to get Fourier coefficients to draw any shape using DFT?

This post then follows up with the relationship between the sampling, epicycle drawing, upsampling, and the normalized sinc function:

Absolute convergence of periodic sinc interpolation

The conclusions are:

Suppose a continuous band limited function (with band limit $L$) is sampled, then reconstructed from the samples.

1) If L < N/2 the function can be unambiguously reconstructed.  The results from an infinite series sinc reconstruction (Whittaker-Shannon) and taking a DFT, zero padding at the Nyquist and then the inverse DFT to upsample, will produce the same path.

2) If L = N/2, the Nyquist bin must be split between N/2 and -N/2 for the sinc reconstruction to match the DFT reconstruction.  It may or not may not match the original function, but it will match at all the sample points.

3) If L > N/2, and N is even, the Nyquist bin must still be split for the sinc reconstruction to match the DFT reconstruction, but neither will match the original function, yet both will still match at all the sample points.

The relationship between the sinc function and the Direchlet Kernel is this:

1) The sinc function is the limit of the Dirichlet kernel as the sample count goes to infinity

2) For odd N, the Dirichlet kernel is an infinite sum of sinc functions.  For even N, it is an adjusted one.

See the posts for details and discussion.

What I forgot to overtly mention is that the Dirichlet Kernel formula is directly derived from, and equivalent to, doing the FFT, zero pad at Nyquist (split Nyquist on even N) then inverse DFT for resampling or coefficients for Fourier Series.  Your link in the comment is about the reverse case, but since the forward DFT and inverse DFT are mirror image twins (especially when $1/\sqrt{N}$ normalized), so what gets said in one direction has to be true in the other.  This means if you zero pad a signal, and have an even number of points, one of the outside points needs to be halved and the other half appended at the other end before zero padding outward, leaving your data centered.

I already explained the relationship between the Sinc and the Dirichlet, just above.  The math can be found in the second link in the section of my answer that has $y_m = Y(w)$ in it.  Towards the bottom.


When can the $\mathcal Z$-transform be inverted? When not? What are the conditions that must be satisfied to be able to invert the $\mathcal Z$-transform?

If I understand your question correctly then you're asking under which conditions the inverse $\mathcal{Z}$-transform of a given function $F(z)$ exists. Since the inverse $\mathcal{Z}$-transform is just a Laurent expansion of $F(z)$ around $z=0$, the inverse transform exists if $F(z)$ can be developed in a Laurent series about $z=0$. This is the case if $F(z)$ is analytic in an annulus centered at $z=0$.


Inverse Laplace transform Using Inversion Formula Use the complex inversion formula to calculate the inverse Laplace transform $f(t)$ of the following Laplace transform:
  $$F_L (s) = \frac{1}{(s+2)(s^2 +4)}.$$
  When the region of convergence is:
  \begin{align}(1)& \quad Re(s)<-2;\\(2)&\quad -2<Re(s)<0;\\(3)&\quad Re(s)>0.\end{align}

Attempt:

Here is an explanation of the complex inversion formula. Plugging the function into the formula:

$$f(t) = \frac{1}{j2\pi}  \int^{\sigma + j \infty}_{\sigma - j \infty} \frac{e^{st}}{(s+2)(s^2+4)}\ ds \tag{1}$$

So, how do I need to choose $\sigma$?
And how do I evaluate this for each of the three regions?

P. S. I tried to solve this without the complex inversion formula, just to see what the answer should look like. I started out by expanding using partial fractions as:
\begin{align}
\frac{1}{(s+2)(s^2 +4)}&= \frac{1}{8} \left( \frac{1}{s+2} + \frac{1}{s^2+4} \right)\\
&=\frac{1}{8} \left( \frac{1}{s+2} + \frac{1}{(s+2j)(s-2j)} \right)\\
&=\frac{1}{8} \left( \frac{1}{s+2} + \frac{j}{4(s+2j)} - \frac{j}{4(s-2j)} \right).
\end{align}

Looking at a Laplace transform table, $\frac{1}{s-a} \leftrightarrow e^{at},$ so

$$f(t) = \frac{1}{8} \left( e^{-2t} + \frac{j}{4} \left(e^{-2jt} + e^{2jt}\right) \right).$$

Is this correct?
If so, how can I get to this using the complex inversion formula?

In engineering practice, the complex inversion integral is hardly ever used. As an engineer, you will almost exclusively need to invert rational functions, and this can be done by partial fraction expansion and elementary inversions. So first I'll show you how to obtain the inverse Laplace transform by partial fraction expansion, then I'll explain the evaluation of the inversion integral using Cauchy's residue theorem.

You have an error in the partial fraction expansion. Furthermore, you don't need to split up the complex pole pair. I would rewrite the Laplace transform like this:

$$F(s)=\frac{1}{(s+2)(s^2+4)}=\frac{A}{s+2}+\frac{Bs+C}{s^2+4}\tag{1}$$

with $A=\frac18$, $B=-\frac18$, and $C=\frac14$. The terms on the right-hand side of $(1)$ are elementary Laplace transforms. Now you just have to consider the different regions of convergence (ROC):

$$\begin{align}\frac{1}{s+2}&\Longleftrightarrow e^{-2t}u(t),&\quad \text{Re}\{s\}>-2\\\frac{1}{s+2}&\Longleftrightarrow -e^{-2t}u(-t),&\quad \text{Re}\{s\}<-2\\
\frac{s}{s^2+4}&\Longleftrightarrow \cos(2t)u(t),&\quad\text{Re}\{s\}>0\\
\frac{s}{s^2+4}&\Longleftrightarrow -\cos(2t)u(-t),&\quad\text{Re}\{s\}<0\\
\frac{1}{s^2+4}&\Longleftrightarrow \frac12\sin(2t)u(t),&\quad\text{Re}\{s\}>0\\
\frac{1}{s^2+4}&\Longleftrightarrow -\frac12\sin(2t)u(-t),&\quad\text{Re}\{s\}<0
\end{align}$$

So for the ROC $\text{Re}\{s\}<-2$ you get the anti-causal signal
$$f(t)=\frac18\left[-e^{-2t}+\cos(2t)-\sin(2t)\right]u(-t)\tag{2}$$

For the ROC $-2<\text{Re}\{s\}<0$ you get the two-sided signal

$$f(t)=\frac18\left[e^{-2t}u(t)+(\cos(2t)-\sin(2t))u(-t)\right]\tag{3}$$

And, finally, for the ROC $\text{Re}\{s\}>0$ you get the causal signal
$$f(t)=\frac18\left[e^{-2t}-\cos(2t)+\sin(2t)\right]u(t)\tag{4}$$

If you need to use the inversion formula then it is very helpful to know Cauchy's residue theorem, which says that

$$\frac{1}{2\pi j}\oint_Cf(s)ds=\sum_kR_k\tag{5}$$

where $f(s)$ is analytic with finitely many poles, $C$ is a positively oriented closed curve, and $R_k$ are the residues of the poles inside $C$. It can be shown that the inversion integral equals a contour integral if the curve $C$ is chosen appropriately:

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=
\frac{1}{2\pi j}\oint_CF(s)e^{st}ds\tag{6}$$

In the case of a rational function $F(s)$ the curve $C$ is chosen as a Bromwich contour, as shown here in Fig.2. The straight line part of the curve is the actual integration path we're interested in. The contribution from the circular part of $C$ approaches zero. Depending on the chosen ROC, we have to choose the position of the straight line (i.e., the value of $\sigma)$ differently. For ROC $\text{Re}\{s\}<-2$ (i.e., the anti-causal solution), the straight line is anywhere to the left of the left-most pole, and so the Bromwich contour enclosing all poles is negatively oriented, which results in a sign change:

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=-\sum_kR_k,\quad\sigma<-2,\quad t<0\tag{7}$$

where $R_k$ are the residues corresponding to the poles of the function $f(s)=F(s)e^{st}$.
For ROC $\text{Re}\{s\}>0$ (i.e., the causal solution), the straight line is anywhere to the right of the right-most pole, and the Bromwich contour enclosing all poles is positively oriented. Consequently, we have

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=\sum_kR_k,\quad\sigma>0,\quad t>0\tag{8}$$

For the two-side solution with ROC $-2<\text{Re}\{s\}<0$ we need to choose two curves with the straight line inside the ROC. One encloses the pole to its left at $s=-2$, so the curve is positively oriented, and the other one encloses the two poles at $s=2j$ and $s=-2j$ to the right of the straight line, so it is negatively oriented (which adds a negative sign to the corresponding residues). Let $R_1$ be the residue corresponding to the pole at $s=-2$, and let $R_2$ and $R_3$ be the residues of the two poles at $\pm 2j$, respectively. The inversion integral is then given by:

$$\frac{1}{2\pi j}\int_{\sigma-j\infty}^{\sigma+j\infty}F(s)e^{st}ds=\begin{cases}R_1,&t>0\\-R_2-R_3,&t<0\end{cases},\quad-2<\sigma<0\tag{9}$$

What remains is the computation of the residues. The residue at pole $p_k$ is given by
$$R_k=\lim_{s\rightarrow p_k}(s-p_k)f(s)\tag{10}$$

With $f(s)=F(s)e^{st}$ we get for $p_1=-2$

$$R_1=\lim_{s\rightarrow -2}(s+2)F(s)e^{st}=\lim_{s\rightarrow -2}\frac{e^{st}}{s^2+4}=\frac{e^{-2t}}{8}\tag{11}$$

In a similar manner the other residues are obtained:

$$\begin{align}R_2&=-\frac{e^{2jt}}{8}\frac{1+j}{2}\\
R_3&=-\frac{e^{-2jt}}{8}\frac{1-j}{2}\end{align}\tag{12}$$

Now $(7)-(9)$ can be evaluated, and the results are of course the same as $(2)-(4)$.


How to prove Line Spectral Pair coefficient properties? I am trying to understand the definition of LSP coefficients from http://www.dspcsp.com/pdf/lsp.pdf.
We take the polynomial $a(x)$ of degree $M$ from the denominator of the LPC system function, and represent it as a sum of a palindromic and an anti-palindromic polynomial:

$a(x) = 0.5(p(x) + q(x))$,
where

$p(x) = a(x) + x^{M+1}a(x^{-1})$

$q(x) = a(x) - x^{M+1}a(x^{-1})$

Then, the textbook says,

  it is not hard to show that if all the zeros of a(x) are inside the
  unit circle, then all the zeros of $p(x)$ and of $q(x)$ are on the
  unit circle. Furthermore, the zeros of $p(x)$ and $q(x)$ are
  intertwined, i.e., between every two zeros of $p(x)$ there is a zero
  of $q(x)$ and vice versa.

How can this be proved?

UPD.

I found an article by Soong and Juang, "Line spectrum pair and speech data compression", with a proof. It depends on the equality
$ |z_iz - 1|^2 - |z-z_i|^2 = (1-|z|^2)(1-|z_i|^2) $,
where neither $z_i$ nor $z$ are assumed to be real.
But then the equality does not hold:
$ |z_iz - 1|^2 - |z-z_i|^2 = (z_iz - 1)(\overline{z_iz} - 1) - (z - z_i)(\overline{z} - \overline{z_i}) = |z_i|^2|z|^2 + 1 - zz_i -\overline{zz_i}  - |z|^2 - |z_i|^2 + \overline{z}z_i + z\overline{z_i} = (1-|z|^2)(1-|z_i|^2) - (z - \overline{z})(z_i - \overline{z_i})$

Ok, so the plan of the proof is as follows:
First, I am going to fix the Soong and Juang's proof that if all roots of $A$ lie inside the unit circle, then all roots of $P$ and $Q$ lie on the unit circle. (I am going to follow their notation, too.)
Then, I am going to prove that all roots of the LPC polynomial are indeed inside the unit circle.
Finally, I'll prove that the roots of $P$ and $Q$ are intertwined.
The first part is easy. We define
$$
H(z) \triangleq z^{-(M+1)}\frac{A(z^{-1})}{A(z)}
$$
Since A is a polynomial with real coefficients, complex conjugation is a permutation on the set of its roots, therefore H can be factored as
$$
H(z) = z^{-1} \prod_{i=1}^M {\frac{z_iz-1}{z-\overline{z_i}}},
$$
where $z_i$ are the roots of $A_M(z)$.
Now, since
$ |z_iz - 1|^2 - |z-\overline{z_i}|^2 = (z_iz - 1)(\overline{z_iz} - 1) - (z - \overline{z_i})(\overline{z} - z_i) = |z_i|^2|z|^2 + 1 - zz_i -\overline{zz_i}  - |z|^2 - |z_i|^2 + zz_i + \overline{z}\overline{z_i} = (1-|z|^2)(1-|z_i|^2)$,
Therefore,
$$
|H(z)| \text{ is } \begin{cases}
       >1 &\quad\text{if } |z| < 1\\
       <1 &\quad\text{if } |z| > 1\\
       =1 &\quad\text{if } |z| = 1\\
     \end{cases}
$$
Because the solution to P(z) = 0 or Q(z) = 0 requires that H(z) = ± 1, we conclude that P(z) and Q(z) can only have zeros on the unit circle.
OK, now how do we prove that all zeros of $A(z)$ lie inside the unit circle? We are going to prove that the polynomials $A_M(z)$ satisfy a recursive relation:
$$
A_{M+1}(z) = A_M(z) - \kappa_M z^{-(M+1)}A_M(z^{-1}),
$$
where $|\kappa_M| < 1$. ($\kappa_M$ are called reflection coefficients.) Then we'll prove recursively that all roots of $A_M(z)$ are in the unit circle.
$A_M(z) = 1 + a_{M,1} z^{-1} + \ldots + a_{M,M} z^{-M}$ is a polynomial whose coefficients minimize the forward prediction error:
$$
E_{f,M} = \sum_k{\left(x(k) + \sum_{m=1}^M {x(k-m) a_{M,m}}\right)^2}
$$
If we write down the optimality condition for $a_{M,m}$, we get the Wiener-Hopf equations:
$$
R_M \vec{a_M} = - \vec{r_{f,M}},
$$
where
$$
\vec{a_M} \triangleq [a_1, a_2, \ldots, a_M]^T,\\
\vec{r_{f,M}} \triangleq [r(1), r(2), \ldots, r(M)]^T,\\
r(l) \triangleq \sum_k{x(k) x(k-l)}, \\
R_M \triangleq \begin{pmatrix}
r(0) & r(1) & \ldots & r(M-1) \\
r(1) & r(0) & \ldots & r(M-2) \\
\vdots & \vdots & \ddots & \vdots \\
r(M-1) & r(M-2) & \ldots & r(0) \\
\end{pmatrix}
$$
Similarly, the solution to the problem of backwards linear prediction
$$
\text{minimize  }
E_{b,M} = \sum_k{\left(x(k-M) + \sum_{m=1}^M {x(k-m+1) b_{M,m}}\right)^2}
$$
must satisfy the equation
$$
R_M \vec{b_M} = - \vec{r_{b,M}},
$$
where
$$
\vec{r_{b,M}} \triangleq [r(M), r(M-1), \ldots, r(1)]^T
$$
Define the coidentity matrix $J_M$ as
$$
J_M = \begin{pmatrix}
0 & 0 & \ldots & 0 & 1 \\
0 & 0 & \ldots & 1 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 1 & \ldots & 0 & 0 \\
1 & 0 & \ldots & 0 & 0 \\
\end{pmatrix}
$$
It is easy to see that
$$
J_MR_M = R_MJ_M,
$$
therefore
$$
R_M \vec{a_M} = -\vec{r_{f,M}} = -J_M\vec{r_{b,M}} = J_M R_M \vec{b_M} =
R_M J_M \vec{b_M}.
$$
Since, as I have shown, the matrix $R_M$ is nonsingular, it follows that $\vec{a_M} = J_M \vec{b_M}$.
Furthermore,
$$
E_{f,M} = r(0) + \vec{r_{f,M}^T} \vec{a_M} = r(0) + \vec{r_{f,M}^T}J_M J_M \vec{a_M} = r(0) + \vec{r_{b,M}^T} \vec{b_M} = E_{b,M} = E_M
$$
Consider the following expression:
$$
\begin{pmatrix}
R_M & \vec{r_{b,M}} \\
\vec{r_{b,M}}^T & r(0) \\
\end{pmatrix}
\begin{pmatrix}
1 \\
\vec{a_{M-1}} \\
0\\
\end{pmatrix} =
\begin{pmatrix}
E_{M-1}\\
0_{(M-1)\times 1} \\
K_M\\
\end{pmatrix}
$$
where
$$
K_M = r(M) + \vec{r_{b,M-1}^T}\vec{a_{M-1}} = r(M) + \vec{r_{b,M-1}^T}J_M\vec{b_{M-1}} = r(M) + \vec{r_{f,M-1}^T}\vec{b_{M-1}}
$$
We define the reflection coefficient as
$$
\kappa_M = \frac{K_M}{E_{M-1}}
$$
From backward linear prediction, we have:
$$
\begin{pmatrix}
r(0) & \vec{r_{f,M}^T} \\
\vec{r_{f,M}} & R_M \\
\end{pmatrix}
\begin{pmatrix}
0 \\
\vec{b_{M-1}} \\
1\\
\end{pmatrix} =
\begin{pmatrix}
K_M\\
0_{(M-1)\times 1} \\
E_{M-1}\\
\end{pmatrix}
$$
Multiplying both sides of this equation by $\kappa_M$ and subtracting it from the previous matrix equation, we get:
$$
R_{M+1}\begin{pmatrix}
1\\
\vec{a_{M-1}} - \kappa_M \vec{b_{M-1}} \\
-\kappa_M\\
\end{pmatrix} =
\begin{pmatrix}
E_{M-1} (1 - \kappa_M^2) \\
0_{M \times 1} \\
\end{pmatrix}
$$
From which it follows that
$$
\vec{a_M} = \begin{pmatrix}
\vec{a_{M-1}} - \kappa_M \vec{b_{M-1}} \\
-\kappa_M\\
\end{pmatrix}
$$
and $E_M = E_{M-1}(1-\kappa_M^2)$. Since $E_M > 0$ and $E_{M-1} > 0$, we get $|\kappa_M| < 1$.
Now let us prove by induction that all roots of $A_M$ lie inside the unit circle.
Base: For $A_1$, $\kappa_1$ is the only root.
Induction step: Suppose we proved that all roots of $A_M$ lie in the unit circle.
$$
A_{M+1}(z) = A_M(z) (1 - \kappa_M H(z)),
$$
So, if $A_{M+1}(z) = 0$, then either $A_M(z) = 0$ or $|H(z)| = 1/\kappa_M > 1$, therefore $|z| < 1$.
The last thing to prove: zeros of P and Q are intertwined. Zeros of P are solutions to $H(z) = -1$ and zeros of Q are solutions to $H(z) = 1$, so it is sufficient to prove that the phase of H monotonously decreases as $z = e^{jw}$ travels the unit circle.
H(z) can be factored as
$$
H(z) = z^{-(M+1)}\prod_{i=1}^M {\frac{1 - \overline{z_i^{-1}}z^{-1}}{1 - z_iz^{-1}}}
$$
We are going to prove that the phase of each multiplier is a monotonously decreasing function of w. An easy way to do this is by using the law of tangents.

$$
\frac{\tan{(\beta - \alpha)/2}}{\tan{(\beta + \alpha)/2}} = \frac{1/r - 1}{1/r + 1}
$$
$$
\angle \frac{1 - \overline{z_i^{-1}}z^{-1}}{1 - z_iz^{-1}} = - 2 \arctan{(\beta - \alpha)/2} = - 2 \arctan{\left(\frac{1/r - 1}{1/r + 1} \tan (90° - (w_i - w)/2)\right)}
$$
is a monotonously decreasing function of w, q.e.d.


What is the scaling function and wavelet function at wavelet analysis? I'm trying to looking the meaning and functionality about scaling function and wavelet function at wavelet analysis. I have googling already. But I can't find and understand the meaning.

What does those affect to analysis?
Does anyone help me to understand roughly?

Wavelets are functions which are form by two resulting coefficients. The detail coefficient and the scale coefficient. From wikipedia Wavelets are defined by the wavelet function ψ(t) (i.e. the mother wavelet) and scaling function φ(t) (also called father wavelet) in the time domain.

The wavelet function is in effect a band-pass filter and scaling it for each level halves its bandwidth.

The scaling and detail basically divide the signal into two applying a high-pass filter resulting into the detail coefficients - (which is the highest level of the transform) and a low-pass filter which results in the scaling coefficients - (which is the lowest level of the transform). This coefficients are calculated for each frequency of the discrete wavelet.

used to for example to extract the deterministic part of a signal and the stockastic part of the signal. They are also used through different procedures to denoise via thresholding the coefficients.


What is the correct gain of an RRC Filter? Breaking my brain all morning with this reading previous questions and googling ...
I have made an RRC filter from the equation on wikipedia. It works fine and I compared it to commpy library in python.
But what I dont understand is the frequency amplitudes.

Should the RRC filter have normalised gain, i.e. 0 dB by dividing the
impulse response by sum(impulse response)?
Or should I be dividing
the FFT result by fftsize?
Or should I infact be doing neither of
these two?

I browse the internet and many images have lots of different gains, none of which I can recreate ... but I can recreate the impulse and general frequency shape. The frequency response is just not shifted vertically correct (in dB scale) to show same gain as other images. What am I missing :(
#Define all parameters needed for RRC
alpha = 0.35
Fs = Fs
span = 10
sps = 10
num_taps = span*sps
Rs = Fs / sps

#taps must be odd for impulse to be symmetrical on time zero
if not num_taps % 2: num_taps += 1

#time vector
t = np.arange(-((num_taps-1)/2) , (((num_taps-1)/2)+1) , 1) * (1/Fs)

#Root Raised Cosine and deal with special cases
h = np.zeros(len(t))
for element in enumerate(t):
    i,x = element #i holds array index, x holds value of t
    if x == 0:
        h[i] =  (1 - alpha + (4 * alpha / np.pi))

    elif (alpha != 0 and x == (1 / Rs) / (4 * alpha)) or (alpha != 0 and x == -(1 / Rs) / (4 * alpha)) :
        h[i] = alpha / np.sqrt(2) * ( ((1 + 2 / np.pi) * np.sin(np.pi / (4 * alpha))) + ((1 - 2 / np.pi) * np.cos(np.pi / (4 * alpha))))

    else: #not special cases
        h[i] = (np.sin(np.pi * x * Rs * (1 - alpha) ) + \
                4 * alpha * x * Rs  * np.cos(np.pi * x * Rs * (1 + alpha))) / \
                (np.pi * x * Rs * (1 - (4 * alpha * x *Rs) ** 2 ))

h = h # / np.sum(h) #Should I divide to normalise gain?

#Fourier Transform
fftsize = 4096
f = np.linspace(-Fs / 2, Fs / 2, fftsize)

H = np.fft.fft(h,fftsize)  #Should I divide by fftsize here?!
H = 10*np.log10( np.abs(H)**2)
H = np.fft.fftshift(H)

#Make the plots
fig = plt.figure()
ax = fig.add_subplot(211)
ax.title.set_text('Raised Cosine Impulse Response')
ax.plot(t,h)
ax.grid(True)

ax = fig.add_subplot(212)
ax.title.set_text('Raised Cosine Frequency Response')
ax.plot(f,H)
ax.grid(True)

The gain is completely arbitrary and you can scale it as desired for the overall receiver or transmitter design. Where special attention must be paid is with fixed point design where the best practice is to let the filters grow the signal- do not scale the coefficients or the input as that only introduces more quantization noise and degrades SNR. Let the filter grow (extended precision accumulators!) and then scale the output as you would do for any gain adjust to bring the signal level to within the range of the precision used.
Typical design practice for fixed point design which avoids precision errors is to use 2 more bits of quantization for the coefficients, sum the weighted taps in an extended precision accumulator and then scale the output to any level desired. In this case the quantization of the coefficients (if the same or less than the datapath) would cause the implementation to not meet the expected rejection and passband ripple given by the floating point equivalent design. Typically we want rejection to be greater than the quantization noise floor, and a reasonable rule-of-thumb is that the achievable rejection will be 5-6 dB/bit of coefficient quantization.
I detail this further at this post showing the quantization error contributions (such that you can also see if you want to limit the accumulator precision that you can scale after so many taps since the trade between quantization noise and accumulator precision is clearer, but in general the safest/easiest approach without further computation is to just wait until the out and let the accumulator be $log_2(N)+2$ higher in precision than the datapath where N is the number of taps in the filter.
Inter-filter bit width
This is further detailed here within the context of a general approach to FIR filter design:
what is the suitable design Method to the filter?
The graphics below demonstrate the quantization effects on a typical low pass FIR filter. This motivates the "best-practice" to use 2 more bits of quantization for the coefficients over what is used in the datapath, and the noise summation process in the filter motivates using extended precision accumulators and let the filter grow the signal (then truncate at the end) as detailed in the links above:


Sliding DFT: I think my results are accurate, but I can't buy a swap-swap IDFT I've got a sliding DFT implementation that appears to be working (judging from an output plot). I would like to be able to invert this implementation using the standard tricks of swapping the real and imaginary components or conjugating my data. However, when I attempt this my output looks nothing like the original. I also end up with a strong imaginary component where that was not expected. (I'm not entirely sure on the proper range of output for each pass through this filter.) Can you see anything obviously wrong with this code? Thanks for your time.

    [TestMethod]
    public void TestInverseDft()
    {
        var rand = new Random(42);
        var dataSize = 50000;
        var data = new List<Complex>(dataSize);
        for(int i = 0; i < dataSize; i++)
            data.Add(new Complex(rand.NextDouble() * 2.0 - 1.0, 0.0));

        var clone = new List<Complex>(data);

        var n = SlidingDFT(data, 44100, 12, 55, 13);

        Swap(data);
        SlidingDFT(data, 44100, 12, 55, 13);
        Swap(data);

        for(int i = 15000; i < dataSize - 15000; i++)
            Assert.AreEqual(clone[i].Real, data[i].Real / n, 0.00001);
    }

    public static int SlidingDFT(IList<Complex> data, int sampleRate, int binsPerOctave, double minFrequency, int stepsFromMin)
    {
        var Q = 1.0 / (Math.Pow(2.0, 1.0 / binsPerOctave) - 1.0); // denominator * Fk gives size of bin
        var k = (int)(Q + 0.5);

        // frequency at k:
        var Fk = minFrequency * Math.Pow(2.0, (double)stepsFromMin / binsPerOctave);
        var Nk = (int)(Q * sampleRate / Fk + 0.5);
        Complex exp = Complex.Exp(new Complex(0.0, 2.0 * Math.PI * k / Nk));
        var current = data[0];
        for(int x = 0; x < data.Count - 1; x++)
        {
            var lookahead = x + Nk >= data.Count ? new Complex() : data[x + Nk];
            var fx = (lookahead - current) + data[x]; //  * 0.9999999 if unstable?
            current = data[x + 1];
            data[x + 1] = fx * exp;
        }
        return Nk;
    }

    public static void Conjugate(IList<Complex> data)
    {
        for(int i = 0; i < data.Count; i++)
            data[i] = Complex.Conjugate(data[i]);
    }

    public static void Swap(IList<Complex> data)
    {
        for(int i = 0; i < data.Count; i++)
            data[i] = new Complex(data[i].Imaginary, data[i].Real);
    }

I haven't tried to run your code (is that C#?) or looked at it in great detail, but SlidingDFT() looks more like a Goertzel algorithm implementation to me. It should be a red flag to you that you're able to perform the calculations in place. If you were truly doing a sliding DFT, then you would have a data rate expansion by a factor of $N$, the length of the transform; for each input time-domain sample, you generate a full DFT's worth of outputs.

The Goertzel algorithm is an efficient, recursive structure for implementing a single DFT bin at an arbitrary center frequency, which can be useful for applications like DTMF decoding. However, you can't "invert" it in the sense that you're trying to; it is more of a filtering operation, not a reversible transform.

It doesn't seem to make sense to me to even try to invert a sliding DFT in the way you've shown; the outputs of a sliding DFT don't "slide" sample-by-sample the same way as the time-domain inputs do. I'd read up some more on how the DFT and its sliding version works.


What is plotted against the DFT of a windowed sampled signal? I'm given a continuous-time analog signal $x_a(t) = \cos(2\pi f_1t)+\sin(2\pi f_2t)$, for some frequency $f_1, f_2$. I'm asked to sample $x_a(t)$ at $F_s=1024\textrm{ Hz}$, apply a 128-point Hamming Window, and take the DFT(Discrete Fourier Transform) of the resulting windowed, sampled signal.

The expression for a sampled signal is:

$$x[n] = x_a\left(\dfrac{n}{F_s}\right)$$

To apply the hamming window, we multiply $x[n]$ by the 128-point window function (courtesy of ML):

$$x_w[n]=x[n].*\mathrm{hamming}[128]$$

We then take the DFT of the windowed function $x_w[n]$:

$$X_w[k] = \mathrm{DFT}\left\{x_w[n]\right\}$$

I've read that the DTFT(Discrete-time Fourier transform) is a continuous spectrum of $x[n]$ and that the DFT of $x[n]$ (in a nutshell) is a sampling of that spectrum. When we take the $N$-point DFT of $x[n]$, we are taking $N$ samples of the DTFT of $x[n]$ where frequencies are unique (i.e., we sample within $2\pi$). That is:

\begin{align}
\mathrm{DTFT}\left\{x[n]\right\} &= X(\omega)\\
\mathrm{DFT}\left\{x[n]\right\} &= X\left(\omega = \frac{2\pi k}{N}\right) = X[k],\quad \textrm{ for } k=0,1,...,N-1.
\end{align}

My questions:

If we plot $X[k]$, we are still plotting the sample of the spectrum  against frequency, right? Is the frequency, given $k$, denoted by $2\pi k/N $ radians?
My professor asked me to plot the FFT versus $k/F_s$ where $k=0,\ldots,127$. What does $k/F_s$ represent?
A homework problem asks (given application of the 128-point Hamming window, above) to find the "continuous-time frequency spacing $\Delta F$ between DFT samples"? What does this mean? Isn't this just $2\pi/N$?

My professor told me that $2\pi/N$ is the discrete-time frequency spacing in radians and that I have to multiply that by $Fs/(2 \pi)$. I'm not sure how to make sense of $Fs/N$.

As you've noted, for a length-$N$ time-domain sequence $x[n]$, the $\mathrm{DTFT}$ (aka Fourier transform) is given by the equation:

$$
\mathrm{DTFT}\left\{x[n]\right\}=X(e^{j\omega})=X(\omega), \quad 0\leq n\leq N - 1.\tag{1}
$$

And the $N$-point $\mathrm{DFT}$ is given by uniformly sampling equation $(1)$ to obtain $N$ equally spaced frequencies $\omega_k = 2\pi k/N$:

$$
\mathrm{DFT}\left\{x[n]\right\}=X(e^{j\omega_k})=X\left(\frac{2\pi k}{N}\right)=X[k], \quad 0\leq k\leq N - 1.$$

If you plot $X[k]$, you're plotting the $N$ equally spaced samples from the $\mathrm{DFT}$ equation not against the continuous frequency $\omega$ but against the index number $k$ which in frequency corresponds to $2\pi k/N$.
The $\mathrm{FFTs}$ are algorithms to compute the $\mathrm{DFT}$. I don't know what's the goal in plotting the $\mathrm{DFT}$ against $k/F_s$. If you check the units you see that $k$ is just a bin number, unit-less, but $F_s$ is in Samples per second. So, that division has units in Seconds per samples...
The spacing between the discrete frequencies is indeed $2\pi/N$ radians. By the continuous-time frequency spacing, I think he meant the corresponding frequency spacing in $\rm Hz$ and this is the frequency resolution:
$$f_k=k\cdot \frac{F_s}{N}$$
From here you can guess what's the frequency spacing since two consecutive bins from the $\mathrm{DFT}$ samples are at $\Delta k=1$.


Is this calculation of inverse z-transform proper I wonder whether my calculation of inverse z transform are correct. My IIR system is described as follows in Z-domain  $H(z) = \frac{z^{-2}}{1-0.5z^{-2}}$

After using partial fraction decomposition I obtained $H(z)$ in Z-domain as follows $H(z)=\frac{\frac{1}{2}}{1-\sqrt{\frac{1}{2}}z^{-1}} + \frac{\frac{1}{2}}{1+\sqrt{\frac{1}{2}}z^{-1}}$.

After applying inverse Z transform I obtained such form: $h[n]=\frac{1}{2}(\sqrt{\frac{1}{2}}^{n}-\sqrt{\frac{1}{2}}^{n})u[n]$. What confused me is that this is 0 as the insides of parentheses cancel out. Maybe this is some dumb mistake I am making (it is quite late)

HINT:

Write $H(z)$ as

$$H(z)=G(z^2)\tag{1}$$

with $$G(z)=\frac{z^{-1}}{1-0.5z^{-1}}\tag{2}$$

Determine the inverse transform $g[n]$ of $G(z)$. Figure out what replacing $z$ by $z^2$ means in the time domain. From this, the inverse transform of $H(z)$ is easily obtained from $g[n]$.

EDIT: Your partial fraction expansion lacks the factor $z^{-2}$, and the major mistake is that you transform both terms to the same sequence, even though the poles have opposite signs. If you do it right, the two terms of the sequence only cancel each other for odd $n$, and you can combine them for even $n$. The method shown as a hint in this answer is probably more straightforward because it avoids partial fraction expansion.


What is an Image Prior? I am having trouble in understanding the concept of image-priors. What is it exactly? What is the advantage of using them? Where are they used?

It is "prior information" on your set of images, that you can use in your image processing problems to enhance results, ease the choice of processing parameters, resolve indeterminacies, publish fancy papers no one ever use, etc.

For instance, you may know that the image, albeit noisy, should contain only 4 colors. Or that pixels follow a specific distributions. These priors, or their approximations, can be put into math form and can be merged into the processing (filtering, deconvolution, segmentation), and reduce the set of feasible solutions, generally through optimization algorithms.


difficulties implementing windowed sinc interpolation (C++) I'm currently trying to implement a windowed sinc interpolation. I've already written some code for that, but it only seems to work for cases where phaseInc <= 1.0 (phaseInc = outSampleRate / inSampleRate). By "not working" I mean that something is causing a lot of distortion or aliasing. I hope that the comments in the code do explain what I'm trying to do:

// SINC_FILT_THRESH = 0.8, SINC_WINDOW_SIZE = 16

float sincStep = phaseInc > SINC_FILT_THRESH ? SINC_FILT_THRESH / phaseInc : 1.00f;
/* ^^ this should in theory detect the transition between phaseInc <= 1.0 and phaseInc > 1.0
 * but with a little headroom for the transition band. "sincStep" then shall be used later
 * to make the sinc function wider to narrow down the rectangle in frequency domain
 * The code below seems to work as long as sincStep = 1.0 */

int i = 0;
do {
    float sampleSum = 0.0f;
    float kernelSum = 0.0f;
    for (int wi = -SINC_WINDOW_SIZE + 1; wi <= SINC_WINDOW_SIZE; wi++) {
        float sincIndex = float(wi) * sincStep - phase;
        float windowIndex = (float(wi) - phase) / SINC_WINDOW_SIZE;

        float s = sinc(sincIndex);
        /* ^^ this is a normalized sinc function with zeroes at whole numbers */

        float w = windowFunc(windowIndex);
        /* ^^ windowFunc = 0.5 + (0.5 * cos(pi * t))
         * I know there is stuff like blackman window, but I'd rather like to start with
         * something simple */

        float kernel = s * w;
        sampleSum += kernel * fetchBuffer[i + wi + SINC_WINDOW_SIZE - 1];
        kernelSum += kernel;
    }
    phase += phaseInc;
    int istep = static_cast<int>(phase);
    phase -= static_cast<float>(istep);
    i += istep;
    /* ^^ This simply advanced phase and i by the appropriate amount.
     * I used it with a simple linear interpolator and didn't have problems with that */

    *outData++ = sampleSum / kernelSum;
    /* ^^ Store the result in the output buffer with unity gain */
} while (--numSamples > 0);

So, as I've mentioned it only seems to work right in cases where sincStep = 1.0 and I don't know why. My intention is to widen the sinc function to lower the cutoff. Doesn't seem to work quite as I want.

What I do know is that there are no memory access violation problems (according to fsanitize) and my fetchBuffer always has enough data (so no out of bounds thing).

I'm sorry if I'm missing something quote obvious but I really don't want to just copy someone elses code. I'd like to actually understand it. Most examples however just seem to focus on the equations and not on the actual thing with all the possibly faulty boiler plate code.

Thanks in advance!

EDIT:
I want to use that interpolation for resampling different waves at different pitch for music. I thought I could add this for clarification for my problem:

So, sleeping a night before posting the question on stackexchange won't help but trying out random things AFTER posting the question results in magic results.

It seems like the following line was faulty:

float sincIndex = float(wi) * sincStep - phase;

The solution was to change it to this:

float sincIndex = (float(wi) - phase) * sincStep;


State of the art algorithms for estimating noise variance Suppose I have a signal $X[n]= f[n] + W[n]$, where $W$ is uncorrelated Gaussian white noise satisfying $E\{W[n]W[k]\} = \sigma^2 \delta[n-k]$.

What are the state-of-the-art algorithms for determining the noise variance $\sigma^2$ from $X$? (i.e., without oracle information?)

If you hope that $f[n]$ can be sparser in some linear representation $L$, then robust statistics can be used to estimate variance bounds on coefficients of $L(X[n])$. One common method is orthogonal wavelet denoising, where a common standard deviation estimate is:
$$ \mathrm{median} \frac{|c_i|}{0.67449}
$$
where $c_i$ denotes the highest scale subband. The magic $0.67449$ scale factor relates the standard deviation to the robust Median absolute deviation through the error function erf ("MAD equals the half-normal distribution median").

How robust is that? I tested it on a classical "crane" image, with several noise realizations, and here is the $\sigma$ versus  $\hat{\sigma}$ estimates with dispersion.

Not to bad: somehow linear, albeit with an offset at zero, and an increasing  dispersion. This can suffice for some uses. Meanwhile, for a better interval of estimates, I usually use several representations (wavelets, filter banks) to capture unknown sparsity optimality, with a certain amount of redundancy, and combine them into a range of estimates. Since redundancy kills independence, a calibration is required to get thre proper factor, that easily done with transforming white noise with the set of transformations you have. This is what we use with redundant Optimization of Synthesis Oversampled Complex Filter Banks or for denoising with Noise Covariance Properties in Dual-Tree Wavelet Decompositions.

Although the estimate is excellent with simulated Gaussian noise, this is of course less efficient in practice, as noise estimated from one realization is rarely Gaussian (or something easy to handle).

Additional references for images:
- Noise Estimation from a Single Image, with Matlab code

  In order to work well, many computer vision algorithms require that
  their parameters be adjusted according to the image noise level,
  making it an important quantity to estimate. We show how to estimate
  an upper bound on the noise level from a single image based on a
  piecewise smooth image prior model and measured CCD camera response
  functions. We also learn the space of noise level functions how noise
  level changes with respect to brightness and use Bayesian MAP
  inference to infer the noise level function from a single image. We
  illustrate the utility of this noise estimation for two algorithms:
  edge detection and feature preserving smoothing through bilateral
  filtering. For a variety of different noise levels, we obtain good
  results for both these algorithms with no user-specified inputs.

Single-Image Noise Level Estimation for Blind Denoising, with Matlab code

  Noise level is an important parameter to many image processing
  applications. For example, the performance of an image denoising
  algorithm can be much degraded due to the poor noise level estimation.
  Most existing denoising algorithms simply assume the noise level is
  known that largely prevents them from practical use. Moreover, even
  with the given true noise level, these denoising algorithms still
  cannot achieve the best performance, especially for scenes with rich
  texture. In this paper, we propose a patch-based noise level
  estimation algorithm and suggest that the noise level parameter should
  be tuned according to the scene complexity. Our approach includes the
  process of selecting low-rank patches without high frequency
  components from a single noisy image. The selection is based on the
  gradients of the patches and their statistics. Then, the noise level
  is estimated from the selected patches using principal component
  analysis. Because the true noise level does not always provide the
  best performance for nonblind denoising algorithms, we further tune
  the noise level parameter for nonblind denoising. Experiments
  demonstrate that both the accuracy and stability are superior to the
  state of the art noise level estimation algorithm for various scenes
  and noise levels.


How can I cut exactly the first harmonic (H1) from a FFT of an image I have 4D( 2D + slices along z axis + time frames) gray-scale image for the heart beating on different moments.
I do like to take Fourier Transform along time axis(for each slice separately), and analyze the fundamental Harmonic (also called H1 component, where H stands for Hilbert Space) so I can determine pixel regions corresponding to ROI which show strongest response to cardiac frequency.
I'm using python for this purpose, and I tried to do that with the following code, but I'm not sure that this is the correct way to do it, because I don't know how to determine the cut-frequency to keep only the fundamental Harmonic.

This link to the image which I'm dealing with

import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt

img = nib.load('patient057_4d.nii.gz')
f = np.fft.fft2(img)
#  Move the DC component of the FFT output to the center of the spectrum
fshift = np.fft.fftshift(f)
fshift_orig = fshift.copy()
# logarithmic transformation
 magnitude_spectrum = 20*np.log(np.abs(fshift))
# Create mask
rows, cols = img.shape
crow, ccol = int(rows/2), int(cols/2)
# Use mask to remove low frequency components
dist1 = 20
dist2 = 10
fshift[crow-dist1:crow+dist1, ccol-dist1:ccol+dist1] = 0
#fshift[crow-dist2:crow+dist2, ccol-dist2:ccol+dist2] = fshift_orig[crow-dist2:crow+dist2, ccol-dist2:ccol+dist2]

# logarithmic transformation
magnitude_spectrum1 = 20*np.log(np.abs(fshift))
f_ishift = np.fft.ifftshift(fshift)
# inverse Fourier transform
img_back = np.fft.ifft2(f_ishift)
# get rid of imaginary part by abs
img_back = np.abs(img_back)
plt.figure(num = 'Im_Back')
plt.imshow(abs(fshift[:,:,2,2]).astype('uint8'),cmap='gray')
plt.show()

```

The solution was to take Fourier transform 3D for each slice, then to chose only the 2nd component of the Transform to transform it back to the spatial space, and that's it.
The benefit of this is to detect if something is moving along the third axis(time in my case).

for sl in range(img.shape[2]):
   #-----Fourier--H1-----------------------------------------
   # ff1[:, :, 1] H1 compnent 1, if 0 then DC
   ff1 = FFT.fftn(img[:,:,sl,:])
   fh = np.absolute(FFT.ifftn(ff1[:, :, 1]))

   #-----Fourier--H1-----------------------------------------
```


Can linear transforms of vector space be seen as LTI systems? Linear transforms of vector spaces has the linearity (i.e., homogeneity and superimposability), which is shared by LTI system. There are also shared concepts such as eigenvectors/eigenvalue...

I am wondering, can linear transforms (of vector spaces) be seen as LTI systems? If not, why? If yes, in particular, are there counterparts in vector spaces for the concept of Dirac $\delta$ and $h(t)$ in LTI systems, where $\delta$ represents all possible inputs (e.g., from frequency space point of view), and thus $h(t)$ characterizes the LTI system? If these two counterparts exists in Linear transforms, what's the counterpart for convolution?

A little background: I had this question recently as I am studying linear algebra (starting from some basics) ...actually I first read some books on DSP/DIP a couple months ago, and then I find out that I need to learn some basics of linear algebra (as well as other branches of math), otherwise most of the DSP/DIP literature are foreign to me...in other words, I am also new to the field of DSP/DIP...so I beg your pardon that my question may be a bit fuzzy.

Thanks.

The question of abstract nature so a rigorous answer could require more elaboration, however on the surface one can at least observe the following.

An LTI system is also a linear transform, which maps its input signal (vector) $x(t)$ or ($ x[n]$) into an output signal (vector) $y(t)$ (or $y[n]$), where both the input and output belong to the same signal space, one kind of vector space.

For such transforms and vector spaces the time invariance property can be checked easily such as $y(t-d) = \mathcal{T} \{x(t-d) \}$, but for other vector spaces such as the n-tuples or matrix spaces etc, the operation of time invariance is not trivial to check, at least intuitively. The shift of a vector (an element of the set) should be defined first.

May be the following links could help further:

link1

link2


Exponential decaying step response in LTI System I'm attempting to better understand the relationship between step responses, impulse responses, and convolutions. Say that I have a system where if I apply a constant input, my output decays from a certain value to an equilibrium value. This gives me the step response of a system that I believe to be LTI based on the system's differential equation.

PDE:

$$\frac{\partial^2u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}-\frac{u}{r^2}=\frac{1}{H_A k}\frac{\partial u}{\partial t}+\frac{1}{H_A k}\frac{r}{2}\frac{\partial \epsilon}{\partial t} $$

ODE in S Domain

$$\frac{\partial ^2 \bar{u}}{\partial \hat{r}^2}+ \frac{1}{\hat{r}}\frac{\partial \bar{u}}{\partial \hat{r}}-\frac{\bar{u}}{\hat{r}^2}-s\bar{u}=\frac{s\bar{\epsilon}\hat{r}}{2} $$

My specific step response is fairly complex and I went more in-depth into this specific system in another post, but a simpler version that still captures my issue could be given as:

$$
Y(t) = 1+e^{-0.1t}
$$

In an LTI System, it's my understanding that the derivative of the step response yields the impulse response, which can then be used to obtain the output of the system for any arbitrary input through convolution.

$$
h(t) = \frac{dY(t)}{dt} = -0.1e^{-0.1t}
$$

So now, if I have this impulse response, which is negative, and I attempt to convolve this with a unit step function to recover my original Y(t), I get something that is both the wrong shape and the wrong sign. Can anyone help point out my mistake?

You're forgetting the discontinuity at $t=0$. Your step response is actually

$$Y(t)=(1+e^{-0.1t})u(t)\tag{1}$$

where $u(t)$ is the unit step function. Differentiating $(1)$ gives

$$h(t)=-0.1e^{-0.1t}u(t)+(1+e^{-0.1t})\delta(t)=-0.1e^{-0.1t}u(t)+2\delta(t)\tag{2}$$

where $\delta(t)$ is the Dirac delta impulse. Convolving $(2)$ with a unit step gives the original step response:

$$\begin{align}\int_{-\infty}^{\infty}h(\tau)u(t-\tau)d\tau&=-0.1u(t)\int_{0}^{t}e^{-0.1\tau}d\tau+2\int_{-\infty}^{t}\delta(\tau)d\tau\\
&=u(t)(e^{-0.1t}-1)+2u(t)\\&=(1+e^{-0.1t})u(t)\end{align}$$


PyAudio recording is taking much longer than planned I'm trying to record sounds in chunks using PyAudio for 15 seconds, chunk size being 1024 and sampling rate is 1024 *40. I understand the number of chunks would be (time of recording) * (sampling rate) / (chunk size).

However, when I timed the time it takes for the for loop to finish, it's much longer than 15 seconds, around 2 than 3 times longer than expected. I wonder why? Granted I have a lag in each loop to display the FFT of each chunk, but that's only 0.001 second.

Is this because Windows is not a real time system?

import time
import image
import numpy as np
import os
import pyaudio
import matplotlib.pyplot as plt
from scipy.fftpack import fft

CHUNK = 1024
FORMAT = pyaudio.paInt16
CHANNELS = 1
RATE = 1024*40
RECORD_SECONDS = 15

p = pyaudio.PyAudio()
stream = p.open(format=FORMAT,
                channels=CHANNELS,
                rate=RATE,
                input=True,
                frames_per_buffer=CHUNK)

input('Press key to start recording:\n')
print('2 sec delay started')
time.sleep(2)
print("* recording")

n = CHUNK
k=np.arange(n)
T = n/RATE #reciprocal of freq resolution
frq = k/T #k* freq resolution for plotting
window=np.hanning(n) #apply Hanning window to minimize spectrum leakage
num_frames=int(RATE*RECORD_SECONDS/CHUNK)

for i in range(0, num_frames):
    data = stream.read(CHUNK)
    decoded = np.frombuffer(data, dtype=np.int16)
    #decoded= decoded/2**15
    windowed=window*decoded
    fft_decode=fft(windowed)/(len(decoded)/2)
    mags=np.absolute(fft_decode)

    plt.ylim(top=200)
    plt.xlabel('Freq (Hz)')
    plt.ylabel('|Y(freq)|')
    plt.plot(frq[range(int(n/2))], mags[range(int(n/2))],'b')
    plt.pause(.001)
    plt.gcf().clear()

print("* done recording")
plt.close()

First off, a small thing:

The sampling rate usually has to be one of the standard sampling rates.
Closest to your 1024*40 would be 44100.

The hardware has to support the sampling rate. Some drivers will resample to the given rate, but that will be somewhat slow as it is in software. Try using a proper sampling rate.

The real bottle neck is the plotting.

Your are calculating the FFT for every block, and plotting it.

The FFT is blazing fast, no worries there.

Plotting on the other hand is horribly slow.  The plotting routine is not optimized for speed, and you are also rescaling and formatting the plot on every go around.

Comment out the lines that do the plotting.  Leave the FFT in.  You'll see a tremendous change in speed.

Edit:
Just noticed.  Besides plot being slower than molasses, you also have a pause in there.  It is only 1 millisecond, but I don't know how precise it really is.  Many timers under Windows use 15 millisecond steps, so that pause could be zero, or 1 millisecond, or 15 milliseconds.

In any case, it is not doing you any good.  Take it out.


Beamforming but only for total power? Is there a method to speed up beam-forming if you only care about the total power in a given direction, or total power in frequency groups?  I've googled this, and can't find anything.  My intuition says that wouldn't make sense given that you need the power in each frequency to properly cancel out signals to create your beam.

I am not sure regarding total power but you can certenly focuse on ccertain frequencies of interest.

Beamforming is a process in which different channels (a number of microphones for example) gather their information to cancel the uncorrelated component between them. Let $y(f)=[Y_1, Y_2, \ldots, Y_M]^T$ be signals received from $M$ channels S.T. $y(f)=x(f)+v(f)$, where $v(f)$ is noise. We construct:

$$Z(f)=\sum_{m=1}^{M}H^*_m(f)Y_m(f)=X_d(f)+V_r(f)$$

Where $X_d$ is our reconstructed desired signal and $V_r$ is the residual noise. Given a directional beam $d(f,cos\theta_d)=[1,e^{-j2\pi f\tau_0\cos\theta},\ldots,e^{-j2(M-1)\pi f\tau_0\cos\theta}]^T$, the white noise gain (WNG) of such a filter is given by:
$$\mathcal{W}[h(f)]=\frac{|h^H(f)d(f,cos\theta_d)|^2}{h^H(f)h(f)}$$
Note that all of these are functions of $f$, meaning they can be applied in only a single or couple of specific frequencies of interest.

To find such a filter, simply solve for the maximum on $\mathcal{W}[h(f)]$ (the higher WNG is the better you reconstruction of the signal is) which yields the common delay and sum beamformer:

$$h_{DS}(f,cos\theta_d)=\frac{d(f,cos\theta_d)}{M}$$

There are many types of beamformers that can be produced such as the maximum directivity factor (DF) beamformer. I find this book is a good source. All of them can be applied as a function of $f$.


Covariance vs Autocorrelation I'm trying to figure out if there is a direct relationship between these concepts. Strictly from the definitions, they appear to be different concepts in general. The more I think about it, however, the more I think they are very similar.

Let $X,Y$ be WSS random vectors. The covariance, $C_{XY}$, is given by $$C_{XY}=E\left[(X-\mu_x)(Y-\mu_y)^H\right]$$ where $H$ stands for the Hermitian of the vector.

Let $Z$ be a WSS random vector. The autocorrelation function, $R_{XX}$, is given by $$R_{ZZ}(\tau)=E\left[\left(Z(n)-\mu_z\right)\left(Z(n+\tau)-\mu_z\right)^H\right]$$

Edit Note There is a correction to this definition as applied to signal-processing, see Matt's Answer below.

The covariance does not involve a concept of time, it assumes each element of the random vector is a different realization of some random generator. The autocorrelation assumes a random vector is the time evolution of some initial random generator. Yet in the end, they are both the same mathematical entity, a sequence of numbers. If you let $X=Y=Z$, then it appears $$C_{XY}=R_{ZZ}$$ Is there something more subtle that I am missing?

According to your definition of autocorrelation, the autocorrelation is simply the covariance of the two random variables $Z(n)$ and $Z(n+\tau)$. This function is also called autocovariance.

As an aside, in signal processing, the autocorrelation is usually defined as

$$R_{XX}(t_1,t_2)=E\{X(t_1)X^*(t_2)\}$$

i.e., without subtracting the mean. The autocovariance is given by

$$C_{XX}(t_1,t_2)=E\{[X(t_1)-\mu_X(t_1)][X^*(t_2)-\mu^*_X(t_2)]\}$$

These two functions are related by

$$C_{XX}(t_1,t_2)=R_{XX}(t_1,t_2)-\mu_X(t_1)\mu^*_X(t_2)$$


Amplitude spectrum (transfer function) of signal? I have one question related to finding amplitude spectrum (transfer function) of signal knowing that output signal is time derivative of input signal. I have the answer graph but I don't understand the reason.

HINT:
If $X(j\omega)$ is the Fourier transform of a signal $x(t)$, try to remember (or look up) how $X(j\omega)$ is related to the Fourier transform of the derivative $x'(t)$.
Some people are more familiar with this in the Laplace transform domain: the time derivative corresponds to which operation in the Laplace transform domain? Replace $s$ with $j\omega$, and you'll understand the given figure.


Eliminating inertia effect from the acceleration signal I am struggling with a problem, and would be happy if you could assist me. I put an accelerometer sensor inside a spherical enclosure, and release it to the water. While it was translating in the flow, I recognized that the enclosure was also wobbling (like a kicked pendulum), which I think that it is likely caused by the inertia of the enclosure. I then added some weight inside the enclosure and observed that the peak was more pronounced. Left panel is the spectra of the acceleration magnitude for three different tests without additional weight (inertia is only described by the sensor + enclosure + battery).

Since the sensor is like a black box (I do not know the mass distribution inside so that I can try to match it with the geometrical center of the enclosure), I would like to find a way to exclude the effect of the wobbling in the signal without disturbing it because the wobbling might also affect the lower frequency part of the signal, as seen in green line in page 118.

Based on the comment below, I have added the spectra after notch filter. Notch filter does its work very well, and attenuated the signal in the region (2 - 5 Hz), without affecting the rest. However, the oscillations (around 3.5 Hz) likely influences the lower frequency range of the spectra. Is there any quantity that I can derive to see the effect of these oscillations in the lower frequency range or any function that I can try to model the effect of inertia? Any ideas and suggestions would be appreciated.

It sounds like you want to remove the frequency component at about 2.5Hz in the plot, but leave the rest of the plot as untouched as possible.

The usual first approach to that problem is to just use a notch filter:

$$
H(z) = K \displaystyle\frac{1 - 2 \cos \omega_0 z^{-1} + z^{-2}}{1 - 2r \cos\omega_0 z^{-1} + r^2 z^{-2}}
$$

where $\omega_0$ is the (normalized) frequency you want to remove and $r \lt 1$.

There are more complex approaches that use, for example, Kalman filtering... but that might require more assumptions about the dynamics of your system that you are prepared to make given your "black box" comments.


How to determine precision needed for sin approximation used for sound synthesis? I am doing some procedural sound synthesis in Java. I want to have a sine wave as one of the possible basic sounds. When experimenting with that, I have found the default Java Math.sin implementation causes significant performance issues to me. I would like to replace it with either a table lookup or a simple polynomial (Taylor or Chebyshev) series.
My question is: how to decide what table size or what polynomial degree to use so that the resulting artifacts are inaudible?
I can (and will) perform some basic experiments, but I do not trust my ears and speakers that much.
The sound produced will be a tone with a frequency in range 100 Hz - 10 kHz. I want the result to be indistinguishable from a real sine wave by a human. My sound output will be 16 bit at 44.1 kHz or 48 kHz, if that is important.

how many entries in the lookup table do I need, assuming the table will be used for -pi/2 to pi/2 range?

If I want to use a polynomial series (Taylor or Chebyshev), what degree does the polynomial have to be?

This is very good question.
Fortunately for you, I have a very good answer.
Assuming standard CD quality, your sound level has just over 4 significant digits accuracy, thus this is the level you need to be truly indistiguishable.
Is this level necessary for auditory purposes?  Let's assume so.
Your fastest solution, by far, is going to be linear interpolation on a segmented domain.  You need two tables.  The first hold the sine value of the center of each interval.  The second holds the first derivative, aka cosine.
To calculate your sine value, multiply your angle by the spacing you are using (so a unit of one corresponds a segment).  The integer portion will give you the index into your table and use the fractional portion (-0.5 to 0.5) for the interpolation.
Now the question becomes "How many entries do I need in my table?"
I'll let you take a crack at that.

I have a better answer (I think, you need to test it).
$$ \sin( x + d ) = \sin( x ) \cos( d ) + \cos( x ) \sin( d ) $$
Back to two tables.  One for the broad range ($x$) , and one for the fine range ($d$).
If the range of $d$ is small enough, then you can use (Taylor or find the Remez instead):
$$ \cos(d) \approx 1 - x^2/2 + x^4/24 $$
$$ \sin(d) \approx x - x^3/6 + x^5/120 $$
But that will require more computation.
For uber accuracy, which you don't need, you could do the interpolation thing on the fine table value.

I was curious, so here you go:

import numpy as np

#==========================================================
def main():

        N_coarse = 128
        N_fine   = 128

#---- The coarse table

        sc = np.zeros( N_coarse )   # Sine Coarse
        cc = np.zeros( N_coarse )   # Cosine Coarse

        theSlice = np.pi * 0.5 / N_coarse
        theAngle = 0.0

        for n in range( N_coarse ):
          sc[n] = np.sin( theAngle )
          cc[n] = np.cos( theAngle )
          theAngle += theSlice

#---- The fine table

        sf = np.zeros( N_fine )   # Sine Fine
        cf = np.zeros( N_fine )   # Cosine Fine

        theSlice /= N_fine
        theAngle = 0.0

        for n in range( N_fine ):
          sf[n] = np.sin( theAngle )
          cf[n] = np.cos( theAngle )
          theAngle += theSlice

#---- The test

        theFactor = N_coarse * 2.0 / np.pi

        for a in range( 157 ):
          theAngle = a * 0.01
          n = theAngle * theFactor

          nc = np.floor( n )
          nf = np.floor( ( n - nc ) * N_fine )

          sine = sc[nc] * cf[nf] + cc[nc] * sf[nf]

          print sine, np.sin( theAngle ), sine - np.sin( theAngle )

#==========================================================
main()

This old fart got it right on the first try!
Here are the first few rows of the output, the rest is comparable:

0.0 0.0 0.0
0.00997070990742 0.00999983333417 -2.91234267487e-05
0.0199404285515 0.0199986666933 -5.82381418187e-05
0.0299081647675 0.0299955002025 -8.73354349791e-05
0.0399687249608 0.0399893341866 -2.06092258625e-05
0.0499294807897 0.0499791692707 -4.96884810173e-05
0.059885272753 0.0599640064794 -7.87337263996e-05
0.0699307504776 0.0699428473375 -1.20968599355e-05
0.0798735744039 0.0799146939692 -4.11195652507e-05
0.089808457497 0.089878549198 -7.00917010057e-05
0.0998298073783 0.0998334166468 -3.60926852212e-06
0.109745746461 0.109778300837 -3.25543758203e-05
0.119650774894 0.119712207289 -6.14323952304e-05
0.129543907942 0.12963414262 -9.02346778173e-05


How do I add AWGN to an I and Q representation of a signal? I have a wireless communication system that I am simulating in Matlab. I am performing some watermarking through slightly adjusting the phase of the transmitted signal. My simulation takes the original I (inphase) and Q (quadrature) values and adds in the watermark. I then have to simulate the resulting bit error rate after being transmitted. For now I just need to add varying amounts of thermal noise to the signal.

Since I have the signal represented as its I and Q channel it would be easiest to add AWGN(additive white Gaussian noise) to the I and Q directly. One thought was to add noise to both channels independently, but my intuition tells me that this isn't the same as adding it to the signal as a whole.

So how can I add noise to it when it is in this form?

Yes, you can add AWGN of variance $\sigma^2$ separately to each of the two terms, because the sum of two Gaussians is also a Gaussian and their variances add up. This will have the same effect as adding an AWGN of variance $2\sigma^2$ to the original signal. Here's some more explanation if you're interested.

An analytic signal $x(t)=a(t)\sin\left(2\pi f t + \varphi(t)\right)$ can be written in its in-phase and quadrature components as

$$x(t)=I(t)\sin(2\pi ft) + Q(t)\cos(2\pi ft)$$

where $I(t)=a(t)\cos(\varphi(t))$ and $Q(t)=a(t)\sin(\varphi(t))$. If you wish to add AWGN to your original signal as $x(t)+u(t)$, where $u(t)\sim\mathcal{N}(\mu,\sigma^2)$, then you can add AWGN to each of the terms as

$$y_1(t)=\left[I(t)\sin(2\pi ft) + v(t)\right] + \left[Q(t)\cos(2\pi ft) + w(t)\right]$$

where $v(t), w(t)\sim\mathcal{N}(\mu/2,\sigma^2/2)$

Also note that because the in-phase and quadrature terms are additive, the AWGN can also be simply added to either of the two terms in the $IQ$ representation of $x(t)$ above. In otherwords,

$$y_2=I(t)\sin(2\pi ft) + \left[Q(t)\cos(2\pi ft) + u(t)\right]$$
$$y_3=\left[I(t)\sin(2\pi ft) +u(t)\right]+ Q(t)\cos(2\pi ft)$$

are statistically equivalent to $y_1$, although I prefer using $y_1$ because I don't have to keep track of which component has noise added to it.


Is the power or root-power scale correct for FFT analysis of an audio file? I wrote Python scripts to create a spectrum analysis of audio files. The goal is to compare multiple spectrums to get the dB difference at certain frequencies.

When I started using the Sonic Visualizer tool for my first analyses, then switched to my automated solutions using Python with SciPy and NumPy, I noticed a large difference in the values of the spectrum analysis.

It seems Sonic Visualizer uses the power scale, while most other audio tools use the root-power scale for the spectrum analysis.

Now I am a little bit confused about the two quantities.

Shall one use the power scale or the root-power scale for spectrum analysis of audio files?

If you use $dB$ properly there should be no difference. You should use $20 \cdot log_{10}$ for root power quantities and $10 \cdot log_{10}$ for power quantities. If you do this, the results should be the same.


Proof for determining Fourier coefficients While determining Fourier coefficients we have this equation
$$\int^{T}_{0} x(t) e^{-jn\omega_0t} dt = \sum^{+\infty}_{k\ =\ -\infty} a_k [\int^{T}_{0} e^{j(k-n)\omega_0t}dt]$$

I want to ask that how can we get final equation of $a_n$ as:
$$a_n= 1/T \int^{T}_{0} x(t)e^{-jn\omega_0t}dt$$  from the above equation.
please any one help to explain this

For the periodic signal $x(t)$ with fundamental period $T=1/F_0$, proceeding from where you are:

$$
\Rightarrow \sum^{+\infty}_{k\ =\ -\infty} a_k \left[\int^{T}_{0} e^{j(k-n)\omega_0t}dt\right] = \sum^{+\infty}_{k\ =\ -\infty} a_k \left[\frac{e^{j(k-n)\omega_0t}}{j(k-n)\omega_0}\right]^{T}_{0}
$$

For $k \neq n$ you got $0$ (you can prove that)
And for $k=n$ you see that $\int^{T}_{0}dt = T$

You therefore have:

$$\int^{T}_{0} x(t)e^{-jn\omega_0t} dt= a_kT = a_nT$$

You can proceed from there.

EDIT:

We do not assume the existence of the two possible solutions. The right-hand side of the equation evaluates to that. It can be written in compact form using the unit impulse which has two distinct values depending on the index. (see the Kronecker delta in digital signal processing). Think $\delta_{kn}$.
We still need the summation, but we can start with the integral. It's only that the integral is zero for all values of $k$ except at $k=n$ where it is $T$. So you remain with one term in your summation.

All in all, your initial equation equals: $\displaystyle \sum^{+\infty}_{k\ =\ -\infty} a_kT\delta_{kn}$


Signal approximation using linear combination of functions How I can approximate the signal $x(t)=0.001\,t^3 \exp(-0.1t)$ in the interval $[0,100]$ using a linear combination of the following functions:

$f_1(t)=A_1$

$f_2(t)=A_2\cos(0.05t)$

$f_3(t)=A_3\cos(0.1t)$

$f_4(t)=A_4\cos(0.2t+1)\exp(-0.2t)$

$f_5(t)=A_5\,t^3$

I try to write a matlab code. Is this correct?

t=[0:100];
x_t=0.001*(t.^3).*exp(-0.1*t); %signal given for aproximation
f1_t=x_t.^0;
f2_t=cos(0.05*x_t);
f3_t=cos(0.1*x_t);
f4_t=cos(0.2*x_t+1).*exp(-0.2*x_t);
f5_t=x_t.^3;
M=[f1_t' f2_t' f3_t' f4_t' f5_t']; %matrix with linear components
A=M\x_t'; %matrix with coefficients
f_t=(A(1)*f1_t)+(A(2)*f2_t)+(A(3)*f3_t)+(A(4)*f4_t)+(A(5)*f5_t);
figure(1),plot(f_t);

The method you chose is a standard least squares approximation, i.e. you minimize the sum of squared errors on the chosen grid. You obtain your coefficients by solving an overdetermined system of linear equations in a least squares sense. I think that's a very sane approach.


Is this sampling frequency right? I designed a Butterworth second order bandstop filter on Matlab with the following parameters:

Notch frequency: 150 Hz (the frequency I want to remove)
Bandwidth = 8 Hz, that is
Fc1=146 Hz
Fc2=154 Hz
Fs (sampling frequency) = 10,000 Hz

According to Nyquist Criterion, 0.5xFs > Fm
where Fm is my signal’s highest frequency. Let's say that my signal frequency ranges from 0 to 3 kHz, so my sampling frequency should be 6 kHz or higher. As you saw, I chose a samplig frequency of 10 kHz.

MATLAB gave me the filter coefficients for my difference equation, that is

y[n]=b0x[n]+b1x[n−1]+b2x[n−2]−a1y[n−1]−a2y[n−2] (IIR Filter).

I implemented this equation in a STM32F767 Nucleo Board. The input signal goes to the ADC, the signal is processed and then go to the DAC for output. I used a timer interrupt every 100 us for this. I mean, the process takes place in the timer interrupt loop.

**

The problem is this:
My filter works, but it's no filtering the frequency I chose to filter. I wanted to remove the 150Hz, but it's removing a frequency of 298~300Hz.

**

What could be happening here? Does the sampling frequency have to do with this? Should I choose a lower sampling frequency? If it's needed I can attach my C code as well as the MATLAB design. But this is the basic idea:

struct difference_equation
{
  // coefficients
  double b0;
  double b1;
  double b2;
  double a1;
  double a2;
  // input history
  double x1;
  double x2;
  // output history
  double y1;
  double y2;

  double process(double x0)
  {
    double y0 = b0 * x0 + b1 * x1 + b2 * x2
      - a1 * y1 - a2 * y2;

    x2 = x1;
    x1 = x0;

    y2 = y1;
    y1 = y0;

    return y0;
  }
}

The STM32 is a 32bit micro-controller and all of your processing seems to be taking place over a double precision data type.

Doubles are 64bit "numbers" but not in the straightforward sense that 64bit integers are. Therefore, mathematical operations on doubles can take some additional steps (see page 47 onwards) on CPUs that are not equipped with special hardware to handle floating point operations. If the CPU does not have an FPU then mathematical operations are handled in software which means that they take more CPU cycles to complete (For example, you can try double precision multiplications on an Arduino (like Uno), it's going to be fun).

The STM32F767 does have (see page 13) a Floating Point Unit, which means that it can potentially perform double precision operations faster than when these operations are implemented in software and take up CPU cycles. Notice that the FPU is better than not having it but it still has its own limits, for indicative timings see page 16 on this document

BUT, this has to be enabled (see page 20).

Therefore, my guess is that your interrupt routine is striving to fire at the intended 10kHz but it doesn't complete before $\frac{1}{10000}$ seconds. So, your effective sampling frequency is not 10kHz and the whole system is transposed to a different frequency range.

You might want to check if the FPU is enabled at double precision or move to fixed point numbers which every CPU will happily crunch. With fixed point numbers, you have another problem which is that your filter might go out of spec because of the rounding that is inevitable due to the reduction in precision. But these are aspects that you have to manage in your solution.

Hope this helps.


Why ideal filter cannot realised practically? I know the difference between ideal and practical filterI just wanted to know why it not realisable?

Most typical examples of the ideal filters are the classical brickwall frequency selective filters (lowpass, highpass, bandpass etc) which are defined as ideal because of their impossible to realize ideally selective frequency responses.

Their frequency responses include exactly flat passbands, exactly flat stopbands having absolute zero gain, and zero width transition bands. These features make it impossible to realize them using either existing physical devices or finite time computable algorithms...

Assuming an LTI system model for such filters, then their impulse responses tend to be of infinite length and typically be acausal; another reason for denying their physical realization. Yet causality alone is only a concern for real-time processing of time domain signals. It does not prevent realization of such filters in offline procesing where the whole signal data with its past and future are already available. Yet even in such an offline processing mode the fact that they require an infinite length of data is the reason for preventing their realization.

Another point of view into their unrealizability is that input/output relationship of such ideal filters cannot be described by finite order differential / difference equations. In principle any filter whose input/output relation can not be described by a finite order differential / difference equation is said to be unrealizable; i.e., there is no finite time computable algorithm to implement it...


Why is this DFT of a real symmetric signal resulting in complex valued coefficients? I am trying to understand exactly how sampling the DTFT to get the DFT works. The signal I'm trying to analyze is $x(n)$ seen below.

$$x(n) = \delta(n\pm2) + 2\delta(n\pm1) + 3\delta(n)$$

Taking the DTFT, we have
\begin{align}
X(\omega)&= \sum\limits_{i=-\infty}^\infty x(n)e^{-j\omega}\\
&= \left(e^{j2ω} + e^{-j2ω}\right) + 2\left(e^{jω} + e^{-jω}\right) + 3\\
&= 2\cos(2ω) + 2(2\cosω) + 3\\
&= 3 + 4\cosω + 2\cos(2ω)
\end{align}

I next implemented this in MATLAB:

x = [1 2 3 2 1];
N = size(x,2);
w = -pi:0.01:pi;
X_DTFT_computational = freqz(x,1,w);
X_DTFT_analytical = 3 + 4*cos(w) + 2*cos(2*w);

This results in the following graphs:

Next, I compute the DFT in two ways:
1. I apply an FFT to the original signal $x(n)$.
2. I sample the DTFT.

The code that achieves this is as follows:

% FFT of x
X_DFT_computational = fftshift(fft(x))

% DFT = Sampled DTFT
X_DFT_analytical(1) = 3 + 4*cos(-4*pi/N) + 2*cos(2*(-4*pi/N));
X_DFT_analytical(2) = 3 + 4*cos(-2*pi/N) + 2*cos(2*(-2*pi/N));
X_DFT_analytical(3) = 3 + 4*cos(0)       + 2*cos(2*(0));
X_DFT_analytical(4) = 3 + 4*cos(2*pi/N)  + 2*cos(2*(2*pi/N));
X_DFT_analytical(5) = 3 + 4*cos(4*pi/N)  + 2*cos(2*(4*pi/N));

Plotting the DFT (while showing the DTFT for comparison), I have

The magnitude of the FFT does give me the correct result.
But just plotting the FFT alone gives me a complex signal (the graph shows the real part).

Why is this signal giving a complex FFT when $x(n)$ is both real and symmetric?

You defined the signal vector as x = [1 2 3 2 1]. Since the DFT is defined by

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

the command fft(x) computes the DFT of the signal

$$x[n]=\delta[n]+2\delta[n-1]+3\delta[n-2]+2\delta[n-3]+\delta[n-4]$$

This signal is not symmetrical with respect to $n=0$, and it is not equal to the signal you computed the DTFT of.

If you want to compute the DFT of the signal

$$x[n]=\delta[n+3]+2\delta[n+1]+3\delta[n]+2\delta[n-1]+\delta[n-3]$$

you have to periodically continue it in the interval $n\in [0,N-1]$ (with $N=5$):

x = [3 2 1 1 2]

This results in a real-valued and symmetrical DFT, as expected:

X = fft(x);
X =

   9.00000   2.61803   0.38197   0.38197   2.61803


Effective Nyquist frequency / signal resolution for overlapping time averaged signal samples Say that I have a continuous signal. If I sample the signal at sampling rate f_sample, then the highest signal frequency that I can resolve without worrying about aliasing, etc. is the Nyquist rate f_Nyquist = 0.5 * f_sample, so arguably I could interpret this as meaning that the time resolution of my signal is dt_min = 1 / f_Nyquist = 2 / f_sample. Fine.

However, the above is for instantaneous sampling. How does the picture change if each of my samples involves averaging over a given finite time period t_width, with the possibility that successive time periods overlap by a time t_overlap? (See picture. Equation for t_overlap in the picture ignores cases where overlap extends beyond that of nearest neighbors, but it's easy to generalize.)

I.e., in the case that t_width = t_overlap = 0, I have f_Nyquist = 0.5 * f_sample, dt_min = 2 / f_sample, as above. For the case that t_width, t_overlap != 0, is there still a simple statement that I can make regarding f_Nyquist, or more importantly some generalization of f_Nyquist which is relevant to the quantity that I am actually interested in, namely the time resolution of my signal t_min? Is there a simple equation for something that is interpretable as t_min that could be used in such a situation?

Extra:

The above is the main thing I am interested in. However, just to add one more generalization in case anyone happens to know it: If the averaging weight function over t_width is not uniform, but is instead something more complicated (e.g., Gaussian), is there still a relationship or theorem that can help me figure out what t_min is?

Thanks for any help!

P.S.

If my interpretation of the Nyquist rate is wrong or anything, let me know and I'll correct it in the post.
I think this is the right StackExchange site for this question, but if there is another that anyone knows of where an answer is more likely, please suggest it.

The process that you are describing can be broken down into two steps

Filter with a moving average filter with length $t_{width}$
Sampling with an interval $delta\_t$

A moving average filter is a low pass filter, but not a very good with lots of side lobes. It will help reduce aliasing in the sample process, but it won't eliminate it with any type of hard "threshold". The longer $t_{width}$ is and the the smaller $delta\_t$ is the less aliasing you get.

I think you approaching this backwards: Your "time resolution" is simply determined by the highest frequency in the signal. You need to determine this first and work backwards from there.

In practice, there is always aliasing. No real world signal is actually bandwidth limited in the mathematical sense. So you need to determine the frequency range "of interest", figure out how much energy the signal has outside this region, determine how much noise or aliasing you can tolerate, and design a suitable anti-aliasing filter to meet your specific requirements.


Matlab - Using findpeaks with fvtool Is there any way I could utilise findpeaks in conjunction with fvtool? I would like to know where on the frequency axis the maximum and minimum of my frequency response happens.

fvtool is interactive so it doesn't let you recover an data in a script. Fortunately you can use the freqz function to get the same frequency response. Something like this might work:

[h,w] = freqz(b,a);  % h is frequency response, w normalized frequency
mag = 20*log10(abs(h));
[minval,minidx] = min(mag);
minfreq = w(minidx);
[maxval,maxidx] = max(mag);
maxfreq = w(mcxidx);


For Direct-Sequence Spread Spectrum, can we reduce the transmit power of the signal to maintain a fixed BER? How to reduce the transmit power of the signal to maintain a fixed BER in CDMA ? I am thinking noise power spectral density, signal energy , and signal rate.

BER depends on your E_b/N0. You need to keep that constant if you want constant BER.
Now, you don't say how your comparing your DSSS system, but it would make sense to assume that your symbol rate constant, thereby increasing your bandwidth by the spreading factor, and keeping the duration of a symbol constant.
In that case, Energy is still the product of power and duration, and the latter stayed constant.


Tapped Delay Line Filter Distinction Is there a distinction between Tapped-Delay-Line Filters and the more general FIR (Finite Impulse Response) filters. From what I understand the Tapped-Delay-Line is a type of FIR filter since it's impulse response settles to zero in a finite duration...but I'm having trouble determining what might make these two things separate in the world of DSP.

A tapped-delay-line is really a methodology or implementation. An FIR filter is a more general description of a purely feedforward digital network.

While your impulse response is finite (no feedback components) you can consider the tapped-delay-line filter to be a particular FIR filter whose coefficients are non-zero at tap positions and zero otherwise.

Should you introduce feedback components, I think you could still describe the system as a tapped-delay line filter, but it would cease to be an FIR filter at this point and become an IIR filter.


Need help cleaning noise I have been recording sound sample from mouse vocalisations and noticed that a background noise appeared in the recording room. the noise is clearly visible in the spectrograms of the recordings as horizontal lines at frequencies 20, 40, 60, 80, and 100 kHz (ultrasonic range). I have tried to record in an acoustically isolated box, lined with sound absorbing material, and it helps, however, the box retains unwanted odours which effect the behaviour of the mice, rendering it unusable.

My question is if there is a method to isolate and remove the noise without affecting the rest of the sounds. I want to study the vocalisations, with an emphasis on the frequency in which the mice communicate under different social encounters, and having that constant noise in the background may affect the results of the study.

Vocal communications were recorded using a 1/4 inch microphone, connected to a preamplifier and an amplifier (Bruel & Kjaer) from mice. Vocalisations were sampled at 250 kHz with a CED Micro 1401-3 (Cambridge Electronic Design Limited, Sunnyvale, CA).

The samples were recorded using spike2, a "DC removal" filter was added to the recording. the files were then transferred to Matlab for splicing further analyses.

I know that there is a strict policy here against attaching files to questions, therefor I'm adding a link to the audio sample and the spectrogram

link to an audio sample and a spectrogram via google drive.

I took a look at your noise spectrum and it's pretty bad

You have strong lines at 20kHz and harmonics thereof. But also some spurs around  15 kHz, 18 kHz, 25 kHz and 44 kHz. The stuff at 20 you could potentially reduce with a comb filter, but your best bet (by far) is to clean up your measurement/recording setup.

There is no reason for this noise to be there, and you should try to identify the culprit(s) and eliminate them. Look at interconnects, other devices in the room especially lights or communication equipment, electric interference, power supplies, grounding issues etc.

This will be very hard to clean up, so having a good recording in the first place would be much preferable.


The difference of downsampling an image and smoothing an image? I am, at the moment, trying to read up on some simple computer vision elements, in which i have become a bit comfused on the terms downsampling and smoothing, and whether there is a difference between those two terms.

I know that downsampling means to reduce the resolution (DPI) of an image to a lower value. Which will result in a smaller image compared to the original.

But what about smoothing? Or convolution? Would it not result in the same thing? Or am i misunderstanding something?

Downsampling is the process of properly discarding every $(M-1)$ samples out of $M$ to reduce the signal sample rate by $M$. To avoid aliasing the signal being downsampled is initially lowpass filtered at a proper cutoff frequency before compression. The effect is a change of sample rate. Also the spectral positions of low frequencies are stretched to higher frequencies towards $\omega = \pi$. As a consequence, any spectra that passes beyond $\omega = \pi$ becomes aliased into mistreated bands, hence the requirement of a pre-processing lowpass anti-aliasing filter.

Smoothing is a more general thing. Whereas it may naturally refer to lowpass filtering in general, it may also mean more complicated things such as Kalman estimators which are also defined as smoothing operators. In the former case smoothing refers to getting rid of high frequency signal components, whereas in the latter case it refers to getting rid of randomness (usually of high frequency).


Edge following using Hough transform I am trying to improve an edge following algorithm developed by some students who did a project at my work. The algorithm is supposed to make a robot follow a line with the use of a camera. Their approach was to detect edges with Canny's algorithm and then extract lines from this image using the Hough transform. Three suggested lines will then be presented to the user, who may choose which line the robot should try to follow.

What I want feedback about is data association using the Hough transform. Is there any good criteria for matching of lines between frames? The used solution right now is to measure the angle and distance between followed line and extracted lines and to choose the one that is most similar.

One other problem I see is that they did not filter the signal, i.e. detected lines between frames. I have an idea about using Kalman filter to estimate the current parameters of the edge by using the extracted lines so that it is dynamic and less error prone. What do you think about that?

I have thought about making a SLAM algorithm instead but since we only want a proof of concept that will be a huge remake and requires much more data than I currently have.

1) Matching lines across multiple views is a common research problem and is reasonably well studied:
If you have the end points for example : http://cmp.felk.cvut.cz/~werner/software/lmatch/lmatch_memo.pdf

Even tough this deals with more geometric properties, if your scenes are well conditioned, you could as well use the gradient/intensity information to match the lines.

Here is a full source code on multiview line matching from monocular setups:
https://github.com/manhofer/Line3Dpp

You might also like to start with that.

2) Using a KALMAN filter on the lines would be a good idea. There are also more sophisticated algorithms. For example, in each frame you might want to minimize cost function to obtain the subpixel line, which has better resemblence to the line in the previous frame.

You can also track the camera pose and from there hypothesize the lines.

You try to avoid SLAM, but maybe, triangulating the lines and obtaining 3D structure would be the best in reducing the errors.


Digital FM Demodulation I have created a Common Lisp Library to receive samples of an FM signal from BladeRF1 SDR. I want to write a subroutine to implement FM demodulation. The input should be in the format of IQ data.

What is the algorithm for demodulating an FM signal?

Assuming you understand how the FM signal is theoretically generated, here is a simplified strategy to start with. In an FM Signal (https://en.wikipedia.org/wiki/Frequency_modulation), the frequency of the carrier is modulated according to the message signal (mostly audio). The instantaneous frequency $f_c(t)$

$ f_c(t) = k*m(t) + f_0$

where $k$ is a modulation constant. The phase of the carrier as a function of time is

$\phi(t) = \int f_c(t) = k t m(t) + f_0 t + \phi_0$.

Theoretically you FM signal have infinite BW but as per the Carson's rule most of the BW is concentrated around $f_c$. You have to make sure that the IQ samples had sufficient BW around $f_c$.
Once you extract the phase information from IQ samples, you can take first differential $\phi_n - \phi_{n-1}$ to get the instantaneous frequency information $f[n]$. The message signal $m[n]$ will be a proportionate of $f[n]$. With the knowledge of $k$ you should be able to get back (demodulate) the message signal. There are a lot of assumptions here though. You should be aware of $k$, Bandwidth of message signal, $f_c$ the center frequency etc.


Help in understanding the IQ raw file format used by SDRs/GnuRadio I apologise as I feel like I should know this off the top of my head. I will try to explain what I believe I know and then ask a simple question.

The usual file format used in software defined radios should be just numbers in binary format, representing in-phase and quadrature values of the sampled waveform. Those values are usually 32 bits, so every 2*32 bits you have a single sample of the waveform in that moment. When you play back this format, the DAC reproduces the waveform associated with those IQ states at the rate you ask it to. If you are using 10 MHz as your sample rate, the DAC will produce a waveform with the IQ values you ask it to and consume those values 10 million times a second.

The only way to entirely corrupt this file format is to cut 32 bits somewhere, as that would invert I and Q for the rest of the file. If you remove any multiple of 2*32 bits from the middle of the file you would only corrupt that single sample (and maybe the one before/after it), but the integrity of the file as a whole will not be compromised.

I'm asking this because I'm interfacing with the USB driver of a SDR and I'm confused as to what is actually going on.

As @Marcus_Müller points out, dropping samples will corrupt whatever digital symbol or analog waveform you were operating on at that instant.

You'll also introduce discontinuities in the time domain, which means you'll get broadband noise in the frequency domain at that time.

Probably the worst effect is, that if you drop a half sample, you end up swapping I & Q and the spectrum gets inverted.  Mathematically, if:

$$ x(t) e^{j\omega t} = I(t) + jQ(t) $$

then

$$ Q(t) + jI(t) = x(t)e^{-j\omega t +j\frac{\pi}{2}}$$

So you get spectral inversion about DC.  If your signal of interest was not centered and symmetric about DC, then your downstream DSP will yield junk.