Define comb function in terms of exponential function

The frequency domain representation of the discrete time downsampling operation (i.e.  $X_M(e^{j\omega})$ , the $DTFT\{x[Mn]\}$, in terms of $X(e^{j\omega})$, the $DTFT\{x[n]\}$ ) can be most easily observed by the following approach. We will begin from the end, $x[Mn]$ and proceed to its pre-compressed form by first deriving the frequency domain relation between the compressed signal $x[Mn]$ and its expanded form (pre-compressed) $x_s[n]$ and then we shall show the spectral relation between an input signal and its sampled form $x_s[n]$ (which is also the expanded signal at the input of the compressor)

First consider a periodic pulse train ,  $p_M[n] \triangleq \sum_{k=-\infty}^{\infty}{ \delta[n-kM]}$ whose period is $M$ and defined as :

$$p_M[n]=
\begin{cases}
1 ~~~~, ~~~ \text{n = 0, M, 2M, ..., kM , ...} \\
0 ~~~~, ~~~ \text{otherwise} \\
\end{cases} \tag{1}$$

Then the signal $x[n]$  will be first multiplied by this perodic pulse train $p_M[n]$, to produce $x_s[n]$ which will be selecting the required samples as follows:

$$x_s[n]=x[n]\cdot p_M[n] =
\begin{cases}
x[n] ~~~~, ~~~ \text{n = 0, M, 2M, ..., kM , ...} \\
0 ~~~~, ~~~ \text{otherwise} \tag{2}\\
\end{cases}$$

Now, we shall take every $M^{th}$ sample of $x_s[n]$ and discard the rest, which is called as the compression operation as graphically represented by:

$$ x_s[n] \longrightarrow \boxed{\downarrow M} \longrightarrow x_s[Mn]=x[Mn]=x_M[n] \tag{3}$$

Note that eventhough the conversion between $x[n]$ and $x[Mn]$ is not invertable, the conversion between $x_s[n]$ and $x_s[Mn]$ is so, as should be clear:

$$ x_s[Mn]=x_{sM}[n] \longrightarrow \boxed{\uparrow M} \longrightarrow x_s[n] \tag{4}$$

i.e. the inverse is the expander operator and we shall exploit it for easily showing the spectral relation between $x_s[n]$ and $x_M[n]$ as it can be shown that:

$$X_{s}(e^{j\omega}) = X_{sM}(e^{j\omega M}) \Longrightarrow  X_{sM}(e^{j\omega}) = X_{s}(e^{j\omega/ M}) \tag{5}$$

Based on this we can argue that (noting the equation (3)), $$ \boxed{ DTFT\{x[Mn]\} \triangleq X_M(e^{j\omega}) = X_{s}(e^{j\omega/ M}) } \tag{6}$$ where $X_{s}(e^{j\omega})$ is the DTFT of the sampled signal $x_s[n]$.

What comes next is the relation between the DTFTs of the signals $x_s[n]$ and $x[n]$ to complete the picture, which will be a straight forward observation if we consider representing the periodic impulse train by another form namelmy by its DFS (Discrete Fourier Series) representations (as is possible for any periodic discrete time signal with period M) :

$$ p_M[n] \triangleq \sum_{k=-\infty}^{\infty}{ \delta[n-kM]} = \frac{1}{M}  \sum_{k=0}^{M-1}{e^{j \frac{2\pi}{M}kn} } \tag{7}$$

Then we shall express $x_s[n]$ as :
$$ x_s[n] \triangleq x[n]\cdot p_M[n] = x[n] \big( \frac{1}{M}  \sum_{k=0}^{M-1}{e^{j \frac{2\pi}{M}kn} }\big) = \frac{1}{M}  \sum_{k=0}^{M-1}{x[n]e^{j \frac{2\pi}{M}kn} } \tag{8}$$

What comes next is using the linearity and modulation properties of DTFT operator and argue that:
$$DTFT\{x_s[n]\} = DTFT\{\frac{1}{M}  \sum_{k=0}^{M-1}{x[n]e^{j \frac{2\pi}{M}kn} } \} = \frac{1}{M}  \sum_{k=0}^{M-1}{ DTFT \{ x[n]e^{j \frac{2\pi}{M}kn} \} } = $$ $$ \boxed{X_s(e^{j\omega})=\frac{1}{M}  \sum_{k=0}^{M-1}{ X(e^{j \big(\omega -\frac{2\pi}{M}k \big)} ) } } \tag{9}$$

The final step is the mergal of the steps (6) and (9) to produce the result:
$$ \boxed{ X_M(e^{j\omega}) = X_s(e^{j\omega/M})=\frac{1}{M}  \sum_{k=0}^{M-1}{ X(e^{j \big(\frac{\omega}{M} -\frac{2\pi}{M}k \big)} ) } = \frac{1}{M}  \sum_{k=0}^{M-1}{ X(e^{j \big(\frac{\omega-2\pi k}{M} \big)} ) } } \tag{10} $$

I hope it was a clear way of expressing the spectrum of downsampled signal $x_M[n]=x[Mn]$ to that of the original signal $x[n]$

Also note that a more complete understanding of the spectral relation would only be possible if you fully explore the relation by a graphical plot as well, which would clearly demonstrate the resulting frequency stretching towards the $\omega  =\pi$ boundary, by a factor of $M$ after which a consideration of aliasing becomes clear...


This plot depicts how to convert your digital signal back to the analog one, using $\mathrm{sinc}$ functions. The nice property of these functions used in this process, is that maximum of each function always occurs at minimums of the other, shifted function:

Now the process of the D/A conversion that is depicted on your plot is basically taking a $\mathrm{sinc}$ function, and multiplying it by one sample that is at the maximum of $\mathrm{sinc}$. You do it for each and every sample of your digital signal, and then sum all of the outputs. For example, for the above plot, the outcome would be following:

Now maybe some more complicated example of D/A conversion. On the plot below, you can see the actual samples (in blue) and "analog" $\mathrm{sinc}$ interpolated signal (in red).

Plotting each function separately you can observe what will become an analog signal later:

I hope that it is more intuitive for you now. In order to complete this example, here is some MATLAB code I used:

clc, clear all, close all
fs = 32;
dt = 1/fs;

% Generate a samples of the signal
ts = linspace(0, 1, 10);
s = 2*sin(2*pi*ts) + sin(2*pi*3*ts);

td = linspace(0, (length(s)-1)*dt,     length(s));
t  = linspace(0, (length(s)-1)*dt, 100*length(s));

sincM = repmat( t, length(td), 1 ) - repmat( td', 1, length(t) );
SS =  s* sinc( sincM/dt );

% Sampled signal
stem(td, s, 'b', 'linewidth', 2)
hold on
% Reconstructed signal
plot(t, SS, 'r', 'linewidth', 2)
grid on

legend({'sampled signal', 'interpolated signal'})

figure
stem(td, s, 'b', 'linewidth', 2)
hold on
grid on
% Each sinc function
basis = sinc(sincM/dt)';
sbasis = bsxfun(@times, basis, s);
plot(t, sbasis)


Determining a system's causality using its impulse response

Clearly, for negative values of $t$, the system needs to know the future in order to determine its output. Hence, the system can't be causal.
Since the system is also time-varying (show it!), its response to an impulse doesn't say much about its general behavior, unlike it would be the case for a linear time-invariant (LTI) system. So the given system's response to an impulse being zero for $t<0$ doesn't say anything at all about causality or any other general properties of the system.
Takeaway: don't use a system's impulse response for drawing conclusions about its properties before having verified that the system is LTI and can in fact be characterized by its impulse response.


How can resonance be added to a fractional order Butterworth low pass

Given the formula for second order sections:

$$
s_k = \omega_c e^{\frac{j(2k+n-1)\pi}{2n}}\qquad k = 1,2,3,\ldots, n \qquad s = i\omega
$$

$$
H(s)=\frac{G_0}{\prod_{k=1}^n (s-s_k)/\omega_c}
$$

Use k = 1 for the resonance section, add $Q$ in the normal way i.e. calculate $H_{k=1}^2 = H_{k=1}(s) \overline{H_{k=1}(s)}$ where $s = i\omega$ and the denominator will be of the form $w^4 + cw^2 + 1$ where $c$ is a constant calculated from $s_k$. We add $Q$ to the first order term (note that this is a quadratic in $\omega^2$). Using any positive real $n$ in $s_k$ gives valid values for $c$.

Let $R^2$ be a modified $H_{k=1}^2$ with $Q$ added in to the first order term of the quadratic in $w^2$ then to add resonance to the original equation in the question multiply it by:

$$
\frac {R^2}{H_{k=1}^2}
$$

This effectively removes the $k = 1$ second order section and adds a resonant section back in. It turns out once you do this, $n$ can be any positive real value when used in $c$ but note this wouldn't work for fractional roots of the sections (you would have to use fractional calculus if you wanted the roots of the fractional equation).

Without resonance:

$$
c = 2 \cos^2(\frac{\pi (1+n) } {2n}) - 2
$$

With resonance:

$$
c = \frac{2}{Q^2}\cos^2(\frac{\pi (1+n) } {2n}) - 2
$$


Mid/Side encoding/decoding in frequency domain instead of time domain

Yes, it is possible. The discrete Fourier transform (DFT, of which the FFT is an efficient implementation) is a linear transformation, i.e.,

$$
a\cdot\text{DFT}\{x\}+b\cdot\text{DFT}\{y\}=\text{DFT}\{a\cdot x+b\cdot y\}.
$$

The same holds for the inverse transform.

As a consequence, you can directly convert the equations for the mid and side signals from the time domain to the frequency domain. Simply replace the subscript $t$ with $f$:

$$
M_f = \frac{L_f + R_f}{2}\\
S_f = \frac{L_f - R_f}{2}\\
\vdots
$$

From $M_f$, $S_f$ you can then compute the inverse discrete Fourier transform to obtain $M_t$, $S_t$.


A DFT produces a finite number of results (frequency bins), so the complex result vector can be indexed by 0..(N-1).

A DTFT produces a result over a continuous range of frequencies (more than a finite number of results) so must be specified by a real number as opposed to an integer.

As MBaz says, you can sample the DTFT at a fixed integer number of points.


Why does my I/Q data from radar are not centered?

The IQ data from an FMCW radar would have the frequency difference from the delay between the transmitted and reflected signal given the transmitted signal is a changing frequency with time; so the result is what should be expected. Given the ramp rate of frequency the round trip delay from target can be determined from the frequency of this result.

The result here appears to have a stronger higher frequency AND lower frequency tone as the result of two different targets. The lower frequency that is superimposed would be due to the closer target.


There's a lot of confusion in your question. You refer to "a channel $H(s)$", then to "my signal $H(s)$", so for the sake of clarity, let's properly define:

Input signal - continuous time $f(t)$, impulse sampled $f[n]=f(n\Delta T)$, Fourier transform of continuous time signal $F(s)$, Z-transform of impulse sampled signal $F(z)=F(e^{s\Delta T})$.
Continuous time channel - impulse response $h(t)$, transfer function (Laplace transform of impulse response) $H(s)$.
Discrete time channel - discrete impulse response $h[n]$, transfer function (Z-transform of discrete impulse response) $F(z)$.

Note the definitions above follow a convention:

time-domain signals/functions are lower case ($f, h$); frequency-domain signals/functions are upper case ($F, H$).
arguments of continuous signals/functions are enclosed in parentheses $(t)$; arguments of discrete signals/functions are enclosed in brackets $[n]$.

Then, you write:

  (...) z-transform for my channel but am unsure how I can go about the convolution with my sampled signal which is just an array of numbers

It seems you do not understand the basics of LTI systems. Note that $H(s)$ and $H(z)$ are frequency-domain and $f(t)$ and $f[n]$ are time-domain.

Convolution makes sense when the two convolved functions are in the same exact domain (continuous/discrete time/frequency), although for the usual purposes the two convolved functions should both be in the continuous/discrete time domain.

If you're dealing with a discrete-time signal ("signal which is just an array of numbers"), then you might want to convolve $f[n]$ with $h[n]$. If you already have $H(z)$, then an inverse Z-transform on $H(z)$ would produce $h[n]$.


Augmentation for EEG signal classification using Deep Learning

...if we are using augmentation for EEG classification and say we generate spectrograms and feed them as images to a CNN

If you are approaching this as an image based problem, it could possibly work. There are not enough details about the problem to provide a more accurate answer on this.

In any case, EEG records brain activity typically as a set of signals. Later, each one of those signals would lead to a time-frequency spectrogram representation. This means, one "image" per EEG channel which does not seem very practical from an augmentation point of view.

  ...can we use Augmentation to increase the training dataset?

A better approach to "augmentation" would be to generate surrogate data via suitable models.

This is closer to "augmentation" but still does not solve all of the problems you might be facing.

Jansen's model (for example) will generate realistic EEG waveforms with similar dynamics but only on a per channel basis (out of the box). If you are interested in replicating an EEG recording as a whole (that is, not only the dynamics of each channel but also their correlations), you would need to couple instances of Jansen's model and at that point you are still left with the question of "which coefficients are supposed to be associated with this task"?

So, if you are interested in creating a classifier that discriminates between "This is most likely an EEG signal rather than an electrode just dangling in the air", augmenting your dataset will work as expected. But if you are interested in creating a classifier that discriminates between "The subject is now more likely day dreaming rather than sleeping", then you need either very complex models or a large amount of data (under a similar condition) to infer some kind of a rule to exploit it later.

Hope this helps.


Analog vs. digital communications in PCM and noise

This is quite a broad topic.  Let's narrow the scope to toll quality telephone audio.  Further assume both analog and digital channels are bandlimited to the same passband.  I assume you don't want to figure out the difference between a pair of wires and a PCM channel (i.e. a pair of wires has a very slow decay in frequency response where as a digital channel is hard band limited).

Ok a for sort of apples to apples comparison of a roughly 4Khz voice band channel, you can look for sampling artifacts in the channel.  The sampling frequency for telephone audio is 8Khz.  The true passband for telephone quality audio is 200Hz-3600Hz (still roughly speaking here).

Look for quantization noise.  Quantization noise is not random and can be made nearly periodic by choosing the right test signal.  The PCM encoding for standard CODEC ditigal telephone channels uses a type of compression that follows a log scale (u-law or A-law).  The analog reconstruction filter does not completely eliminate the quantization noise.  Because of the compression scheme, the quantization steps increase with increasing signal amplitude.  If you create a test signal that is a ramp running from -Vp to +Vp (where Vp is the maximum signal amplitude the channel can handle), you will see discontinuities that are closely spaced near 0V and spaced farther appart near +Vp and -Vp.  Use a frequency of 1KHz so the signal frequency is harmonically related to the sampling frequency.  This will place the quantization errors nearly at the same location for each period of the ramp signal.

You can also use a spectrum analyzer and look for spurs in the frequency content of a test signal.  In this case, use a 1Khz sine wave.  Once again the quantization errors will remain at nearly the same position on each period of the sinewave.  The quantization errors will appear as distortion of the sine wave and show up as harmonics in the signal spectrum.  These harmonics will be greatest when the sampling frequency is harmonically related to the signal frequency so 1KHz and 2KHz are good choices.  As the test tone drifts away from intger divisions of the sampling rate, the quantization noise spreads and begins to look like a spectral contribution from random noise.  In telephone systems, a 1KHz test tone is generally not used for this reason.  The normal test tone frequency is 1004Hz.

Note that modern CODECS with high compression rates such as G.729 have a difficult time passing periodic signals.  These CODECs are designed in such a way as to interpret all signals as human voice or noise, so the reproduced periodic signal will be highly distored or modulated.


Difference between cv::subtract and operator - in cv::Mat OpenCV


Discrete Time Fourier Transform (DTFT) cross correlation property

Your first formula in your question is generally wrong, that's why you can't prove it. The correct formula is

$$\sum_{n=-\infty}^{\infty}x[n]y^*[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})Y^*(e^{j\omega})d\omega\tag{1}$$

which is just Parseval's theorem.

If $x[n]$ and $y[n]$ are real-valued, $(1)$ can be written as

$$\sum_{n=-\infty}^{\infty}x[n]y[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})Y(e^{-j\omega})d\omega\tag{2}$$

which is the form you've correctly obtained.

In the special case that $y[n]$ is not only real-valued but also even, i.e., $y[n]=y[-n]$, then $Y(e^{j\omega})=Y(e^{-j\omega})$ holds, and $(2)$ turns into the first formula in your question.


How does sampling rate impact Discrete-Time Kalman Filter state space modeling assumptions?

The answer to your question is : maybe. When one is at the discrete time model, it is generally assumed that you solved the state transition equations from the continuous time differential equations where the physics have already been appropriately modeled, which is where your answer lies. Time discretization is based on sampling the continuous time state transition matrix, thus the sampling intervals aren't related to Nyquist, they correspond to the fidelity of the modeling. Shorter intervals tend to be more accurate approximations. The time intervals don't have to be uniform, but uniform intervals provide a less complicated implementation.

While shorter time intervals correspond to a heuristic where shorter has better approximation, you raise a valid concern, sampling too closely, may indeed introduce correlation in the process noise.  This can be solved with an augmented state.

These are all considerations. Kalman Filters need to be tuned.  The measurement noise has been neglected in this discussion.

So, the answer to your question is that it really depends.


how to get the location long and lat out of the raw measurements?

I only know the vastly simplified form: all satellite nav systems work by the satellites telling you where they are and what time it is.  In GPS terms, that time is the "pseudorange", and it's all you need (in theory) to locate yourself.

With no satellites in view, you have no clue where or when you are.  With one satellite in view, you can place yourself on a 3-sphere* in time and space**.  With two satellites you can place yourself on a 2-sphere, etc.  When you have four satellites in view you can locate yourself to a point, and you're done***.

I suggest you look around and get a book on satellite navigation.  I'd love to recommend one, but first, I'm not current, and second, StackExchange discourages product recommendations.  But they're out there, and these days there should be one that covers GPS, GNSS and Galileo.

* I'm using the mathematical terminology here: a 3-sphere is the locus of points equidistant from a point in 4-dimensional space, a 2-sphere is the surface of a ball in 3-D space, a circle on a plane is a 1-sphere, and a point is a 0-sphere.

**  I.e., you know that if it's the same time as the satellite's time you're at the satellite's position, if it's 1ms of a second after the satellite's time you're on a 2-sphere 300,000 meters in diameter around the satellite's position, etc.

*** Except for noise, of course, and I can't recall if, even in the perfect case, you end up with just one possible location or a finite number of them.


Differences between two closed-loop systems

There are more than just one way to approach the solution:

1) You can find the steady state error by finding the equation of the denominator of the transfer function.

In system 1:

$$D(s) = a(s)p(s) + b(s)q(s)$$

The order of the system is determined by this order of $D(s)$. Solving for $D(s)$ will give you the poles of the system. If the poles 'solutions of $D(s)$' all are negative. That means if all the values of the solution to $D(s)$ are in the left half plane $<0$. This means the system is stable.

2) The error in time domain is defined as "input - output".
That is:

$$e(t) = y_r(t) - y(t)$$

Now, to compute the steady state error, that is time at infinity 'after a long time'.
The equation will be: $$ e_{ss} = \lim_{t \to \infty}  y_r(t) - y(t) $$
That is simply taking the limit as t goes to infinity.

Now, using the final value theorem to compute the steady state error in the $s$ domain the equation will be:

$$E_{ss} = \lim_{ s \to 0} E$$

That is in your terms

$$E_{ss} = \lim_{s \to 0} s\cdot e_{ss1} \quad(1)$$

$$E_{ss} = \lim_{s \to 0} s\cdot e_{ss2} \quad(2)$$

Therefore; you can note that the result depends on the 'Type of the system and the Input'.
Type of the system: is how many poles at zero do you have?
That is writing:

$$a(s)p(s) + b(s)q(s) = s^n (s+a_0)(s\cdot a_1) \cdot ... \cdot(s\cdot a_{n-1})$$

the value n represents the order of the system.

Type zero means no poles to zero.
Type I means one pole at zero, etc.

Now if the input is a unit step that is $$y_r(s) = \dfrac{1}{s}$$
this will cancel out with the $s$ as specified in equations $(1)$ and $(2)$. The equation will be:

$$ e_{ss1} = \lim_{s \to 0} \dfrac{r(s)a(s)}{a(s)p(s)+b(s)q(s)} $$

Now, looking at the type of system:

If it is of type one or higher, it is clear now that the $e_{ss1}$ goes to
infinity yielding unstable system.
If type zero, it will have a certain value, yet still stable.

These are the techniques, and the answer highly depends on the type of the system or even the characteristic equation that is the equation of the denominator of the transfer function.


Time domain maximum from frequency domain data?

Suppose that Alice has a vector $\mathrm x \in \mathbb R^n$. She computes the DFT of $\mathrm x$

$$\mathrm y := \mathrm F \mathrm x \in \mathbb C^n$$

where $\mathrm F \in \mathbb C^{n \times n}$ is a Fourier matrix. Alice then tells Bob what $\mathrm y$ is. Since the inverse of the Fourier matrix is $\mathrm F^{-1} = \frac 1n \, \mathrm F^*$, Bob can recover $\mathrm x$ via

$$\mathrm x = \frac 1n \, \mathrm F^* \mathrm y$$

and then compute $\| \mathrm x \|_{\infty}$ to find the maximum absolute value of the entries of $\mathrm x$. What if computing matrix inverses and Hermitian transposes is not allowed? Bob can then write $\mathrm F$ and $\mathrm y$ as follows

$$\mathrm F = \mathrm F_{\text{re}} + i \,\mathrm F_{\text{im}} \qquad \qquad \qquad \mathrm y = \mathrm y_{\text{re}} + i \,\mathrm y_{\text{im}}$$

and, since $\mathrm x \in \mathbb R^n$, the equation $\mathrm F \mathrm x = \mathrm y$ yields two equations over the reals, namely, $\mathrm F_{\text{re}} \, \mathrm x = \mathrm y_{\text{re}}$ and $\mathrm F_{\text{im}} \, \mathrm x = \mathrm y_{\text{im}}$. Bob can then solve the following linear program in $t \in \mathbb R$ and $\mathrm x \in \mathbb R^n$

$$\begin{array}{ll} \text{minimize} & t\\ \text{subject to} & - t 1_n\leq \mathrm x \leq t 1_n\\ & \begin{bmatrix} \mathrm F_{\text{re}}\\ \mathrm F_{\text{im}}\end{bmatrix} \mathrm x = \begin{bmatrix} \mathrm y_{\text{re}}\\ \mathrm y_{\text{im}}\end{bmatrix}\end{array}$$

which can be rewritten as follows

$$\begin{array}{ll} \text{minimize} & \begin{bmatrix} 1\\ 0_n\end{bmatrix}^{\top} \begin{bmatrix} t\\ \mathrm x \end{bmatrix}\\ \text{subject to} & \begin{bmatrix} -1_n & \mathrm I_n\\  -1_n & -\mathrm I_n\end{bmatrix} \begin{bmatrix} t\\ \mathrm x \end{bmatrix} \leq \begin{bmatrix} 0_n\\ 0_n\end{bmatrix}\\ & \begin{bmatrix} 0_n & \mathrm F_{\text{re}}\\ 0_n & \mathrm F_{\text{im}}\end{bmatrix} \begin{bmatrix} t\\ \mathrm x \end{bmatrix} = \begin{bmatrix} \mathrm y_{\text{re}}\\ \mathrm y_{\text{im}}\end{bmatrix}\end{array}$$

and not only recover $\mathrm x$ but also obtain $t = \| \mathrm x \|_{\infty}$. However, is solving a linear program cheaper than computing a Hermitian transpose?

MATLAB code

The following MATLAB script

n = 8;

% build n x n Fourier matrix
F = dftmtx(n);

% -----
% Alice
% -----

% build vector x
x = randn(n,1);

% compute DFT of x
y = F * x;

% ---
% Bob
% ---

% solve linear program
c = eye(n+1,1);
A_in = [-ones(n,1), eye(n);
        -ones(n,1),-eye(n)];
b_in = zeros(2*n,1);
A_eq = [zeros(n,1), real(F);
        zeros(n,1), imag(F)];
b_eq = [real(y); imag(y)];
solution = linprog(c, A_in, b_in, A_eq, b_eq);

% extract t and x
t = solution(1);
x_rec = solution(2:n+1);

% check results
disp('t = '); disp(t);
disp('Infinity norm of x = '); disp(norm(x,inf));
disp('Reconstruction error = '); disp(x_rec - x);

produces the output

Optimization terminated.
t =
2.2023

Infinity norm of x =
2.2023

Reconstruction error =
1.0e-013 *

0.0910
0.0711
0.0167
-0.1077
0.1049
0.0322
0.1130
0.2776

The original vector is

>> x

x =

-1.1878
-2.2023
0.9863
-0.5186
0.3274
0.2341
0.0215
-1.0039


Why do we get ripples in the frequency spectrum of a chirp, and how do I remove it?

Another thing the OP can do is use the "complex" chirp. What you really have to do is (at two different times and both times sweeping all frequencies from -Nyquist through DC to +Nyquist) pass the real part, $x_r(t)=\cos(\pi \beta t^2)$ and then pass the imaginary part $x_i(t)=\sin(\pi \beta t^2)$ of the complex chirp through the material under test and synchronously measure the two real results ($y_r(t)$ and $y_i(t)$), then assemble them as the real and imaginary parts of the response.

Let
$$ \begin{align}
x(t) & =  e^{j \pi \beta t^2} \\
& = \cos(\pi \beta t^2) + j \sin(\pi \beta t^2) \\
& = x_r(t) \qquad + j \ x_i(t) \\
\end{align} $$

Then
$$ X(f) \ \triangleq \ \mathscr{F}\left\{ x(t) \right\} \ = \ \sqrt{\frac{j}{\beta}} e^{-j\frac{\pi}{\beta} f^2} $$

You will see that $ |X(f)| = \sqrt{\frac{1}{\beta}} $ for all $f$.

The respective outputs that are synchronously measured with the two driving signals are $y_r(t)$ and $y_i(t)$:

$$ \begin{align}
y_r(t) & \triangleq h(t) \circledast x_r(t) \\
 & = h(t) \circledast \cos(\pi \beta t^2) \\
\\
y_i(t) & \triangleq h(t) \circledast x_i(t) \\
 & = h(t) \circledast \sin(\pi \beta t^2) \\
\end{align} $$

Then this aggregate complex output, $y(t)$, can be assembled (in the mind of the computer) from the responses of the cosine sweep (the real part) and the sine sweep (the imaginary part):

$$ \begin{align}
y(t) & =  y_r(t) + j \ y_i(t) \\
& = h(t) \circledast x_r(t) + j \ h(t) \circledast x_i(t) \\
& = h(t) \circledast \left(x_r(t) + j \ x_i(t)\right) \\
& = h(t) \circledast x(t) \\
& = \int_{-\infty}^{+\infty} h(u) \, x(t-u) \ du \\
& = \int_{-\infty}^{+\infty} h(u) \, e^{j \pi \beta (t-u)^2} \ du \\
& = \int_{-\infty}^{+\infty} h(u) \, e^{j \pi \beta (t^2-2tu+u^2)} \ du \\
& = e^{j \pi \beta t^2} \int_{-\infty}^{+\infty} \left(h(u) \, e^{j \pi \beta u^2}\right) \, e^{-j 2 \pi \beta t u} \ du \\
& = e^{j \pi \beta t^2} \int_{-\infty}^{+\infty} \tilde{h}(u) \, e^{-j 2 \pi (\beta t) u} \ du \\
& = e^{j \pi \beta t^2} \tilde{H}(\beta t) \\
\end{align} $$

where
$$ \tilde{H}(f) \triangleq \int_{-\infty}^{+\infty} \tilde{h}(t) e^{-j 2 \pi f t} \ dt = \mathscr{F}\left\{\tilde{h}(t) \right\} $$

and
$$ \tilde{h}(t) \triangleq h(t) e^{j \pi \beta t^2} $$

So this is what you get: first consider making this adjustment (dechirping) in the measured output:
$$ y(t) e^{-j \pi \beta t^2} = \tilde{H}(\beta t) $$

Now if the sweep rate $\beta$ is low enough, then $e^{j \pi \beta t^2} \approx 1$ for as long as $h(t)$ is not close to 0. Then $\tilde{h}(t) \approx h(t)$ and $\tilde{H}(f) \approx H(f)$ and $\tilde{H}(\beta t) \approx H(\beta t)$ and what you have measured in the output is another chirp, but with relative magnitude and phase changed by $\tilde{H}(\beta t) \approx H(\beta t)$.

Even if the sweep rate is not so low, you can compute either the true frequency response $H(f)$ by first computing

$$ \begin{align}
\tilde{h}(t) & =  \mathscr{F}^{-1}\left\{ \tilde{H}(f) \right\} \\
& = \int_{-\infty}^{+\infty} \tilde{H}(f) e^{j 2 \pi f t} df \\
& = h(t) e^{j \pi \beta t^2} \\
\end{align} $$

Or from that result you can get the true impulse response $h(t)$ and then $H(f)$ if that's what you want. It should be the same $H(f)$ you would get if you divided $Y(f) = \mathscr{F}\left\{ y(t) \right\} $ by $X(f) = \mathscr{F}\left\{ x(t) \right\} = \sqrt{\frac{j}{\beta}} e^{-j\frac{\pi}{\beta} f^2}$. And there is no division by zero, because there are no ripples in $X(f)$ because the magnitude is a constant $\sqrt{\frac{1}{\beta}}$.

That is how to remove ripples (and any danger of division by zero) in the spectrum of a linearly-swept chirp.


How can I compare two pictures from the camera and tell if there are enough differences to detect movement?

It seems to me what you are looking for is background subtraction technique. With noisy images and changing lighting conditions it could be nontrivial. The current state of the art technique for this is the low-rank matrix representation, but it require not two but many (~dozen) images. Some heavy duty math follow:
Each image considered as  a vector of pixels, vectors combined into matrix and this matrix is decomposed into low-rank matrix and remnant. Low-rank matrix columns are backgrounds and remnant is moving objects. There are some open-sourced implementations, but only for factorization itself, not complete image pipeline IIRC

Here is a paper and code for matrix factorization
http://www.ece.rice.edu/~aew2/sparcs.html

Here is overview from CS blog and link to other code:

https://sites.google.com/site/igorcarron2/matrixfactorizations

http://nuit-blanche.blogspot.com/search/label/MF

Survey of some other techniques:
http://www.vis.uni-stuttgart.de/uploads/tx_vispublications/Brutzer2011-2.pdf


Digital filters with more zeros than poles

If you consider the transfer function of a causal IIR filter

$$H(z)=\frac{B(z)}{A(z)}=\frac{\sum_{m=0}^M b_mz^{-m}}{\sum_{n=0}^N a_nz^{-n}},\quad a_0=1$$

then you always get the same number of poles and zeros, regardless of the choice of $M$ and $N$ (as already pointed out by Robert). However, what is meant by a system with "more zeros than poles", is a system with a numerator degree $M$ greater than the degree $N$ of the denominator. In this case "poles" refers to the poles away from the origin. It is only these poles that have to be implemented. So a system with some poles at the origin is cheaper to implement than a system with all its poles away from the origin. So your argument that you "might as well use those poles ... if you have them" does not really hold, because poles at the origin come without any implementation costs.

The remaining question is why such filters with more zeros than poles (away from the origin) can be useful. Let me just give two examples:

In implementation as well as in the design process, poles can give a lot of trouble (instability, noise enhancement, local optima in the design process, etc.). On the other hand, FIR filters with good frequency selective properties must often have a very large degree (resulting in a large delay and in high implementation costs). Adding just a few poles can give us the best of both worlds: better filter behavior at a greatly reduced (numerator) degree $M$, and only little trouble because there are only very few poles.
For frequency selective filters with a small desired phase distortion in the passbands, IIR systems with $M>N$ are very useful, because the poles ideally contribute only to the passbands, whereas the zeros must contribute to the passbands (for phase equalization) as well as to the stopbands (zeros on or close to the unit circle). Consequently, we need poles only at angles corresponding to passband frequencies, but we need zeros everywhere.

The figure below shows such a system (lowpass filter with approximately linear passband phase response, $M=12, N=6$):


How to perform a block-wise frequency down-shift?

As an aside, are there other caveats I should be aware when performing
  a down-shift, that I haven't accounted for here?

Yes there's one and let me put one of the most important caveats when performing frequency modifications. An acoustic sound signal (in the form of a music, speech etc) will physically exist in a so called harmonic structure (for more indepth discussions see sound synthesis, sound analysis etc.) whose immediate implication is that the perception of the frequency content is in octaves rather than linear frequency decades (an octave is the span of frequencies between $f_0$ and $2 f_0$).

Now, eventhough for most technical applications (such as in those communication algorithms) a band of frequencies can be shifted to any desired frequency without a loss of information (which forms the essence of the modulation theory) when it comes to audio perception this does not hold.

It's easy to see that the harmonic structure of the audio band will be distorted when a complete band is shifted to other frequencies in a linear manner. As a solid example consider a primitive audio waveform composed of $440$ Hz, $880$ Hz and $1320$ Hz (these are $f_0$, $2f_0$ and $3f_0$ of an absurd instrument, more natural instruments will have more complicated structures including inharmonicities etc.)

Now if you linearly shift this waveform base from $f_0$ to $1.5 f_0$ then the new frequency distribution will be $1.5f_0$ , $2.5f_0$ and $3.5f_0$ which are no more harmonics. This creates an important problem in natural pitch shifting applications.


Discrete Hilbert transformation of a spatial signal

You need the Fourier transform of the impulse response of the Hilbert transformer

$$h(t)=\frac{1}{\pi t}\Longleftrightarrow H(j\omega)=-j\textrm{ sign}(\omega)$$

So all positive frequencies are multiplied by $-j$ and all negative frequencies are multiplied by $j$. This corresponds to a phase shift of $-\pi/2$.

EDIT:

I'm referring to your updated question including the code. I can't see if the computation of the wavenumber is OK, but one problem I definitely see is in applying the phase shift. You change A[i].r and then you re-use it (as if it were the orignal value) to compute A[i].i. Furthermore, I think the signs are wrong. If signfun=1, A[i].r=A[i].i and similarly for the computation of A[i].i. But the most important thing is to use a temporary variably to store the old value of A[i].r for computing A[i].i.


Weird ringing at signal start and end for default Matlab lowpass filter

You are using a linear phase filter. These are non causal and hence you get pre-ringing at your transients. Use a minimum phase or IIR filter. Do NOT use filtfilt, since that makes an IIR filter linear phase again and you end up with the same problem
The artifact at the beginning is the step response of the filter. The filter needs access to "past" samples and since it doesn't know what these are, it assumes they are 0. In your case, you have a DC jump and it looks like they are more like -0.05. There multiple work-arounds. In your case it's probable the easiest to pre-pend the signal with a few hundred constant samples that are the average of the first few hundred sample of your signal, filter the whole thing and cut the pre-pended section off again.


Time-Frequency Analysis of Big Data - Data Size Reduction: averaging the most appropriate method?

So, first of all, CSV seems to me the least suitable format imaginable for this amount of data. It needs to be parsed, is memory hungry, and wastes precision, and isn't linearly addressable (ie. to get to the 99999. element, you need to parse the preceeding 99998 elements).

So I'd recommend keeping the structure from these files, but converting them to binary files, eg. numpy arrays, HDF5, whatever, but keep it binary, parse-free. Convert once and keep that data as what you load – it will reduce your file loading times immensely, and your software (whatever you're using) possibly doesn't even need to load all the file at once into memory, because when the data format has a fixed bit width per element (e.g. always 4-Byte integers and 8-Byte floating point numbers), it's trivial for the software to just load the part of the file into memory that you need. In fact, operating systems do that for you ("mmaped files").

Reducing your text file format to e.g. an array of single-precision floats or 32bit integers, your total data size would be a mere 8GB – and that fits multiple times in the RAM of every modern workstation, so I'd argue it's not really "big data" anymore! (Just to take the fear of handling so many numbers away)

I've ranted about this before, and people have answered that "CSV files are human-readable". I dare anyone who says so to read a 100k-column, 500 row CSV file and justify their impression of the data as being representable for the whole set.

  I would like to know what is the best method to make sure that those 1-4 columns represent the total data set in each file?

As always in engineering, there's no "best" without you defining a measure for good. If you want to highlight the differences among different frequencies, pick the column e.g. that has the highest sample variance. If you want to represent the whole measurement "fair", don't pick a single column, but average all columns. It's impossible to answer that question without being you – it's you that wants to show something. What you have is data – meaning is something that you humanly attribute (and, yes, you must be aware that the choice of data reduction method is your biased decision, and you should definitely communicate that).

  Is averaging simply a good method in this case?

It's the first thing I'd try to see if I can still see something interesting.
That's a pure coincidence – I might as well just pick a random one!

  What are the alternatives?

Impossible to tell, again – picking one column based on some metric, calculating variance, third statistical moment, postulate a certain functional shape would be handy and fit that shape to your data set … The sky is the limit, and what you want to find out is your only guideline.


Calculating the peak power of a signal on python

Peak power is simply the highest square of any of the samples you encounter.
Average power is simply the average square.


GPS C/A Code Time Delay Using Ultrasound

Instead of using "square" for your carrier generation, use "cos".

Also the -1 square in your receiver will not generate a quadrature signal from + square as you have used, copied below:

This does not produce a sine and cosine:

x2 = square(2*pi*fc*t); %sin
x3 = -1*square(2*pi*fc*t);%cos

Here as well sin and cos directly would be more appropriate, but the following is how I would approach your receiver, using exp instead of sines and cosines since $e^{-j\omega t} = cos(\omega t)-jsin(\omega t)$:

modulated =  send();  %generate incoming signal at carrier fc
rx = modulated * exp(-j*2*pi*fc*t)   %downconvert to baseband
%the above will have a high frequency image at -2fc that must be filtered
coeff= firls(31, [0 .2 .5 1], [1 1 0 0]);
rxfilt= filter(coeff,1,rx)
corrOut = ifft(fft(rxfilt).*conj(fft(CA_local));

The $e^{-j\omega t}$ will shift the positve half of your modulated signal to baseband (0) and the negative half to $-2f_c$, so the example filter shown is to filter out the image. (Note there was no consideration given for actual filter design, this was quickly "pulled out of the hat" just to demonstrate what will work).


Estimating signal strength from PCM-encoded signal

In general, it is very hard to estimate the power of a signal after it has gone through an analog front-end and sampled. The front-end filters and mixes the signal, which can amplify or attenuate it. It's also common to have an automatic gain control (AGC) stage, which adjusts the signal's amplitude to match the ADC's input range and maximize dynamic range. In many cases, the digital back-end will also normalize the signal amplitude to a pre-defined range, such as -1 to 1.

If you understand every piece in the system very well, you may be able to estimate the signal's power from its samples. Another alternative is to use a scope or a spectrum analyzer to measure the signal as close as possible to the receiver's input.


price of dsp kit TMS320c6678(used for signal processing)

Yes. If you buy at least 1000 units, each will cost you 120


How to derive low pass filter $\frac{N}{N+S}$

I assume that you want to see that the transfer-function magnitude is inverse proportional to the frequency. In this case, you replace $S$ with $i\omega$ and you get that the transfer function magnitude decays with the frequency$|G(\omega)|=\left|\frac{n}{ n+i\omega}\right|$


what does "fixed multipath channel model" means?

Your "h" is a valid multipath channel impulse response.  It is a special case, though, in that it represents a multipath delay that is an exact multiple of your sample period.  Usually the delay will not be an exact multiple of the sample period.

I am not an expert on multipath channels, so hopefully someone else can chime in, but I think that the easiest way to generalize your multipath channel is to make the gain of the multipath complex rather than real.  The phase of the gain would alter the phase of your signal, simulating a non-integer delay.

If you really wanted to be hardcore about it you could use this answer on how to make fractional delay filters to create a filter that delays the signal by as much or as little as you want, and then add it to [1 0 0 0 0...] (the LOS channel) to get the complete multipath channel.


Applying a window function to a speech signal

I don't really understand what "apply a suitable windowing function" "prior to the STFT" means. Most STFT implementation incorporate apodizing windows. And the best is always related of some cost function, objective or processing.

Note that there are people called Mr Hanning, but they are not (AFAIK) at the origin of the raised-cosine or von Hann window (Julius von Hann, Austrian meteorologist). You can find a mention at page 97 from the book  The measurement of power spectra from the point of view of communications engineering (archive.org), 1958, by Blackman and Tuckey. Noteworthy, they spell hamming and  hanning in small letters.

Hamming or Hann windows are very standard. $50\%$ or $75\%$ overlap are common. In  the recent works I know of, especially for blind source separation,  people often use at least two window sizes: a shorter one for transients, and a longer one.

This is a first step toward using different windows depending on the scale of observation. Some have used a multiscale context, such as in the
ERBlet non-stationary Gabor filterbank, ERB for Equivalent Rectanguar Bandwidth, with a constant-Q transform:

Some have been using asymmetric windows, better suited to the skewed shape of speech onsets (see eg Using asymmetric windows in automatic speech recognition).

The choice of the best windows, to me, even for speech, is not  clearly settled.


Determine filter given transfer function

No $$H(z) = 1 + \frac{\alpha}{1-\alpha z^{-1}}$$ is not a lowpass filter. It's a lowpass filter in parallel with an all pass filter (kind of low-boost shelving filter). If you want a lowpass filter which rejects the high frequencies then you may take the following:

$$H_{lp}(z) = \frac{\alpha}{1-\alpha z^{-1}}$$ which is a lowpass filter for the range of parameter $\alpha$. And therefore $1 - H_{lp}(z)$ will be the complementary highpass filter in simplest terms.

So you can see it by the following code:

freqz(alfa, [1 -alfa])

Note that the gain of the filter will be dependent on $\alpha$! Note also that as Matt.L commented you might have $1 + \alpha$ in the numerator instead...


How range of the signal increases when modulation is done?

I think that you are interpreting the Friis equation to mean that RF signals carry for less distance than lower frequency RF signals. The “path loss” term “predicts this”.

The first point which is still confused of why we modulate should be fairly obvious, because we don’t want a 5000 meter antenna. One also uses modulation with amplification.  One also needs to be aware that frequency spectrum is regulated. In the US the FCC is not tolerant of unlicensed transmission so no signal regardless of carrier frequency will carry if you are in jail.

The Friis equation covers lossless spherical spreading just like the homogeneous wave equation.

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

The wave equation predicts that the signal amplitude from an ideal point source drops off as $1/r$ and if that ideal point source has managed to radiate a signal there is no frequency dependence on the electromagnetic field.

The Friis equation doesn’t predict electromagnetic field behavior. It predicts the voltage at the output of an antenna. Path loss refers to the path loss with respect  to an antenna with a fixed aperture where voltage does drop off with the carrier frequency.

For someone like myself who started out as a SONAR engineer, the term “path loss” is an unfortunate choice.

Terrestrial RF propagation almost never conforms to the physical assumptions which the Friis equation is based on. It does have `relevance because it is actually useful as a first order heuristic.  Higher frequencies are more susceptible to scattering and shading than low frequencies.
Low frequencies penetrate buildings deeper. High frequencies tend be favored in line-of-sight applications.

There are other reasons why we want to use a modulated carrier. There isn’t a deep theoretical contradiction in play.

Essentially your conception of Electromagnetic attenuation is a bit too simple. Waveguides for instance have a characteristic known as the cut off frequency. A field has to have a frequency above the cut off frequency to propagate. Photons are electromagnetic waves and they travel billions of miles. The wavelength is not in itself the cause of attenuation, it is the interactions with matter which have the various effects


DSP interview question: use of the identity in development of a significant transform

This is related to Chirp Z-transform (CZT) (refer to the Bluestein's algorithm). Using this identity, the CZT can be expressed in terms of a convolution. Hence, it can be efficiently implemented using FFT.


Discrete time signals connected in series

You simply need to put $w[k]$ into $y[k]$:

$$w[k]=0.25x[k-2] \Rightarrow \begin{align}w[k-1]&=0.25x[k-3]\\
 w[k-2]&=0.25x[k-4]
\end{align}$$

\begin{align}
y[k] &= 0.5w[k-1] + 0.5w[k-2]\\
&=0.5\cdot0.25x[k-3] + 0.5\cdot0.25x[k-4]\\
&= 0.125(x[k-3]+x[k-4])
\end{align}

You can directly read the impulse response from it.


I know this is an old thread, but for reference for future enthusiasts in this field:

The equation for the full PCP algorithm as described by Fujishima is as follows:

$$ PCP(p)=\sum_{ls.t.M(l)=p}^{l}||X(l)||^2$$

$$M(l)=\mathrm{round}\left[12\log_2 \left(\dfrac{f_\text{s} \ l}{N f_\text{ref}} \right) \right]\operatorname{mod} 12 $$

The confusion that OP has is in regards to the 2nd equation:

$N$ is how many samples are used to calculate the FFT, or typically known as the FFT size.  Note that a larger $N$ (more samples, longer time) would mean higher resolution of the FFT.

$f_\text{s}$ is the sampling frequency of the audio file you are processing.  Typically this is  44 100Hz.

$f_\text{ref}$ is the reference frequency for the pitch class you are trying to match.  This is typically a set of 12 frequencies, representing the fundamental frequencies of the 12 semitones of a classical keyboard.  For your reference, they are
[16.35, 17.32, 18.35, 19.45, 20.60, 21.83, 23.12, 24.50, 25.96, 27.50, 29.14, 30.87], starting from C0 to B0.

To understand this equation better, the collective term of $f_\text{s} l/N$ is equivalent to converting the FFT bins into the actual frequency that bin represents.  $l$ here is just the bin counts for the FFT output.  So for every bin of the FFT, the associate frequency can be calculated with this term.

Because every octave is effectively a double in frequency to the previous octave (if C0 is 16.35Hz, then C1 = 2*C0 = 32.70  C2 = 2*C1).  So to reverse this, we take the $\log_2$ of the calculated frequency from the previous step.

Now to the first equation, note the summation if for all $l$, such that $M(l) = p$.  This means for all of the different pitch classes p (0-11), we must only use bins in the FFT such that $M(l) = \textrm{the pitch classes}$.  To find such conditions, we sweep the fref for all frequencies in the FFT so we find the frequencies that contributes the most to the individual pitch classes.

Note that the $M(l)$ equation should not give you a negative number.  The reason for this is the FFT creates both positive and negative frequencies.  But for all real signals (such as audio signals), the negative frequencies is just a reflection of the positive frequencies.  You only really need to deal with the positive frequencies, which resides in the first half of the output from the FFT.  Or you can use a rFFT (such as the one in python - numpy) algorithm that specifically deals with real signals, and will only return to you the positive frequencies.

If you were to implement this algorithm, you would do the following:

Pick a pitch class profile (lets say you want to look at $PCP(0)$ to find all notes that falls under the musical C)
Now, according to the 1st equation, you need to find all $l$ such that $M(l) = 0$
You go to the 2nd equation, going through all the bins of your FFT (which is l).  You also set your $f_\text{ref}$ to be the reference frequency for the PCP bin you're after (at this case, since we are looking for C, we can set the $f_\text{ref}$ to 16.35Hz).
Find all values of $l$, such that $M(l) = 0$.  Go back to the 1st equation, now with the subset list of $l$, take the magnitude of the FFT at those specific bins ($l$), to find the final sum for pitch class 0.


Why would FFT interpolation using zero-padding undershoot the true frequency of a single tone sinusoid?

The FFT of a strictly real input results in a conjugate mirrored result. For a sinewave, both a positive and a negative frequency peak will appear.  If the input is not integer periodic over the full aperture or is zero-padded, the two response peaks will be Sinc shaped, with lots of extended ripples in the frequency domain.  The positive and negative frequency Sinc images overlap and thus can constructively or destructively interfere, depending on the phase of the input sinewave and/or the window.  For the lowest few FFT bins (and those very near N/2), this overlap can be quite substantial, and thus so is the potential peak shift of the summation.


For a trigonometric Fourier series:

$$
\tilde{x}(t) = a_0 + \sum\limits_{n=1}^{\infty} a_n \cos(n \omega_0 t) + \sum\limits_{n=1}^{\infty} b_n \sin(n \omega_0 t)
$$
with
$$
a_0 = \frac{1}{T} \int_{-T/2}^{T/2}  x(t) \ dt
$$

$$
a_n = \frac{2}{T} \int_{-T/2}^{T/2}  x(t) \cos (n \omega_0 t) \ dt, \quad n = 1, 2, \ldots
$$

$$
b_n = \frac{2}{T} \int_{-T/2}^{T/2}  x(t) \sin (n \omega_0 t) \ dt, \quad n = 1, 2, \ldots
$$
and $ \omega_0 = \frac{2\pi}{T} $ .

the Fourier series of an odd square wave will be a sum of only sine terms, and the Fourier series of an even square wave will be a sum of only cosine terms.

This is because the product of an odd and an even function is an odd function, which when integrated over a symmetric interval (i.e. from $-T/2$ to $T/2$), will return zero. For an odd square wave, this means that all the $a_n$ will be zero, and for an even square wave, all the $b_n$ will be zero.

You can save yourself time and potential mistakes by exploiting this fact and only computing the coefficients that are nonzero.


Construct IFFT from only the major frequencies of an input signal?

Without having the original data it is difficult to confirm, but the DC offset is due to bin 0 in the DFT. Bin 0 corresponds to the DC term or average offset value for the sequence. The plot in the DFT result shows a large value for bin 0 while the time domain waveform appears to be closer to 0 average. I would need to see the original data to determine how the DFT bin 0 occurred.

Update: Upon reviewing the original data from the OP, the average of this data is 0.4123 which is consistent with the result achieved.


The input-otput relation for frequency selective channel:

$y_i=\sqrt{\frac{P}{M_T}}H_is_i+n_i$

More information can be found here.


You could apply a 2d FFT. You will get a frequency spectrum as z values on a 2d plane:

 ← That's log(abs(fftshift(fft2(image))))) where image is from the question, but with white borders trimmed

From the z values you can determine the highest peak and calculate its position on the 2d plane. With some geometric calculation, you can derive wavelength, direction and phaseshift from this.


Difference between Gabor filter implementations

Both equations are related by Euler's formula: http://en.wikipedia.org/wiki/Euler%27s_formula
They are mathematically equivalent, it's well explained on the Wikipedia page.

In a practical point of view, $\cos(x)+i\sin(x)$ allows for easy separation of real and imaginary. On the other hand, $\exp(ix)$ makes it easier to deal with Gaussian since $\exp(x)\cdot\exp(y)=\exp(x+y)$.

Hope this helps!


The linspace function as used is going from -fs/2 to +fs/2 for 3000 samples so the +fs/2 is getting counted twice. In contrast the FFT and freqz, which goes from 0 to N-1 do not duplicate the two end points (in this case they go from DC to 1 bin less than fs where fs is the sampling rate). Therefore the sample locations are not exactly the same leading to the difference observable in the two freqz methods.

Further you can get samples of the DTFT by zero padding as another option: fft(x, 3000).

Instead of linspace which will work just fine if you choose the proper start and stop, I like to do:

t = [0: length(x)-1]*1/fs


Zeroing out bins / attenuating them in discrete Fourier domain is universally a bad idea, due to the undesirable time-domain effects of that.

Instead, use a audio processing program to apply an adjustable equalizer to the song of choice. Let's walk you through Free software:

Get audacity;
load your song in that, open "Effects"->"Equalizer"
start with one  of the presets to see how things work, and then
sit down with your grandpa and figure out what frequency ranges to dampen.
Save that preset!
You can later apply it to all your music collection using batch processing.

You should also, in that Equalizer dialogue, "manage / export" your freshly designed filter frequency response. With a little work, you could, for example, build a digital equalizer (in music, your granddad probably won't mind the latency that brings, which would be a bad thing in live music) using e.g. a Raspberry Pi and a USB soundcard, which you could put in between you granddad's CD player and his amplifier, or so. (of course, it could also act as media player/server/center/whatevs.)


Algorithm to Count Zeros Outside Unit Circle for FIR Filter

Here is one answer, if someone can improve on this I will select it as the "right" answer (also comments very welcome on obvious flaws with this approach):

Given Cauchy's argument principle, the number of zeros outside the unit circle is given by the number of encirclements of the origin for the frequency response of the filter as plotted on a complex plane. Each encirclement causes the unwrapped phase to go further than $\pm \pi$ in any interval $\pi$ in normalized radian frequency.

Thus the algorithmic approach would be to unwrap the phase, remove the initial phase at $\omega=0$ and then count the number of crossings of the horizontal lines given by $\pm \pi$, using the Bentley-Ottman algorithm to efficiently count the crossings: https://en.wikipedia.org/wiki/Bentley%E2%80%93Ottmann_algorithm.

Further details on algorithmic phase unwrapping given here:

https://www.ljmu.ac.uk/~/media/files/ljmu/about-us/faculties-and-schools/fet/geri/onedimensionalphaseunwrapping_finalpdf.pdf


solution on the time domain becomes “periodic” after the inverse fourier transform

Multiplying a kernel and signal spectrum in Fourier domain lead to a circular convolution and not a linear convolution, so in order to it become linear convolution you must zero pad your signal and kernel before taking the Fourier transform (up to M+N-1 where M is the signal's length and N is the kernels's length). (if you compare the blue and orange signals, the beginning of orange signal has higher values because of of higher value at the end of blue signal and it has lower values in the end because of lower values in the beginning of the blue signal).


Filter ECG signal with Wavelet and Python

I will try to answer the question after working on this school project.

First, why filter an ECG using wavelets? I had a raw signal, full of noise. The filtering of the signal using the wavelet method makes it possible to capture spatial and temporal information very important for an unusual detection.

Here is a tutorial that inspired me, it perfectly describes the role of coefficients:: https://medium.com/@andrewtan_36013/electrocardiograms-qrs-detection-using-wavelet-analysis-a1070505efee


Differential BPSK vs BPSK what are the advantages and disadvantages?

Differential BPSK

As the information is encoded in the phase difference of two consecutive symbols, phase estimation can be omitted in systems where the carrier phase can be assumed to be constant during two symbols. Specifically, it can be used in a self-coherent receiver which mixes one symbol with its predecessor in order to reconstruct the phase difference between the two symbols. In such a scheme, no local oscillator is required which simplifies the system (e.g. frequency synchronization is not needed).

Non-differential (coherent) BPSK

In differential BSPK, two symbols are "added" and as a consequence the noise of the two symbols adds up. This results in a 3 dB sensitivity gain of coherent BPSK. In other words, coherent BPSK is more robust against noise.


FFT/PSD/IFFT analysis on single axis piezoelectric accelerometer signals for curb impacts

Filtering by zeroing out bins is not a recommended approach as it will introduce significantly more time domain ringing, as detailed here
Why is it a bad idea to filter by zeroing out FFT bins?
Consider using the firls function available in MATLAB / Octave or Python scipy.signal to design optimized multiband filters around your frequencies of interest.
However the frequency content will be driven by repetition in the data; if the actual application will repeat (or not) at unknown and variable intervals, any such filtering techniques will not be useful.


How to filter out everything but a single frequency in the time domain?

Depending on the exact characteristics of the signal and implementation requirements (if any) I can think of a number of approaches to extract content at a "single" frequency. The most common:

Apply an FFT to calculate the spectrum of the signal and extract the frequency of interest.
Apply Goertzel's algorithm to effectively extract signal content at a single spectral line.
Apply a bandpass filter around the frequency of interest.
Apply an inverted notch filter around the frequency of interest.

Depending on the amplitude and frequency content of the signal, the additional noise present and other implementation requirements (e.g. should the signal be processed in real time or offline, should the implementation be digital or analog, etc.) a specific approach might be preferred.


Fourier Transform of the Hilbert Transform of cos(t) (using Fourier time-shifting property)

I forgot about the property of a dirac delta function: $f(t)\delta (t)=f(t)|_{t=0} \delta (t)=f(0)\delta (t)$

$Y(\omega)=F\{x(t-\frac{\pi}{2})\}=X(\omega)e^{-j\omega\frac{\pi}{2}}=(\pi \delta(\omega-1)+\pi \delta(\omega+1))e^{-j\omega\frac{\pi}{2}}=\pi \delta(\omega-1)e^{-j(1)\frac{\pi}{2}}+\pi \delta(\omega+1)e^{-j(-1)\frac{\pi}{2}}=-j\pi \delta(\omega-1)+j\pi \delta(\omega+1)$


How the Grayscale images are converted back to RGB?

It is not possible to get color back. A grayscale image can be converted to RGB image, by copying the grayscale data to all three RGB channels,but it will still be grayscale.

Note that in general, RGB is not simply averaged to grayscale, but a weighted average is used to match eye sensitivity.


Generate Complex White Gaussian Noise in MATLAB

If you want a Circular Complex Gaussian Noise (Independent):

vComplexNoise = sqrt(noiseVar / 2) * (randn(1, numSamples) + (1i * randn(1, numSamples)))

For correlated noise you'll need to define the Co Variance Matrix and use Cholesky Decomposition.

Update

Following @Stanley Pawlukiewicz advise, run the following code:

numSamples = 100000;
noiseVar   = 4;

mA = sqrt(noiseVar / 2) * (randn(numSamples, 1) + (1i * randn(numSamples, 1)));

var(mA)

You should see result which is very close to noiseVar on the screen.


Reason of Fast Fourier and Inverse Fast Fourier

FFT and IFFT are algorithms that implement the (inverse) discrete Fourier transform (DFT). These transforms convert a signal into another representation, namely the frequency domain, and back. This other representation allows us to analyze certain properties of the signal and also to change these properites which would be hard to do in the original representation.

Consider a digital audio signal: it's a sequence of samples that represent amplitude over time. This representation is useful to study the temporal properties of the signal and to play it back to make it audible. We can hear that the audio signal contains low and high sounds (frequencies) but from looking at the signal it's impossible to tell what frequencies exactly. After applying the DFT to the signal we can analyse what frequencies it contains and how strong they are. We can even amplifiy or attenuate certain frequencies and then apply the inverse DFT (IDFT) to make the signal audible again. In this way we have built a filter. The application of IDFT is necessary in this case, because the frequency domain representation of the audio signal is useless for playback.

Image processing deals with two-dimensional signals and therefore two-dimensional (I)DFT. The concept of frequencies with images is not so illustrative as with audio signals. Edges, i. e. sudden changes of intensity are represented by high frequencies whereas smooth gradients are represented by low frequencies. This allows us to do edge detection, for example. If high frequencies are attenuated this will result in a smoothed image after application of IDFT. Again, for just looking at the image, the frequency domain representation is useless but it allows analysis and filtering.


Change pitch of a wav file without changing the Sample rate?

If you're looking for a library solution, there's a few, for instance librubberband, binaries available in many Linux distros.

If you need a wheel to reinvent, well, approaches may vary, from granular stretching (more applicable to time stretch, but these two effects always come hand in hand), to FFT uppitch (probably somewhat similar in the frequency domain), I doubt there's a standard consensus on how to implement this effect.


as the comments to your question mention you have not done the first step: system classification

Is the system linear? We have to input three different signal in the system and compare the outputs: the two signals $x_1$ and $x_2$ and the sum signal $x_{1+2} = x_1 + x_2$

\begin{align}
 y_1 &= (-1)^t x_1 \\
 y_2 &= (-1)^t x_2 \\
 y_{1+2} &= (-1)^t (x_1 + x_2) \\
\end{align}

and even at this point it is easy to see that $y_1 + y_2 = y_{1+2}$ that means the system is linear

Is the system time variant? it's usually very easy to say yes here when your system is a function of the time (or n in your case) $z(t) = (-1)^t $ Mathematically you have to input a single signal $x(t)$ at two points in time and check if the system gives you different outputs at different times:

\begin{align}
 y_1 &= (-1)^t x(t) \\
 y_2 &= (-1)^t x(t-\Delta t) \\
\end{align}

to compare 1 and 2 you have to time shift the output of the first system to make sure they overlap the way you want.

\begin{align}
 y_1 &= (-1)^t x(t) \\
 y_{1-\Delta t} &= (-1)^{t - \Delta t} x(t - \Delta t) \\
 y_2 &= (-1)^t x(t-\Delta t) \\
\end{align}

the comparison shows that those two are not the same, which means your system changes over time

\begin{align}
(-1)^{t - \Delta t} x(t - \Delta t) &= (-1)^t x(t-\Delta t) \\
(-1)^{t - \Delta t} &= (-1)^t  \\
(-1)^{ - \Delta t} &= 1  \\
\end{align}

this is not true in general, that means your system is linear and time variant and you can not use the description of a LTI System. But there's a definition of the impulse response of a time variant system that you could use:

Time-varying "impulse response"

\begin{align}
y(t) = \int h(\tau, t) x(t - \tau) d\tau
\end{align}

since your system is memory less and kausal (i don't think there's a trivial mathematical proof for those two properties but it is easy to see), your time variant impulse response is rather simple:

\begin{align}
h(\tau, t) = \delta(\tau) (-1)^t
\end{align}

which leads to the collapse of the integral in the time variant convolution (just check the rules for integration over the dirac function if you don't remember them):

\begin{align}
y(t) = \int \delta(\tau) (-1)^t x(t - \tau) d\tau \\
y(t) = (-1)^t x(t)
\end{align}

And the last line ( and basically your first line) is the best description for linear memory less time variant systems - a simple multiplication with a time variant function or in general just:

\begin{align}
y(t) = x(t) z(t) \\
\end{align}

bye bye stackexchange you managed to annoy me enough in my first two answers for me to decide not to come back.


$$H(s) = \frac{R C s}{R C s + 1}$$

If $R C s$ is much lower than 1 (i.e., $R C s \ll 1$, in Math), then you can make an approximation: $$H(s) \simeq \frac{R C s}{1}$$

Basically, the effect of the $R C s$ in the denominator becomes insignificant.

So, just pick a radian frequency well below $1 / R C$ and compute the magnitude of $H(s)$ at that point -- mark that on your plot, and draw your line rising at 20dB/decade from there.


The Fourier transform can be generalized for functions that are not absolutely integrable. We can define a Fourier transform for functions with a constant envelope (e.g., sine, cosine, complex exponential), and even for functions with polynomial growth (but not with exponential growth). In these cases we must be prepared to deal with generalized functions in the expression of the Fourier transform, such as the Dirac delta impulse or its derivatives. This is also true for the Fourier transform of the step function.

The multiplication property of the Fourier transform remains true, albeit with certain restrictions, so we can generally compute the Fourier transform of the convolution of two functions by multiplying their Fourier transforms, provided that the convolution exists.


Scipy butter filter - remove DC offset w/ high sample rate

Why do I have ringing and such an aggressive response?

Because you designed a really aggressive filter, which is borderline unstable. All the poles are less than 1/1000 away from the unit circle.

How does the low cut off frequency and high sample rate affect the filter response?

The lower the cutoff, the higher the order, and the higher the sample rate, the more it will ring (in terms of samples).

Why did Scipy go from the b,a approach in the butter() filter to the "sos"?

Because in transfer function representation the filter would actually BE unstable, i.e. the numerical precision of 64-bit floating point would NOT be enough to run this filter as a high order polynomial.

I"m curious why they switch to this new sos appraoch. My understanding is that it is more accurate and/or prevents errors in the b,a approach.

NO, NO, NO. Transfer function representation is the worst and should be avoided in favor of poles, zeros or second order section. Once you are in transfer function, you often can't get back to poles and zeros. It requires finding the roots of a high order polynomial which is a numerically ill defined problem.

Also, you can notice (if you print(w) to see my freq axis) that my frequencies seem to run from 0 to 2*pi. Is that normal or should it be from 0 to pi?

If you use a function, it's always a good idea to read the documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.sosfreqz.html#scipy.signal.sosfreqz
Since you pass in the sample rate as an argument, the function already returns the frequency axis. By your attempt to normalize it again, you get complete non-sense. The dead give away here is the scale factor $10^9$ on your frequency axis. Your current axis tops out at 1.25 GHz.

If question #4 is is "yes this should be 0 to 2*pi" than should line #23 in script:

It should simply be ax1.plot(w, h_db)
The question you didn't ask:

What I should be doing ?

What @Ben said: use a first order filter AND either subtract the mean from the signal or seed the filter state with the signal mean.


Effect of sampling rate on BER simulation

Another hypothesis: as you increase the sample rate you also increase the bandwidth. Depending on how you generate the noise, this probably means that the spectral density of your noise decreases. If the spectral content of your signal stays constant, that means that the in-band SNR goes up.


What's the difference between a DSP and microcontroller?

a Digital Signal Processor is one that has, in its instruction set, some instructions and addressing modes that are optimized for processing digital signals.

usually these optimizations can be shown around what is needed to perform the dot-product needed for an FIR filter.

$$ y[n] = \sum\limits_{i=0}^{L-1} h[i]\,x[n-i] $$

to do this in, say, $L$ instructions, a DSP must be able to do in one instruction:

multiply $h[i]$ and $x[n-i]$ together.
accumulate that product into an existing sum.
fetch the next $h[i+1]$ and $x[n-i-1]$ in anticipation of the next multiply-accumulate.  since these are two numbers to fetch, a DSP will use something called a Harvard architecture that has at least two separate memory spaces for $h[i]$ and $x[n]$ so the DSP can fetch these two numbers simultaneously.
addressing $x[n]$ must be in a circular queue.  a DSP will perform the modulo (or "wrap around") arithmetic on the index or address of $x[n]$ necessary without additional instructions.
the result $y[n]$ will eventually go to an output DAC or fixed-point stream and there is some way to saturate the value of $y[n]$ against some $\pm$ maximum without additional instructions.  if the DSP is a fixed-point DSP, then this accumulator register will have width in bits that is the sum of the bitwidth for $h[i]$ and the bitwidth for $x[n]$ plus a few more bits on the left as "guard bits".

a general purpose CPU can do all of these, but not likely in a single instruction cycle and things like modulo arithmetic and saturation will need their own specific instructions in a general-purpose CPU.

a DSP may also have some instructions and an addressing mode that facilitates the operations of the Fast Fourier Transform (FFT).  this may include instructions necessary to perform an FFT "butterfly" in as few as four instruction cycles (or two instruction cycles if SIMD is operational).


If I understand correctly, you are mixing a 1500 Hz sinusoid with an MSK signal, which uses 1200 Hz to transmit 0 (or 1) and 1800 Hz to transmit 1 (or 0). The spectrum of the MSK signal, though, is quite large and includes energy at 1500 Hz. So, when mixing, you can expect to find a DC value.


Constant used to determine the minimum samples in Compressed Sensing

The constant C is usually left unspecified as it can be very hard to calculate.  It is independent of the signal dimension, and results in Compressed Sensing are generally only concerned with big-O results, so it falls out when we say something like $m = \mathcal{O}(s \ln{(N/s)})$.  The important part in this statement is the strong dependence on the sparsity level, and that the signal dimension $N$ only logarithmically affects the number of measurements.

To give an example of where such a constant would come from:  Let $V_d$ be the volume of the unit $\ell_2$ ball in $\mathbb{R}^d$.  Some simple cases are that $V_2 = \pi r^2$ and $V_3 = \frac{4}{3} \pi r^3$.  Note that both of these results are of the form $V_d = C r^d$, and this pattern generalizes for finite $d$.  In CS, it is $r^d$ that we are interested in, as it captures the dependence on the signal dimension.

So $C$ will vary depending on the algorithm as different inequalities may be used in the proofs.  I don't know of any publications that specifically calculate C.


How does sampling rate of $x[n]$ relate to sampling rate in frequency domain after DFT transformation?

A better way to think about it is the DFT is agnostic about the sampling rate.  What matters is the number of sample points.  If you look at the DFT definition, $F_s$ does not appear in it, but $N$ does twice.

The sampling rate can be used to figure out the duration of the sampling frame given the number of samples, or used to figure out how many sample points for a given duration.  Once the DFT is taken on the sample, the bin number corresponds to the number of cycles per frame.  The number of cycles per frame can be divided by the number of seconds per frame (the frame duration), or whatever unit of time is being used, to get the number of cycles per second (or time unit).

I think the equation the r b-j gave is better expressed as:

$$ f = k \cdot \frac{F_s}{N} $$

Where the units are:

$$ \frac{cycles}{timeunit} = \frac{cycles}{frame} \cdot \frac{\frac{samples}{timeunit}}{\frac{samples}{frame}} $$

The latter factor can be considered a conversion factor for the frequency.

When $ k = N/2 $:

$$ f = \frac{N}{2} \cdot \frac{F_s}{N} = \frac{F_s}{2} $$

Which confirms r b-j's statement about the halfway bin corresponding to the Nyquist frequency.

Hope this helps.

Ced


No, the signal
$$
h[n] = \frac{\sin(\omega_c \, n)}{\pi \, n}, \qquad -\infty < n < \infty
$$

is not absolutely summable, as the sum

$$
\sum_{n=-\infty}^{\infty} \Big|h[n]\Big| = \sum_{n=-\infty}^{\infty} \left |\frac{\sin(\omega_c \, n)}{\pi \, n} \right|
$$

diverges to infinity. Hence the absolute sum does not converge and therefore it's considered as an unstable system's impulse response. Yet, its square is absolutely summable:

$$
\sum_{n=-\infty}^{\infty} \Big|h[n]\Big|^2 = \sum_{n=-\infty}^{\infty} \left |\frac{\sin(\omega_c \, n)}{\pi \, n} \right|^2
$$
does not diverge. This kind of condition is stated as mean square convergence; i.e., a typical condition for discontinuous frequency responses that has finite jumps in its frequency response which should be continuous for all $\omega$ according to the Fourier theorem that would require uniform convergence at all frequency points $\omega$.


It seems to me that the energy carried by a sound wave should be
  proportional to frequency squared.

Why? What makes you think that this is the case ?

  But in DSP they are all considered to have the same signal energy
  (integral of x(t) squared)

The concept of associating power or energy with a a digital signal is somewhat tricky. Digital signals are just a bunch of numbers so assigning actual physical properties to them doesn't really work unless you carefully define the context.

Mostly the concept is used assuming that the digital signal is an accurate representation of an analog signal that has actual physical properties. If this relationship is well defined and constant you can compare the sum-of-squares of two different digital signals and draw some conclusions on the energy difference of the actual analog situations. However, that really requires "all things being equal" and there are actually a lot of things that need to be equal for this to hold.

STEP 1: basic physics of power

Power, intensity or energy are typically defined or calculated as the product of two field quantities (not one). In electricity power is defined as the product of of voltage and current. Simple example: if your signals were defined as "voltage over a one Ohm resistor", than you can indeed calculate the power by summing the squares and it would be the same power for all of your three examples since the current is proportional to the voltage.

If it were the voltage over an ideal inductor or capacitor then the answer would be wrong. Average power over one of these would be zero in all cases. For an ideal inductor voltage and current are 90 degrees out of phase and the average power is therefore zero.

STEP 2: physics of sound

The two field quantities that make up a sound wave are sound pressure and particle velocity. First you need to define which one your digital signal actually represents. Most of the time it's pressure but many second order microphones (cardioid, dipoles) will create a signal that's proportional to the particle velocity. If you are far enough away from the sound source than pressure and velocity are in-phase and proportional to each other and the proportionality factor is the free field impedance of air (density times speed of sound). In this case you can actually calculate the intensity at the microphone position by simply squaring and scaling. If you are close to a sound source that doesn't work anymore and you actually need to measure  both signals.

As others have said: per it's physical definition sound can't exist at 0 Hz. Sound pressure is a variation of pressure around a steady-state average and if there is no variation there is no sound.

Step 3: human perception

Sound energy and perceived loudness are only loosely related. Human perception is very dependent on frequency and its quite non-linear. Building good loudness models is quite difficult. Simple example: 1 Watt of sound at 1kHz radiated into a typical residential room would extremely loud and quite painful. At 50 kHz it would be utterly inaudible, even though it's the exact same amount of physical sound power.


Fourier transform of a rect*half triangle

There is a Fourier transform property called as "differenciation in frequency domain" which is as follows:

If the Fourier transform of $x(t)$ is $X(\omega)$, then the Fourier transform of $tx(t)$ is as below: $$\mathcal{F}(tx(t))  = j\frac{d}{d\omega} (X(\omega))$$

Where $\omega = 2\pi f$,

I think you would be able to easily derive the answer you mentioned in your question now.

Hope that helps


The Fourier transform property in your first equation is for continuous time but the image is a discrete time signal. $\partial_x$ isn't exactly differentiation, it's a discrete time filter which does finite differencing.

Since you are only taking derivatives in one direction, let's think about a 1-D signal for simplicity (or, equivalently, let's assume your image has only one row). Let the signal be $f(x)$ and its continuous time Fourier transform $F(u)$. Then we have the transform pair $\frac{df(x)}{dx} \leftrightarrow j 2\pi u F(u)$. In discrete time, if $x[n] \leftrightarrow X[k]$, you can't really say $\{x[n]-x[n-1]\} \leftrightarrow j2\pi k X[k]$.

Instead, we can think of $\partial_x$ as a high pass filter. It then agrees with the intuitive meaning of the continuous time Fourier transform property too (where multiplying by $u$ attenuates lower frequencies more than higher frequencies). Try this code, to see the high pass response:

d = np.zeros(10000)
d[0] = -1.0
d[1] = 1.0
D = np.fft.fft(d)
plt.plot( np.abs(D) )


Designing a multi-band FIR filter on MATLAB

According to your filter equations, your polynomials for $b1$, $b2$ and $a$ should only contain 1, zeros and -1. Your don't need to evaluate the last coefficient (num_exp, and others) at the complex exponential - That is done inside the freqz()
Take a look at this example here
So you should have:
a=[1, -1]
b1 = [1, 0, 0, ..., 0, -1]
b2 = [1, 0, 0, ..., 0, -1]

So $a$, $b1$, $b2$ contain the coefficients of the $z^{-n}$ terms.


Fourier Series Representation of Continuous-Time Periodic Signals

When you are given a continuous-time periodic signal $x(t)$ and you want to find out the corresponding CTFS (continuous-time Fourier series) coefficients $a_k$ associated with $x(t)$, then you use the analysis equation; i.e, analyse $x(t)$ to find out $a_k$ use

$$ \boxed{ a_k  = \frac{1}{T} \int_{0}^{T} x(t) e^{-j \frac{2 \pi}{T} k t} dt  }\tag{1}$$

On the other hand, when you are given a set of CTFS coefficients $a_k$ and you want to obtain the corresponding continuous-time periodic signal $x(t)$, then you use the synthesis equation; i.e., sum up exponentials weighted by $a_k$ and create $x(t)$ as in

$$ \boxed{ x(t) = \sum_{k=-\infty}^{\infty} a_k e^{j \frac{2\pi}{T} k t } } \tag{2}$$


How to derive the results that averaging $N$ signals yields a $\sqrt{N}$-fold increase in signal-to-noise ratio?

I will show how to calculate the SNR for the case of $N=2$ measurements; it is easy to extend the result to general $N$. Assume a signal $s(t)$ has power $S$, and the noise $n(t)$ has variance $\sigma^2$ and zero mean. Then, the signal $s(t)+n(t)$ has SNR equal to $S/\sigma^2$.

Now assume you observe $s(t)$ twice, each time with different, uncorrelated noise, $n_1(t)$ and $n_2(t)$, each with zero mean and variance $\sigma^2$. You average the two observations to get $$\frac{2s(t)+n_1(t)+n_2(t)}2 = s(t)+\frac{n_1(t)}2+\frac{n_2(t)}2.$$

The variance of both $n_1(t)/2$ and $n_2(t)/2$ is $\sigma^2/4$, so the total noise variance is $\sigma^2/2$ and the SNR is $2S/\sigma^2$, for an improvement of 2.

For some reason I don't understand, the wikipedia page on signal averaging defines the SNR as $S/\sigma$. Under that definition, the improvement will indeed be $\sqrt{N}$.


Can someone explain me how to generatee a doppler spectrum?

When we talk about Doppler, we typically refer to time-varying systems.

A popular model in system theory is to assume systems to be LTI (linear time-invariant). This means that beside being linear, their input-output transformation does not change in time, it does not matter when we apply a certain input signal $x(t)$, we'll always get the same $y(t)$. Such systems are completely characterized by an impulse response $h(\tau)$ which specifies exactly this transformation: $y(t) = h(t) * x(t)$, where $*$ denotes convolution.

In practice, this assumption is often not true. Wireless communications are a good example. While the channel is typically (approximately) linear, its characteristics change in time, e.g., due to movement of transmitter, receiver or interacting objects such as cars, trees, people, etc. For this reason, we need to extend the impulse response by another time variable that tracks the change of $h(\tau)$ in time. We often simply write $h(\tau,t)$. It may seem weird to have two time variables but there is solid theory behind this. As an engineering perspective, we can think of $\tau$ as the "fast time" (at which signals travel, orders of $\mu$s) and $t$ as the "slow time" (at which the channel changes, order of ms) and treat them more or less as independent.

The rate at which the channel changes is called Doppler frequency, say $\alpha$. Mathematically we can think of it as the Fourier transform of the time variable $t$. Given that we have two time variables and we can Fourier both of them, this leads to four possible representations of the input-output relation: $h(\tau,t)$ (time-varying impulse response), $H(f,t)$ (time-varying transfer function), $h(\tau,\alpha)$ (Delay-Doppler spread function), and $H(f,\alpha)$ (name it as you like). They are also referred to as the (first set of) Bello functions to drop a keyword for further reading.

Since this is the title of your question: the Doppler spectrum measures how much power is associated to a certain Doppler frequency, it is the power spectral density of $h(\tau,t)$ over the time variable $t$.

Regarding the implementation of a certain Doppler spectrum into your simulations: in general what you would do is to generate the channels for different time points independently and then add a filter over time that introduces the proper amount of correlation. If you're using Matlab, have a look at [MATLAB Fading channel page], it explains it quite well.


How do I calculate the SNR of a noisy signal without having the noise?

If the noise is additive, and your signals are in sync (i.e. same number of samples, and temporal aligned). You need to subtract your original signal from your Noisy signal - which will give you the additive noise. You can then calculate the power of this and compare it to the clean version of your signal to compute the SNR

In Matlab this should just be:

Noise = NoisySignal - OriginalSignal; %//Should give a vector or Matrix depending on what stereo or mono.
SNR = snr(OriginalSignal,Noise) %// should give the single SNR reading in dB.

However, in a case where the signals are not in sync, it becomes much harder, and in general, we consider it a mix of different elements that could have caused the distortion (or lose in sync), and usually a combination of SAR (Signal to Artefact Ratio), SIR (Signal to Interference Ratio) and SDR (Signal to Distortion Ratio) are used instead. The SIR is very similar to SNR in this case.


Gender classification using Periocular Region

You can't if you have nothing to start labelling them with.

As you describe it, there's simply no info available (not even about statistics) about the gender of your reference images, so there's nothing you can train.

What you can do is use a different technique (as you say "biometrics say..." I assume there's a non-ML method of classification) and run it over the dataset, and use the resulting labels to train your neural network. Discussion of the sensibility of that is left to you as an exercise for your report!


Regarding your integral, it is not quite correct I think. You used partial integration to solve it, which works. Let's consider $\int t {\rm e}^{-\jmath \omega t}{\rm d}t$ and call $f(t) = t$ and $g'(t) = {\rm e}^{-\jmath \omega t}{\rm d}t$ so that $g(t) = -\frac{1}{\jmath \omega}{\rm e}^{-\jmath \omega t}$. Then using $$\int f(t) g'(t){\rm d}t = f(t)g(t)-\int f'(t) g(t) {\rm d}t,$$ we obtain $$\int_0^1 t {\rm e}^{-\jmath \omega t}{\rm d}t = \left[-\frac{t}{\jmath \omega}{\rm e}^{-\jmath \omega t}\right]_0^1 + \frac{1}{\jmath \omega}\int_0^1{\rm e}^{-\jmath \omega t} {\rm d}t.$$ It's close to what you did, but you used $g'(t)$ in the first term instead of $g(t)$.

  ended up finding but still there's e(−ω∗j) stuck which confuses me as i believe the plot should be 2 dimensional only (toward x and ω, not toward j).

Note that $j$ is not a variable! It's the imaginary unit $\sqrt{-1}$. Hence, your expression only depends on one variable, which is $\omega$. Matlab will recognize $j$ as the imaginary unit by default, so you should be all set for plotting it.


Stitching together audio buffers with potential discontinuity at the boundary

In 1998 I made a "zero level adjustment click removal" system that shifted the zero level to fix the discontinuity, and this "zero level" was gradually moved back to 0, by decreasing it by 1/2048 of its value every sample. It worked well for low-frequency sounds, but when a high-frequency sound ended abruptly at a large sample value, that created a low-frequency click. I ended up removing the system and simply adding more channels in order to fade out sounds. For this I used a one-pole lowpass filter on all volume envelopes, and when a sound was killed the volume envelope was set to zero and some time was allowed for the fadeout before making the channel available to other sounds. If too many channels were reserved for fading out, the one with the smallest amplitude value was freed to use.

But the thought kept bugging me that if the zero level adjustment would be done only to a lowpass filtered portion of the output, then it might sound better. Purely high-frequency sounds don't really need click removal.


Calculating power and energy of a signal in matlab

The Matlab function int needs a symbolic expression as an integrand, and a symbolic variable as an integration variable. In your example s and t are both vectors, and you can't use int in this case. If you want to compute the power or energy of discrete time signals, then you need to use the corresponding definitions:

$$E_x=\sum_{n=-\infty}^{\infty}|x_n|^2\\
P_x=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{n=-N}^N|x_n|^2$$

In practice you will need to consider finite length signals unless you have an analytic formula for $x_n$ and some software package that can deal with the analytic evaluation of infinite sums. Obviously, finite length signals always have finite energy and $P_x=0$ according to above definition.


Recovering a signal that was upsampled by non-integer factor

The Sampling operation (both upsample and downsample) depends on two very critical conditions:

1- The existance and applicability of ideal frequency selective filters

2- The operated signal being strictly bandlimited

For most practical systems these two critical conditions are only approximately met. Hence the computational results by those practical systems are only approximate and not exact.

Note that in case of simple but crude interpolation (upsampling) by zero-order hold, which requires no filter at all as being just a sample replication, you can exactly recover the original samples due to the simplicity of the operation.

Unless some advanced algorithms are utilized, therefore, some deviation (but controllable) is expected at the output of any practical sample rate conversion systems.


Finding causal impulse response given the imaginary part of the frequency response

You're right about the relation between the imaginary part of the frequency response and the odd part of the impulse response:
$$\textrm{DTFT}\big\{h_o[n]\big\}=jH_I(e^{j\omega})\tag{1}$$
So from the given $H_I(e^{j\omega})$ you can obtain $h_o[n]$. Now note that the odd part of $h[n]$ is defined by
$$h_o[n]=\frac12\big(h[n]-h[-n]\big)\tag{2}$$
The important extra information that you have is that the system is causal. This means that $h[n]=0$ for $n<0$. Consequently, from $(2)$ we obtain
$$h[n]=\begin{cases}2h_o[n],&n>0\\0,& n<0\end{cases}\tag{3}$$
The only ambiguity we have is the value of $h[0]$, which cannot be obtained from the imaginary part of the frequency response.
Note that the same solution can be obtained by using the fact that for causal signals, the real and imaginary parts of the frequency response are related via the Hilbert transform. But the above method is probably even more intuitive and doesn't require any knowledge of the Hilbert transform.


Gradient Domain Reconstruction - Scaling Problem

The Poisson equation may suffer from ill-posedness if the boundary condition is not suﬃcient to yield a unique solution to the problem. The most often used boundary condition is Dirichlet boundary condition, in which the values of $I(x,y)$ are specified for those points $(x,y)$ on the boundary of $I$.

In your case, after you stretch your solution to another scale, you got the correct solution. It is possible the solution to your Poisson equation is not unique, and your algorithm just found out the correct one based on some restrained part (such as minimal norm of least square) in the optimization. It is much much more possible that it is already sufficient to determine $I(x,y)$ up to an arbitrary additive constant, and this constant is determined by the solution of Laplacian equation: $$\nabla^{2}\hat I= 0$$


The notes made me wonder as well, and I cannot suspect the author made a mistake. The first graph with periodicity is generic. The second one corresponds to dyadic iteration (wavelets), often displayed in $[0,\pi]$, except when dealing with complex wavelets where  $[-\pi,\pi]$ is used. I agree, $[0,2\pi]$ could have been used instead, but I have never seen the $2^k$ partitioning mirrored at the end.

You can see from the rest of the paper that the author uses $[0,\pi]$.


Reconstruction of bandpass filtered signal from decimated version of itself

Zero stuffing (which is what you propose) will cause a lot of harmonics in the upsampled signal.  They will be proportional to the strength of the original signal.  Your brick wall filter (a low pass will do, no need for a bandpass) will greatly reduce them, but won't get rid of them entirely.

Depending on your application, the result of zero stuffing and filtering may be fine.  If not, you could minimize the harmonics with a couple of methods.

You could duplicate samples (original sample+ (n-1) copies.)

You could interpolate (calculate values along the line between sample und sample+1.)

Duplicating is easy and will cost about as much processor time as zero stuffing.
Interpolating would be smoother but would cost more time as you must calculate individual values.


2D deconvolution of recorded electron beam data

Your method of deconvolving seems to work.  Below is some Gnu Octave (more or less Matlab) code that implements your deconvolution.

Feed it two files that represent your target and your beam and it will perform a convolution of the two, then use your method to reverse the convolution.
The result shows the original two images, the convolved image, and the recovered shape of the beam.
The attached picture shows the results from my test.

Top left is the target (50x50px)

Lower left is the beam (14x16px)

Top right is the convolution of the two (the map of your electron counts)

Lower right is the recovered beam shape

The padding is a little off, so there's some wrap-around in the recovered image.  The result is clear enough, though.

I would check the images that you used for beam and target, and also check the convolved image to see if it is what you expected.  Part of the problem might also be the "same" that you used in the convolution.  That returns an image the same size as one of the convolved images, which will cause you to lose a lot of information.  If I use "same" then I don't get a recovered shape at all.  The convolved image is larger than the larger of the two images, so you have to use "full" then pad the target image to do the deconvolution.

ETA:
I got so busy looking at the details of the deconvolution that I forgot to mention your thoughts on the effect of scanning the beam across the target.

I agree that this is equivalent to convolving the beam shape with the target shape, so that I think you are on the right track to determining the beam shape - if I weren't convinced of that I wouldn't have bothered to check the deconvolution.

There's a few things that you will need to keep an eye on when you go to actually do the scan to determine the beam shape:

Scan area - you will need to overscan the target.  Take the radius of the target as R and the maximum radius (estimated) of the beam as r, then you will need to scan a square of at least R*2+2*(r*2) - actually more, but that is the absolute minimum.  Since you don't know the beam shape (and therefore its true size,) it may be better to go overboard here rather than have to repeat the measurements
The image you use to do the deconvolution must match (as perfectly as possible) the true shape and size of your target.  The position of the target shape within the image isn't as important.
When you pad your target shape image, center the image in the padding.  As I've done it, the padding is all on the top and left side.  This puts the right and bottom edges of the target shape right at the edges of the image.  Rounding errors can (and do, as you can see) cause part of the image to wrap around.  Since you will be creating the images, you could just create your image with the appropriate border instead of padding it in code.

target=imread("target.png");
beam=imread("beam.png");
convolved=real(fftconv2(target,beam,"full"));
widthc=size(convolved,2);
heightc=size(convolved,1);
widtht=size(target,2);
heightt=size(target,1);
padw=int32(widthc-widtht);
padh=int32(heightc-heightt);
paddedtarget=padarray(target,[padh,padw],0,'pre');
recovered=real(ifft2(fft2(convolved)./fft2(paddedtarget)));
subplot(2,2,1)
imshow(single(target));
subplot(2,2,2)
imshow(single(convolved));
subplot(2,2,3)
imshow(single(beam));
subplot(2,2,4)
imshow(single(recovered));


You can't really use Q5.27 in 32-bit storage. You need at least one bit for the sign, so I guess you'd call it Q4.27 since the sign bit is implicit.  The conversion from float is exactly the same as for Q0.27 except you can convert any numbers between -16 to just less than 16.

int(0.9 * 2^27)  = 0x07333333
int(-0.9 * 2^27) = 0xf8cccccd

As you can see you need the sign bit to represent a negative number.

-16 and 16 arise from:

0x7fffffff / (2^27) = 15.9999999925  (biggest +ve number)
0x80000000 / (2^27) = -16            (biggest -ve number)


Amplitude of frequency bin in FFT doesn't match time-domain amplitude

You seem to know how to handle the DFT/FFT scaling so that its output matches the
sampled signal's power/amplitude levels.

Yet there is a discrepancy between your expectation and the obsevred DFT result? There are a few reasons. Your assumptions on the sampling rate or the signal frequency could be wrong. Indeed as it seems from your figures, the signal is not a pure sine wave. In that case the amplitude looses its crisp definition and spectral leakage and main lobe width smearing becomes an issue. For details please search the site for spectral estimation, spectral display etc...

For your convenience, I have produced an ideal case to your description with MATLAB / OCTAVE code. I'm getting exactly what is expected. So this suggests me that your signal is not a pure sine wave. It's either modulated, or your sampling rate and frequency information is not consistent.

clc; clear all; close all;

T = 0.0225;     % observaton interval in seconds
f0 = 750;      % sine frequency in Hz.
Fs = 357E3;     % sampling frequency in Hz.

t = 0:1/Fs:T;   % analog time

N = length(t);  % number of samples taken...

A = 1.5;        % amplitiude of the sine
x = A*sin(2*pi*f0*t);   % sampled analog signal...

figure,plot(t,x);
title('Analog current waveform');

M = N;          % DFT length (mostly equal to sample size)
X = (2/N)*fft(x,M);

f = linspace(-Fs/2,Fs/2,M); % analog frequenc range in Hz.
figure,stem( f , abs(fftshift(X)));
axis([-1500 1500 -0.1 2]);
title('FFT magnitude of the obtained samples');

The resulting figures are as follows: