SNR calculation with Noise Spectral Density I am having difficulty in calculating the SNR for the following scenario:

In a cellular network

Transmit power=36dBm
Noise spectral density=-174dBm/Hz (AWGN)
Bandwidth=10MHz

I have some other parameters such as pathloss and Rayleigh fading, but I think my biggest problem is the conversion of dB to power.

The comment above already gave you a roadmap to get to an answer, but I'll flesh it out a bit:

SNR is the ratio of signal power to noise power. To calculate the signal power at your receiver, you'll need to take the path loss and fading model into account, as you noted.
The noise power at the receiver is described by a (flat) noise power spectral density and receiver bandwidth. To calculate the total noise power over that bandwidth, you simply multiply the amount of power per Hertz times the width of the band. When power is specified in logarithmic units (dB), you need to first convert it to a linear scale. Recall the relationship between mW and dBm:
$$
P|_{dBm} = 10 \log_{10}\left(\frac{P}{1\text{ mW}}\right)
$$

That is, to convert a linear-scale power quantity to dBm, take ten times the base-10 logarithm of the ratio of that power quantity to 1 milliwatt. For example:

30 dBm: 1 W
0 dBm: 1 mW
-30 dBm: 1 uW
-60 dBm: 1 nW
-90 dBm: 1 pW

To go the other way, invert the calculation:
$$
P = 10^{\frac{P|_{dBm}}{10}} * 1\text{ mW}
$$
So, to calculate the total noise power at your receiver, you would convert the noise power spectral density to linear units using the above equation:
$$
S_n = 10^{\frac{-174}{10}} \frac{\text{mW}}{\text{Hz}} = 3.981 * 10^{-18} \frac{\text{mW}}{\text{Hz}}
$$
then multiply by the bandwidth to get the total amount of noise power:
$$
P_n = BS_n = \left(10 * 10^{6}\text{ Hz}\right)\left(3.981 * 10^{-18} \frac{\text{mW}}{\text{Hz}}\right)
$$
$$
P_n = 3.981 * 10^{-11} \text{ mW}
$$

Then, use your model to predict the received signal power, take the ratio, and you've got your SNR.


Numerical Calculation of impulse response out of frequency response Currently I work on my Masterthesis which deals with structural dynamic simulation. As a part of it I have to filter a Signal with a FIR-Filter. From that Filter I know the Frequency Response (picture).

The code is written in Python and I'm not sure how to apply the IFFT since I have just a real input.

Thank you for your answer robert bristow-johnson!

In python there is function called irfft which stands for inverse real fft. So a ifft for real input. But it seems that something goes wrong with the time signal. I'm aware of that I have no phase in this case. But I guess for Filters in Signal processing it is pretty common to just have the frequency response and calculate the impulse response out of it? So there should be a right way to do it...
Unfortunately I'm no expert in Signal processing and at the moment I'm pretty much stuck at this problem.

UPDATE: I think the correct way is to add just zero phase (so +0j) to the frequency response. For the ifft the first bin has to be 0 valued, from 1 - N/2 i add my frequencies as I have it from the response, and from N/2 - N I add the conjugate complex numbers (as the phase is zero the values are just in  descending order).
Hope this will do the trick.

I would suggest the frequency sampling design method, which results in a linear phase FIR filter. You sample the desired magnitude on an equidistant grid, add a linear phase term, and then you compute the impulse response by taking an inverse DFT. The resulting impulse response is usually windowed. I've explained the details in this answer, which also includes some Matlab/Octave code. Note that windowing is not done in that simple code fragment.


Using neural networks for rotation invariant OCR I have video with text which can be rotated in any direction. If I trained a neural network on images of each letter at all possible degrees of rotation, would it be possible to classify the characters accurately? Alternatively would another technique like HAAR, or KNN be better?

Aha, this is a quite funny story that I've heard about CV. Are you a bio guy? Any way, here are my suggestions.

If you are a bio guy and just want to finish this project ( I mean successfully identify each insect in a video frame ), go for barcode, QR etc. They are labels though their contents are not directly readable by eyes. However, you will have a much better chance to handle all kind of problems you might encountered when you use manual labels, e.g. occlusion, error correction, lighting etc. There are libraries can decode these labels, and you simply call them.

In this way, you will save a lot of energy to not working on something you are in lack of experience, but use the state-of-the-art technologies.

If you are not a bio guy or these manual labels are not something can be changed, you should first locate these labels and then recognize them. Donot mess up these two things. The first task is detection, and the second one is recognition. In other words, you need to do two things:
a). Train a text/none-text classifier. In this classifier you can use those rotational invariant features.
b). Train a character classifier ( e.g. using NN ) to decide the class of a character candidate you found from (a).

By the way, I donot believe you need to train a NN accepting character at an arbitrary angle. Instead, all you need is to train a NN accepting normal characters, and normalize a candidate using image moments or other techniques before you feed it into the classifier. Or simply rotate your candidate to a certain degree (e.g. for every 10 degree) and then feed them into your classifier and pick the result with the highest confidence.


Matching guitar signal with tablature, near real-time algorithms? I'm looking for candidate algorithms that might work for a use-case comparable to Rocksmith's analysis system. To compare a clean electric guitar signal with a tablature reference.

A higher latency is acceptable, as well as compromises in accuracy, but the goal is to have near real-time information on whether what you actually played, is within a threshold of what you were expected to play.

The output would be a probability value per expected strum (or pluck or tap). Which is then matched against a configurable threshold.

My current line of thinking is to take multiple algorithms and predetermine which algorithm should be used per segment of the song.

One for chords (compare lowest note & 12-bin PCPs?)
One for sustained / single notes (FFT cross-correlation?)
One for fast or polyphonic riffs (Supervised training? Classification models? Don't know.)
One for bends... etc.

My two question are:

Is a (context-sensitive) combination of sub-algorithms required for such a use-case?
Which algorithms would make good candidates for the individual tasks, given this use-case?

Bonus question:
do any "master" algorithms already exist to do such switching / tweaking of sub-algorithms?

Extra details,

What is considered acceptable processing time depends on the reference. For example if you're instructed to play a sustained single note for 3 bars, it's perfectly acceptable to do a slow analysis of the plucking envelope.

However if you're instructed to play a solo, alternating picking and hammering, 1/16th notes and bends. You're working with shorter samples, more human error and residual noises while playing. But we still need to attribute a binary decision if it was "close enough" for each note within, say ~200ms after each note was supposed to be released.

The guitar signal comes directly from a dedicated USB guitar audio interface, making it single instrument and a good quality signal.

The tablature is highly detailed. Assume a modern guitar pro format. Making relatively accurate training samples possible, but they will not contain the characteristics of the exact guitar, pickups, noise, etc. that the guitar signal would have.

Assuming that the tablature is well defined; let's say it's a MIDI file complete with pitch bend commands.  you can make a synthesizer where each note is actually a tracking comb filter tuned to the midi note (and offset by the pitch bend command) with notches at each integer harmonic (including the fundamental).  if those tracking comb filters processes a sound file playing that MIDI file and the two are properly aligned in time, the tracking comb filters should null out most of the input audio, even polyphonic.

Then magnitude-square-LPF the result and a minimum value means a pretty good match.

Just an idea.  No equations, but I have posted about pitch detection before. And I thought I did about tuned comb filters, too (you need a precision delay line to tune the comb filter to whatever fine-tuned pitch that is specified).


Filter Banks and FFT, are they similar conceptually? I'm using Filter Banks (Narrow Band Pass) created using the functions in Matlab. However,  I wanted to know if there are fundamental limitations like that of FFT which also affect the Filter Banks. Is there any co-relation between the two?

Basically, I don't know if they share the same problem of time and frequency being inversely related. If they do share this problem, this means that we cannot have higher time and higher frequency resolution at the same time.

The answer is yes, the discrete Fourier transform (DFT), of which the fast Fourier transform (FFT) is an efficient implementation, can be interpreted quite literally as a filter bank. Recall the formula for the DFT:

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

The "filterbank" interpretation for the above equation goes something like this:

In order to get the $k$-th DFT output $X[k]$, we start with the input signal $x[n]$.
Multiply $x[n]$ by a complex exponential function $e^{\frac{-j2\pi k n}{N}}$. By the frequency-shifting property of the DFT, this has the effect of shifting the content at frequency $\frac{2\pi k}{N}$ down to zero frequency.
Convolve the result with a filter that has the impulse response:
$$
h[n] = \begin{cases}
1, && 0 \le n \le N-1 \\
0, && \text{otherwise}
\end{cases}
$$
This filter (known as a moving average or a "boxcar") has a lowpass characteristic.

Thus, we can view the DFT as placing a bank of uniformly-spaced critically-sampled filters across the signal's spectrum. Each filter has the same frequency response as an $N$-sample moving average filter, centered at frequencies $e^{\frac{j2\pi k n}{N}}, k = 0, \ldots, N-1$. I described this similarly in an answer to a previous question here.


How to obtain digital IIR filter coefficients for known z-plane zero-pole plot? My clue is just to expand brackets and obtain coefficients of numerator and denominator of my transfer function. Am I right?

This might look trivial, but with the pole-zero diagram alone one can get the transfer function with a scale (gain) ambiguity.

That is, all transfer functions $(Kb,a)$, where $K$ is a constant, have the same pole-zero representation. For this reason, we usually use pole-zero-gain terminology.


Audio processing FFT of chunks of audio Im new to Matlab. I need to split the time dimensions of an audio track into frames of 10mS and for every 10mS frame have to compute the FFT to generate the spectrograph.
I need help on this.

http://www.mathworks.com/help/toolbox/signal/ref/spectrogram.html

This pretty standard fare. For audio I would recommend a hanning() window and an overlap of half a window length so

x = randn(44100,1); % some random noise sample of a second
S = spectrogram(x,hanning(512),256); % short term FFT

Assuming a sample rate of 44.1 kHz this gives you frames about 11.6 ms length. You can get to exactly 10ms but that requires FFTs that are not a power of 2.


Camera calibration using checkerboard instead of 3D calibration object What is the purpose and advantage of using a checkerboard than a 3D calibration object for camera calibration ?
With a 3D calibration object (e.g. a transparent cuboid like a table), I can take one image with a camera, then  find the corners of the cuboid, which gives me 8 2D-3D correspondences.
Since a projective camera matrix $P$ has 11 DoF, I have enough equations to obtain $P$.
Now with a checkboard given one image, I obtain the internal camera matrix and external kamera matrix.
However, the external camera matrix describes the change of basis from the checkerboard coordinate system to the camera coordinate system, and thus changes from image to image.
Thus is the idea of the checkerboard to take multiple images to obtain a more accurate estimation of the internal camera parameters, which should be same across all taken images ?
So for image i, we have 2D points $p_{i}$ with corresponding 3D points $P_{i}$ (on the checkerboard). Let $K_{i}$ be the internal camera matrix, $E_{i}$ be the external camera matrix and $M = \begin{pmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0  \\ 0 & 0 & 1 & 0\end{pmatrix}$  the  projective camera model.
So I am assuming that the optimization problem that is being solved is
$\min_{K,E_1,\ldots,E_m} \sum_{j=1}^{m}\sum_{i=1}^{n} || (K M E_{j}) P_{i}- p_{i,j}||^2$,
so that we obtain one internal camera matrix that fits to all images ?
Accordingly, am I right that it is necessary to do both, the calibration with the checkerboard to get an accurate estimate of $K$ and then the calibration with the 3D calibration object to obtain the extrinsic camera matrix, so that the 3D calibration object defines the world coordinate system ?

Now with a checkboard given one image, I obtain the internal camera matrix and external kamera matrix.

No, you don't. You get a single coordinate transformation matrix that contains only the ratio between focal length and object distance. You cannot split that into an internal camera matrix (focal length) and external matrix (object distance) without further information.
If you want that extra information from only a single image, then the points of your calibration target must not be all in a single plane. This is why people use cuboid calibration targets. Alternatively, you can take at least two images of checkerboards from different perspectives, which also allows you split focal length and object distance.

Accordingly, am I right that it is necessary to do both, the calibration with the checkerboard to get an accurate estimate of K and then the calibration with the 3D calibration object to obtain the extrinsic camera matrix, so that the 3D calibration object defines the world coordinate system?

No, this is not necessary. A single image of a cuboid calibration target is enough to determine both the internal and external camera matrix and a single image of a checkerboard gets you nothing.


Denoising short, non-uniformly spaced, bandlmited sequences Consider a bandlimited signal $x(t)$ with bandwidth $BW$.  Samples of this signal are observed at non-uniformly spaced sample points in the presence of noise:

$y(t_i)=x(t_i) + n_i \qquad i\in\{0,1,\cdots,N-1\}$

where $N$ is on the order of 10-30 samples.

The samples are "densely-spaced" in the sense that:

$\frac{1}{t_i - t_{i-1}} > 2\times BW$

The noise samples can be modeled as a realization of independent, identically distributed random variables of zero mean and some variance $\sigma^2_n$.  A variant of the noise model would be the inclusion of occasional outliers.

My goal is to recover estimates $\hat{x}(t_i)$ of the desired sequence which explicitly exploit the bandlimited nature of $x(t)$.

Savitzky-Golay filters can be generalized for non-uniformly spaced samples, but it's not clear how the bandwidth can be taken into account.  (For uniformly spaced samples, such an interpretation is available, see "What Is a Savitzky-Golay Filter?" by R. Schafer).  I guess I could use the frequency response interpretation of the SG filter for a worst case effective, constant sample frequency, but I'm wondering if other approaches may be better.

An extension would be interpolation/resampling of the data on a uniform grid.

Any suggestions?

Not a lot of references that include nonuniform, band limited, and filtering

D. Bonacci and B. Lacaze, "Lowpass/bandpass signal reconstruction and digital filtering from nonuniform samples," 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), South Brisbane, QLD, 2015, pp. 3626-3630.
doi: 10.1109/ICASSP.2015.7178647
Abstract: This paper considers the problem of non uniform sampling in the case of finite energy functions and random processes, not necessarily approaching to zero as time goes to infinity. The proposed method allows to perform exact signal reconstruction, spectral estimation or linear filtering directly from the non-uniform samples. The method can be applied to either lowpass, or bandpass signals.
keywords: {band-pass filters;digital filters;linear phase filters;low-pass filters;random processes;signal reconstruction;signal sampling;band-pass signal reconstruction;digital filter;finite energy function;linear filter;low-pass signal reconstruction;nonuniform sampling problem;random process;spectral estimation;Baseband;Estimation;Fourier series;Interpolation;Nonuniform sampling;Random processes;Time-frequency analysis;Nonuniform filtering;Periodic nonuniform sampling;Sampling theory;Signal reconstruction;Spectral estimation},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=7178647&isnumber=7177909

H. Johansson and P. Lowenborg, "Reconstruction of nonuniform sampled bandlimited signals using digital fractional filters," ISCAS 2001. The 2001 IEEE International Symposium on Circuits and Systems (Cat. No.01CH37196), Sydney, NSW, 2001, pp. 593-596 vol. 2.
doi: 10.1109/ISCAS.2001.921140
Abstract: This paper considers the problem of reconstructing nonuniformly sampled bandlimited signals using a synthesis system composed of digital fractional delay filters. The overall system can be viewed as a generalization of time-interleaved ADC systems. By generalizing these systems, it is possible to eliminate the errors that are introduced in practice due to time-skew errors
keywords: {analogue-digital conversion;bandlimited signals;digital filters;filtering theory;signal reconstruction;digital fractional filters;nonuniform sampled bandlimited signals;signal reconstruction;synthesis system;time-interleaved ADC systems;time-skew errors;Analog-digital conversion;Delay;Digital filters;Electronic mail;Filter bank;Frequency;Interleaved codes;Nonuniform sampling;Sampling methods;Signal synthesis},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=921140&isnumber=19923

Y. C. Eldar and A. V. Oppenheim, "Filter bank interpolation and reconstruction from generalized and recurrent nonuniform samples," 2000 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings (Cat. No.00CH37100), Istanbul, 2000, pp. 324-327 vol.1.
doi: 10.1109/ICASSP.2000.861961
Abstract: This paper introduces a filter bank interpretation of various sampling strategies, which leads to efficient interpolation and reconstruction methods. An identity, referred to as the interpolation identity, is used to obtain particularly efficient discrete-time systems for interpolation to uniform Nyquist samples, either for further processing in that form or for conversion to continuous-time. The interpolation identity also leads to a new class of sampling theorems including an extension of Papoulis' (1997) generalized sampling expansion
keywords: {channel bank filters;discrete time filters;filtering theory;interpolation;signal reconstruction;signal sampling;discrete-time systems;filter bank interpolation;generalized sampling expansion;interpolation identity;recurrent nonuniform samples;sampling theorems;signal reconstruction;uniform Nyquist samples;Data acquisition;Digital signal processing;Filter bank;Fourier transforms;Interpolation;Iterative algorithms;Laboratories;Sampling methods;Signal processing;Signal sampling},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=861961&isnumber=18684


Finding causal impulse response given the imaginary part of the frequency response I understand that I would need to calculate inverse Discrete Time Fourier Transform (iDTFT) to find $h(n)$. Since $h(n)$ is real, iDTFT of the imaginary part of $H(e^{j\omega})$ gives the odd part of h(n) and  iDTFT of the real part of $H(e^{j\omega})$ part would give the even part of $h(n)$ (Correct me if I am wrong here). But, how do I find the even part of $h(n)$ when I only have imaginary part of frequency response and find all the possible $h(n)$?

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


Question About $ Q $ Matrix (Model Process Covariance) in Kalman Filter I am implementing my own discrete Kalman filter to estimate velocity from acceleration and position measurements (using Matlab ).

I think I managed to deal with the $R$ matrix (measurements noise covariance matrix) in this way:

M = [x;x2dot];
R = cov(M);

What I am not sure about is the matrix $Q$.

In all the examples I found on the web and on this website also, the values inside that matrix are really really small. Moreover, I found here this statement:

  if you select an overly large Q, then it doesn’t seem like the Kalman
  filter would be well-behaved.

The fact is that, in my model, to obtain acceptable results I used values in the order of $1e6$ on the diagonal of the $Q$ matrix.

In this way I got the estimated displacement equal to the measured one and the estimated acceleration equal to the measured one. The velocity is in the right range of values but it still looks noisy.

So the question is: is there a limit for the values to be used in $Q$?

And also, are there guidelines to choose those values?

EDIT

Last question, does the $Q$ matrix have to be diagonal or it can be full also?

I'm the author of the textbook linked to above. This is a new account and thus I am not allowed to reply to that answer.

Anyway, the Gaussian chapter covers the definition of a covariance matrix. In general though your Q matrix will be full, not a diagonal, because there is correlation between the state variables. For example, suppose you are tracking position and velocity in x-axis for a car. Process noise (which Q is modelling) might be things like wind, bumps in the road, and so on. If there is a change in velocity due to this noise, there will also be a change in position. They are correlated, and so the off-diagonal elements will be non-zero for those terms.

The Kalman math chapter tells you how to compute the Q matrix for various scenarios. Admittedly there is more work to be done on that section. It's one of the more opaque subjects as most sources either gloss over it or just throw an integral at you and expect you to figure out the rest.

Finally, the companion FilterPy project contains code to compute the full Q matrix for you given common scenarios. The code is in common/discretization. For analytically difficult situations there is a routine to find the result numerically following an idea from Golub and Van Loan; the textbook is lacking on this topic, however.


Alternate Octave Commands for Matlab can you tell me how to write video data to file..
what is the alternate of these 4 commands in octave

videoreader('bg.mp4');
readframe(v);
videowriter('bg1.avi');
writevideo(newv,k);

I am not an Octave user, but the Video package appears to do these things.

Good luck!


Trouble with fixed-point arithmetic 2nd order IIR implementation in c DSP newcomer here!

I am tinkering around with a TI DSP and am trying to implement a second order IIR filter in C. Input is 16 bit 2's complement, as is the output, the accumulator is 32 bit wide.
I have tried several structures: Direct form I, Transposed direct form II, Gold & Rader structure. However, I can't seem to get anything to work and I am not sure whether it is a typo in the code of all of my filter implementations or whether it is because of me misunderstanding the principle of fixed-point arithmetic. I assume - or better: hope - that it's the same mistake that is affecting all implementations.

Here's my code for the transposed df II implementation:

// state variables
int w[3] = {0};

short filter_dfii(short value) {
    w[0] = ((value * sos_filter_coeffs_b[0] + w[1]) >> 15 ) & 0xffff;
    w[1] = value * sos_filter_coeffs_b[1]
       - w[0] * sos_filter_coeffs_a[1]
       + w[2];
    w[2] = value * sos_filter_coeffs_b[2]
       - w[0] * sos_filter_coeffs_a[2];

    return (short) w[0];
}

My df I implementation looks like this:

int last_value[2] = {0};

short filter_dfi(short value) {
    w[0] =
        ((value * sos_filter_coeffs_b[0]
        + last_value[0] * sos_filter_coeffs_b[1]
        + last_value[1] * sos_filter_coeffs_b[2]
        - w[1] * sos_filter_coeffs_a[1]
        - w[2] * sos_filter_coeffs_a[2])
        >> 15) & 0xffff;
    last_value[1] = last_value[0];
    last_value[0] = value;
    w[2] = w[1];
    w[1] = w[0];
    return (short) w[0];
}

The filter coefficients are:

int sos_filter_coefficients_a = {32767, -51150, 21015};
int sos_filter_coefficients_b = {658, 1316, 658};

How did I get those coefficients? I utilized the Matlab 'butter' function and specified a 2nd order Butterworth lowpass with a cutoff frequency of 0.05*fs or 2400 kHz in my specific case. Then I multiplied each coefficient with 32767 and rounded it to an integer value. If I call Matlab freqz function with those quantized coefficients, I get a bode plot with the expected result, so I believe quantization error is not the issue here.

However, the result is not a lowpass but just garbage. Unfortunately, I lack the proper testing equipment to describe the result scientifically, but it sounds like digital oscillation that is slightly modulated by the actual input signal and is present even when there's no input at all.

So, in my opinion, this can't be an unstable filter since the poles are well within the unit circle. It can't be overflow limit cycle either because then it shouldn't be happening when there's only low input level or no input at all.. right?

If I alter the coefficients to

int sos_filter_coefficients_a = {32767, 0, 0};
int sos_filter_coefficients_b = {32767, 0, 0};

input gets returned basically unprocessed as expected, so my issues most likely aren't caused by any external factors - there must be something wrong with my coefficients, with the function itself, or both.

Any help is much appreciated!

The " & 0xffff" operation is questionable. It works for positive numbers but for negative numbers it will mask out the sign bit as well, so you will turn everything into positive numbers.

For the DF1, it would be better to delay the shift and mask until you return the actual output and do the conversion to short. This way you can keep the state variable in full double precision.


Graph of $[u(\tau)-u(\tau-4)]\cdot[u(t-\tau)-u(t-\tau-4)]$ I'm trying to visualize what would the following signal be like
$$[u(\tau)-u(\tau-4)]\cdot[u(t-\tau)-u(t-\tau-4)]$$
Consider that $\tau$ is the independent variable here, $t$ is the shift variable.

Note: I'm doing this to calculate the continuous-time convolution

As this signal (function) exists inside an integration operator, and also assuming that there are no impulses residing at any discontinuities introduced by $u(t)$ step functions, then we can simply disregard the isolated discontinuity points and rely on the following:

$$[u(\tau)-u(\tau-4)][u(t-\tau)-u(t-4-\tau)] =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{for t < 0}}\\
&u(\tau)-u(\tau-t) &\scriptstyle{\text{for 0 < t < 4}}\\
&u(\tau-t+4)-u(\tau-4) &\scriptstyle{\text{for 4 <t < 8}}\\
&0 &\scriptstyle{\text{for t > 8}}\\
\end{align}
\end{cases}$$

Where the intervals are interpreted with the $t$ viewed as a parameter for an integrand argument function of a dummy variable $\tau$.

This convolution is reminiscent of convolving a rectangular pulse of duration 4 units with itself, whose output result should be a triangular pulse of duration 8 units, extending from $t=0$ to $t=8$ with an apex at $t=4$. Assuming that this reminiscence is correct, I would define the convolution result with this more apparent form as (defining the pulse as $p(t) = u(t) - u(t-4)$) :
$$y(t) = p(t)*p(t) =  \int_{-\infty}^{\infty} {p(\tau)p(t-\tau)d\tau} =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{for t < 0}}\\
&\int_{\tau=0}^{\tau=t}{1d\tau} = t &\scriptstyle{\text{for 0 < t < 4}}\\
&\int_{\tau=t-4}^{\tau=4}{1d\tau} = 8-t &\scriptstyle{\text{for 4 <t < 8}}\\
&0 &\scriptstyle{\text{for t > 8}}\\
\end{align}
\end{cases}$$

Where the limits of the integrals are deduced from the integrand arguments given at the first set of definition above: For example, for the interval $0<t<4$ the limits of integral are choosen based on the integrand $u(\tau) - u(\tau -t)$. Looking at this integrand one can see that for this to be nonzero (for a given fixed t) we need both $\tau > 0$, (from the first $u(\tau)$) and also $\tau < t$ from the second $u(\tau-t)$. So combining these two conditions requires that $0<\tau<t$ must be the limits of integration. The other cases are treated similarly.

On the other hand the operation is most easily understood and visualied by the help of an application of "flip & drag" technique to evaluate the intervals and to deduce the resulting integrand functions, instead of the way I presented as a piecewise definition form.


Blending artifacts in OpenCV image stitching I am using OpenCV to blend a set of pre-warped images. As input I have some 4-channel images (*.png or *.tif) from where I can extract a bgr image and an alpha mask with the region related to the image (white) and the background (black). Both image and mask are the inputs of the Blender module cv::detail::Blender::blend.

When I use feather (alpha) blending the result is ok, however, I would like to avoid ghosting effects. When I use multi-band, some artifacts are appearing on the edges of the images:

The problem is similar to the one raised here, and solved here. The thing is, if the solution is creating a binary mask (that I already extract from the alpha channel), it does not work for me. If I add padding to the ovelap between both images, it takes pixels from the background and messes up even more the result.

I guess that probably it has to do with the functions pyrUp and pyrDown, because maybe the blurring to create the Gaussian and Laplacian pyramids is applied to the whole image, and not only to the positive alpha region. In any case, I don't know how to fix the problem using these functions, and I cannot find another efficient solution.

When I use another implementation of multiresolution blending, it works perfectly, however, I am very interested in integrating the multi-band implementation of OpenCV. Any idea of how to fix this issue?

You need a edge detection algotithm witch will filter the critical edges and then create out of them an binary image. That image will give you an region to blur and the rest should tay the same. This should avoid the water stains or ghost efects.


Better understanding of downsampling (decimation) and upsampling (interpolation) Although some questions were asked about this topic, I have not seen any that answers all the basic questions, that is why I took the liberty to ask more about this.

I suggest to limit, (or at lest to start by), 1D signals. Of course if someone wants to extend they are welcome.

For example, in some lectures it is mentioned that

an ideal 1D interpolator is an upsampler (by the wanted factor, e.g. 2) in cascade with an ideal lowpass filter i.e. an sinc.
Decimation would be to opposite. A convolution with a filter (sinc?) then a downsampler (by the wanted factor e.g. 2).

In Matlab, given a sine wave at a given frequency for example, how would you proceed to downsample it to a lower frequency? And could you make the link with aliasing ?

I have started to do the following, "manually" coded:

    %Working with audio frequencies, for example,
    %to listen the results with audioread() function afterwards
Fs=8000;%sampling frequency
t=0:(1/Fs):1;
f=4000;
x=sin(2*pi*f*t);
x_dwnsamp=zeros(length(x));
%decimate x to see what happends wrt aliasing
decimation_step=2;
for i=1:decimation_step:length(x)
   x_dwnsamp(i)=x(i);
end
figure;
subplot(1,2,1);
plot(abs(fft(x)));
subplot(1,2,2);
plot(abs(fft(x_dwnsamp)));

This was an attempt to resample the signal x which is just at the Nyquist freq. (4KHz) fora sampling freq. of 8KHz.

There are few problems in the code and question.

In your question 'given a sine wave at a given frequency for example, how would you proceed to downsample it to a lower frequency' this statement is not right. We are not reducing the frequency of sinewave, we are reducing the sampling frequency with which we sampled it.

Second, the code just inserts zeros alternatively. It doesn't downsample (means remove those samples altogether). This means, it is still sampled with original sampling frequency (8kHz) but now with alternate zeros. And what is the effect of alternate zeros in signal? It creates a replica of spectrum.
Also, in the code, the signal frequency is exactly half of sampling frequency, which causes sampling to happen at zero crossings. You can plot the signal $x$ and see.

In the code, I changed $f$ to 3000, and can see 2 peaks at 3000Hz and 5000 ($fs -  f$)Hz. For the downsampled signal, apart from original signal, there is 2 peaks corresponding to replicated spectrum.


Find a stable transfer function $G(z)$ such that $|G(z)| = |H(z)|$ Consider the following causal IIR transfer function:
  $$ H(z) = \frac{2z^3 - 4z^2 + 9}{(z-3)(z^2+z+0.5)} $$

  Is $H(z)$ a stable function? If it is not stable, find a stable transfer function $G(z)$ such that $|G(z)| = |H(z)|$.

Now obviously the transfer function is not stable due to the pole at $z = 3$.
And I would have just thought that the solution was to flip this pole inside the unit circle, so $z = 1/3$.

And so  $G(z)$ should be $  = \frac{2z^3 - 4z^2 + 9}{(z-1/3)(z^2+z+0.5)}$ right?

The solutions however go $  G(z)= \frac{2z^3 - 4z^2 + 9}{(1-3z)(z^2+z+0.5)}$, which seems fine, except for the next step, where they go $G(z) = H(z)A(z)$ where $A(z) = \frac{1-3z}{z-3}$

That part makes no sense to me, since $$H(z)A(z) =  \frac{(1-3z)(2z^3 - 4z^2 + 9)}{(z-3)^2(z^2+z+0.5)} \neq G(z)$$

  Can someone please explain a) what am I doing wrong; and b) why we introduce $A(z)$?

Your idea of replacing the pole $p$ outside the unit circle (i.e., $|p|>1$) with a pole $1/p^*$ inside the unit circle is correct. However, as mentioned in this answer, you need to be careful with the scaling. In order to obtain the correct scaling, a factor $(z-p)$ in the denominator must be replaced by a factor $\pm |p|(z-1/p^*)$. In that way the magnitude of the transfer function will be unchanged on the unit circle, i.e., for $|z|=1$. This should also be added to the problem formulation:

  ... such that $|G(z)|=|H(z)|$ for $|z|=1$.

Clearly, the new stable transfer function $G(z)$ can be written as the original transfer function $H(z)$ cascaded with an allpass filter $A(z)$, since $|A(z)|=1$ for $|z|=1$:

$$G(z)=H(z)A(z)\quad\Longrightarrow\quad |G(z)|=|H(z)|,\quad |z|=1\tag{1}$$

For the given problem, the transfer function $A(z)$ is given by

$$A(z)=\frac{z-3}{1-3z}\tag{2}$$


How does DCT decorrelate images? I've read in multiple places that DCT decorrelates Toeplitz matrices and images usually have Toeplitz structure. Can you explain with an example how DCT decorrelates a Toeplitz matrix?
Example for DFT:
DFT decorrelates circular matrices. This is how I was able to understand that.
Suppose $X$ is a matrix whose correlation matrix is not diagonal. We want to find a transformation $Y=AX$ such that correlation matrix of $Y$ is diagonal.
$$\mathbb{E}[YY^T] = \mathbb{E}[AXX^TA^T] = A\mathbb{E}[XX^T]A^T$$
Let the Eigen Value Decomposition be
$$\mathbb{E}[XX^T]=U\Lambda U^T$$
Then
$$A=U^T \qquad\text{gives}\qquad\mathbb{E}[YY^T]=U^TU\Lambda U^TU = \Lambda$$ which is diagonal.
So given a matrix, its eigenvector matrix decorrelates it.
Consider a circular matrix
$$A =
\begin{bmatrix}
a & b & c \\
c & a & b \\
b & c & a \\
\end{bmatrix}$$
A $3 \times 3$ DFT matrix is given by
$$\begin{bmatrix}
1 & 1 & 1 \\
1 & w & w^2 \\
1 & w^2 & w \\
\end{bmatrix}\qquad\text{and}\qquad w^3 = 1 $$
We can easily see that all the columns (or rows since it is symmetric) of the above matrix are eigenvectors of the considered circular matrix $A$. Thus DFT decorrelates circular matrices.
Is it possible to show in a similar way that DCT decorrelates a Toeplitz matrix?
PS: The answers here, here and here didn't solve my doubt.

[EDIT] In 1991, Nasir Ahmed wrote: "How I Came Up with the Discrete Cosine Transform". Interesting to read, on how he was inspired by  Chebyshev polynomials, and on how he didn't get funding, for a tool at the heart of JPEG and MP3.
Natural images are not very stationary, but locally, their covariance is often modeled by a first- or second-order process. Being able to "concentrate" those matrices is one key to efficient compression (not the only one, however).
In different words: one aims at projecting, or representing, the pixels in another system of vectors (a basis, a frame), such that the information content of the picture is better represented, or sparser. Classically, one preserves energy (hence orthogonal or close to orthogonal transformations. Two main options:

compute or learn the best adaptive representation from each image. PRO: more efficient in compaction; CON: computationally expensive and storage lost because should be sent for decoding
use a known and fixed transform (or a set of known transforms). PRO: fast; CON: less optimal.

The DCT is option 2, with the performance of option 1: a fixed transform with close to adaptive performance, for most standard images. It worked, people still use JPEG and MP3, at the heart of which one finds the DCT.
The goal is thus to find the eigenvectors of such (covariance matrices, with less computational cost than matrix diagonalisation. The original DCT (DCT-II) from Discrete Cosine Transform (1974) aimed at finding basis vectors that could approximately diagonalize the first-order Markov Toeplitz  matrix:
$$ \left[\begin{array}{c}
1&\rho& \cdots &\rho^{N-1}\\
\rho& 1&\cdots &\rho^{N-2}\\
\vdots& \vdots&\ddots &\vdots\\
\rho^{N-1}& \rho^{N-2}&\ddots &1\end{array}\right] $$

The original motivation for defining the DCT was that its basis set
provided a good approximation to the eigenvectors of the class of
Toeplitz matrices that constitutes the autocovariance matrix of a
first-order stationary Markov process, with the result that it had a
better performance than the discrete Fourier transform (DFT) and some
other transforms [l], 3, 4 with respect to such kinds of
processes. In fact, as shown in 4, the DCT is asymptotically
equivalent to the Karhunen-Loeve transform (KLT) of a first-order
stationary Markov process as $\rho$ tends to 1, where $\rho$ is the
correlation coefficient.

Looking at the basis functions for $\rho=0.9$

the match is pretty good, and it turns out that the processing performance is close as well to that of the Karhunen-Loève transform.
In their 1995 paper Diagonalizing Properties of the Discrete Cosine Transforms, the authors provide an extension of the above property.

In this paper, we obtain the eight types of DCT’s as the complete
orthonormal set of eigenvectors generated by a general form of
matrices in the same way as the discrete Fourier transform (DFT) can
be obtained as the eigenvectors of an arbitrary circulant matrix.
These matrices can be decomposed as the sum of a symmetric Toeplitz
matrix plus a Hankel or close to Hankel matrix scaled by some constant
factors.

In 1991, Nasir Ahmed wrote: "How I Came Up with the Discrete Cosine Transform". Interesting to read, some excerpts:

What intrigued me was that the KLT was indeed the optimal transform on
the basis of the meansquare- error criterion and the first-order
Markov process model, and yet there was no efficient algorithm
available to compute it. As such, the focus of my research was to
determine whether it would be possible to come up with a good
approximation to the KLT that could be computed efficiently.

to study a “cosine transform” using Chebyshev polynomials

The motivation for looking into such “cosine functions”
was that they closely resembled KLT basis
functions for a range of values of the correlation coefficient
$\rho$ (in the covariance matrix).


references for techniques to reduce recursive filter startup transients I am trying to track down references (published articles or reports ideally) that establish the basis for some of the techniques being used to minimize the effect of recursive digital filter startup transients on the sampled signal. My specific application is the digital filtering of acceleration signals that have already been acquired (post-processing). The examples below are based on Python and the scipy.signal module.

An example of the techniques I am referring to is illustrated in scipy.signal.filtfilt (and presumably the Matlab version). In this function, there are options to specify a padding type and padding length presumably to get the filter initialized and dissipate the transient before running into the "real" data. Another technique and a reference to F. Gustaffson, “Determining the initial states in forward-backward filtering”, Transactions on Signal Processing, Vol. 46, pp. 988-992, 1996, is also implemented in the filter.

Another example is scipy.signal.lfilter_zi that describes a technique in the notes to set the initial filter delays from the initial value of the signal. (Interestingly, this is also implemented within the filtfilt function referenced above, in addition to the padding.) The notes describe the derivation of the technique, but there are no references in the source code for this function. There is also no discussion of the pros and cons of this technique. I have found a reference that seems to describe this technique by Fletcher (1973). Can someone confirm that this is the correct, general reference?

Another source I have found for a more general technique is a report out of MIT by Chornoboy (1990). I have not seen this technique implemented in any libraries and I have not found a comparison of this technique with others, other than what is in the report. This report mentions the Fletcher report.

Finally, I have perused a number of books on digital signal processing. The book I own ("Digital Filters" by Hamming) specifically states that a general-purpose technique that is independent of the input signal has not been found. Other books do not mention the filter startup transient or, if it is mentioned, practical methods to minimize the effect are not discussed.

Is anyone aware of a single resource that discusses each of these techniques in the context of post-processing (filtering) digital signals? Short of that, are there any references for the (obvious) technique of padding the signal and discussing the various techniques of padding (i.e. odd, even, constant)?

For a given fequency response magnitude, a minimum phase filter implementation will have the shortest possible transient duration.

A minimum phase filter will have most of its energy packed at the beginning. Hence it's also called as a maximum energy filter. So it will produce steady state outputs as soon as possible (among the set of filters having the same frequency response magnitude)

However, a minimum phase filter cannot be linear phase simultaneously. Because linear phase filter would have a symmetric impulse response but that would also imply transfer function zeros outside of unit circle. Consequently it cannot be a minimum phase filter at the same time as for a minimum phase filter all the poles and zeros of transfer function are inside the unit circle.

Alternatively, if you have some reliable estimation of the initial samples, before $n=0$,  of the input signal, then you can initialize the filter according to those values instead of a full blind zero initialization.


Keypoint detection with color images Are there any known methods to use existing keypoint detectors (which generally work on grayscale images) on color images? of course, one could always add the results of the detector on all three channels to get $3N$ keypoints, but I was wondering if there exist other established methods.

Specifically, I am looking to use a FAST like detector.

Note, that I am not looking for a color based descriptor, but rather for a color based detector.

Short answer is no, there are none, or maybe none that made it until becoming a de facto standard like FAST or SIFT.
This is still an open area of research by the way if you plan to work in this field.

Mostly, the issues are linked to the fragility of color representation in digital images. The common RGB color space is not really relevant (distances between colors in RGB do not match the perceptual intuition of color distance), so you probably have to go Lab or YCbCr first, then work on the luminance channel. Also, be aware of the white balance that can vary a lot between different cameras, or with a given camera but varying light conditions (sunlight, artificial light...).

Furthermore, in most cases the detection and/or description have good reasons to be color agnostic. For example, when you watch a black-and-white movie, your brain can still recognize objects although they have lost their color.


What is the intuition behind convolution? I have been using convolution for finding outputs of various systems .I know how to use it.But I still don't know what does convolution exactly  means? how one can define convolution ?

No doubt there are many links and expositions on wikipedia and elsewhere about convolution, its relation to linear systems, and its link to frequency analysis, however I do not think that is what you are after.

Forget convolution for a moment. Just think of nothing but averages:

If I give you a series of say, $10$ numbers, and ask you to take the average, you can do that very easily.

Let us say that I give you another set of $10$ numbers again, and also ask you to take the average. You can again, do that very easily.

I keep giving you a set of $10$ numbers, and you keep giving me averages. What a life!

Now imagine I tell you a secret: Those sets of $10$ numbers I have been giving you, are actually from a taken from a long series of numbers. In fact I have a list of $100$ numbers. The first set of $10$ numbers I gave you were the first $10$ numbers of that list. (First sample to 10th sample). You averaged them.

The second list of $10$ numbers were actually the 2nd sample, to the 11th sample. You averaged them.

The third list of $10$ numbers I gave you were actually the 3rd sample, to the 12th sample. You averaged them.

After all this, I ask you to plot all the averages you gave me, one after the other.

Congratulations! You just showed me a convolution result!

(Disclaimer: There are some border conditions that I have skipped over, but this is basically all a convolution is).

So in summary, think of a convolution as just a sliding average.

Now, remember how I said "average" those $10$ numbers I gave you? Well you can also do a weighted average. (That is, instead of multiplying all the $10$ numbers by $1$ before averaging them, we multiply them by different numbers, before averaging them. This is a weighted average). This will still be a convolution though. We call those numbers "weights", and sometimes they are referred to as the "filter kernel".

Some advanced systems actually change their weights as they run, and they are called adaptive filters.

Eitherway, this should provide you with a solid intuition: Just remember that convolution can be thought of - is - a moving average of your data. Hope that helped!


Type of filter, given difference equation? Given a difference equation, how do we tell if it is an IIR filter or FIR filter?
For example,

y(n) = x(n-3)+y(n-1) . Is it FIR or IIR? Can you please give me a way to figure this out?
Thanks! :)

The difference equation of IIR filter is recursive, which means the current output depends upon the previous output (the output depends upon the infinite past). Your example is IIR.

On the other hand, the difference equation of FIR filter depends only upon the input. If the input is finite in duration, then the output is also finite in duration.


how fast can you sweep sinewave I am asking specifically about discrete sound signal, non dithered 16 bit 48000 Hz pcm/wav format.

My question is, how fast can you change sinewave frequency before it distort? Imagine a frequency sweep,it starts at 3000 Hz and ends at 6000 Hz, 3000 Hz full cycle sinewave at 48k sampling frequency is 16 samples long, 6000 Hz is half that, so 8 samples long...

If we do the sweep from 3k to 6K long duration,lets say one second,it will cleanly sweep from one point to another,but what is the precise limit,the threshold of how fast/short we can do this sweep before it becomes distorted?

Obviously we know its at least as long as single cycle of the lowest frequency we want to sweep,so its 16 minimum becose remember 16 samples is minimum lenght of single 3000hz cycle,but is it precisely 16 or is it longer? maybe one full cycle of starting frequency sample length(16) + same for final highest frequency(6000hz = 8) so 16 + 8?

I believe its also called "sweep rate",and sinewave frequency sweep is called "chirp",I want to know precisely the mathematics behind this to figure out where exactly this threshold is.

Also does this change depending on loudness/amplitude level? if we do the 3000 Hz to 6000 Hz sweep -6 db instead of maximum amplitude,would it make difference? considering signed 16 bit is 32768 amplitude levels,-6 db would be 16384.

I want some kind o mathematical formula that would precisely tell me exactly how many samples do I need to make shortest sweep possible at specific loudness and frequency points.So for example 3000 -> 6000 full amplitude sweep? atleast XXX samples needed! 375hz to 12000hz -12db? that be XXXXX samples minimum,you get the idea...

remember,we are talking about sweeping sinewave,not waveforms like squarewave that contain harmonics

edit
Definition of distortion,when the sinewave stops being sinewave,its simple as that.For example,if you module amplitude of sinewave slowly and only one way (up or down ) then it will always be just pure clean non-distorted sinewave,but if you start changing its amplitude too fast,so its going up and then changes direction downwards per single cycle of that sine,harmonic distortion will grow from that sine,so you will have multiple tones at same time for example 1000hz sine will get 2000hz,3000hz,4000hz... overtones/partials/harmonics all from amplitude changes

same things happens with freqency change instead of amplitude change,you do it too fast,and new sinewaves(hamonics) will appear at same moment on top of the original moduleted sine.

To ask "how fast can you sweep a sinewave?" is to ask "How fast can you hit $\frac{F_s}{2}$?".

Because, sweeping a sinewave is equivalent to modulating it.

A sine wave looks like:
$$x(n) = \sin(n 2 \pi f / F_s)$$

Let's focus on that instantaneous phase. I have put $n$ first, on purpose, to show that phase, in discrete time, is a line that departs from $0$ at a slope of $\frac{2 \pi f}{F_s}$. Even more descriptively, $2 \pi \frac{f}{F_s}$.

The sharpest slope for a given $F_s$ that can exist is at $\frac{2 \pi \frac{F_s}{2}}{F_s} = \frac{2 \pi F_s}{2 F_s} = \pi$.

Anything even slightly higher than that rate will start "catching" the sinusoid at its "past", or in other words, create aliasing.

So, how fast can you sweep a sinusoid?

Consider a cartesian plane. It's $x$ axis corresponds to time and it's $y$ axis corresponds to rate of phase change. We mark a red line on the $y$ axis at $\pi$. That's our limit. We can't cross that. The time axis changes in terms of $T_s = \frac{1}{F_s}$. So, theoretically, we could go $0 \ldots \frac{F_s}{2}$ in just one step, in just one $T_s$.

Is that going to look like a sweep?

It is going to look like a sweep, as much as a sinusoid looks like a sinusoid when it is sampled at some $\frac{F_s}{2}$. That's two samples per period and it looks exactly like a triangle and nothing like a sinusoid, but when you are sampling at $F_s$, that's the best sinusoid you can "see" at $\frac{F_s}{2}$.

This cartesian plane representation though, helps convey another point too. That is, how many frequencies did the sweep end up going through anyway? In the extreme example given above, instantaneous phase was so rapidly modified that effectively just two frequencies showed up, 0 and $\frac{F_s}{2}$. Sometimes, when you are trying to obtain the frequency response of a system, you might "miss" a part of it because your sweep is so fast that it "jumps" over the mode you should have excited.

Therefore, a more meaningful way to determine your sweep rate is to accept an average power per Hz (or frequency band) that you want your chirp to deliver and then find the rate of phase change that "allows enough samples" to deliver that power.

Hope this helps.


Understanding my Butterworth filter My input is x and my filtered output is y. The sampling rate of x is Fs. I create a low-pass Butterworth filter in MATLAB as follows.

[b,a]=butter(order,cutoff*(2/Fs));
y = filter(b,a,x);

Letting order vary (cutoff fixed at 20Hz) yields these outputs y.

Letting cutoff vary (order fixed at 6) yields these outputs y.

I have a few very rookie questions about this sort of filter.

Why does the filter constrain y to equal zero at the first sample?
Why is there an initial ripple in y?
Why does the time-delay of y increase with order and cutoff?
The time step between successive samples of x is 1/Fs. Is the time-step between successive samples of y also 1/Fs?

FYI, my goal in filtering x is to remove high frequency noise above 20 Hz or so.

The filter doesn't constrain the first sample to be zero. The first sample just has a relatively small value, because it takes the filter some time to react to the input signal. This is because initially the delay elements of the filter contain zeros (i.e. the initial state is zero). This means that for the filter it looks like there has been a long input signal consisting only of zeros, and suddenly the input jumps to a value of around $-20$. This is also why you see the ripples in the output signal. What you see is basically the (negative) step response of the filter. And any frequency selective filter with a decent (frequency domain) performance will have some overshoot in its step response.

Note that the function filter() accepts an initial state as an input argument. So you could compute the filter state for a constant input signal (say about 100 samples or more) equal to the first value of the actual input signal. The filter function will give you the final state after having processed that constant input signal, which you can use as initial state for filtering the actual data. This will greatly reduce any transients in the filter output:

[dum,s] = filter(b,a,x(1)*ones(100,1));    % get initial state
y = filter(b,a,x,s);                       % filter with initial state s

The group delay of the filter is frequency dependent, so it's hard to make general statements about the delay of a specific input signal. What is true is that the average group delay increases with increasing filter order.


Getting started with wavelet-based real-time anomaly detection for temporal signals using Python I have an anomaly detection problem I wanted to consider addressing in software (python) with wavelets.   I’m new to wavelets as well as the python’s spicy.signal wavelet package, let alone with the numerous additional packages floating around out there (e.g., PyWavelets).

Here’s the problem formulation:  Consider a desired low-pass signal subject to various anomalies and a small amount of additive white Gaussian noise (AWGN):

$$
r(t) = d(t)x(t) + n(t) + a(t)
$$

where $x(t)$ is the desired low-pass signal, $n(t)$ is the AWGN, and $a(t)$ are the additive anomalies and $d(t)$ are the multiplicative anomalies.  The anomalies are roughly modeled as a linear combination of scaled and shifted versions of two types of effects: impulsive spikes, exponentially decaying transients.  Multiplicative square pulse-type signal drop-outs directly affect the desired signal:

$$
a(t) = s(t) + e(t).
$$

The spikes can be written as:
$$
s(t) = \sum^{N_s-1}_{k=0} s_k (t)\\
s_k(t)  =  S_k \delta(t-t_k)
$$
where $N_s$ is the number of spikes, $\delta(t)$  is the Dirac delta function, and $S_k$ and $t_k$ are the amplitude and delay associated with the $k$th spike

The exponential transients are written as:
$$
 \begin{align*}
 e(t) & = \sum^{N_e-1}_{k=0}e_k(t)\qquad
 e_k(t)  = \left \{  \begin{array}{lll}
         E_k \cdot \exp({-[t-t_k]/\tau_k} )&& t\ge t_k\\
         0 && t< t_k
        \end{array}
       \right .
 \end{align*}
$$
where it is assumed that the transients are well-separated in time and $\tau_k$ are known not to exceed some maximum.

Dropouts are crudely modeled as product of scaled and shifted square pulses:
    \begin{align*}
 d(t) & =\prod^{N_d-1}_{k=0}d_k(t)\\
 d_k(t) & = \left \{  \begin{array}{lll}
               \epsilon_k && t\in[t_k, t_k+T_{d_k}]\\
               1 &&  t\notin[t_k, t_k+T_{d_k}]
        \end{array}
    \right . \qquad \epsilon_k\in[0,1)
 \end{align*}
where it is assumed that the drop-outs do not overlap in time and that $T_{d_k}$ are less than some maximum.

The goal is to flag the time intervals affected by the anomalies which are further assumed to arise one at a time (e.g., a spike will not occur in the middle of a dropout).  The exponential transients are considered to be non-negligible for a limited period of time (say up to 90% attenuation from the start of the transient).

I don’t want to rely heavily on a precise parametric model for detecting the anomalies since the above anomaly modeling is only a crude approximation.  It seems like wavelets may be a good method for localizing these effects in the received signal.  The desired signal’s bandwidth is far lower than that of the anomalies which overlap and extend beyond the desired signal bandwidth.  It is necessary for the detection procedure to be at least near real-time.

I don’t know much about wavelets aside from a few very basic tutorials and I haven’t played with any python-based wavelet packages.

My questions are as follows:

1) Are wavelets indeed well-suited to this type of detection problem?  If not, what other approaches would make sense to try?

2) If wavelets do make sense here, what type of wavelet transform should be considered (e.g., continuous or discrete)?

3) What type of mother wavelets should be tried here?

4) Are the wavelet routines in scipy.signal adequate for a preliminary, meaningful study?  If so, which are the relevant functions?  If not, which package would be preferable?  PyWavelets?

5) Can anyone point me to resources that have addressed directly similar problems?

I suggest that you look at the Wikipedia article:

https://en.wikipedia.org/wiki/Change_detection

and first focus on the statistical aspects of the problem.  The references cited reflect most of the foundational material of the topic.  The article is of high quality. You should also think about the systems aspect of your problem like, the architecture of the system in the sense of issues like having a set of essentially special purpose detectors, or something that allows some interactions.  Your problem statement says that only a single anomaly is present at a time, and this should be reflected in the organization of system.  Questions such as Do you need automatic gain control (AGC) or do you need a constant false alarm (CFAR) behavior do have  indirect effects.

The major departure from the usual change detection problem, is $a(t)$. Typically, it is in $a(t)$ that one wishes to find changes.  I don't think that this is necessarily a problem but $a(t)$ will propagate through your detectors.

Change detectors work on changes from something typical and $a(t)$ seems like a good place to look for "typical". This is probably where your likelihood ratio test or score metric is going to reside.

I'm guessing that your essential problem is that $a(t)$ is from a remote sensor, and your anomalies are associated with a transmission channel.  If so, transmission glitches are not typically subtle low SNR events that necessitate near optimal techniques needed to find them. They are obvious and are often managed by simple heuristics. Much that follows will be irrelevant if this is the case.

Finally, Wavelets. The answer is, you are going to determine this for yourself, by trial and error and thinking through what makes sense. One heuristic is that time-bandwith matching has a small loss relative to something like a matched filter.

The choice of basis functions is specific to your application.  People use Wavelets. People use STFT, People use Wigner.  People use Cyclostationarity. People use running histograms.

It all depends on how parsimonious a basis set represents your signals of interest, but you can't neglect $a(t)$. Those basis function have a direct bearing on characterizing "typical", which is why I strongly recommend you think about the statistical organization of your system, and then work on selecting your basis functions.


Relation between frequency and phase in PLL? Apparently both frequency and phase are independent properties but what is the relation between frequency and phase in PLL?so that when one is locked ,other is also locked

no, frequency and phase are not independent quantities.  phase is the integral of frequency.  and in the VCO (or NCO) inside a PLL, there is an integrator.  the frequency is proportional to the input of the VCO or NCO, but the phase (which is what goes into the phase discriminator) is the integral of that.

with feedback, a PLL is a servo-control mechanism.  but because of the inherent integrator in the loop, the "controller" is either a PI or an I.  not much of a D in that controller inside a PLL.


Sampling from a distribution with arithmetic coding I'm struggling to implement a sampling scheme where arithmetic decoding is used to generate the samples.

This is based on excercise 6.3 from David MacKay's book, where you compare rejection sampling with arithmetic decoding sampling. For each symbol you want to generate with rejection sampling you have to generate an (uniformly distributed) integer from $1$ to $2^{32}$, rescale it to be in $\left(0,1\right)$ and then if your random number is less than 0.99 emit a 0 (or 1 accordingly).

To generate symbols by arithmetic coding, you feed standard random bits into an arithmetic decoder equipped with the appropriate model.

(Note: I think this works because when we compress data we go from symbols generated under P to symbols under Q - P,Q are probability distributions - and decompression works the opposite way. If we choose Q to be an easy distribution to sample from then by decoding symbols from Q we can generate symbols from P.)

I've implemented the methods in the code below:

clear all;
pdist = [0.99 0.01] ;
counts = [99 1];
N=100;
U = randi(2^32, [1 N]);
U = U/2^32;
rejection_samples(1:N) = 0;

for i=1:N
    if U(i) < pdist(1) ;
       rejection_samples(i) = 0;
    else
       rejection_samples(i) = 1;
    end
end

seq = randsrc(1, N, [0 1;0.5 0.5]);
arith_samples = arithdeco(seq, counts, N);

I'm confused becuase the arith_samples variable sometimes contains a 2 (but I guess this just maps onto (0,1) as the output never seems to contain a 0 - just 1 and 2).

The answers say that the rejection sampling method takes 32 bits per generated symbol, but that arithmetic coding only takes 0.081 = $H_2\left(0.01\right)$. I can't really see that from my code.

Is there some way I could measure 'bits used'?

Also this seems a bit too good to be true - surely there are better ways of generating random numbers than just rejection sampling?

The relevant chapter is found here:

http://www.cs.toronto.edu/~mackay/itprnn/ps/105.124.pdf

While this answer may have come a bit too late for OP, here it is.

Let's start by replacing the phrase 'bits used' with 'queries to an unbiased bit generator'. Then assume that the only random number generator you have at your disposal is an unbiased bit generator, that generates $0$ or $1$ with probabilities $p_0=p_1=0.5$. e.g. something like the following (well, this is probably not the best way to implement an unbiased bit generator in octave/matlab, but suffices for our purpose):

function b = unbiased_bit_generator() b=randi([0, 1], 1,1); end;

Now your goal is to use this function as less as possible (1 use = 1 query) to generate a given number of iid biased bits with a specified bias.

In case of rejection sampling, you would use $\texttt{unbiased_bit_generator()}$ 32 times to get a 32-bit integer $i$, then check if $\frac{i}{2^{32}}<p_0$ to generate one biased bit. Thus you require $32$ queries/biased-bit. Here is a biased bit generator based on 32-bit rejection sampling:

function b=biased_bit_generator1(p0)
    for i=1:32
        bin_array(i)=unbiased_bit_generator();
    end;
    str=num2str(bin_array);
    str=str(str~=' ');
    if bin2dec(str)/2^32<p0
        b=0;
    else
        b=1;
    end;
end;

In case of arithmetic decoder based sampling, you feed the decoder with a sequence $s$ of $k$ unbiased bits (generated using $\texttt{unbiased_bit_generator()}$ $k$ times) and get a sequence $d$ of of length $n$ biased bits. Thus you needed $\frac{k}{n}$ queries/biased-bit.  Because arithmetic coding is a variable-length coding technique, two encoded sequences of same length may correspond to two message sequences of different length. Thus, if you were to end up with  sequence $s_k'\neq s_k$ after using the unbiased bit generator $k$ times and feed it to the arithmetic decoder, then you may get as output a sequence $c_{n'}$ s.t. $n\neq n'$.  In this case, you needed $\frac{k}{n'}$ queries/biased-bit.

The number 0.081 specified in the book is the expected number of queries/biased-bit where the expectation is over all encoded input sequences to the decoder. As you keep feeding your decoder with more and more unbiased bits and keep generating more and more biased bits the queries/biased-bit will approach this number.

  Also this seems a bit too good to be true

I assume you are talking about 0.081 queries/biased-bit. I think this statement is equivalent to saying that Shannon's source coding theorem is too good to be true ;)


Compression of Multivariate Phase Signals I have a data matrix $X$ whose columns are composed of instantaneous phases of recordings from different electrodes. In other words, each column keeps a periodic sawtooth signal with a range from $[-\pi,\pi]$. I am interested in compressing this dataset by PCA, but the application of PCA on phase signals sounds counter-intuitive to me. Is there any principled way of achieving this?

Your signal is easily described by two parameters (the initial value and the frequency), the sawtooth may be linearized, by stitching each linear segment to the end of the other, and then you can use a linear regression.
The function that stitch the linear segments of a sawtooth of amplitude $2 \pi$ is the unwrap, there is a function with the same name in octave/matlab if you want to use them.
Here is a code snipped that you could use to get the frequency of your columns in python
import numpy as np
omega = np.arange(2, 8).reshape(1, -1) # a row vector
t = np.linspace(0, 10, 1001).reshape(-1, 1) # a column vector
# this is a matrix as you described
matrix = np.angle(np.exp(1j * omega * t))
linearized = np.unwrap(matrix, axis=0);
# this function will get the frequency of your underlying signal
omega_ = np.mean(np.diff(linearized, axis=0), axis=0)

Then for the phase you could use
phase = np.mean(linearized - omega_ * t)


audio signal time delay estimation(TDE) in indiviual and separate mics common TDE problems algorithms(such as GCC) often require audio signals in different mics concurrently, however, is it possible to estimate audio signal time of arrivals(ToA) in different mics separately in the noisy environment in order to estimate TDoA&TDE?

Also, since memory space limitation, cross-correlation related algorithms may not be feasible, starting point detection may be a possible solution(in noisy environment)

single mic, i'm not sure.  it seems to me that all you can do is mark different events (perhaps with envelope following or perhaps by comparing the signal or parameters measured from the signal to templates) and just use the time difference of the onset times of each marked event.

for two mics, there is cross-correlation.  but some people have also suggested an adaptive LMS filter (one of the mics is the "desired signal", $d[n]$) of some sort and computing some measure of group delay from the FIR coefficients that result in the LMS.  something in me makes me think is the Julius Smith/Perry Cook thing, but it might be someone else.


How to apply a Butterworth filter to data of varying sample rate? I am trying to apply a Butterworth bandpass filter to accelerometer data of my smartphone. However, the accelerometer samples I receive do not come at regular intervals. Sample frequency varies between 1 Hz and 128 Hz in as little as a second.

I read the source code of a few Java implementations of Butterworth, but all of them appeared to assume a constant sample rate. Is there a way to implement Butterworth that allows for a variable sample rate? Or should I interpolate the accelerometer data to obtain a constant sample rate?

If you answer with math, please include a layman's explanation. I'm a bit weak on the theoretical side of signal processing.

I can figure out two options:

Interpolate data
Remember that filters are dynamical systems and solve the corresponding ordinary differential equation (ODE)

In the second method, you basically use the retrieved time step as the time step in the ODE solver. Assuming you have a second order system:

$$\ddot{y}+\frac{\omega_0}{Q}\dot{y}+y\omega_0^2=x\omega_0^2$$

For butterworth, $Q=\frac{1}{\sqrt{2}}$ Introduce two variables $u=y$ and $v=\dot{u}$, and you get two coupled first order ODE:

$$\left\{\begin{aligned}
     \dot{u}&=v\\
     \dot{v}&=x\omega_0^2-u\omega_0^2-\frac{\omega_0}{Q}v
    \end{aligned}\right.$$

Assuming you do not want to have a too high cutoff frequency

$$\left\{\begin{aligned}
   u&\leftarrow u+v\Delta t\\
   v&\leftarrow v+ \left(x\omega_0^2-u\omega_0^2-\frac{\omega_0}{Q}v\right)\Delta t\end{aligned}\right.$$

Where $\Delta t$ is the time difference between two steps. If you want $\omega_0 > \frac{0.25\text{ Hz}}{2\pi}$, you need to consider the Euler backward method, or if you need more accuracy, the trapezoid method.

Addendum

You may not want the overshoots and undershoots of the Butterworth filter. Then, set $Q$ a bit lower. Probably, you want a critically damped system, which do not have any such effects. In that case, set $Q=\frac{1}{2}$. An interesting side-effect of choosing the ODE approach is that you get a filtered derivative for free (it is the $v$ component), if you were to trig on changes in acceleration.


Is the Laplace transform a special case of Fourier transform? (Not the other way around) Always had a thought about why Laplace transform reveals the transient properties of the system?
My doubt is based on the following fact,
Fourier transform is given as

\begin{equation}
\mathscr{F}\left\lbrace f(t)\right\rbrace = \int_{-\infty}^\infty f(t) e^{ -j \omega t}  dt
\end{equation}

Where Mathematically and intuitively we believe that the angular frequency $\omega$
takes only real value.

What if, instead of taking real angular frequencies, if the variable
$\omega$
assumes a complex angular frequency in the form
$\beta - j \alpha$
, then,

$$
j \omega t = j (\beta - j  \alpha) t = (\alpha + j \beta ) t = s t
$$

While taking Fourier transform w.r.t $\omega$, the quantity  $\beta$ will be real angular frequency in radians per second and $\alpha$ will be the $\textbf{imaginary angular }$ frequency in radians per second.

\begin{equation}
\int_{-\infty}^\infty f(t) e^{ -j \omega t}  dt = \int_{-\infty}^\infty f(t) e^{ - s t}  dt = \mathscr{L}\left\lbrace f(t)\right\rbrace
\end{equation}

Hence is it mathematically correct to consider bilateral Laplace transform as a special case of Fourier transform (not the other way around) when $\omega$ takes a complex angular form $\beta - j  \alpha$ ?  I believe the fact that $\omega$ can take complex values is the reason why we get transient properties of the system when using Laplace transform.

The Fourier Transform is the Laplace Transform with the complex variable s restricted to be the imaginary axis on the s plane. For this reason the Fourier Transform only exists when the imaginary axis is within the region of convergence. The variable s is called a "complex frequency" as it is the frequency variable that can take on real ($\sigma$) and imaginary ($\omega$) components. That said, I would view the Fourier Transform as a subset of the Laplace Transform, or the Laplace Transform as an expansion on the Fourier Transform that provides a lot more functionality and can exist when the Fourier Transform can't.

This is also the reason that the frequency response for a system with a general transfer function $H(s)$ is given as $H(j\omega)$.

When a system is restricted to $s= j\omega$ as the input, then the input is restricted to be only sinusoids or signals given by $e^{st}$ with $s = j\omega$ which maintain a constant magnitude with time. By allowing s to have real and imaginary components as in $s = \sigma + j\omega$ then we also allow the input to grow or shrink with time, depending on which point in the s plane is used as the input to the system.


How to efficiently calculate points of intersection of a straight line and a contour? Ex.
I have an image, let's say 100x100 pixels with some shape on it. And a set of straight lines that are passing through origin located at some point of image, for example origin located at point (50, 50). Number of lines determined by step - angle. If step = 15 degrees, then we have 180 / 15 lines. I want to determine points of interesection in efficient way, but all that I came up for now is printing a line on matrix and check for non-zero elements if they overlap non-zero elements of a shape(and check for diagonal connectivity).

Calculate the straight lines angles with respect to the x-axis

Calculate the angles of each and every point on the contour

Search for points on the contour with same angles as the lines with respect to the origin.

These are the points of line-contour crossing.


Plane Settings of the Matched $z$-transform Method I've come across that the matched $z$-transform maps poles of the $s$-plane design to locations in the $z$-plane. My question is, what is the $s$-plane and what does this mean? I'm aware that the variable $s$ can often be used to denote the complex plane but the fact that this has been used in conjunction with the $z$ variable has confused me slightly.

The $s$-plane is the complex plane associated with the Laplace transform, i.e., with transfer functions of continuous-time systems, whereas the $z$-plane is the complex plane associated with the $\mathcal{Z}$-transform, i.e., with transfer functions of discrete-time systems. The matched $Z$-transform is one way of transforming continuous-time systems to discrete-time systems.

Take a look at this answer for more information about transformations from continuous-time to discrete-time.


periodicity coefficient I wonder if an efficient method exist to compute how much a signal is periodic, it should be ~1.0 when the signal is totally periodic (like a sinusoïdal signal) and ~0.0 when totally random, like a white noise.

Edit : This question came with a big misunderstanding, it has been modified to fit the answer.

I suggest using Spectral Flatness, aka Wiener Entropy. It is defined as a ratio of geometric and arithmetic mean of the magnitude spectra $X(k)$:
$$\Xi=\dfrac{\sqrt[K]{\prod_{k=0}^{K} X(k)}}{\frac{1}{K}\sum_{k=0}^{K}X(k)} $$
For signals which have flat spectra, its value tends towards $1$, whereas for tonal signals it is close to $0$. In your particular application, you might want to consider $1-\Xi$ as a measure of "tonalness" in your signal.

Equivalently, in order to avoid underflows, one can replace the geometric mean implementation with $$GM=\exp\left( \frac{1}{K} \sum_{k=0}^{K}\ln X(k) \right) $$


Intuition for how a Fourier transform handles signals where the frequency varies with time? I created a discrete signal sampled at 44.1kHz that is 330Hz for 1 second, 360Hz for 0.5s, and 390Hz for 0.25s.  I did an FFT transformation of this, and got the magnitude graph below:

I don't understand how this magnitude information (combined with the phase information not shown) can reconstitute the original signal.  I know it does -- I'm not asking for a proof.  I'm hoping someone can supply the intuition for how that strange jagged FFT image ends up reconstructing the original signal.

Each frequency burst is equivalent to a rectangular window on an infinite sinusoid.  A rectangular window in the time domain is the same as circular convolution with a Sinc function in the frequency domain.  Thus, you end up with a spike with side humps (the spike for the sinusoid, and the humps due to the convolution with the transform of the rectangular window) for each frequency.

Then, since the FFT operation is linear, the 3 sine bursts or, equivalently, the 3 spiky hump groups, just sum additively in both the time and the frequency time domain, respectively.

You are correct in that the phase information is important for reconstruction.  The FFT phase results includes the information about the location of each (non-overlapping in your case) rectangular window for each sine burst.  Twist (or untwist) each spiky hump group's phase in the frequency domain, and you will move the burst's location (the rectangular window) in the time domain.


Application of the time-shifting property in case of Fourier-Transform of cosine Time-shifting property: $x[n-n_d] \xrightarrow{\mathscr{F}} e^{-j\omega n_d} X(e^{j\omega}) $
Fourier-Transform of cosine-signal: $\cos(\omega_0n) \xrightarrow{\mathscr{F}} \frac{1}{2}(\delta(\omega - \omega_0) + \delta(\omega+\omega_0)) $

Combining 1. & 2. together, I am getting:
$\cos(\omega_0n - \frac{\pi}{2}) \xrightarrow{\mathscr{F}} e^{-j\frac{\pi}{2}}\frac{1}{2}(\delta(\omega - \omega_0) + \delta(\omega+\omega_0)) $, but instead the Fourier-Transform of

$$\cos(\omega_0n - \tfrac{\pi}{2})$$

is

$$\cos(\omega_0n - \tfrac{\pi}{2}) \xrightarrow{\mathscr{F}} \tfrac{1}{2}\delta(\omega - \omega_0)e^{-j\frac{\pi}{2}} + \tfrac{1}{2}\delta(\omega+\omega_0)e^{j\frac{\pi}{2}}$$

Can anyone tell what I'm doing wrong here?

Two things

Time shift equation is wrong. Should be $e^{-j \omega n_d}$.
You need to express the cosine as a time shift $cos(\omega (n - \pi /2 /\omega))$ Time shift is not the same as phase shift and the  time shift is actually pi/2 divided by the frequency


Causality of the given system I am given the following discrete system

$$y[n]=x[-n]$$

where

$$x[n]=n+3$$.

Now from what I understand, since

$$y[n]=x[-n]$$

therefor the system is causal as the output is depending on the past values of the input. Right?

You are correct that $y[n]$ only depends on past input values for $n \gt 0$, but when $n \lt 0$, that is not the case --- then $y[n]$ is dependent on future input values, regardless of what $x$ is.

To see this:
$$
y[-2] = x[-(-2)] = x[2]
$$
so the output at time $n=-2$ is dependent on the input at time $n=2$.


Lower bound of weighted average of sequence Can anyone prove that
$$\mathrm{avg}\left(\frac {a_i}{\left(1+a_i\right)^2}\right) \ge \frac{\mathrm  {avg}({a_i})}{\left(1+\mathrm{avg}{(a_i)}\right)^2}$$
for a sequence of positive valued elements $a_i, \quad\text{with}\quad i= 1,2,\ldots,N$?

You can recognize the function $$f: x\mapsto \frac{x}{(1+x)^2}\,.$$
The question can be rephrased as: how do compare the average of the function values $\overline{f(a_i)}$ with the function of the average value $f(\overline{a_i})$? Graphically, this can be illustrated by the following graph (spoiler: from Behold! Jensen's inequality). Here, the average is represented by $\mathbb{E}$, the statistical expectation, akin to the classical average.

The intuition here is the following: for a convex function, as the above, the average of the $y=f(x)$ ($\mathbb{E}[f(x)]$) is higher than the function of the averages of $x$  (or $f(\mathbb{E}[x])$).
In the discrete form, Jensen's inequality can be proven by induction Proof 1 (finite form).
This inequality is reversed for a concave function. The reverse inequality is the one you are interested in. Now, function $f(x)$ is not concave on $[0,\infty[$, graph from RechnerOnline:

With $a_1=0$ and $a_1=1$, we get $$\overline{f(a_i)}-f(\overline{a_i}) = 1/8 - 2/9 \le 0 \,.$$
With $a_1=1$ and $a_1=9$, we get $$\overline{f(a_i)}-f(\overline{a_i}) = 17/100 - 5/36 \ge 0 \,.$$
Hence, the inequality has a sign change. It would be valide for instance on some interval $[\alpha,+\infty)$, such that $f$ is concave.


Matlab FIR filter numerator coefficients to python I created a equiripple FIR filter in Matlab using filterDesigner and exported it to an ASCII file. The file contains the numerator values as a list of 51 numbers tagged under Numerator: (this is all the data the file has)
Now I read this 51 values into a python array named num and try to recreate the filter in python using scipy.signal.lfilter
from scipy import signal

x_filt = signal.lfilter(num, np.ones_like(num), x)  # x is my input signal

This is my original signal

Here's the filtered signal from Matlab (Fpass-800Hz, Fstop-1200Hz)

But this is what I get from the python code using scipy.signal.lfilter

Can somebody please show me what I am doing wrong here?

EDIT: (After getting the answer) Here's the correct python code
x_filt = signal.lfilter(num, 1, x)

You don't want np.ones_like(num) in the denominator.  That's making a denominator of all ones, with lots of possibly unstable poles.  Your denominator in Python needs to be [1,0,0,0,0,0,...]. Or just [1] if possible.


Dominant frequencies of QRS complex I read in several sources that the dominant frequencies of the QRS complex in the ECG signal are between more or less 10-40 Hz. What I don't understand is that normally heart beats occur every other second (for a normal heart rate of 60 Beats per minute) - that is a frequency of 1 Hz. In each heart beat there is one QRS complex. So why is the dominant frequency not 1 Hz but 10-40 Hz ?

The QRS complex is a "spike" that lasts about 100 ms. There is also the P wave that comes before Q which indicates atrial contraction, and the ST segment after. The QRS complex is associated with the ventricles contracting. Blood is being pushed into the major arteries. This part of the cycle is called systole. QRS is like a spark plug if you would permit a simplistic analogy. It represents timing not blood flow. The waveform that represents blood flow is the instantaneous pressures in each  artery and vein and those waveforms are much smoother. The pulse on your wrist is a pressure waveform.

See here for a good tutorial.

While an average heart rate is about one a second, there is a lot of beat to beat variability even when an individual is resting. It's not really periodic in the Fourier decomposition sense.   Signals that are not periodic have a continuous frequency spectrum and the majority of the energy for QRS is in the 40 to 100 Hz band. Blood pressure waveforms have most of their energy at the lower frequencies.


Making an EQ from cascading filters I hope this isn't a dublicate . I checked but maybe didn't see everything.

Here goes :

I've been scratching my head over some EQ related questions. I read a lot of papers but my limited math and DSP knowledge prevents me from fully understanding them, so I thought I'd ask the experts.

I hghave a good library (http://www.kvraudio.com/forum/viewtopic.php?t=249926 in case you wonder) of standalone filters. There are biquad filters and higher order filters. From that, I'd like to create an equalizer. This is C++ but I don't think it's relevant;. My questions refer to the signal processing side.

1) Can I just casacade the filters :

input signal -> processed through f1 -> ... -> processed through fn -> output signal ?

2) Given that filters are not perfect, does the order matters ?

3) Do you know how I can calculate the transfer curve of the EQ, from each filter's transfer curve ?

4) Finally, what kind of filter is best ? From what I read, I think it's better to use high order filters instead of quadratic, so I think I should use Butterworth or Chebyshev filters. Am I right ?

I'm currently a bit lost so thanks for ANY insight ! What I would like is to maximize the sound quality 9i.e : make a transparent eq, nothing too fancy) , even if the CPU cost is high, or even better, provide options like optimise for speed/optimise for quality

P.S. : If you've got what it takes to create this, and you're interested in contract work, please feel free to contact me.

1) Can I just casacade the filters :

  input signal -> processed through f1 -> ... -> processed through fn -> output signal ?

Yep, that's how you do it.

  2) Given that filters are not perfect, does the order matters ?

Probably not.  You're probably working with ≥32-bit floats on a computer processor?  So you don't really need to worry about clipping or noise floor issues.

Hardware DSPs are often in fixed-point formats, so they can clip and then it matters because you don't want to boost a frequency to the point of clipping and then attenuate it with a subsequent filter, or cut a frequency down into the noise floor and then boost it back up, dragging the noise up with it.

  3) Do you know how I can calculate the transfer curve of the EQ, from each filter's transfer curve ?

If you mean "frequency response", just multiply them together.  Or convert them to dB and then add them together. dB is logarithmic, so it converts multiplication into addition.

  4) Finally, what kind of filter is best ? From what I read, I think it's better to use high order filters instead of quadratic, so I think I should use Butterworth or Chebyshev filters. Am I right ?

No, you should always use biquad filters and combine them in cascade to make higher-order filters.  If you search for "second-order sections" you can find instructions for making higher-order Butterworth and Chebyshev filters this way.  This minimizes numerical error from some intermediate numbers being extremely large while others are extremely small.


How to determine a quantization table for image compression with DCT? The quantization table is given on the JPEG Wiki page. How was this table derived? I'm currently doing a project where I have to come up with custom quantization table for an image, but I don't know how to begin. I figure if I learn how the JPEG table was derived, I might come up with a semi-original idea.

The standard JPEG quantization table is a actually an 8 x 8 matrix of weights whose elements divide the corresponding elements of the 2D-DCT coefficient matrix of an 8 x 8 image block.

The division effectively reduces the magnitude of the DCT coefficients and most of them are forced to be zero, due to the integer divison rounding.

The standard table is determined to yield a minimum distortion under maximum compression applied to the most typical image (i.e., exhaustive average of large image data set) and typical viewer visual characteristics.

As a result of this quantization, high frequency coefficients are reduced as much as possible.

Note that the JPEG quality parameter provides a means of scaling the standard quantization table weights, but no control over individual elements.

You can also design and use a specific quantization table for your needs, according to the rate-distortion criteria you have in your mind. But this requires you to embed Q table information within the compressed data bitstream.


Help in viewing the output of a neural network I used nftool in Matlab 2012 and trained a network. I gave the training inputs as x=[250:1] and targets as t=[250:1]. I used 10 hidden layers. I trained the network and got the results and exported the data to workspace. I finally got net as the variable in the workspace. I want to realize this network in hardware. So, I want the exact weights and bias values. I searched the net and typed net.IW and got 10 values, net.LW and got yet another 10 values and net.b and got 11 values. My network diagram is as shown below:

If IW is the Input Weight and LW is the Layer Weight, I realized the network diagram as shown below:

I got 10 values for IW, which would fit into 10 hidden layers, and 10 values for LW, which will fit into 10 hidden layers according to my realization. I got 11 values for bias, which will fit in to my network, since there are 11 (b) blocks in the network.
But, i am missing one IW and one LW value.

I want to know whether there is a mistake in my realization of the output of the network or I missed any of the values. Please help.

If I understand correctly your IW should be [1*10] as it takes 1-dim input to 10-dim hidden layer
than
LW should be [10*1] as it takes 10-dim hidden to 1-dim output,
b should be 10 for the hidden layer and 1 for the output layer

Why do you think your'e missing something?


Is there any modern information about a 12-hour periodicity in sonar noise? I have a copy of Principles of underwater sound by Urick, 3rd edition. In the chapter on The noise background of the sea: ambient noise level, on pages 215/216, it says
"The analysis of data in the frequency band 20 to 100 Hz revealed periodicity of 12 and 24 hours, with a maximum of noise at noon and midnight, local time, at all locations, independent of longitude. Although the change in level was small (normally 1.5 to 5 dB), it was undoubtedly real, but the origins of these changes in noise, which are apparently free from biological activity and synchronized to the mean solar day at diverse locations, remains a mystery."
The book was published in 1967, and my copy printed in 1983. I'd be fascinated to find out if there has been any new research into this. I did a little googling but couldn't find any references to this issue.

Some further googling seems to have shed some light: it appears that shipping traffic is highest at these hours, as demonstrated by Figure 3 in
http://cetus.ucsd.edu/Publications/Publications/McKennaIslandsSymp2009.pdf.
Why that is, I'm not sure. Also, this publication isn't about global shipping patterns, but only that in a certain region on the west coast of the USA - but it certainly goes some way to explaining the noise detected by sonar.
I know this was a very niche question, but hopefully this is helpful to someone who may stumble on my question in the future.


Can anyone explain how does CZT (Chirp Z Transform) really help in 'spectral zooming'? I found some explanation alongwith the Matlab Code here,
but I can't figure out, without a good example, why would this result in spectral zooming? (or develop some time intuition like viewing it as time resolution vs. frequency resolution thing). Any help (hopefully with an example) will be appreciated.

The CZT allows for a fairly general evaluation of the Z transform - the more general evaluation path looks like a spiral, so it has a radial component step size as well an angular step size.For  spectral zooming, you're only using a subset of this. You're evaluating around the unit circle and only for a small set of frequencies.

The Zoom FFT can be implement in a more traditional manner:

Frequency shift by multiplying time domain signal by a complex exponential - baseband around the center frequency.
Low-pass filter and decimate - to reduce the sampling frequency.
Apply a (shorter) FFT

Note that the FFT in step 3 is shorter because of the reduced sampling rate. To get a frequency resolution of $\nabla f$ you need a time series of length $T=1/\nabla f$ and at a lower sampling frequency you can get that time length signal in fewer samples. This series of steps applies a series of operations to the signal.

Now another way of looking at those series of steps is to apply them to the filter rather than the signal. Now you're shifting a low-pass filter in frequency so it becomes a complex band-pass filter and the decimation step causes duplicates of your band-pass signal to appear around 0 Hz (and throughout the entire frequency spectrum). At this point you can apply a shorter FFT. To make the two series of steps exactly equivalent there are some restrictions on the sampling frequency - If you shift the filter to a center frequency of $f_c$ then you have to be careful of the decimation factor if you want the frequency content at $f_c$ to be mapped exactly to 0 Hz.

The CZT for Zoom FFT is essentially doing these same steps. If you examine the CZT from this point of view you can work the mathematical equivalence.


How can I tell if an image is part of another image? I have one image that may or may not be part of another image (possibly resized).

How can I tell if that's the case, and pointing out the region of the first image that contains the second?
For example, I'd like to tell that the second image is an enlarged portion of the first.

This problem can be treated as signal/image detection or image registration. From the image registration perspective, there is a good tutorial that is exactly what you need: Registering an Image Using Normalized Cross-Correlation


IFFT-based OFDM system Please explain, mathematically, how can a bank of modulators be replaced with the IFFT in an OFDM system?

Say you want to transmit complex symbols $s_0,s_2,\ldots,s_{N-1}$ using frequency-division multiplexing. In the "bank of modulators" approach, you  would assign a carrier frequency $f_k$ to each symbol and create the signal
$$x(t) = \sum_{k=0}^{N-1} s_k \exp(j2\pi f_kt).$$

Contrast this with the definition of the inverse discrete Fourier transform:
$$x[n] = \sum_{k=0}^{N-1} s_i \exp(j2\pi nk/N).$$

So, if you interpret the symbols as the magnitude and phase of each carrier, and choose the carriers as $f_k=k/N$, then you'll see that you can generate $x(t)$ simply by converting $x[n]$ to an analog signal.

Note that:

I've explained how to transmit one block of symbols; extending this to multiple blocks has a slight complication and requires inserting a cyclic prefix to $x[n]$.
I'm assuming a rectangular pulse shape, which results in sinc-shaped carriers in frequency. The beauty of this method is that the carriers are orthogonal even if they overlap. In practice, this means that you can recover the symbols $s_k$ by performing the DFT of $x(t)$.
It is perfectly possible to use other pulse shapes with this method.


How signal moments work? We know that moments extract important features of a signal that has no redundancy. For example Zernike moment gives us rotation invariant data of a signal like image or Fourier gives us frequency of a signal to us.

But what is the theory behind it?
How a moment performs this?
How basis function of a moment is designed?

I think it has roots in many different math theory but I can't find out a good source to know.

I dont think there is a common framework for defining all those operators. Each of them were developed by different people, at several moments of the history. They are being integrated day by day onto a common discipline called DSP.

As clarification, an operator (like a moment or any other) and the "redundancy" of a signal are very different and unrelated concepts.

IF you want a deeper answer, list all the moments you are considering, and all the different math theories you are reviewing for each of them.


Why NRZ-L and NRZ-I line coding have DC component? I can find the effect of having a DC component in the Internet which is "When the AC waveform has a DC component, the average voltage would be equal to the DC voltage instead".

But in my text book, it says NRZ-L and NRZ-I both have a DC component although NRZ-L and NRZ-I can have the average voltage $0$.

Why is it?

Assuming bi-polar encoding, and assuming that the number of transmitted 1's is equal to the number of 0's, then over very long time intervals the positive pulses "cancel out" the negative pulses and the average voltage is zero.

However, over smaller time intervals, there's likely to be more 1's than 0's, or vice versa. This is what introduces a DC component to the signal.

Intuitively, you may think about it this way: if you filter the line-encoded signal with a very narrow low-pass filter, you'll see a non-zero, varying output. The DC value is the expected value of the filter's output.

In order to reduce or eliminate DC you need a line encoding that forces the number of positive and negative pulses to be the same, even over very short time spans. One example is the Alternate Mark Inversion (AMI) encoding.