I'll assume that your data is in the form of discrete samples equispaced in time (as given by a standard analog-to-digital converter).

If the nominal signal is a clean 60Hz sine wave, I would suggest filtering out the 60Hz with a digital notch filter. This will essentially subtract the correct signal, leaving you with the noise. Then, compare the magnitude of the filter output to a threshold (which you will choose empirically) in order to detect the spikes. Note that the filter may or may not introduce a delay between the original signal and the filter output (esp, if you are doing this in realtime), which you will need to compensate in your time-of-occurrence estimate. This time compensation may also apply to the methods below.

If the nominal signal is not a clean sine wave, you can probably track the nominal frequency characteristics of the signal, e.g. using short time fourier transform, and watch for undesirable deviations from this.

You may also be able to perform a linear prediction on the signal, and detect spike noise events based on the magnitude of the prediction error signal.

These are all just suggestions, many other approaches likely exist.


Adaptive Band pass filter for extracting respiratory waveform from accelerometer data

Under the extremely controlled conditions required to measure respiratory rate using this method, you can try to use the autocorrelation (which is more robust in the presence of noise) to determine the "dominant" frequency in your signal.

  2) The red line is spectra after filtering. Why has the amplitude of certain frequencies increased? I am using the 'spec' function in R: uses the hanning window for computing.

If that is not a mere rescaling of the waveform as it gets plotted over the same figure (?), then it would be safe to assume that the filter included some overall gain.

Hope this helps.


How to calculate the mean/center frequency of the spectrum?

[EDIT] If just want to compute the location of the maximum amplitude, you just have to compute the maximum. You can work in a more generic context, however. As you remarked, the notion of center frequency is not evident. The notion of centroid could be more pertinent.

I assume that your signal is real (since you are talking about a maximum, not well defined in the complex domain). Since its amplitude spectrum is symmetric about $0$, this possible center of mass is useless, since it is the same for all real signals (the zero-frequency). More likely, one is looking at a more useful quantity, based on a one-sided spectrum:

$$\overline{\omega}_W = \int_{\mathbb{R}^+} \omega W(X(\omega))d\omega$$

where $W(X(\omega))$ is a weighting function that (may) depend on the spectrum $X$. The "mean" $\overline{\omega}_W $ will depend on the choice of $W$. A code was given in  the answer to Finding the right measure to compare sound signals in the frequency domain.

You can base it on an absolute value, or something else (positive).  But it is important to have a normalized weight: it should, somehow, sum to one, like for a traditional center of mass. The reason is scaling homogeneity. The mean should not change if you multiply every $X(w)$ by the same quantity.

For orthogonal transforms, based on the energy conservation:

$$ \int [x(t)|^2dt = \int [X(\omega)|^2d\omega $$

it is very common to use a mean/center definition based on an energetic weight, ie
a square of the absolute spectrum. So  every often, the mean frequency is defined as:

$$\overline{\omega}_2 = \frac{\int_{\mathbb{R}^+} \omega|X(\omega)|^2d\omega}{\int_{\mathbb{R}^+}|X(\omega)|^2d\omega}\,.$$

This choice is very sensitive when dealing with time/frequency location uncertainty related to the Weyl-Pauli-Heisenberg inequalities.

But you can use other $L_p$ norms, with:

$$\overline{\omega}_p = \frac{\int_{\mathbb{R}^+} \omega|X(\omega)|^pd\omega}{\int_{\mathbb{R}^+}|X(\omega)|^pd\omega}\,.$$

As $p\to\infty$, this estimator converges to the location of a maximum, as demonstrated in Metric Spaces: Why $L_\infty$ selects the maximum value.

When $p=1$, this is called the spectral centroid, or "brightness" (of sound), or simply "centre of gravity of amplitudes". It is described with other moments in Chapter 6.1 Spectral shape description. When $p=2$, this is sometimes called the power spectral centroid.

I have not found yet useful references for other powers (still searching).


Scaling of a rect function,drawing of signals

$f(t)$ is just time-stretched $x(t)$ by a factor of 2. Meaning $f(t)$ will take the same value at $t$, which $x(t)$ was taking at $\frac{t}{2}$. So, if the original $rect(t)$ was extending from $-1$ to $1$, then the time-stretched $rect$ will extend from $-2$ to $2$. It can not go from $-3$ to $3$. One thing to keep in mind is that all  time operations (shifting and scaling) will be applied on the independent variable $t$ in $x(t)$. So,
$$f(t) = x(\frac{t}{2}) = (\frac{t}{2} + 2).rect(\frac{t}{4}) + 2\delta(\frac{t}{2} - 2)$$


Yeah, notation is not ideal.

It is not - he assumes that each of the $M$ antenna elements is connected to its own RF chain, i.e., there are also $M$ receivers available. If you have fewer receivers you need to modify your $A$, it needs to contain the response your $K$ receivers observe given a wave from a certain direction.

He doesn't put it but yeah, $F$ should have $N$ columns (and so should $X$ and $W$) that correspond to $N$ subsequent snapshots (i.e., samples in time) of the received RF data.

The received signal, i.e., the (basedband version of the) RF data at the ouput of the $M$ antennas observed in $N$ subsequent snapshots is the $M \times N$ matrix $X$. The matrix $F$ comprises only the waveforms itself, you could think of it as the signal "before" the antenna array. But that intuition is also slightly misleading. Point is that each antenna is observing a superposition of $D$ wavefronts corresponding to $D$ far-field sources. This can be written as $A F$.

Note that this data model is a gross oversimplification of the actual physics. It ignores near-field sources/reflections, multipath, diffuse components and many things more. It still works well enough in some applications. But the physical interpretation we can put into it has its limits.


Convert frequency response to impulse response

Here is a MATLAB example. Inefficient code but hopefully easy to read

%% Example: first order bandpass from 200 to 2000 Hz sampled at 44.1 kHz
n = 1024; % FFT length
fs = 44100; % sample rate
f0 = [200 2000]; % bandpass range
[b,a] = butter(1,2*f0/fs); % these are the difference equation coefficients
b0 = b(1); b1 = b(2); b2 = b(3);
a0 = a(1); a1 = a(2); a2 = a(3);  % for easier readability

% Step 1: build normalized frequency vector
om = 2*pi*(0:n/2)'/n;

% step 2: calculate the frequency response at every frequency
z = exp(j*om); % this is the "z" of the z-transform
% calculate the Z transform ,
H = (b0+ b1*z.^-1 + b2*z.^-2)./(a0+ a1*z.^-1 + a2*z.^-2);

% Step 3: make it conjugate symmetric
H = [H; conj(H(end-1:-1:2))];

% Step 4: inverse FFT
h = real(ifft(H));

% Step 5: verify against direct evaluation of the time domain difference equation
d0 = zeros(n,1); d0(1) = 1; % unit impulse
href = filter(b,a,d0);
err = 10*log10( sum((h-href).^2)./sum(href.^2));
fprintf('Error = %6.2f dB\n',err);


Channel model: LOS Component and Rician K-factor

For $0 < K < \infty$, the channel is a combination of both a deterministic component (i.e., LOS) and a fading component.

As the $K$-factor is the ratio of the energy in the deterministic Line-of-Sight (LOS) component to the energy in the aggregation of the random scattered paths (i.e., the fading component), higher  means that the channel is more deterministic.


Amplification of a speech signal in the recorded audio

I wouldn't recommend amplifying the low signal because that will increase noise significantly. As was pointed out above, it's better to downmix to mono.
Open your recording in a DAW (e.g. Audacity), delete the quiet channel, then copy the loud signal into the deleted channel. Voila!


Unexpected peaks in power density following downsampling and filtering

2 comments:

First, scipy signal decimate already low-pass filters the signal. It first filters the signal and then subsamples. So there is no need to pre-filter your signal before decimation. In your case, the signal is filtered at 125 Hz, which is half the resulting sampling frequency. (decimate doc here)

Second, a Butterworth filter is not an ideal filter. This filter will have a 3dB attenuation at the cutoff frequency, no matter the order of the filter. So you are expected to have components remaining between 100 Hz and 125 Hz.

If you want nothing detectable above 100 Hz, then design a filter with stop-band at 100 Hz and sufficient attenuation. This will result in a cutoff frequency below 100 Hz.

Depending on your requirements, an elliptic or Chebyshev filter will result in a narrower transition than a Butterworth of the same order.


In the first formula, $x[n]$ is written as if it were a sampled version of the continuous-time function

$$x_c(t)=A\cos(2\pi ft+\phi)\tag{1}$$

If you sample $(1)$ at $t=nT$, $n\in\mathbb{Z}$, you get the first formula. $T$ is the sampling period, or the inverse of the sampling frequency $f_s=1/T$.

If you compare the two formulas in your question you see that $fT=k/N$. Assuming that $k$ and $N$ are integers, then the second formula is a special case of the first one. The discrete-time signal obtained by sampling $(1)$ is not necessarily periodic. However if $fT$ is rational, as in your second formula, then $x[n]$ is periodic with period $N$.


Changing camera and image properties to calibrate a camera

You should not experience any major drawbacks.
Generally such techniques utilize adaptive thresholds and rely on gradient orientations (especially normalized ones), which are intrinsically invariant to uniform illumination and contrast changes. In other words, brightness changes will not effect your subpixel positions, if you are using an appropriate corner detection algorithm.


What is the physical meaning of the cross-spectrum of two signals measured at two locations and different moment?

I figured this out myself. The answer to this is that for the steady state, phase from the cross-spectrum is a contribution of both time and space. And if two signals are measured at the same time or location, the contribution of the other factor to phase will be eliminated. The conclusion is made based on formula deduction, which I will not discuss further here.


Help understanding FFT Phase: code and results provided

You phase isn't zero because your FFT width isn't an exact multiple of your sinewave's period.


I think, this will act like low pass filter with $\pi/2$ instead of $\pi/4$ cutoff frequency


Oversampling and decimation: What filter to use?

Should I, for example, do a low-pass FIR filter to squeeze out as close to 50Hz of signal bandwidth as possible?

Yes, that is exactly what you should do.  That is an extremely low data rate, so even with a wimpy processor I would think that should be able to do a pretty good filter.  Especially since you only need to calculate $\frac{1}{8}$ of the filter outputs.


1 - Too broad, please be more specific
2 - Most functions are supported, you can check the simulink blocks to see if they support code generation. That being said, even if a function is supported, it does not mean that it makes sense to have it in your code. You can also convert Matlab functions to C/C++ but there are restrictions.
3 - What you describe is a model-based design where your simulation model and your executable are based on the same source. I use model-based designs in power systems, I use a controller for simulations and I use the same controller for code generation.
It works well, reduces development time, helps eliminating bugs. However, it's not a silver bullet, it's not as easy as Mathworks wants you to think. You still need some coding skills.


Which Approach Is Better for Decomposing an Image into High Frequency and Low Frequency Components?

The bilateral filter is just that: a filter. It does not decompose the image into anything. It just attenuates "medium" frequency content in a clever manner. In other words, it smooths things that are mostly smooth further, and leaves sharp edges sharp. What you are left with has more low and high frequency content (on a relative basis) than what you started with.

The Fourier transform, on the other hand, actually does reveal the entire spectrum for an image. Unlike a filter, it does not change the information content of an image. It just rearranges it into the frequency, rather than spatial, domain. The bilateral filter does not reveal the spectral content of an image, so really only the Fourier transform (or one of its related transforms, like cosine and wavelet) can be used to get the high and low frequency components of an image.


system response: time vs frequency. Why do I get different magnitudes?

What you are trying to do is to compare the amplitude of your signal with an amplitude of FFT peak, expecting it is the same. The thing is that amplitude of your signal is varying and FFT is taking the "average" over whole duration.

In such case you should use time-frequency representation of your signal, such us spectrogram. It will allow you to observe how amplitude of particular frequency is changing over time. Since you work in MATLAB you can try following code:

win = hamming(256)
noverlap = 0
nfft = 512
X = abs(spectrogram(spost, win, noverlap, nfft))
% Scale by window and multiply by 2 to restore energy from negative freqs
X = 2*X/sum(win)
surf(X)

It will produce something like:

The maximum value of displacement is around $0.053$. Results can be tweaked even further by tweaking even more and using other methods.


Seeing spurious peaks in acquired vibration data; peak frequency changes with sampling rate

Can you please advise what could be possible reason for such spurious FFT peaks?

Most likely it's your measurement setup. Good  measurements are hard to make and require careful setup, characterization, and execution.

Determine the noise floor (spectral shape and level) of your data acquisition system, just short the input and grab few spectra. If that looks good
Determine the noise floor (spectral shape and level) of you sensor: attach the sensor to something heavy that's not moving and grab few spectra. If that looks good
Build a "calibration" setup. Attach your sensor to a well characterized, calibrated or programmable vibration source (a shaker, or maybe a electro mechanical driver) and measure different frequencies at different levels and compare extend the expected results. Determine the frequency and dynamic limits of your setup.

In each of the steps above set targets, requirements, and expectations BEFORE you make the measurement. These should come from the data sheets, general physical considerations and the specifics of your application.
Don't even look at any real measurement data until you have proven that your setup does indeed measure vibration with sufficient accuracy and signal to noise ratio over the required frequency range and dynamic range. If that all checks out and you still see unexpected data you can ask again.


What is the difference between `image normalization`, `contrast stretching` and `histogram equalization`?

I believe that they all do the same task of distributing the frequency of pixel intensities to a bigger range

They do.

  but can these terms be used interchangeably?

no, they can't. because the way they increase the dynamic range of the image diffears.

All this methods do a transformation of the pixel intensity $u$ of and image into a new pixel intensity $v$ $$v = f(u) $$ in order to improve the dynamic range of the image.

for constrast stretching (that as far as I know is a more general case of image normalization) the transformation is a piecewise linear function

$$ f(u) = \begin{cases}
      \alpha u & 0 \leq u < a \\
      \beta(u - a) + v_a& a \leq u < b \\
      \gamma(u - b) + v_b & b \leq u < L
   \end{cases}$$

Where L is the total number of gray levels.

On the other hand, Histogram equalization is a more powerful technique for image enhancement that aims to give the histogram of a given image the desired uniform shape. This is achieve by the transformation
$$v_{aux} = F_u(u) = \sum_{i=0}^u p_u(i) $$

$$v = int\left( \frac{v_{aux} - v_{aux_{min}}}{1 - v_{aux_{min}}} (L-1) + 0.5 \right) $$

where $p_u(i) $ is the normalized value of the histogram at the  grey level $i$


well i'm assuming you mean "conventional" DACs and not $\Sigma \Delta$ DACs.

in a conventional DAC (like an R-2R ladder or something), there are the micro errors that occur between neighboring DAC codes.  e.g. non-monotonicity.  i think the DSP solution to that is adding a teeny amount dither noise to the value that is output to the DAC.

there is a more macro nonlinearity that one can see in the DAC code to voltage curve.  the curve should be a straight line, but isn't exactly.  to linearize that, you need to compute the inverse curve and put that into a lookup table and process the output to your DAC through the lookup table before outputting it.

all this requires getting to know the specific non-linear faults of your conventional DAC first.  to do that requires an ADC (that the DAC output is connected to) that you trust or many repeated measurements with analog dither and various DC offsets added to the ADC input to smooth out the crap in the ADC.


In-Phase and Quadrature components of low-pass equivalent of a band-pass signal?

i will call the bandpass signal:

$$\begin{align}
 y(t) &= \Re \Big\{ (x(t) + j\hat{x}(t)) \, e^{j \omega_0 t} \Big\} \\
\\
      &= x(t) \cos(\omega_0 t) - \hat{x}(t) \sin(\omega_0 t)  \\
\end{align}$$

the complex-valued $x(t) + j \hat{x}(t)$ is the low-pass equivalent to the bandpass signal.

the real part $x(t)$ is the in-phase component and the imaginary part $\hat{x}(t)$ is the quadrature component.

the in-phase and quadrature components may have nothing to do with each other and then the spectrum of the low-pass equivalent is arbitrary with no symmetry implied.

if $\hat{x}(t)=0$, the low-pass equivalent is purely real and spectrum of it has hilbert symmetry with the real and imaginary parts.

if the in-phase and quadrature components are related by the Hilbert Transform,

$$\begin{align}
 \hat{x}(t) &= \mathscr{H} \Big\{ x(t) \Big\} \\
\\
            &= \lim_{\epsilon \to 0} \ \int\limits_{-\infty}^{t-\epsilon} x(u) \frac{1}{t-u} du \ + \ \int\limits_{t+\epsilon}^{+\infty} x(u) \frac{1}{t-u} du \\
\end{align}$$

then the spectrum of the low-pass equivalent is zero for all negative frequencies and has only a positive-frequency component.


Why Does the Rudin Osher Fatemi (ROF) Method Use Variational Methods for Image Denoising When Denoising Problems Are Not Boundary Value Problems?

This turned out to be easier than I thought. Image processing is a boundary value problem, and the boundary is the set of pixels along the edge of the image. The common notation for this boundary in the Vese Guyader book is $\partial(\Omega)$, so that is not the most obvious notation for an image. Of course $\partial$ is often the boundary operator in differential topology too. I think for the Rudin-Osher-Fatemi method, there is a Neumann boundary condition, meaning that the first derivative at the boundary pixels is not changing.
So the answer was very easy after all.


ANSI S1.11.2004 Octave-Band Filter implementation

That's actually rather simple.
A little bit of trial and error will do the trick. The standard specifies upper and lower limit for the transfer function of the filter. So you can just try a normal Butterworth bandpass and see what happens.
The picture below shows a 3rd octave design at 1kHz and a sample rate of 48 kHz. A 2nd order bandpass isn't quite a good enough for class 1 limits, but a 3rd order does the job just fine. The 4th (and higher) order works as well.
Once you have a design that works at 1kHz, you can apply the same design to every other center frequency as well. At very high frequencies you will see some bilinear frequency distortion. Depending on your application you want to adjust for this or not.
Things to tweak are filter order, filter type and exact location of the corner frequencies. In the example the "nominal" corner frequencies worked just fine.


what is the best technique to detect a given signal : energy detection, autocorrelation, and matched filter based sensing, CRC?

There is no best rule, it depends on what we know about the signal.

Energy Detector: The energy Detector makes a decision on presence or absence of a signal based on sum of squared samples. This rule comes from the fact that the signal to be detected is inherently assumed to be random following a wide sense stationary process with a known PDF, the likelihood ratio test then simply determines this rule.

Matched filtering and autocorrelation: The matched filtering is nothing but the correlation of the. Received samples with the known signal itslef, so they are inherently the same, the difference compared to energy Detector is that here the signal is assumed to be a known deterministic quantity. The matched filter is the optimal detector in low SINR scenarios.

CRC data: A check on CRC is an agreement between the receiver and tarsmitter in whether an error has occurred or not while transmitting data on a wired/optical/Wireless channels. By itself it is not a signal detection scheme. It only detects whether an error has occurred or not. Once an error is detected the receiver can then request the trasnmitter to send the data again by sending a NACK (not acknowledge) signal to the trasnmitter (there are other ways to handle it as well, like soft combining etc.)


Spectrogram power/magnitude normalization in analogy with image intensity normalization

there isn’t any one solution to this problem.  the application governs the solution.

There are a few issues associated with mapping a DFT bin from a STFT calculation to a pixel on a display.

In many applications, the spectra tends to fall off with increasing frequency.  The low frequencies tend to dominate in the signal. Spectral leakage can cause these low frequencies to mask the high frequencies. Picking an appropriate window is important. Whitening your signal prior to taking DFTs is usually a good idea and can help picking a proper window. The window influences your overlap fraction.

One can also do some short term averages (of power) in time of sequential spectra to help improve SNR for some classes of signals. This also compresses spectra in time which saves display bits.

One issue is gray scale or colors? This influences how many pixel bits you have.  Color gives you more bits to map to but one needs to tune a color map.

Color is useful for overlays on gray scale like cursors.

Using the log of intensity (dB) does a very good job of mapping the dynamic range out a FFT to the range of a display pixel.

A fixed point FFT can be easier to map to pixels over a floating point FFT.

One reason that there isn’t a standard display technique is that display technology keeps advancing.  I first did this sort of display with a total of 480 bits across.

One trick is to set all the pixels below a threshold (like a median) to zero. This cuts down a lot of clutter. Scene complexity can cut either way.

I would recommend that making your spectra flatter over frequency is worth the effort.

One should also think about the interactions a user will have with the STFT image. Things like zoom, drill down, harmonic cursors, measurements like time differences and frequencies, ......


Compressive Sensing Incoherence Principle

Say you have measurements $y=Ax$, with a basis $\Psi$, and $A\in\mathbb R^{M\times N}$. Then we can write $y=A\Psi\hat x=\Phi\hat x$.

I then believe direct computation of the coherence requires us to solve
\begin{equation}
\mu = \max_{i<j} \frac{| \langle \Phi_i,\Phi_j\rangle |}{\| \Phi_i\|\| \Phi_j\|}
\end{equation}
where $\Phi_i$ is the $i$-th column of $\Phi$, which I think is of order $\mathcal O(N^3)$.

Now, assuming the above is correct you can compute $\Psi_i =\text{ifft}(e_i)$, where $e_i$ is a unit vector, construct $\Phi=A\Psi$, and then with a double for loop find $\mu$.


How to find the group delay of IIR filters in Python/SciPy

I only need the total group delay, not spectrum of group delay.

Group delay is a spectrum, so this doesn't make sense.  The group delay is the derivative of the phase response of the filter, so in Python it can be calculated as

from scipy import signal
from numpy import pi, diff, unwrap, angle

w, h = signal.freqs(b, a)
group_delay = -diff(unwrap(angle(h))) / diff(w)

The units of phase are radians, and the units of frequency are radians/sec, so the units of group delay are in seconds.  Examples at 1 and 2

The group delay of a Bessel filter is maximally flat, so when you say "the" group delay, you probably mean the maximum passband group delay.  Like in this filter:

The passband has about 0.032 seconds of group delay.  Your Butterworth filter does not attempt to have a flat group delay in the passband, however, so it's up to you to decide what "the" group delay is:

Do you want the group delay at 0 Hz, or the maximum group delay at the cutoff frequency?  Or maybe at a specific frequency like the fundamental clock frequency of your digital signal?


Convolution: generating a continuous-time signal from a discrete-time sequence

I believe (2) and (4) are same. One can write $z(t)$ as

$$z(t) = \sum_k z[k]\delta(t-kT) = \sum_k \left(\sum_l g[l]h[k-l]\right)\delta(t-kT)$$

But $h[k-l] = h((k-l)T)$. Thus,
$$z(t) = \sum_l g[l] \left(\sum_k h((k-l)T) \delta(t-kT) \right) = \sum_l g[l] \left(\sum_k h(kT-lT) \delta(t-kT) \right) $$
or
$$z(t) = \sum_l g[l] h(t-lT)$$
which, by change of variables, is same as (2).


Calculating covariance matrix for MVDR beamforming

MVDR is a narrowband beamformer. For broadband signals it is usually applied for each frequency bin. That means that $\mathbf{R}_{xx}$ is frequency dependent. In other words, for each time you should have $M$ matrices, each one is $3\times 3$.

Now, since you usually cannot compute $\mathbf{R}_{xx}$ exactly, you perform covariance estimation $\tilde{\mathbf{R}}_{xx}$. For more information you can look here for example. In short, as you mention you perform averaging of the covariance matrix in order to improve the estimation. In order to adapt you estimation to time-varying statistics, usually exponential averaging is performed, i.e.

$$
\tilde{\mathbf{R}}_{xx}[n] = \alpha\tilde{\mathbf{R}}_{xx}[n-1] + (1-\alpha)\mathbf{x}[n]\mathbf{x}^H[n]
$$

for some $0<\alpha<1$ (usually closer to 1).


Fourier series - time and frequency domain confusion

The coefficients of the Fourier series that you have computed are, in effect, the spectrum of the periodic signal consisting of the sum of signals $f(t)$ delayed in time or advanced in time by integer multiples of $1$ second.  Mathematically, the Fourier transform of a periodic function has impulses in it with the impulse amplitudes being the Fourier series coefficients. The impulses occur at intervals of the period, in this instance $1$ second,
and a "plot" of the spectrum typically shows a bunch of vertical arrows (a common representation of impulses) with the amplitude marked just above the arrowhead. In your case, assuming that you have computed the complex form of the Fourier series,
you would get a series of such arrows extending from $n=0$ out to $\pm \infty$.
Note that your formula
$$
 \hat f_n = \frac{-(-1)^n}{2\pi ^2 n^2} + \frac{1}{2\pi ^2 n^2}
$$
is saying that $\hat{f}_n = 0$ when $n$ is even. Also, putting $n=0$ in the formula
makes $\hat{f}_0$ undefined which is incorrect since $\hat{f}_0$ is the
"DC value" which is $\frac{1}{4}$ by inspection.. So please do check your work.
In particular, note that while integrating by parts as you are likely
to be doing, the antiderivative
$$\int e^{-j2\pi nt}\,\mathrm dt = \frac{e^{-j2\pi nt}}{-j2\pi n}$$
cannot be used when $n=0$: $\int 1 \,\mathrm dt = t$.

Note: if you computed the sine/cosine form of the Fourier series, then you
should have found that the coefficients of the $\sin(2\pi nt)$ terms are all
$0$ and your spectrum will extend from $n=0$ to $n=\infty$ only.


Assuming MATLAB functions to be used, passband frequency 0.15 means $\omega_p = 0.15 \pi$ whereas stopband frequency 0.2 means $\omega_s = 0.2 \pi$ and since $\omega_s > \omega_p$ these settings qualify for a proper lowpass filter.

Also passband ripple being $3dB$ and stopband attenuation being $60 dB$ are valid choices (although more exact specifications requires the knowledge application requirements)

The equiripple option designs a filter in which the designed filter magnitude error is distributed accross the pass and/or stop bands in equal amounts so as to reduce the filter length to a minimum while still meeting the maximum of the allowed peak error criteria.


Is $y(t) = \cos(t) + x(t)$ a time-invariant system?

You've answered the question yourself already. So, does the response to a shifted input equal the shifted response to the original (non-shifted) input? If so, then the system is time-invariant. Otherwise, it's a time-varying system.


SBC was specifically created for this purpose. See http://en.wikipedia.org/wiki/SBC_(codec). I don't know what the current legal situation is.

Other alternatives depend on your application: what is your target bandwidth and what type of artifacts and errors can you tolerate and which ones not.


This is tricky. There are two things to consider:

Can the system be implemented as depicted?
Can the difference equation described by the system (i.e., the input/output relation) be implemented?

The answer to the first question is indeed 'no'. You can't implement the system with the structure given in the figure. The answer to the second question is 'yes'. I'll show you what I mean.

Defining two (time dependent) auxiliary variables $v[n]$ and $w[n]$ at the outputs of the two adders at the top left of the figure, respectively, the system is described by the following set of difference equations:

$$\begin{align}v[n]&=x[n]+G(x[n]-y[n])\\
w[n]&=v[n]+Aw[n-1]+Cw[n-2]\\
y[n]&=w[n]+Bw[n-1]+w[n-2]\end{align}\tag{1}$$

As you've noticed, these equation cannot be implemented directly, because in order to compute $y[n]$ you need $v[n]$, but in order to compute $v[n]$ you need $y[n]$. However, these equations can be 'unrolled' and rewritten in such a way that they can be implemented. This can be shown by using the $\mathcal{Z}$-transform of $(1)$ to convert difference equations into algebraic equations:

$$\begin{align}V(z)&=X(z)(1+G)-GY(z)\\
W(z)&=V(z)+W(z)(Az^{-1}+Cz^{-2})\\
Y(z)&=W(z)(1+Bz^{-1}+z^{-2})\end{align}\tag{2}$$

Now we can eliminate $V(z)$ and $W(z)$ from these three equations, and with a bit of algebra we arrive at

$$H(z)=\frac{Y(z)}{X(z)}=\frac{(1+G)(1+Bz^{-1}+z^{-2})}{1+G+(GB-A)z^{-1}+(G-C)z^{-2}}\tag{3}$$

For $G\neq -1$, Eq. $(3)$ corresponds to the following difference equation:

$$y[n]=x[n]+Bx[n-1]+x[n-2]+\frac{A-GB}{1+G}y[n-1]+\frac{C-G}{1+G}y[n-2]\tag{4}$$

which can be implemented without any problem. Note that $(4)$ implements the same difference equation as described by the depicted system, but it uses another structure such that there are no delay-free loops.

To help your intuition, throw away the block with the two delay elements including the multipliers and the two adders, and just consider the trivial system described by

$$y[n]=x[n]+G(x[n]-y[n])\tag{5}$$

There is a delay-free loop, but if $G\neq-1$ the system can be implemented without any problem by rewriting $(5)$ as

$$y[n](1+G)=x[n](1+G)\tag{6}$$

which (under the mentioned assumption on $G$) is of course equivalent to $y[n]=x[n]$. This principle also applies to much more complicated systems.


find the gradient direction at a point in a shape

Assuming you are using discrete inputs and that s(0..n-1) (or s(1...n) in Matlab) is an array of x,y points along the boundary of a 2d shape ...

Gradient(i) = (s(i).x, s(i).y) - (s(i-1).x, s(i-1).y)
Amplitude(i) = sqrt(Gradient(i).x^2 + Gradient(i).y^2)
Angle(i) = atan(Gradient(i).y / Gradient(i).x)

To make you calculations relative to (Xc,Yc) substract (Xc,Yc) from each point in s prior to calculations.

If your calculations are continuous, replace the above s(i)-s(i-1) with limit(s)/limit(sigma) and use the symbolic toolkit.


What is the Uncertainty or error due to sampling frequency?

Uncertainty IS NOT determined by the sampling frequency (i.e., 100Hz) only.
Uncertainty is mostly determined by both the Bandwidth of the signal you sampled (as stated by Marcus Muller) and the sampling frequency.
In fact, the uncertainty between sampled peaks can be solved by properly designed interpolation (or DAC) as long as the sampling theorem is met. The main reason is depicted as follows.
Although Nyquist-Shannon sampling theorem state that no theoretical reconstruction loss as long as the sampling frequency is twice larger than the bandwidth. But, YOU NEED TO RECONSTRUCT THE SAMPLED SIGNAL before the temporal analysis. Most of the analyzers forget to reconstruct (interpolate) their sampled signal while analyzing.
This sample code sample a simulated 2 Hz analog signal,x_a, (with sampling rate 1500 Hz) to a discrete-time signal,x_d, with 5 Hz sampling rate (nearly Nyquist rate) and then reconstruct this back,x_r.
%% This code simulate the AD/DA processing discussed in Chapter 4.8.3 [1]
close all; clear all;

%% parameters.
% analog
analog_fps = 1500;
analog_window_time = 3; %sec
t = 0:1/analog_fps: analog_window_time-1/analog_fps;

% digital
digital_fps = 5;
n = downsample(t,analog_fps/digital_fps);

% ADC: Quantizer
X_m = 1; % Range
B = 10;% Bit number.

%% Signal generation
freq_hz = 1; % Hz.
x_a_1 = 0.5*cos(2*pi*freq_hz*t+0.1);

% add a small high-frequency component as asked.
signal_freq = 2; %Hz
x_a_2 = 0.5*cos(2*pi*signal_freq*t+pi/2);
x_a = x_a_1 + x_a_2;

%% ADC
% Sampling
x_s = downsample(x_a,analog_fps/digital_fps);

% Quantizing (abs of input value should not over 1)
% x_d = Quantizing(x_s,B,X_m);  % A For complete ADC, a quantizing should
% be added here.
x_d = x_s; % For basic case, we skip the quantizing here.

%% DAC
% up sample / DAC
x_up = upsample(x_d,analog_fps/digital_fps);

% LPF (Reconstruction Filters)
h = intfilt(analog_fps/digital_fps,4,0.9);
%% Important
% please not the parameter 0.9, ideally should be 1 for Nyquist rate.
% 0.9 here is the ratio of Nyquist.
% Given known limit band signal, a shorter ratio can enhance SNR by oversampling.
%  (i,e, here I filterout the freq larger than 2.5(Nyquist rate) * 0.9 = 2.25Hz)

x_r = filter(h,1,x_up);
x_r(1:floor(mean(grpdelay(h)))) = [];
x_r = [x_r zeros(1,floor(mean(grpdelay(h))))];

%% Display

figure;
plot(t,x_a);
hold on;
plot(n,x_d);
plot(t,x_r);
title('analog signal (1500Hz) v.s. digital signal (5Hz) v.s. Reconstructed signal (1500Hz)');
legend('x_a','x_d','x_r');

If you forget to interpolate, the uncertainty of x_d is very high!! However, as long as you properly interpolate the x_d, you will obtain the x_r where the uncertainty is much much smaller.
To sum up, as long as the band-limited signal is sampled higher than the Nyquist frequency, properly reconstruction (interpolation) can theoretically reduce the uncertainty to a very small amount (which depends on the interpolation filters you designed). The interpolation filter sample I provided above is FIR and zeros phase filter, which means finite impulse response (i.e., no error propogation and phase distortion).
DO NOT FORGET to RECONSTRUCT YOUR SIGNAL, or you will overestimate the uncertainty.


The $S$ matrix is finally of $2\times2$ the size of the initial image. Gathering it is a bigger image is somehow misleading. You effectively get a $2\times 2$ matrix at each pixel of the original image, by gathering values from the $4$ submatrices at corresponding locations.

Structure tensor provides local information. Do not forget to use smoothed directional derivatives as well when computing $I_x$ and $I_y$. You can for instance compute them with Difference of Gaussians (DoG) operators.
More discussion in  How do I implement the Structure tensor (in Matlab)? or Structure tensor tutorial.


Interpret periods from fourier transfrom graph of periodic impulse signals

Assuming you are in the continuous-time domain, given the periodic impulse train with a period of $T_0$

$$ x(t) = \sum_k \delta(t - k T_0)$$

then its CTFT will be another periodic impulse train in the frequency domain
$$ X(\omega) = \frac{2\pi}{T_0} \sum_k \delta( \omega - k\omega_0)$$

Where the period $\omega_0$ of the frequency impulses is related to the period of the time impulses as

$$ \omega_0= \frac{2\pi}{T_0} $$

hence the time period will be given as:

$$ T_0 = \frac{2\pi}{\omega_0} $$

In your case, if $T_0 = 30$  then $\omega_0$ will be
$$ \omega_0 = \frac{2\pi}{30} = \frac{\pi}{15} $$

Conversely, given the frequency period as $\frac{2\pi}{30}$ then the corresponding time period will be $T_0 =  30$.

Note that there will be errors in your computation depending on how you simulate the impulse train or CTFT.


Finding linear convolution of two time series

There is a small mistake in your solution. Here $f_k$ and $g_k$ are infinite duration signal, ie. $n \to \infty$. You made a mistake that you took  n in $h_n$ and number of coefficients in $f_k$ or $g_k$ as same, but it is different.

$$h_n = \sum_{k=0}^{\infty} f_k g_{n-k} = \sum_{k=0}^{\infty} a^k b^{n-k} = b^n \sum_{k=0}^{\infty} \left(\frac{a}{b} \right)^k$$

$$h_n = b^n \cdot \frac{1}{1 - \frac{a}{b}}$$

$$h_n = \frac{b^{n+1}}{b-a}$$

If you want to finite the number of coefficients in $f_k$ and $g_k$, use other variables other than n.

Then n in the RHS of your equation will replace with the new variable.


Any basic paper on biomedical signal processing to practise with?

I would suggest you go to PLOS Biology and do a search for "Signaling MatLab".  There are many scientific papers there on various biology research topics.  Some even give the Matlab commands used.


The wavelet transform has the desired property of presenting most images as a sparse signal where most of the signal energy is concentrated in few coefficients.

A noiselet is structured to take such sparse representation and generate a noise-like signal from it. The coefficients resulting from the noiselet transform are spread out to almost equal magnitudes across the whole signal.

Why are they useful?

noiselets are are highly 'incoherent' with the signal representation in the wavelets domain. This means that noiselets can be used as a sensing matrix to recover sparse signals with high probability of success. Section 2.2 in this paper provides some mathematical details and an image recovery application.


Perform Convolution in Frequency Domain Using FFTW?

You need to zero pad x and h.

nZeros = length(x) + length(h) -1

Example in MATLAB:

clear all; close all; home

xSamples = 10;
hSamples = 7;

x = randn(1,xSamples);

h = randn(1,hSamples);

nZeros = xSamples + hSamples - 1;

X = fft(x,nZeros);
H = fft(h,nZeros);

x_conv_h = ifft(X.*H);

figure(1)
plot(real(x_conv_h))
hold on; grid on
plot(conv(x,h),'.r')
legend('fft based convolution','convolution')


Use samples of Fourier transform as DFT?

Check my answer on a similar question : H(z) as N-DFT

Only thing to add it that for both the systems, it is possible to compute $H[k]$ by figuring that the numerator and denomenator are as follows :

Numerator = $[1]$, Denominator = $[1, -\frac{1}{2}]$.
Numerator = $[1]$, Denominator = $[1, -2]$.

Take $N_b$-DFT of numerator and $N_b$-DFT of Denominator and point-wise division of these will give $H[k]$. Take IDFT to get $N_b$ samples of $h[n]$. Convolve with $x[n]$ to get output.

This is possible for stable and unstable both systems, because we are truncating $h[n]$, which will inherently make the system stable. The main requirement is that the system is causal.


Baseband signal processing algorithms would refer to any algorithm that is written after you have down converted your RF or IF signals to baseband(signal centered around 0 Hz).

Example : Suppose there is a difference in the RF carrier frequency between the transmitter and receiver . So , after downconversion at the receiver , there will be a frequency offset (equal to Tx Center Frequency - Receive Center Frequency) in the baseband signal . In this case , an algorithm to estimate/correct this frequency offset on the baseband signal can be termed as a baseband signal processing algorithm.

Other baseband algorithms could include estimation/correction of Timing Offset and IQ Impairements to name a few.


Problem with Kalman filtering accelerometer data

If there is any offset in the acceleration measurement, then a Kalman filter won't help without the addition of some velocity and/or position reference measurements.  Double integration of random-walk noise in the acceleration measurement will eventually lead to a reported position of your device outside the solar system.


Can harmonics be created or simulated a la carte

Without re-hashing all of the maths from RBJs response, I made a crummy little plug in that did more/less what you are proposing. It had a stupid simple GUI with one slide for nominal amplitude and several more to control harmonic levels. It was OK, but not as useful as one would hope for audio work. The idea was that given a sinusoid with a nominal input amplitude, it would calculate the polynomial required for the given harmonics.  The bummer is that audio doesn’t work that way, at all.  The sound doesn’t have any more harmonic ‘richness’ if you can call it that, it just sounds distorted in different ways.
To wit, my peripheral understanding of products advertising harmonic distortion use crossover filters to divide the signal into different bands, on which distortion is applied. By narrowing the bandwidth, you get something more closely approximating a sinusoid. Again though, my impression was that they were fine, but not exactly what I had in mind.
My experience was that as simple as ‘adding harmonics’ seems, it’s a pretty tough thing to implement. My apologies for not providing a more rigorous analysis, but figured some first hand experience might be relevant.


The desired frequency-response specifications

If you have to consider the whole frequency spectrum range in $[0,2\pi]$ for the design of the discrete-time filter, without assuming any type of symmetry, then you are considering the most general case of the filter and its coefficients will be complex valued. And nothing more can be said about them, unless you impose further constraints on the impulse and frequency responses of the filter.


What happens at the edges in fft2 filtering?

Quick answer: FFT2 assumes periodic boundary conditions. This means that performing the FFT2 of an image is equivalent to have your original image periodically repeated in the horizontal and vertical directions, compute the DFT, and then only retain the central region. The weird artifacts you are seeing are probably the result of juxtaposing opposite edges that are not related to each other.

There are ways to work around this by, e.g., smoothing these discontinuities by using MATLAB’s edgetaper function before filtering.


The open-loop transfer function

$$G(s)=\frac{2s-11}{s(s^2+s+1)}\tag{1}$$

given in your question has the following Nyquist plot:

which is different from the one in your question.

Note that in any case it's important to know how the trace is closed, i.e., how the points $\omega=0^-$ and $\omega=0^+$ are connected. In this case they are connected by a semi-circle with infinite radius in the left half-plane (the dashed curve in the figure below):

[generated with WolframAlpha].
This is a consequence of the pole of $G(s)$ at $s=0$.

The Nyquist criterion is actually simpler than you made it look in your question. The number of clockwise encirclements ($N$) of the point $-1+j0$ equals the number of right half-plane zeros of $1+G(s)$ ($Z$) minus the number of right half-plane poles ($P$):

$$N=Z-P\tag{2}$$

Since $G(s)$ has no poles in the right half-plane ($P=0$), the number of clockwise encirclements equals the number of zeros in the right half-plane. From the Nyquist plot we see that $N=1$, so there is one zero in the right half-plane, which corresponds to one pole of the closed-loop transfer function in the right half-plane. Consequently, the system is not stable.

This is also easily verified by directly computing the poles of the closed-loop transfer function

$$Q(s)=\frac{G(s)}{1+G(s)}=\frac{2s-11}{s^3+s^2+3s-11}\tag{3}$$

roots([1,1,3,-11])
ans =

  -1.2833 + 2.3182i
  -1.2833 - 2.3182i
   1.5667 + 0.0000i

which shows that there is one pole with a positive real part, as predicted from the Nyquist plot.


Well i don't have the Symbolic Math Toolbox.  but i think i know why, and it's not a really good reason.  usually we electrical engineers mean, when we say "the unit step" function, we mean what math and physics people mean by the "Heaviside step function but, there is a convention issue with what this function is with zero for its argument.

the root definition of the Heaviside unit step is

$$ H(x) \triangleq \int_{-\infty}^x { \delta(t)} \, \mathrm{d}t $$

where $\delta(t)$ is the Dirac unit impulse function.  But this does not define $H(x)$ at $x=0$.  It's infinitesimally thin, but it could be anything finite at zero.

in continuous time, it is usually most convenient to define the Heaviside step function in terms of the sign function which is:

$$\operatorname{sgn}(x) \triangleq \begin{cases} -1, & x < 0 \\ 0, & x = 0 \\ 1, & x > 0 \end{cases}$$

then the Heaviside comes out as

$$\begin{align}
 H(x) &= \tfrac12 \big(1 + \operatorname{sgn}(x) \big) \\
\\
&= \begin{cases} 0, & x < 0 \\ \tfrac12, & x = 0 \\ 1, & x > 0 \end{cases} \\
\end{align}$$

because this, $H(x) + H(-x) = 1$, is always true, they are complementary.

the problem is extending these definitions and desired behavior to discrete $x$.  then we like to define the discrete impulse as the Kronecker delta which DSPer's express more simply as

$$ \delta[n] \triangleq \begin{cases} 0, & n \ne 0 \\ 1, & n = 0 \end{cases}  $$

and the discrete unit step as

$$ H[n] \triangleq \sum_\limits{m=-\infty}^{n} \delta[m] $$

but that makes this discrete unit step come out to be:

$$ H[n] = \begin{cases} 0, & n < 0 \\ 1, & n \ge 0 \end{cases}  $$

that is not quite the same as

$$\begin{align}
 H[n] &= \tfrac12 \big(1 + \operatorname{sgn}(n) \big) \\
\\
 &= \begin{cases} 0, & n < 0 \\ \tfrac12, & n = 0 \\ 1, & n > 0 \end{cases} \\
\end{align}$$

now here is the stupid problem:  MATLAB is not a continuous-time tool.  it is a discrete-time tool.

so what they are doing with the discrete Heaviside unit step function appears to me to be another one of MATLAB's mistaken conventions of definitions.  they fucked up the index origin convention which, among other problems fucked up the fft() convention putting the DC term into X(1) rather than have the index represent frequency.  they fucked up the order of polynomial coefficients in polyval() and polyfit().  and they fucked up the definition of the discrete Heaviside unit step function.

EDIT: actually, they defined the Heaviside function correctly, because it has a continuous argument.  but they (and we) need to define a discrete function for the unit step that does not have $\tfrac12$ in it.


How to extract prosodic cues from a wav file using Python

This project Speech Signal Processing Toolkit (SPTK) provides several features you are looking for. Here is a good wrapper around it https://github.com/r9y9/pysptk
This wrapper is using a slightly different version of SPTK, but installation is straightforward.

pip install pysptk

1) I would use kaldi or another speech recognition package on phones output and a VAD system to compute the voicing rate.

2) For autocorrelation, convert the waveform to mel-cepstrum with the function mcep then compute the autocorrelation with c2acr

3) Here is an example to compute the Pitch http://nbviewer.jupyter.org/github/r9y9/pysptk/blob/master/examples/pysptk%20introduction.ipynb


This is wrong
Because if your carrier frequency is less than your message frequency, then you will have trouble in filtering .


The logo you are looking at is symmetric. For example, you could rotate it by 120 degrees and get the same pattern. If you take a look, many of the "errors" are not really errors. The same pattern appears in your image many times. It is called self-similarity.

You could do several things:

Change descriptors to avoid confusion between features that are similar in an image.
Instead of finding one match for keypoint, find 3 matches and test multiple combinations by fitting a transformation and checking which fits better.
Limit the transformation by rotation angle, mirror, etc. For example, angle must not be larger than 60 degrees


From the perspective of filtering, a common LTI system will weight each components in the input signal: (1) weighting means changing the magnitude values (it is apparent in multiplication of two complex numbers), thus it is often to observe the magnitudes (this does not mean the phase is not important; (2) of course, phase information is also part of the whole signal, and although the phase information may be also changed, but the phase is not apparent; (3) we can consider a FIR filter, and compare the changes of both magnitude and phase.


it's because all of the filter coefficients (which is the impulse response $h[n]$) are positive and there is no "ringing" in the impulse response.

so, assuming the impulse response is scaled so that all of the $h[n]$ add to 1, the output of the gaussian filter will be a weighted average of the input.  averaging is a process that removes high-frequency variation.  so what is left are the low frequencies, like DC.  that's what a low-pass filter does.


pairwise error probability and union bound of space time code

The union bound can indeed result in a BER higher than 1. As you say, it is just an upper bound, and it becomes tighter as the SNR increases. For low SNR, the union bound can be quite loose.


The left-hand side is

$$
\sum_{i \% M=0}\sum_{j\%M=0} x(i/M,j/M)y(i,j)
$$

where the sum can only be done over all i,j which are dividable by M (otherwise x(i/M)=0). Hence, by substitution $i=Ma, j=Mb$ we get

$$
=\sum_{a\in\mathbb{Z}}\sum_{b\in\mathbb{Z}}x(a,b)y(aM,bM)
$$

which equals the right-hand side.


Long burst of high amplitude noise removal

If I understand you correctly, then your issue is connected with transient that is present when switching on your device, because it is a RLC circuit. Most likely these spikes are not changing over time (because values of passive components are constant), if so, then you can:

Create the template of such spike and use the normalised correlation for detection. If it's above a given threshold, then substitute samples with linear/spline or whatever interpolation.
Use the Matched Filtering approach (which is almost the same as above approach, but offers better SNR performance) and remove detected samples accordingly

If your spikes are of random nature, then there is many methods that are allowing to remove any crackles and clicks. You might refer to the literature for such solutions:

Oudre L. - Automatic detection and removal of impulsive noise in audio signals
Knapp M., Bashir R. - Audio Click Removal using LPC
Godsil S., Rayner P. - Digital Audio Restoration
Alvarez E., et al. - Detection of clicks using sinusoidal modelling for the confirmation of the clicks


non linear equalizer vs linear equalizer

This is not a complete and comprehensive answer but provides some examples showing the key features in a linear and non-linear equalizer and also clarifies consideration for feedforward vs feedback structures in equalizers.

Linear equalizers are typically feed-forward in structure, and linear as the output is a linear combination of scaled and delayed versions of the input. The type of linear equalizer varies depending on the algorithm used to determine the scaling in adaptive cases (LMS, RLS, RLS-Gradient, etc) but can also be used in non-adaptive applications (such as to compensate for an analog filter response in a receiver). Linear equalizers are typically the simplest to implement but have the worst performance compared to non-linear equalizers, especially in channels with deep selective frequency nulls in the passband (from frequency selective fading).

Since a typical multipath channel distortion can be modeled as a linear feed-forward FIR system (multiple copies of the transmit signal arrive at different delays), it would be tempting to use a recursive IIR system as the linear equalizer to generate the inverse response.

If it was certain that the channel being equalized would always be a minimum phase system, then a recursive equalizer could be used. A minimum phase channel would have all "echos" trailing the primary or strongest response, and all zeros would be inside the unit circle. Thus the inverse filter as an IIR would have all poles inside the unit circle.

Typical channels are mixed with both leading and trailing echos, which means they have zeros outside the unit circle. An equalizing IIR filter as the inverse filter would end up having poles outside the unit circle and be unstable and therefore recursive feedback structures cannot be used (as a linear equalizer solution).

Decision Feedback Structures however are recursive and frequently used. They are the common choice for time domain equalization for frequency selective channels with deep frequency nulls. Another example of a non-linear equalizer is the MLSE Equalizer which represents the best that can be done but is the most computationally intensive. More often the MLSE equalizer is used in simulation as a benchmark to compare various equalizer algorithms targeted for implementation.

A block diagram of the Decision Feedback Equalizer is shown below. What specifically makes it non-linear is that hard decisions are made on the waveform prior to being input into the feedback structure. The process of making hard decisions is a non-linear process. The output of the feedback structure is no longer dependent on the linear combination of scaled and delayed versions of past outputs (specifically the effective scaling of each sample changes based on the value of that sample when the decision is made - this is non-linear).


Fractional Delay Filters and Cutoff Frequences

Okay, I think I've learned a sufficient amount about interpolation theory to mostly answer this. Interpolation can be viewed in two equivalent ways. The most obvious is, for each point you are trying to compute the interpolated value for, the multiplication/sum of the source samples with the interpolation kernel centered on the interpolation point. Alternately, you could view the signal itself as the convolution of a set of "hidden" coefficients with the very same convolution kernel; in this case to interpolate between samples you multiply/sum each source coefficient by the point on that sample's corresponding interpolation kernel where the interpolation point falls. The former model is what you actually use for interpolating a signal, but the second makes obvious some things that are less clear in the first.

For upsampling, the equivalent low-pass filter cutoff will always be at the Nyquist rate of the source signal (at least for an ideal low-pass filter), meaning that the interpolation kernel will be fixed regardless of upsampling ratio. Furthermore, it means that the coefficients for a given interpolation point will always have fixed distance on the interpolation kernel curve, which is valuable information for optimization (e.g. if you quantize your frequency ratio to the same resolution as the kernel LUT, you can arrange the kernel LUT such that all coefficients for a single interpolated sample can be fetched at once, rather than requiring numerous array accesses per interpolated sample). And as a consequence of these, the prefilter (if needed) will be fixed as well.

For downsampling, the cutoff frequency must be the Nyquist rate for the target sampling rate, meaning the interpolation kernel as well as the prefilter must be dilated proportionately when applied to the source samples. I am not entirely certain, but I think this would mean interpolation methods whose kernels cannot be dilated (e.g. matrix-based B-spline interpolation) cannot directly be used for downsampling, in which case you would have to use such an algorithm to upsample to some integer multiple of the target frequency then low-pass filter and decimate (obviously simultaneous polyphase low-pass filter and decimate would be an optimization).


What's the expression of the noise for a RF signal?

Yes, you can represent your bandpass noise as:

$ n(t) =n_I(t)\cos(2\pi f_c t)-n_Q(t)\sin(2\pi f_c t),$

where $n_I(t)$ and $n_Q(t)$ represent the inphase and quadrature components of the noise signal. In matlab you can generate this sequence as

n = sqrt(Npwr/2) *(randn(1,len)+j*randn(1,len))

The scaling factor is different than in the baseband case because both the inphase and quadrature components have unit power (before scaling), so when added together the total power is equal to 2 and therefore the $\sqrt{2}$ factor is needed in the scaling factor to get unit power.


Compressing A Range Of Decibel Values to One "Representative" Value

The answer provided by Jason R in the comments --

  "dB values" are really just a way of representing values on a logarithmic scale. Arithmetic-averaging values on a logarithmic scale is equivalent to geometrically averaging the original values (in linear scale). This probably isn't what you want. You probably would like a measure of total power in each band, which you would get by just summing the power spectrum bins across each band. You could divide by the number of bins, but if the bands aren't equal-width, that will add some bias to the result.

does a good job of identifying the problem but I'm not sure it is a correct solution (or maybe I'm just reading it wrong). If a summation is to occur, it should occur on the raw magnitudes of the bins (before conversion to dB).

The issue is this

3dB + 3dB != 6dB

to do a summation on the decibel scale is more complex than simply adding the values together. It can be done but if so you should do some reading on the dB scale and the nature of adding relative power. Doing the addition when still in a linear domain is easiest.


Designing butterworth filter with pole placement

I think I found the problem. Everything was fine except for the denormalization step for arbitrary cut-off frequencies wc. I changed the following lines to obtain the correct state-space matrices (and to handle the case n=1):

...
if (n == 1)
    a = [1,1];
else
    p = ...
    ...
end

a = real(a);
a = a./(wc.^(n:-1:0))   % new!
a = fliplr(a);
...

% last_row = (wc.^(N:-1:1)).*last_row;    % removed!

if( n == 1)
    this.x = last_row;
else
    k = ones(1,N-1);
    this.x = diag(k,1) ;
    this.x(end,:) = last_row;
end

...

this.y(N,:) = 1/a(end);    % changed!

...

I used the function ss2tf to compute the polynomial coefficients from the state-space matrices computed above, and the coefficients were correct for different values of wc. This means that the matrices are correct. If the output signal is still different from what you expect, the problem must be in the filtering routine.


How to Improve the Kalman Filter for Tracking the Periodic Motion of a Car?

You basically have 4 models here:

Accelerating to constant speed.
Moving at constant speed.
Decelerating to zero speed.
Standing.

So the basic solution is building the 4 models and switching using Hard Switch between them.
Yet there is a smoother framework to handle smooth transition between them called Interacting Multiple Model (IMM) Kalman Filter.
Using the IMM Framework (Which is basically adding weighing step for the models) with the knowledge of the Input as described in other answers will give you a much better results.


Ok, here's the simple method illustrated in Python:

# Correlation is convolution with one input reversed
corr = fftconvolve(img1, img2[::-1, ::-1])

# Get coordinates of peak in x, y pixels
peak = unravel_index(corr.argmax(), corr.shape)
# peak = (199, 189)
# Calculate offset from center point
print array(peak) - 199 # Offset is [  0 -10]

So it finds the vertical offset (10 pixels), but not the horizontal, as you've shown (bottom left output below).  This is because the cross-correlation output is stronger in the center when more of the input images are overlapping (gets dark toward the edges even if the inputs are identical).  (The found peak is shown with the red cross.)

And here's the method I described in the comments.  I have no idea if there's a mathematically rigorous foundation for this or if it has a name, but it produces the right results.  It may or may not be the same thing as normalized cross-correlation (I haven't bothered to decipher what Wikipedia is saying):

weighting = fftconvolve(ones(shape(img1)), ones(shape(img2)))
corr = corr / weighting # Normalize
peak = unravel_index(corr.argmax(), corr.shape)[::-1]
# peak = (224, 189)
print array(peak)-199 # Offset of [ 25 -10]

So this method works.  It normalizes every element of the output array by the number of overlapping pixels that were used to calculate that element, so that when only the thin edges of the images are overlapped, it still produces a strong output.  Now the offset from the dim square is not swamped by the white bar (lower right of image).

Note that while this works perfectly with your example images, this method doesn't work as well for real images, since the edges of the output tend to vary wildly compared to the middle, and can produce spurious peaks, but if you know the offset is constrained to less than the image width, you can look for the peak only in realistic locations near the center, ignore the edges, and it works.

Wikipedia suggests a method using windowed images which I thought might improve on this, but unless I'm doing it wrong, it doesn't:

h = hamming(shape(img1)[0]) # 1D Hamming window
ham2d = outer(h, h) # 2D Hamming window

h1 = img1 * ham2d
h2 = img2 * ham2d

corr = fftconvolve(h1, h2[::-1,::-1])
peak = unravel_index(corr.argmax(), corr.shape)[::-1]
# peak = (216, 189)
print array(peak) - 199 # Offset is [ 17 -10]

The windowing does too much, and causes it to find an erroneous peak that's in between the previous two results.


Is calculus background a must for learn DSP?

With Google's search engine and Google books, I investigated that someone, without calculus background, can learn DSP. I expected learn from you that for avoid laborious searching. I'm very sorry for did not learn that from you. Regards.


How was this result on discrete Fourier series achieved?

Note that $x[n]$ has period $N$ and the sum you refer to has $NM$ terms, so the sum has the form

$$\sum_{n=0}^{NM-1}x[n]f[n]=x[0]f[0]+x[1]f[1]+\ldots+x[N-1]f[N-1]+x[0]f[N]+\ldots +x[N-1]f[2N-1]+x[0]f[2N]+\ldots=\\=x[0]\big(f[0]+f[N]+f[2N]+\ldots+f[MN-N]\big)+x[1]\big(f[1]+f[N+1]+\ldots+f[MN-N+1]\big)+\ldots+x[N-1]\big(f[N-1]+f[2N-1]+\ldots+f[MN-1]\big)=\\=\sum_{n=0}^{N-1}x[n]\sum_{l=0}^{M-1}f[n+lN]$$

With $f[n]=e^{-jk2\pi n/NM}$ we obtain

$$\frac{1}{NM}\sum_{n=0}^{NM-1}x[n]e^{-jk2\pi n/NM}=\frac{1}{NM}\sum_{n=0}^{N-1}x[n]\sum_{l=0}^{M-1}e^{-jk2\pi (n+lN)/NM}\tag{1}$$

For the other sum over $y[n]$ you can do exactly the same. The final result is obtained as follows. The right-hand side of Eq. $(1)$ can be written as

$$\frac{1}{NM}\sum_{n=0}^{N-1}x[n]\sum_{l=0}^{M-1}e^{-jk2\pi (n+lN)/NM}=\frac{1}{NM}\sum_{n=0}^{N-1}x[n]e^{-j(k/M)2\pi n/N}\sum_{l=0}^{M-1}e^{-jk2\pi l/M}\tag{2}$$

The last sum over $l$ in Eq. $(2)$ is given by

$$\sum_{l=0}^{M-1}e^{-jk2\pi l/M}=\begin{cases}M,&k=mM,\;m\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{3}$$

Consequently, for $k=mM$, Eq. $(2)$ evaluates to

$$\frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{-j(k/M)2\pi n/N}=\frac{1}{N}a_{k/M}\tag{4}$$

where $a_k$ are the discrete Fourier series (DFS) coefficients of $x[n]$.
For other values of $k$ the expression in Eq. $(2)$ equals zero. Exactly the same manipulation can be done with the sum over $y[n]$.


Quantization SNR of sine wave doesn't match 1.761 + 6.02 * Q

Some issues here:

Your SNR formula only applies to full scale sine waves, your sine wave has -6dB amplitude so your SNR will be 6 dB lower
The formula also implies rounding, not truncation, that's another 6 dB
You use a frequency that's a small integer divider of the sample rate, that means you are just repeating the same samples over and over again and don't get enough sample coverage to get a statistically meaningful result.
Your SNR analysis in the frequency domain is needlessly complicated and prone do imprecision and masking errors. Just do it directly in the time domain.

Here is how this would look in Matlab
%% quantization noise of a 16-bit sine wave
fr = 975.3; % something odd
n = 8192;
quant = 2^15;
% make the sine wave
y0 = sin(2*pi*(0:n-1)'/n*fr);
% quantize
yq = round(quant*y0)/quant;
% noise
yNoise = yq-y0;
% SNR
fprintf('SNR = %6.2fdB\n', 10*log10(mean(y0.^2)/mean(yNoise.^2)));

Technically, you would also have to work around the fact that the positive max amplitude of the sine might clip, but for large quantization this makes no meaningful difference.
EDIT
It may be good to review where the formula comes from and what it actually means. It all starts with the quantization noise. If we quantize and round than the quantitation noise is uniformly distributed between $[-0.5 \delta ,-0.5 \delta ]$, where $\delta$ is the quantization step. For truncation it would be uniformly distributed $[0,\delta ]$. For rounding, the resulting noise powers are
$$ P_{round} = \frac{\delta ^2}{12}, P_{trunc} = \frac{\delta ^2}{3}$$
If $B$ is the number of bits, then for a signed signal, we simply have $\delta = 2^{B-1}$ and hence we get the noise level for 16 bit rounding as
$$P_{16} = \frac{2^{-30}}{12} \Rightarrow -101.1 dB $$
That's true for ANY signal that has a reasonable wide distribution of samples. Since a full scale sine wave has a power of $-3dB$ the resulting SNR would be $-98.1 dB$


Biquad filter from Matlab to C using "tf2sos"

What I can see from the code is the following mapping between the array coeffs[] and the filter coefficients $a_i$ and $b_i$ (note that $a_0=1$ and is not stored in coeffs[]):

$b_0=$coeffs[2], $b_1=$coeffs[1], $b_2=$coeffs[0]

$a_1=$coeffs[4], $a_2=$coeffs[3]

For the time being choose the scaling equal to $1$ (scalefactor=0) to check if everything else works. The final scaling per section will depend on the order of the sections, and on the goal of scaling (minimization of the possibility of overflow, or minimizing quantization noise).


Pitch detection of two following music notes (C0 and C0#)

There are many ways to do pitch detection (seee e.g. https://ccrma.stanford.edu/~pdelac/154/m154paper.htm ) and which one is the most suitable depends on your application requirements and constraints. FFT isn't great if you need a lot of accuracy of accuracy (in tune, out of tune?) or good resolution at low frequencies.

You also need to be aware that for many instruments, the harmonic have higher energy than the fundamental, so the peak in the FFT isn't always the pitch.


Determine the Z-Transform for the following sequence: $ |n|(\frac{1}{2})^{|n|} $

It's always good to do a sanity check on such results. E.g., you could try to see that $X(1)$ equals the sum of all time domain samples:

$$X(1)=\sum_{n=-\infty}^{\infty}x[n]\tag{1}$$

With $x[n]=|n|\left(\frac12\right)^{|n|}$ it is clear that the result of $(1)$ must be positive. However, for the solution from the manual you get $X(1)<0$, whereas for your solution you obtain $X(1)>0$. So I think you can be confident that your solution is the correct one.


Matched filter loss due to frequency offset

Short answer; I'll try to come back and improve this later.

Recall that a matched filter is equivalent to a sliding correlator. At each time instant, the sliding correlator multiplies the last $T$ (where $T$ is the pulse duration) of the input signal by the conjugate of the pulse shape and coherently integrates the result over the pulse duration. If there is no frequency offset, then the result after multiplying by the conjugated pulse shape is a constant (DC) value, so you don't observe any loss.

If there is a frequency error, then the result after multiplying by the conjugate of the pulse shape is a sinusoid whose frequency is commensurate with the frequency offset. The scalloping loss is incurred is due to the fact that you're integrating the value of said sinusoid over a time window of length $T$. Therefore, the amount of attenuation can be well-modeled by the frequency response of a boxcar filter of length $T$:

$$
\text{rect}\left(\frac{t}{T}\right) \Leftrightarrow T \text{ sinc}\left(\frac{2\pi f T}{2}\right)
$$

You can make the following observations from this relationship:

As the pulse duration $T$ increases, the frequency response of the correlator rolls off more quickly; for a given frequency offset $f$, you'll experience more attenuation from the ideal peak value.
As the frequency offset $f$ increases, the frequency response of the correlator rolls off more quickly as well.

These combine to tell you what you probably already intuitively know; the longer you coherently integrate a signal (i.e. the longer the pulse duration in your case), the more sensitive the result will be to frequency offset.


What are some good signal and image processing blogs worth following?

There are many for different subjects -

Efg's algorithm collection : http://www.efg2.com/Lab/Library/ImageProcessing/index.html
DSP Forum : http://www.dsprelated.com/
Data compression - http://datacompression.info/
About rendering - http://www.realtimerendering.com/portal.html
For all research papers - http://ieeexplore.ieee.org/Xplore/guesthome.jsp
Resources on Mp3 and Audio - http://www.mp3-tech.org/programmer/docs/index.php
Steve on Image Processing - http://blogs.mathworks.com/steve/
Image Processing and Retrieval http://savvash.blogspot.com/
Accelerated Image Processing - http://visionexperts.blogspot.com/
The Digital Signal Processing Blog - http://centerk.net/dspblog/
Noise & Vibration Measurement Blog - http://blog.prosig.com/
Image Processing with Matlab, Open Blog - http://imageprocessingblog.com/


Realtime fixed point IIR filter implementation in C

Although this seems like a remarkably simple questions, it requires a remarkably complicated answer.
I don't think there is a "one-size" fits all solution. The best choice of algorithm will depend on what noise you can tolerate and the type of low pass (steepness & frequency). For example at 44.1 KHz sample rate a 4th order Butterworth at 10 kHz is fairly straight forward, whereas a lowpass at a 100Hz is a royal pain. In essence it depends on how close your poles are to the unit circle.
Quantization and rounding error of IIR filters are usually transfered to the output weighed by the pole-only transfer function. A 4th order Butterworth 10 kHz low pass filter has a worst case noise amplification of only 5dB, so that is not much of a problem.
However at a 100 Hz low pass (again 4th order BW) the noise gets amplified by a whopping 75 dB. If you use Q15 math, your basic noise floor is maybe at -100dB or so. After the filter, your signal to noise ratio will only be 25 dB.
That's one of the reasons why fixed point IIR filters are fairly complicated. IF you need low cutoff frequencies and half way decent signal to noise ratio, then the basic algorithms will not work. You need to look into double precision math and/or error spectrum shaping or related methods.


As applesoup says in the comments the term
$$
\mathbf{u}(kL+i)e(kL+i)
$$
is a vector, not a single value (some integer).

Why do you think it's a scalar?

To answer your question: no, it's incorrect to state that all the algorithm does is add a constant to all elements of $\hat{w}(k)$.

As per this slide from here

the block LMS error is still a scalar. It's just that the errors across multiple blocks are used to form the updates.

Another way to look at it is that $\Phi(k)$ can be written as:
$$
\mathbf{u}_k \mathbf{e}_k
$$
where $$\mathbf{u}_k =
\left[
\begin{eqnarray*}
&\mathbf{u}(kL+0)
&\mathbf{u}(kL+1)
&\mathbf{u}(kL+2)
&\ldots
&\mathbf{u}(kL+L-1)
\end{eqnarray*}
\right]$$
is an $N_u \times L$ array and
$$\mathbf{e}_k =
\left[
\begin{eqnarray*}
&e(kL+0)\\
&e(kL+1)\\
&e(kL+2)\\
&\vdots\\
&e(kL+L-1)\\
\end{eqnarray*}
\right]$$
is an $L\times 1 $ vector.

As a result $\mathbf{u}_k \mathbf{e}_k$ is an $N_u \times 1$ vector.


Without knowing what you find hard, it is difficult to recommend a specific book. Some are hard to read to. So, oldies but goldies, I suggest you to start with video lectures, and for instance by godfather A. V. Oppenheim

  Digital Signal Processing

and tell us where you have trouble (and why).


Causal and Non-causal filtering. Strange behavior

the output causal filtered signal inside the red box seems to be very uncorrelated with the zero-phase filtered signal.

So let's look at this statement. Suppose the original, unfiltered signal is $x[n]$. Filtering this using filtfilt with a filter with impulse response $h[n]$ gives:
$$
y_{ff}[n] = h[-n] \star h[n] \star x[n]
$$
and filtering with just filter is:
$$
y_f[n] = h[n] \star x[n]
$$
so
$$
y_{ff}[n] = h[-n] \star y_f[n]
$$
Regardless of what $x[n]$ is, the 'correlation' between $y_f$ and $y_{ff}$ will look like the anticausal (perhaps anticausal + constant term) $h[-n]$.


Confusion regarding discrete signal value?

They can very well have decimal values, represented as either fixed point or floating point numbers.


What is the correct gain of an RRC Filter?

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


Can you calculate the absolute/amplitude values of a Fourier transform without complex/imaginary?

Yes.  You can do the calculation.  But, no, it's not a shortcut.

If you want to obtain your DFT magnitudes fast (e.g. using the O(NlogN) FFT algorithm), then you will need the complex math (or its computational equivalent).  This is because the cosine and sine forms of the DFT basis vectors are orthogonal, and you need both to completely be able to represent arbitrary input waveforms.  The sine portions of the basis vectors are represented by the imaginary components of an FFT's complex result.  So the full complex result of an FFT needs to be computed to end up with enough information for a magnitude or amplitude result.

If a slow DFT computation is acceptable, then you could just compute the cosine components using strictly real arithmetic and then separately compute the sine components also using reals.  Then compute the hypotenuse of the real triangle for the amplitude.  No complex numbers needed.  But this method is O(N^2), which would take a lot longer.

For just magnitude results (computed slowly), N Goertzel filters is another possibility, and may save in trig function calls.


You can shift your spectrum in the time domain by multiplying it with a complex exponential. It is the same way you implement an AM modulator.

Let $A(t)$ be an unspecific baseband signal. You multiply in time domain with a complex exponential and the result will be the bandpass signal:
$$
y(t) = A(t)e^{j2\pi ft} = A(t)*\cos(2\pi ft) + j*A(t) \sin(2 \pi ft)
$$
Since cos and sin both contain only a single frequency component, a fourier transform (if you look at the real part) will result in two impulses at $\pm f$ and the spectrum of $A(t)$ shifted to $\pm f$.

The same applies to your problem: Multiply your time domain signal with a complex exponential of the desired frequency offset and you will observe your spectrum to be shifted down.


The formula is premised on the wavelet being analytic, or being nonzero only over non-negative frequencies: ${\hat\psi} (\omega < 0) = 0$. (Note all wavelets also have ${\hat \psi (0)}=0$ per the admissibility criterion). From Daubechies et al,
$$
\begin{align}
\int_0^\infty W_f(a, b) a^{-1} da & = \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{0}^{\infty}{\hat f}(\xi) \overline{ {\hat \psi}(a\xi)} a^{-1} da d\xi  \tag{1a} \\
&= \frac{1}{2\pi} \int_0^\infty \int_0^\infty {\hat f}(\xi) \overline{ {\hat \psi} (a\xi)} e^{jb\xi} a^{-1} da d\xi  \tag{1b} \\
&= \int_0^\infty \overline{ {\hat \psi}(z)} \frac{dz}{z} \cdot \frac{1}{2\pi} \int_0^\infty {\hat f}(\xi) e^{jb\xi} d\xi  \tag{1c} \\
\end{align}
$$
where a rather confusing change of variables was replaced with a clearer one, $z=a\xi,\ dz=\xi da$. Note that this c.o.v. is only valid if integrals share limits, which is enabled by $\psi$'s analyticity; otherwise, $z$ would integrate $-\infty $ to $\infty$, but $z$ includes $a$, which only spans $0$ to $\infty$. An identity was also utilized in splitting integrals.
Assuming $f$ is real, we have ${\hat f(\xi)} = \overline{\hat f (-\xi)}$, thus from $0$ to $\infty$,
$$
\frac{1}{4\pi} \Re e \left\{ \int_0^\infty {\hat f(\xi)} e^{jb\xi} d\xi \right\} = f(b),
$$
and setting $C_\psi = \int_0^\infty \overline{ {\hat \psi (z)}} z^{-1} dz$ thusly derives
$$
\begin{align}
f(b) &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty W_f(a, b) a^{-1} da \right\}  \tag{2a} \\
     &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty f \star \overline{\psi_a(b)} \frac{da}{a} \right\}  \tag{2b} \\
     &= \boxed{ \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty \left< f(t), \psi(t) \right> \frac{da}{a} \right\} }  \tag{2c} \\
\end{align}
$$
Since $W_f$ sweeps $b$ over all $t$, above equivalently holds for $f(t)$.

The paper you linked uses $a^{-3/2}$ -- Indeed, as it assumes an L2 norm on the wavelet, whereas MATLAB uses L1; further info can be found here.
What about the $\delta$? -- The $\delta$ is from MATLAB writing the expression as a special case of the more general inverse, as it shows above near it, with $\left<f, g\right>$, which is derived here.
and the $2 \Re e$? -- Good question - perhaps MATLAB's $C_{\psi, \delta}$ is a bit different, but per above definition we clearly divide by two, not multiply. If in doubt, try both.

Interpretation: it helps to observe that $b$ in all derivations can be replaced with a constant - so the relation holds pointwise; if we integrate along $a$, i.e. a column of the 2D time-frequency representation, we recover the function at that point, $f(b)$.
Notice in the double-integral iCWT we use a dual wavelet to "undo" the inner product with the CWT wavelet, analogous to $e^{+j\omega t}$ in $\mathcal{F}^{-1}$; here, we don't use any inverting kernel - the transform undoes itself when integrated over all scales. This implies some symmetry property on the wavelet $\psi$ in encoding information across scales - though I'm unsure how this exactly works out; the most important is analyticity which permits this in the first place, so one may search there for an answer.

Discretized case with exponentially distributed $a$
From same page of linked paper, the discretized linear case reads
$$
f(b) \approx \frac{1}{2} \Re e \left\{ C_\psi^{-1} \sum_k W_f (a_k, b) \frac{(\Delta a)_k}{a_k} \right\},  \tag{3}
$$
with $(\Delta a)_k$ taking role of $da$. For the linear case, $(\Delta a)_k = \Delta a$, but not for the log case, as the increment size increases for greater $k$. That's the only distinction between log and linear.
In practice, however, it's better to directly use the discretized log-scaled transform - which derives as follows:
$$
\begin{align}
& z = e^a, \ dz = a da \\
& \frac{da}{a^2} = \frac{d \ln (z)}{z} \tag{4}
\end{align}
$$
Suppose $z(a) = 2^{a/32}$. Note that the log-plot of an exponential scale is linear, so the difference $\Delta(\ln (z))$ or $d(\ln (z))$ is constant, and equal to $\ln (z(1))da$. To better see this, let's derive directly for this example, using $n_v$ for "number of voices":
$$
z(a) = 2^{a/n_v} \rightarrow dz(a) = z \frac{\ln(2)}{n_v} da = \ln(a_0)da \tag{5}
$$
When discretizing, $da=1$, and we end up with
$$
\boxed{f[] \approx \frac{\ln (a_0)}{2 C_\psi}  \Re e \left\{\sum_{k=1}^{K} \frac{1}{a_0^k} W_f (a_0^k, .) \right\} } \tag{6}
$$
same as Mallat Eq 4.67 (note his derivation is for double-integral, but same scaling argument applies here).


It depends, but convolutional codes like these are usually decoded using maximum likelihood sequence estimation, typically through the Viterbi decoder.
The design methodology for any decoder doesn't start with an electronics implementation of the encoder (that is something that's highly specific to the speed and technology you want to use), but with a general design of the encoder, usually in shape of the polynomial that your circuitry implements.
Draw the state diagram to the encoder, then draw the trellis diagram to that, and then you gain something you can implement as viterbi decoder.
I can't see a case where someone would ever actually implement this encoder using discrete 74xx series logic gates like you do – I presume this is just a representation of the logic.
You can implement encoding logic very easily in software, running on a microcontroller; if you need 100s of megabits per second in throughput, you'd implement this in an FPGA, but then, this encoder wouldn't, from a code design point of view, be what you use to begin with, probably. So, sure, this is also something you can easily implement using logic gates.
Decoding is always more complex than encoding. Decoding convolutional codes is kind of inherently something that you'll need to test hypotheses while doing; this will be no fun / make no sense to implement in discrete logic components (gets way too large, and thus, too bug-prone very quickly), you'd do that in software, really, or on an FPGA.
Just like for your previous question: you can't understand how this works by implementing it in hardware first without reading any theory. So, wikipedia: convolutional codes, Trellis diagram, Viterbi decoder.


Matlab: How do I implement the Expectation operator for time series

I'm not sure I understand the question. If I'm missing something, leave a comment and I'll try to edit/expand/clarify as needed.

First of all, "mean" and "expected value" are related concepts, but are not the same thing. You find the mean of a set of samples, or observations. If you sampled/observed properly, the mean will be close to the expected value, which is the theoretical, "true mean" of the underlying random process.

When you sample a random process, you end up with a set of random numbers. If you sample the same process again, you end up with a different set of numbers. So, you can think of each sample as a random variable (a black box that produces random numbers); the random process is then a sequence of random variables.

These random variables may be correlated; that means, intuitively, that if one of them has a certain value, you can use it to predict the values of some of the other samples. For example, you may find that if one sample is positive, its neighboring samples tend to be positive too.

In your question you say that you have one set of samples, y(t). In general, you don't have enough information to generate the correlation matrix. However, you could do it like this:

M = length(y);
C = zeros(M,M);
for i = 1:M
    for j = 1:M
        C(i,j) = mean(y(i)*y(j));  % Note that the mean is useless here,
                                   % since y(i)*y(j) is scalar
    end
end

To find a closer approximation to the actual correlation matrix of the underlying process, you need to gather many observations. Let's say that you have stored the observations in the rows of a matrix Y (every row is a different set of observations). Then you could do:

M = size(Y,2);
C = zeros(M,M);
for i = 1:M
    for j = 1:M
        C(i,j) = mean(Y(:,i).*Y(:,j));
    end
end

To test this, let's try sampling a Gaussian process of zero mean and variance 1. Every time I look at the process, I'll gather three samples, and I'll look at the process 1000 times:

Y = randn(1000,3);

If I run the code above with matrix Y, I obtain this correlation matrix:

0.9598    0.0590    0.0238
0.0590    0.9675   -0.0286
0.0238   -0.0286    1.0055

which is expected (the "true" correlation matrix is the identity).

Regarding calculating a time shift in Matlab, it's helpful to arrange vectors in a table as follows (I've assumed $t=1,2,\ldots$ and $y(t)=10t$):

t    | y(t)  |  y(t-1)  | y(t-2)
     |       |          |
1    |  10   |  --      | --
2    |  20   |  10      | --
3    |  30   |  20      | 10
4    |  40   |  30      | 20
...

It's clear from the table that, to find yn=y(t-n), you cut n elements from the end of y, and shift the time vector:

tn = t(n:end);    % shifted time vector
yn = y(1:end-n);  % shifted signal

The following code provides an example, which doesn't assume integer time or integer delay:

t = 0:0.1:3;  % time vector
y = exp(t/5); % y(t)
n = 0.8;      % time delay
tn = t(t>n);  % delayed time vector
yn = y(1:length(tn));  % delayed signal
plot(t,y,'r',tn,yn,'g');  % plot

This produces the following plot:


Are the RLS filter and Kalman filter gradient methods?

Recursive least squares (RLS) filters don't use gradient descent. As their name suggests, they use a least-squares fit to determine the optimum coefficients at each time step. Via clever formulation of the filter structure, one can use the calculations done from time step $n$ to recursively calculate the updated coefficients for time step $n+1$ without having to do the full least-squares fit again.


HINT: Your partial fraction expansion is wrong; you should always check it by recombining the factors. The correct expansion is

$$X(z)=\frac{z^{-1}-\frac12}{\left(1-\frac12 z^{-1}\right)^2}=\frac{\frac32}{\left(1-\frac12 z^{-1}\right)^2}-\frac{2}{\left(1-\frac12 z^{-1}\right)}\tag{1}$$

You should try to derive this result yourself.

BUT: note that you don't even need partial fractions here because if you know (or look up) the $\mathcal{Z}$-transform correspondence

$$\frac{az^{-1}}{\left(1-a z^{-1}\right)^2}\Longleftrightarrow na^nu[n]\tag{2}$$

then you can directly apply it by writing $X(z)$ as

$$X(z)=2\frac{\frac12z^{-1}}{{\left(1-\frac12 z^{-1}\right)^2}}-z\frac{\frac12z^{-1}}{{\left(1-\frac12 z^{-1}\right)^2}}\tag{3}$$

Finally, note that the given result is wrong. It should be

$$x[n]=2n\left(\frac12\right)^nu[n]-(n+1)\left(\frac12\right)^{n+1}u[n+1]\tag{4}$$


what is matlab function is taking as sign of second order IIR filter?

The filter function assumes this form (for your case):

$$
\tilde{H}(z) = \frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}
$$

where the $1$ on the bottom is from the assumption that $a_0=1$, and call it $\tilde{H}(z)$ so not to be confused with your $H(z)$. Since you have written the denominator with the signs of $a_1$ and $a_2$ inverted, then you just have to fix that before using filter.

Suppose someone wrote the transfer function from your question down $\big(H(z)\big)$, handed you the coefficients, and said use MATLAB to find the output. The only thing you need to do before is to invert $a_1$ and $a_2$:
filter(..., [a0, -a1, -a2], ...).

When using besself, it returns the coefficients in the form according to the transfer function $\tilde{H}(z)$, and you can then directly use those coefficients.


The autocorrelation matrix is a Toeplitz symmetric matrix that can be formed as:

Rxx = toeplitz(xcorr(x));

and that should give you a sensible (square) matrix that can be passed to eig.


parseval for a continuos but limited signal

Sure. You can can just integrate in the time domain from $0$ to $T$ since the area outside  $[0, T]$ is 0.

Please note that you still must integrate from $-\infty$ to $+\infty$ in the frequency domain since finite support in the time domain implies infinite support in the frequency domain (and vice versa).

Things are different if the function is periodic.


Interpolation is the process of adding zeros in between your samples and then low-pass filtering to get rid of the aliases that result.  Linear interpolation is just one out of many different possible filters.  Linear interpolation filters can always be represented by FIR filters.  In one-dimension the interpolation by two linear interpolation filter is $[\frac{1}{2}, 1, \frac{1}{2}]$.  There are better low-pass filters.

It is not clear from your quote what they meant by a "wider" filter.  I'm guessing they mean "one with more taps", but it could also mean "one with a wider passband".


You're missing a step: pulse shaping. The sequence [.7+.7i, ...] is a sequence of pulse amplitudes. As such, it should not be considered to be signal; in consequence, it does not have a spectrum.

Say your sequence of amplitudes is $A=a_0, a_1, a_2,\ldots$. Now, you need to specify your symbol rate (how many pulses per second you'll transmit). Call this rate $R_p$ and let $T_p=1/R_p$ be the pulse interval. Finally, you need to specify a pulse $p(t)$ that has this property: $$\int_{-\infty}^\infty p(t-kT_p)p(t-lT_p)dt=\begin{cases}1,\textrm{ if $k=l$}\\0,\textrm{ if $k\neq l$}\end{cases}$$for $k$ and $l$ integers. Now you're ready to create your actual signal: $$s(t)=\sum_i a_ip(t-iT_p).$$ This signal:

Has a spectrum; specifically, its bandwidth is the same as the bandwidth of the pulse $p(t)$.
Can be up-converted to a carrier frequency.
Can be converted to an analog signal and physically transmitted over a medium.

For completeness, I'll explain how to recover the amplitudes $a_i$ from $s(t)$ (ignoring noise for simplicity). The receiver can estimate $a_k$ by doing this operation: $$\hat{a_k}=\int_{-\infty}^\infty s(t)p(t-kTp).$$ You can easily verifty that this is the case, using the properties of $p(t)$. Such an operation can be implemented with a matched filter.


Estimating Plane Pose Without Knowledge of Intrinsic Camera Parameters

I think there is a mistake in those calculations.

How can you assume that $K$ has only a scaling effect? The principal point is integrated in $K$, it's not just a matrix of focal lengths.
How come is this true?

$$det(K)K^{-T} = \left[\begin{matrix} 1&0&0\\0&1&0\\-u_x&-u_y&1 \end{matrix}\right]$$

You don't end up with such equal numbers in the diagonal if you compute $det(K)K^{-T}$.

To me, your question is a bit unclear and I guess some assumptions are fishy.


first of all, if i were the OP, i would fix the notation a bit.  $x(t)$, $y(t)$ are continuous-time signals and $t\in\mathbb{R}$ is a continuous argument.  $x[n] \triangleq x(nT)$ and $y[n] \triangleq y(nT)$ are discrete=time signals and and $n\in\mathbb{Z}$ is a discrete argument, an integer.  square brackets mean that an integer goes into them.  the normal convention is that $T=\frac{1}{f_s}$ is your sampling period.

a simple moving-average filter with equal-weighting in the average is

$$ \begin{align}
   y[n] & = \frac{1}{N} \sum\limits_{i=0}^{N-1} x[n-i] \\
        & = \frac{1}{N} \left(x[0] + \sum\limits_{i=1}^{N} x[n-i] - x[n-N]\right) \\
        & =  \frac{1}{N} \sum\limits_{i=1}^{N} x[n-i] + \frac{1}{N} \left(x[0]  - x[n-N]\right) \\
        & =  \frac{1}{N} \sum\limits_{i=0}^{N-1} x[n-1-i] + \frac{1}{N} \left(x[0]  - x[n-N]\right) \\
        & =  y[n-1] + \frac{1}{N} \left(x[0]  - x[n-N]\right) \\
\end{align} $$

that is representing both the non-recursive implementation and the recursive.  the recursive implementation works as long as you subtract exactly what you earlier added, something that might not happen in a floating-point context.  both are FIR, even though the recursive implementation is something we would call a Truncated IIR,  still an FIR.

to get the frequency response, apply the Z-Transform to either form, and solve for $Y(z)$ in terms of $X(z)$.  then compute the transfer function:

$$ H(z) \triangleq \frac{Y(z)}{X(z)} $$

then substitute $z = e^{j \omega}$  and solve for the $\omega$ such that $\left|H(e^{j \omega})\right|^2 = \frac{1}{2}$.  that $\omega= 2\pi f T$ is your -3.01 dB cut-off frequency.

try doing that math yourself.  we'll be here to help guide you.


Obtain Envelope and Pitch curves of my Singing Voice

A zero-padded FFT is not the best method for accurately estimating musical pitch in (near)real-time.  I suggest researching other pitch detection/estimation algorithms (lots of literature and research papers), such as interpolated autocorrelation, AMDF, YAAPT, cepstrum, harmonic product, phase vocoder, and etc., and see if one of those algorithms better meets your exact requirements.
Once you have a pitch estimate, you can sum the waveform energy over a single or small integer number of pitch periods, as that will include all the harmonics, even if fundamental frequency energy is low or missing.  Using an integer number of periods will reduce the amounts windowing artifacts (scalloping or so-called "leakage") that can appear in a non-integer-period-multiple sized FFT magnitude.
Hilbert methods of measuring envelopes work best for narrow-band signals, but where one does not know or can not extract the exact periodicity.  And a human singing voice is rarely, if ever, narrow-band.
The "attack" portion of a pitched vocalization may not only be rich in harmonics, but temporarily ambiguous in pitch (e.g. not purely periodic, even disregarding envelope).


To remove doppler for equalizer doppler does it need decision directed like carrier recovery?

Doppler offsets are typically small fractions on the symbol rate and as such do not significantly effect the task of the equalizer which can perform under such carrier offsets. The equalized waveform is then passed through timing recovery and carrier recovery synchronization loops.


What is the amplitude range of a spectrogram?

Recall how a spectrogram is computed. You place a window of length $N$ around a sample of interest in your signal $x$ and then compute its DFT. This corresponds to one column of the spectrogram. In symbols, a fixed column is given by a DFT operation:

$$
X(k) = \sum_{n=0}^{N-1} x(n) w(n) e^{\frac{-j 2 \pi k n}{N}}.
$$

I believe you want an upper and lower bound on the magnitude of the spectrogram. An obvious lower bound is $0$. To find an upper bound:

\begin{eqnarray}
|X(k)| &=& \left | \sum_{n=0}^{N-1} x(n) w(n) e^{\frac{-j 2 \pi k n}{N}} \right|\\
&\leq& \sum_{n=0}^{N-1} |x(n)||w(n)|\cdot 1\\
&\leq& \sum_{n=0}^{N-1} 1\\
&=& N.
\end{eqnarray}
where I assumed that the window function has a maximum value of $1$.

This means that the magnitude of any element in the spectrogram matrix must be between $0$ and $N$ where $N$ is the window length.

Note that there are other ways of scaling a spectrogram which will change these calculations. Some people plot squared magnitude instead of magnitude, or $20\;\log_{10}(\cdot)$ of the magnitude. Some people like putting a $\frac{1}{N}$ normalization. Matlab provides 'psd' and 'power' options and scipy provides 'density' and 'spectrum' options which scale things differently.


As an example of what I mean, let's have a look at a simple sinusoid (though I see you're actually assuming double gamma functions; I may rework this later to use those).

This models the signals as:
$$
y_p(t) = \sin(\omega t + \phi_p)
$$
so the only difference between them is the phase $\phi_p = (p-1)\pi/6$.

The top plot is the mean signal $\bar{y}(t)$ and the bottom plot shows the Pearson coefficient between $\bar{y}$ and $y_p$ for $p=1\ldots 6$ (so $n=6$).

As you can see, the Pearson coefficient varies quite widely even in this simple example.

R Code Below

# 26438

Nsigs <- 6
T<- 1024
omega <- 2*pi*0.00982734982
t <- 0:(T-1)

y <- matrix(,nrow = T, ncol = Nsigs)
phi <- c(0, 1*pi/6, 2*pi/6, 3*pi/6, 4*pi/6, 5*pi/6 )

for (idx in 1:Nsigs)
{
  #phi[idx] <- runif(1)*2*pi
  y[,idx] <- sin(omega*t + phi[idx])
}

mean <- rowSums(y) / Nsigs

total <- matrix(,nrow = T, ncol = Nsigs+1)
total[,1] <- mean
total[,2:(Nsigs+1)] = y

pearson <- cor(total,use="complete.obs", method="pearson")

par(mfrow=c(2,1),pty="m")
plot(t,mean, type="l")
title("Mean value")

plot(1:6, pearson[1,2:7], pch=19, col="blue")
title("Pearson coefficient")


What algorithms can automatically determine a 3D scene from one or a few 2D images?

SLAM(Simultaneous Localization and Mapping) algorithms can be used to for 3D reconstruction. They offer solutions for both monocular as well as stereo cameras. With single camera they estimate depth with few images and reconstruct the scene. You can find some of the open source solutions here.

Real time 3D reconstruction can done using ORBSLAM and it is also opensource.Here is the related paper for monocular SLAM.