Help with a $z$-transfrom Problem I have the function $$(1-e^{-2n})u[n-1]$$ where $u[n]$ is the step input.

I want to find the $z$-transform for this function. I know that the transform of $1-e^{-2n}$ will be $$\frac{z}{z-1} - \frac{z}{z-e^{-2}}$$ Now do I shift each element by one unit to the left to get the entire $z$-transform?

You might rewrite your signal as

$$x(t)=(1-e^{-2t})u(t-1)=(1-e^{-2t})u(t)-(1-e^{0})\delta(t)=(1-e^{-2t})u(t)$$

The idea is to express the signal $x(t)$ as the difference of two signals. The first signal has the step function starting at $t=0$, and the second function removes the contribution of the step at the sample $t=0$. In your case, even the function has $x(0)=0$, so it does not matter, if you have the step function starting at $t=0$ or $t=1$. So, you can just apply your solution above and you have it.

To explain my steps in a bit more detail: First note, that
$$u(t)=u(t-1)+\delta(t).$$(we are in a discrete setting with $t\in\mathbb{Z}$). Therefore $$u(t-1)=u(t)-\delta(t).$$ Let's insert this into your signal (we take a general $h(t)$. In your case you'd have $h(t)=1-\exp(-2t)$:

$$x(t)=h(t)u(t-1)=h(t)(u(t)-\delta(t))=h(t)u(t)-h(t)\delta(t)$$

The last part, i.e. $h(t)\delta(t)$ becomes equal to $h(0)\delta(t)$ due to the sifting property of the Dirac. So, in general you have

$$\mathcal{Z}\{h(t)u(t-1)\}=\mathcal{Z}\{h(t)u(t)\}-\mathcal{Z}\{h(0)\delta(t)\}=H(z)-h(0),$$

where $H(z)$ is the Z-transform of $h(t)$, and we assume that $h(t)$ is causal, i.e. $h(t)=h(t)u(t)$. (This property is implicitely assumed in the standard Z-transform).

In your case, you have $h(0)=0$, so the second term vanishes.


How does MATLAB recover picture from magnitude spectrum alone? This is the transformation I did.
The code fft2() the Lena picture than ifft2() it back to the original. Add some abs() so it actually displays it.

According to this post an image cannot be reconstructed without the phase information, yet I see no indication of phase being used anywhere here.

I'm guessing MATLAB still uses it but keeps it out of sight for convinve. If so is it possible to view the phase component of a given picture as well as the magnitude spectrum?

CODE:

[x,map] = imread('PATH'); %import image
i = ind2gray(x,map);   % to greyscale

I = fft2(i); % 2D FFT

I = fftshift(I); % centre

ir = ifft2(I);

figure(1); imshow(abs(i),[]) % original
figure(2); imshow(log(double(abs(I))),[]) % original
figure(3); imshow(abs(ir),[]) % original

I guess I worded this really badly. How do I see the phase plot of the image's transform? I guess that sounds less stupid.

Your code uses the phase for the reconstruction. Have a look at the output of fft2(x); they are complex numbers, i.e. the contain phase and magnitude. Have a look at this code:

%%
[x,map] = imread('http://www.cs.cmu.edu/~chuck/lennapg/lena_std.tif'); %import image
%%
i = rgb2gray(x);   % to greyscale

I = fft2(i); % 2D FFT

I = fftshift(I); % centre
mag = abs(I);
phase = angle(I);
I_magonly = abs(I);
I_phaseonly = I ./ abs(I);

figure(1);
subplot(321); imshow(abs(i),[]); title('Original');
subplot(322); imshow(log(double(abs(I))),[]); title('spectrum magnitude');
subplot(323); imshow(phase,[]); title('spectrum phase');
subplot(324); imshow(abs(ifft2(I)),[]); title('reconstructed');
subplot(325); imshow(abs(ifft2(I_magonly)),[]); title('reconstructed magnitude only');
subplot(326); imshow(abs(ifft2(I_phaseonly)),[]); title('reconstructed phase only');

It loads the Lena image, creates the FFT, and displays phase and magnitude information. Then it tries to reconstruct only based on phase or magnitude. Apparently, phase contains more information (for our eyes) than the magnitude. But you need to both exactly reproduce the image from its Fourier transform.


Alternative to Extended Kalman Filter When Prediction Function Is Not Differentiable I am looking at a tracking problem. It can be modelled similarly to the Extended Kalman Filter:

$$
\begin{array}{rcl}
\mathbf{x}_k &=& \mathbf{f}(\mathbf{x}_{k-1}, \mathbf{u}_k) + \mathbf{w}_k\\
\mathbf{z}_k &=& \mathbf{h}(\mathbf{x}_k) + \mathbf{v}_k
\end{array}
$$

Here $\mathbf{x}$ is what I want to track, $\mathbf{u}$ are the controls, $\mathbf{z}$ are measurements and $\mathbf{w}$ and $\mathbf{v}$ are noise. The function h is the identity in my case. My problem is that the nonlinear function $\mathbf{f}(\mathbf{x},\mathbf{u})$ is not given in closed form (it's not something simple like $\sin(\mathbf{x}^T \mathbf{u}))$. I can only write a function in Java (or another programming language) which evaluates $\mathbf{f}(\mathbf{x},\mathbf{u})$ with an algorithm. Hence, I can't use the Extended Kalman Filter, because I can't build a derivative of $\mathbf{f}(\mathbf{x},\mathbf{u})$.

The only alternative that came to my mind is to use a particle filter, where it should not be a problem that $\mathbf{f}()$ is computed by an algorithm.  However, it will probably be too computationally expensive, because $\mathbf{f}$ is computationally expensive and I would need to call it for each particle.

Does anybody have an idea for an alternative filter?

I will tell you something, even if it is differntiable, use Unscented Kalman Filter for any non linear case.

This flavor of Kalman Filter, based on the Unscented Transform, is almost always superior to the Extended Kalman due to its properties.
The main reason is it is able to better predict the mean and variance (Which all Kalman Filter needs) of the distribution after the non linear function applied.

In most cases it is even simpler to calculate hence more computationally efficient.
Think of it as Poor's man Particle Filter (Though the sampling points are deterministic).


what does it mean to have a decorrelated colour space? I'm working on the problem of colour transfer between images. In the literature there's a common practice that consists in performing the transformations on images in what is called a decorrelated colour space, such as CIE L* a * b* or YUV. The argument that's given most of the time is that a decorrelated colour space allows us to transform each one of the 3 channels of the image without affecting the other 2 channels, but I don't quite see how so. What is a decorrelated colour space in the first place?. As I conceive it, CIE Lab is decorrelated means that there's an independance between channels, if for example the value of l* grows bigger, this doesn't necessarly imply that the values of a* or b* increase too, but so is the case in rgb I guess. So, what am I missing?

edit: to be clear this answer describes why Lab can be described as a decorrelated color space. This does not imply that decorrelation is the main benefit of using Lab (see many answers on why Lab is useful)

If you plot all the RGB colors of a standard RGB image of the natural world you will notice something (see below), the values tend to fall on a diagonal line that isn't one of the axes of the color space.

In fact if you know the red and green values of a color in RGB color space you could make a good guess as to the blue color, it would be the blue that makes the point fall near that clustered line. This is because the values are all very correlated, they are good predictors of each other.

  As I conceive it, CIE Lab is decorrelated means that there's an independence between channels, if for example the value of l* grows bigger, this doesn't necessary imply that the values of a* or b* increase too, but so is the case in rgb I guess. So, what am I missing

So as you can now see, as the values of one RGB dimension increase, the others are also likely to increase.

As it turns out, that clustered diagonal line running through the RGB cube represents the luminance dimension.

In CIE Lab, that luminance dimension we see in the RGB cube as a diagonal running from (0,0,0) to (1,1,1) has been made into one of the axes of the color space. If you know the a* and b* values of a color you still don't know much about the L value. The CIE Lab dimensions are decorrelated from each other because they don't predict each other.


What is the physical meaning of the convolution of two signals? If we convolve 2 signals we get a third signal. What does this third signal represent in relation to the input signals?

There's not particularly any "physical" meaning to the convolution operation. The main use of convolution in engineering is in describing the output of a linear, time-invariant (LTI) system. The input-output behavior of an LTI system can be characterized via its impulse response, and the output of an LTI system for any input signal $x(t)$ can be expressed as the convolution of the input signal with the system's impulse response.

Namely, if the signal $x(t)$ is applied to an LTI system with impulse response $h(t)$, then the output signal is:

$$
y(t) = x(t) * h(t) = \int_{-\infty}^{\infty}x(\tau)h(t - \tau)d\tau
$$

Like I said, there's not much of a physical interpretation, but you can think of a convolution qualitatively as "smearing" the energy present in $x(t)$ out in time in some way, dependent upon the shape of the impulse response $h(t)$. At an engineering level (rigorous mathematicians wouldn't approve), you can get some insight by looking more closely at the structure of the integrand itself. You can think of the output $y(t)$ as the sum of an infinite number of copies of the impulse response, each shifted by a slightly different time delay ($\tau$) and scaled according to the value of the input signal at the value of $t$ that corresponds to the delay: $x(\tau)$.

This sort of interpretation is similar to taking discrete-time convolution (discussed in Atul Ingle's answer) to a limit of an infinitesimally-short sample period, which again isn't fully mathematically sound, but makes for a decently intuitive way to visualize the action for a continuous-time system.


Duality of the continuous-time Fourier transform - derivation and notation Suppose we have the Fourier transform pair $x(t)$ and $X(\omega)$ such that
 $$X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \mathrm{d}t$$

The duality property states that $X(t)$ and $2\pi x(-\omega)$ constitute a Fourier transform pair. I was trying to prove this statement when I ran into the following problem: the Fourier transform of a signal is sometimes denoted by $X(j\omega)$, such that $$X(j\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \mathrm{d}t$$

Using this notation, how do we even state the duality property? Note that evaluating $X(t)$ in this case would effectively replace $j\omega$ (as opposed to just $\omega$ in the previous case) by $t$ which is clearly wrong. Would we have to write $X(jt)$ for consistency?

Returning to the question of deriving the duality property and using the notation I used in the beginning, what is the flaw in the following approach?

Since $$X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \mathrm{d}t$$
we let $t$ be dummy variable $u$, and evaluate the function $X(\cdot)$ at $t$. We then have $$X(t) = \int_{-\infty}^{\infty} x(u) e^{-j ut} \mathrm{d}u$$
With the change of variables $\omega = -u$ we get $$X(t) = -\int_{-\infty}^{\infty} x(-\omega) e^{j\omega t} \mathrm{d}\omega$$ $$X(t) = -\frac{1}{2\pi}\int_{-\infty}^{\infty} 2\pi x(-\omega) e^{j\omega t} \mathrm{d}\omega$$

Which is almost the correct result except for the additional minus sign. Where did I go wrong? Thank you in advance.

I would also like to add that I feel uncomfortable about simply changing $t$ and $\omega$ in the functions' argument because they have different dimensions.

You've answered your first question yourself: if you use the notation $X(j\omega)$ then you'll need to use $X(jt)$, because you simply replace $\omega$ by $t$.

The sign error in your derivation comes from the fact that you forgot to change the sign of the limits of the integral. After substituting $\omega=-u$ you should get

$$X(t) = -\int_{\infty}^{-\infty} x(-\omega) e^{j\omega t} \mathrm{d}\omega = \int_{-\infty}^{\infty} x(-\omega) e^{j\omega t} \mathrm{d}\omega$$


Power of an angle-modulated wave Let $$\phi(t) = A\cos(\omega_c t + km(t)) \tag{1}$$ be an angle-modulated wave(FM or PM). What's the power of $\phi(t)$? Intuitively, it seems that the answer is $$P = \frac{A^2}{2} \tag{2}$$ since the amplitude $A$ remains constant but I couldn't prove that using the definition which gives $$\lim_{T \to +\infty} \frac{1}{T}\int_{-\frac{T}2}^{\frac{T}2} A^2\cos^2(\omega_ct + km(t))dt \tag{3}$$ It seems we can't evaluate $(3)$ in the general case $m(t)$.

Let's consider the complex signal
$$\phi_c(t)=Ae^{j(\omega_ct+km(t))}\tag{1}$$
The power of $(1)$ is readily computed as
$$P_{\phi_c}=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}\big|\phi_c(t)\big|^2dt=A^2\tag{2}$$
IF $(1)$ were the analytic signal corresponding to the real-valued signal $\phi(t)$ then the power of $\phi(t)$ would simply be $P_{\phi_c}/2=A^2/2$.
In general $(1)$ is not analytic. However, in all practical cases of angle modulation, a sufficiently high carrier frequency makes $(1)$ analytic for all practical purposes, i.e., there are no measurable negative frequency components. Or, equivalently, a practical angle modulated signal is a bandpass signal and has no components at or close to DC. If this is the case, the power of $\phi(t)$ is indeed given by $A^2/2$.
In all other (artificial) cases, the power of $\phi(t)$ cannot generally be shown to equal $A^2/2$.


FIR or IIR filter from poles and zeros Is it possible to differentiate between FIR and IIR Filter just by knowing the poles and zeros of the system?

If yes, how?
If not, why?

Thank you.

Axel Mancino's answer is correct for causal filters. In general, FIR filters have poles at either $z=0$ or $|z|\rightarrow\infty$, or both.

Take as an example a fourth-order causal FIR filter:

$$H_1(z)=a+bz^{-1}+cz^{-2}+dz^{-3}+ez^{-4}\tag{1}$$

Clearly, $H_1(z)$ has all its poles at $z=0$. An anti-causal FIR filter such as

$$H_2(z)=az^4+bz^3+cz^2+dz+e\tag{2}$$

has all its poles at $|z|\to\infty$. And, finally, a general non-causal FIR filter has poles at $z=0$ as well as at $|z|\rightarrow\infty$:

$$H_3(z)=az^2+bz+c+dz^{-1}+ez^{-2}\tag{3}$$

In sum, a filter with all its poles at $z=0$ or $|z|\to\infty$ (or both) is an FIR filter. Filters with poles in the region $0<|z|<\infty$ must be IIR.


Calculating SNR for ECG signals I'm reading through the literature for finding the best way to calculate SNR for a ECG signal. While there are many algorithms for improving SNR, I couldn't find any that details any standard procedure for evaluating the SNR itself. The following are my questions.

What is the SNR of a typical ECG signal from the MIT/BIH database? What is considered a good SNR for a typical ECG signal?
http://www.physionet.org/physiobank/database/mitdb/
What would qualify as Signal and Noise for ECG? The assumption that I worked with is to apply a bandpass filter with cut-offs fc=[0.5Hz,45Hz] on the recorded ECG signal to get the Clean ECG signal, then subtract the clean ECG signal from recorded signal to obtain noise and then use the matlab snr function to evaluate  SNR with the clean signal and noise as arguments. But, in some papers it is recommended to reject everything except the QRS complex to get the clean signal. Is there any standard/recommended procedure for calculating SNR for ECG?
http://en.wikipedia.org/wiki/QRS_complex
What would be the best way to calculate ECG SNR in real-time in an embedded platform?

I don't know about real time application normally what I am doing is to calculate $\mathrm{PRD}$ first and calculate $\mathrm{SNR}$ from $\mathrm{PRD}$

$$\mathrm{SNR}=-20\log_{10}(0.01\mathrm{PRD})$$

$\mathrm{PRD}$ -  Percentage Root Mean Square Deviation


Is deep learning killing image processing/computer vision? I'm looking forward to enroll in an MSc in Signal and Image processing, or maybe Computer Vision (I have not decided yet), and this question emerged.

My concern is, since deep learning doesn't need feature extraction and almost no input pre-processing, is it killing image processing (or signal processing in general)?

I'm not an expert in deep learning, but it seems to work very well in recognition and classification tasks taking images directly instead of a feature vector like other techniques.

Is there any case in which a traditional feature extraction + classification approach would be better, making use of image processing techniques, or is this dying because of deep learning?

On the top of this answer, you can see a section of updated links, where artificial intelligence, machine intelligence, deep learning or and database machine learning progressively step of the grounds of traditional signal processing/image analysis/computer vision. Below, variations on the original answer.
For a short version: successes of convolutional neural networks and deep learning have been looked like as a sort of Galilean revolution. For a practical point of view, classical signal processing or computer vision were dead... provided that you have enough or good-enough labeled data, that you care little about evident classification failures (aka deep flaws or deep fakes), that you have infinite energy to run tests without thinking about the carbon footprint, and don't bother causal or rational explanations. For the others, this made us rethink about all what we did before: preprocessing, standard analysis, feature extraction, optimization (cf. my colleague J.-C. Pesquet work on Deep Neural Network Structures Solving Variational Inequalities), invariance, quantification, etc. And really interesting research is emerging from that, hopefully catching up with firmly grounded principles and similar performance.
Updated links:

2021/04/10: Hierarchical Image Peeling: A Flexible Scale-space Filtering Framework
2019/07/19: The Verge: If you can identify what’s in these images, you’re smarter than AI, or do you see a ship wreck, or insects on a dead leaf?

2019/07/16: Preprint: Natural Adversarial Examples

We introduce natural adversarial examples -- real-world, unmodified,
and naturally occurring examples that cause classifier accuracy to
significantly degrade. We curate 7,500 natural adversarial examples
and release them in an ImageNet classifier test set that we call
ImageNet-A. This dataset serves as a new way to measure classifier
robustness. Like l_p adversarial examples, ImageNet-A examples
successfully transfer to unseen or black-box classifiers. For example,
on ImageNet-A a DenseNet-121 obtains around 2% accuracy, an accuracy
drop of approximately 90%. Recovering this accuracy is not simple
because ImageNet-A examples exploit deep flaws in current classifiers
including their over-reliance on color, texture, and background cues.
We observe that popular training techniques for improving robustness
have little effect, but we show that some architectural changes can
enhance robustness to natural adversarial examples. Future research is
required to enable robust generalization to this hard ImageNet test
set.

2019/05/03: Deep learning: the final frontier for signal processing and time series analysis? "In this article, I want to show several areas where signals or time series are vital"
2018/04/23: I just come back from the yearly international conference on acoustics, speech and signal processing, ICASSP 2018. I was amazed by the quantity of papers somewhat relying on deep Learning, Deep Networks, etc. Two pleanaries out of four (by Alex Acero and Yann LeCun) were devoted to such topic. At the same time, most of the researchers I have met were kind of joking about that ("Sorry, my poster is on filter banks, not on Deep Learning", "I am not into that, I have small datasets"), or were wondering about gaining 0.5% on grand challenges, and losing the interested in modeling the physics or statistical priors.
2018/01/14: Can A Deep Net See A Cat?, from "abstract cat", to "best cat" inverted, drawn, etc. and somehow surprizing results on sketches
2017/11/02: added references to scattering transforms/networks
2017/10/21: A Review of Convolutional Neural Networks for Inverse Problems in
Imaging
Deep Learning and Its Applications to Signal and Information Processing, IEEE Signal Processing Magazine, January 2011

Deep learning references "stepping" on standard  signal/image processing can be found at the bottom. Michael Elad just wrote  Deep, Deep Trouble: Deep Learning’s Impact on Image Processing, Mathematics, and Humanity (SIAM News, 2017/05), excerpt:

Then neural networks suddenly came back, and with a vengeance.

This tribune is of interest, as it shows a shift from traditional "image processing", trying to model/understand the data, to a realm of correctness, without so much insight.
This domain is evolving quite fast. This does not mean it evolves in some intentional or constant direction. Neither right nor wrong. But this morning, I heard the following saying (or is it a joke?):

a bad algorithm with a huge set of data can do better than a smart  algorithm with pauce data.

Here was my very short try: deep learning may provide state-of-the-art results, but one does not always understand why, and part of our scientist job remains on explaining why things work, what is the content of a piece of data, etc.
Deep learning used too require (huge) well-tagged databases. Any time you do craftwork on single or singular images (i. e. without a huge database behind), especially in places unlikely to yield "free user-based tagged images" (in the complementary set of the set "funny cats playing games and faces"), you can stick to traditional image processing for a while, and for profit. A recent tweet summarizes that:

(lots of) labeled data (with no missing vars) requirement is a deal
breaker (& unnecessary) for many domains

If they are being killed (which I doubt at a short term notice), they are not dead yet. So any skill you acquire in signal processing, image analysis, computer vision will help you in the future. This is for instance discussed in the blog post:  Have We Forgotten about Geometry in Computer Vision? by Alex Kendall:

Deep learning has revolutionised computer vision. Today, there are not
many problems where the best performing solution is not based on an
end-to-end deep learning model. In particular, convolutional neural
networks are popular as they tend to work fairly well out of the box.
However, these models are largely big black-boxes. There are a lot of
things we don’t understand about them.

A concrete example can be the following: a couple of very dark (eg surveillance) images from the same location, needing to evaluate if one of them contains a specific change that should be detected, is potentially a matter of traditional image processing, more than Deep Learning (as of today).
On the other side, as successful as Deep Learning  is on a large scale, it can lead to misclassification of a small sets of data, which might be harmless "in average" for some applications. Two images that just slightly differ to the human eye could be classified differently via DL. Or random images could be set to a specific class. See for instance Deep neural networks are easily fooled: High confidence predictions for unrecognizable images (Nguyen A, Yosinski J, Clune J. Proc. Computer Vision and Pattern Recognition 2015), or Does Deep Learning Have Deep Flaws?, on adversarial negatives:

The network may misclassify an image after the researchers applied a
certain imperceptible perturbation. The perturbations are found by
adjusting the pixel values to maximize the prediction error.

With all due respect to "Deep Learning", think about "mass production responding to a registered, known, mass-validable  or expected behaviour" versus "singular piece of craft". None is better (yet) in a single index scale. Both may have to coexist for a while.
However, deep learning pervades many  novel areas, as described in references below. Many not-linear, complex features might be revealed by deep learning, that had not been seen before by traditional processing.

Deep learning for image compression
Real-Time Adaptive Image Compression, ICML 2017
Full Resolution Image Compression with Recurrent Neural Networks
End-to-end optimized image compression, ICRL 2017
Deep learning for video compression
Can deep learning be applied to video compression?
Deep learning for denoising, restoration, artifact removal
CAS-CNN: A Deep Convolutional Neural Network for Image Compression Artifact Suppression
Super-Resolution with Deep Convolutional Sufficient Statistics

Luckily, some folks are trying to find mathematical rationale behind deep learning, an example of which are scattering networks or transforms proposed by Stéphane Mallat and co-authors, see ENS site for scattering. Harmonic analysis and non-linear operators, Lipschitz functions, translation/rotation invariance, better for the average signal processing person. See for instance Understanding Deep Convolutional Networks.


How to compute the impulse response from a transfer function Say I have a transfer function, for example:

$$H(z)=\frac{1}{1+0.1z^{-30}}$$

How can I compute the impulse response?

(This is just an example, the important thing is that it is in closed symbolic form!)

If I understand correctly, if I could just rewrite $H(z)$ as some polynomial in $z^-1$ (by using a Taylor expansion), I could just read off the coefficients.

But is there another way?

Your filter is the all-poles IIR, this simplifies things a bit. Normally you can write transfer function in following form:

$H(z)=\dfrac{\sum_{i=0}^{P}b_{i}z^{-i}}{\sum_{j=0}^{Q}a_{j}z^{-j}} $

Going back to the discrete time domain you will get:

$y[n] =\dfrac{1}{a_0}\left( \sum_{i=0}^{P}b_ix[n-i]-\sum_{j=1}^{Q}a_jy[n-j]\right)$

Therefore in your case it is:

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

Calculation of particular values can be done by feeding Kronecker delta to your system ([1 0 0 0 0 0 0 ...] signal). Basically you will observe that something is happening every 30 samples ;) You can easily generate the plot in MATLAB/Octave, just remember that your filter has following coefficients:

b=[1]
a=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1]

By calling impz(b,a) function you will get similar plot (first 100 samples):


Output shape of numpy's rfft function For numpy's rfft function, the length of the output is about half of N, the length of the input sequence. What is the reason behind this? What are the frequency bins?

Check out the numpy docs ( https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft.html )

To sum up, real DFT normally takes a real input of N samples and returns a complex output of N/2+1 values.

What you end up with are the positive frequencies (up to the Nyquist frequency) and the
negative frequencies, which are just complex conjugates of their positive counterparts, are discarded or are not even computed.


How to simulate High Power Amplifier(HPA), like Solid State Power Amplifier (SSPA) I cannot understand how to simulate Solid State Power Amplifier (SSPA). When I directly write the Rap's SSPA model in Matlab, it does not produce BER of OFDM like published papers by other researchers'. Even no change is observed between BER with and without SSPA. But theoretically there should be significant changes between these (for low nonlinear parameter r, Asat). I have seen several matlab code in Matlab fileexchange, but those could not produce BER like published papers.  What is the wrong in my code given below.

Please help me out.

The model works fine.  You need to slow down and check your (or someone else's? ) steps.  I just worked through one pass of this, without pskmod, and only one value was high enough voltage to push the model into a nonlinear region.  Work through it step by step, you'll see that your x values are too low in magnitude to really get out of the linear region.  The x values have to approach A to get distortion.


Defining frequency axis for periodogram I am trying to simulate example 8.2.2 from Hayes, Statistical Digital Signal Processing and Modelling, on MATLAB.

Let $x(n)$ be a WSS process consisting of random phase sinusoid in unit variance.
$$ x(n) = A\sin(n\omega_{0} + \phi) + v(n)$$

$A = 5 $, $\omega_{0} = 0.4\pi$, $N = 256$

Can anyone help me figure out how to generate the frequency vector such that the peak appears at $0.4\pi$? I am able to achieve the spectrum that is shown in the example, but my x-axis can only be called frequency bin as of now.

As you have taken 1024 point FFT, multiply the x-axis with a factor of 2*pi/1024.

For now your peak is at sample number 200.
200*2*pi/1024 approximately equal to 0.4*pi


How to calculate a delay (correlation peak) between two signals with a precision smaller than the sampling period? I would like to calculate a delay between two signals in MATLAB. Normally, I would do this with a cross-correlation (xcorr), and then calculate the position of the correlation peak, but in this case I know that the delay is smaller than one sampling period. Therefore, the method I propose will yield either $0$ or $1$ sample as a delay.

One solution I can come up with, is to construct an interpolated version of the two signals, and then calculate the delay. However, then I still rely on the sampling period of this upsampled signals. Another possible solution is to calculate the cross-correlation peak for different blocks, and then calculate the mean delay. However, I am not sure whether this solution will converge to the right delay. A third solution I can come up with, is to apply a fractional delay to one of the two signals by means of a sinc reconstruction, and then calculate for which fractional delay the correlation peak is the largest. This requires a lot of calculations however, and I still rely on the fractional delays I apply (so I end up with the same problem as with the upsampling).

Therefore, I am wondering whether there is a more elegant or precise solution?

This is basically what @hooman suggests: fit a parabola to the three points near the peak of the sample cross-correlation of the data.

Using the formula for $p$ here:
$$
p = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma}
$$
where $\alpha,\beta,$ and $\gamma$ are the values of the sample cross-correlation just before the peak, at the peak, and just after the peak respectively.

I've assumed that the signal is a white noise sequence and that a simple FIR filter can be used to simulate a fractional delay.

The plot below shows the FIR filter in the top plot, and 100 estimates of the fractional delay using this technique.

R Code Below

#35112
library(signal)

fractional_delay_filter <- function(delay, length)
{
  eps <- 10^-10
  t <- seq(0,length-1)
  t[t == 0 ] <- eps
  filter <- sinc(t-delay,length)
  return(filter)
}

# https://ccrma.stanford.edu/~jos/parshl/Peak_Detection_Steps_3.html
peak_location  <- function(ykm1,yk,ykp1)
{
  return(0.5*(ykm1 - ykp1)/(ykm1 - 2*yk + ykp1))
}

estimate_delay <- function(x,xd)
{
  cc <- ccf(xd,x)
  idx <- which.max(cc$acf)
  return(peak_location(cc$acf[idx-1], cc$acf[idx], cc$acf[idx+1]) + cc$lag[idx])
}

true_delay <- 10.5
d5p1 <- fractional_delay_filter(true_delay,32)
fz <- freqz(d5p1)
#plot(fz)
gd <- grpdelay(d5p1)
#plot(gd)

T <- 1000
Nres <- 100

results <- rep(0,Nres)
for (n in seq(1,Nres))
{
  x <- rnorm(T,1)
  xd <- filter(d5p1, 1, x)

  #plot(x, type='l')
  #lines(xd[1:T],col='red')å

  results[n] <- estimate_delay(x,xd[1:T])
}

par(mfrow=c(2,1))
plot(d5p1, type='l', col='blue')
title('Fractional filter')

plot(results,col='blue')
lines(c(1,Nres), c(true_delay,true_delay), col='green')
title('Estimates and true delay')

The sinc function is from the signal library (which you've said privately you have trouble installing in R). It's defined as:

sinc <- function(omega, N)
{
  eps <- 10^-56
  return( sin(omega*N/2)/(omega/2 + eps) )
}

So changing that sinc function to:

sinc_pjk <- function(x)
{
  return(sin(pi * x) / (pi * x))
}

and the test delay to 10.1 samples, Olli's version is less biased:

<>


How can we prove the correctness of the integration property of the Laplace transform? I was going through an Electrical Engineering textbook for understanding the Laplace transform and came across the following proof for one of the properties of the Unilateral Laplace transform.
Integration property of the unilateral Laplace transform:

In the proof, it is stated that:
$$ e^{-st} \rightarrow 0 \mbox{ as } t \rightarrow \infty \ \ \ (i)$$
and therefore the term:
$$ -\frac{e^{-st}}{s}\int_{0-}^{t}{f(\tau) \ d\tau} = 0 \mbox{, as } t \rightarrow \infty \ \ \ (ii)$$
But my doubt is: Isn't there a case where the real part of $s$, gets cancelled by the function obtained after the integration of $f(\tau)$?
If that happens then we can't guarantee that (ii) would hold true right - i.e. the LHS in (ii) would not be zero right?
So, isn't the proof that has been provided wrong?

You are right that the argument in the proof is not correct, or at least misleading. The fact that $e^{-st}$ becomes zero as $t\to\infty$ is true for any $s$ with $\textrm{Re}\{s\}>0$. However, the integral
$$\int_0^{\infty}f(\tau)d\tau$$
might not exist, so we can't just compute the limit by claiming that the first term becomes zero.
We have to compute
$$\lim_{t\to\infty}-\frac{e^{-st}}{s}\int_0^{t}f(\tau)d\tau\tag{1}$$
The limit $(1)$ only equals zero for $s$ inside the region of convergence (ROC) of $f(t)$, i.e., we don't just require that $e^{-st}\to 0$ for $t\to\infty$, but we require it to decay sufficiently fast to compensate for the growth of the integral.


Does the cv2.waitKey argument determine the frame-rate while capturing a video? Kinda new to image processing, so pardon me if the question is too stupid.

Using the library cv2 in python to capture video from my webcam, which, so far as I can see as a programmer, is basically an array of images. But nothing specifies the frame rate of a video. Suppose I use the following code snippet to capture a video.

while True:
    capture, frame = cap.read()
    cv2.imshow('Being Captured. Hit Escape to Stop',frame)
    if cv2.waitKey(1)&0xFF==27:break

As I understand about the waitKey function, the camera captures a frame, waits for a millisecond, then captures the next one etc. So if a millisecond for each still image, does it mean the frame-rate is 1 KHz? Typical rates for most videos is 25-30 Hz, right? So can I adjust that just by changing the argument of waitKey in inverse proportion? Or, is my thought process wrong somewhere?

It does not determine the framerate, but it has an influence on the framerate.

The first caveat is that although you're only waiting 1ms between frames, the program has other work to do in order to render each frame, so the delay between frames is not purely caused by that 1ms wait.

The next caveat is that (from the doc)

  Since the OS has a minimum time between switching threads, the function will not wait exactly delay ms, it will wait at least delay ms, depending on what else is running on your computer at that time.

So, you may be thinking, okay, I won't get 1/delay framerate, but I will get at most 1/delay framerate.  Alas!  Caveat 3 strikes:

What if you spam keyboard input, e.g. by holding down a key and letting it repeat?  Then the wait function will get cut short.  Again, you may be led astray by thinking, but if I use delay=1 (millisecond) this will be trivial because no keyboard at a rate anywhere near 1 KHz!  But remember that this delay is a minimum amount of time it will wait, in reality the wait is much longer.

In an application I wrote (more complex than this though) using cv2.waitKey(1) after showing each frame, I get a 50% increase in fps by simply holding down a key on my keyboard.

 Hands-off FPS: 64.09
Key-repeat FPS: 90.99

The moral is that if you want stable, controlled FPS, you should pick another route, because there are so many caveats and unpredictabilities here.


Averaging magnitude squared coherence across multiple time series In a previous post, A. Donda had suggested that, in order to calculate the average magnitude squared coherence of more than one pair of time series (e.g. y1 and x1, and y2 and x2), one ought to follow the following steps:

sx1 = pwelch(x1, w, o, n);
sy1 = pwelch(y1, w, o, n);

sxy1 = cpsd(x1, y1, w, o, n);

sx2 = pwelch(x2, w, o, n);
sy2 = pwelch(y2, w, o, n);

sxy2 = cpsd(x2, y2, w, o, n);

sxa = (sx1 + sx2) / 2;
sya = (sy1 + sy2) / 2;

sxya = (sxy1 + sxy2) / 2; % NOTE THIS LINE %
ca = abs(sxya) .^ 2 ./ sxa ./ sya;

...where w,o,n are the window parameters, and the other variables should (hopefully) be self-explanatory.

Let me illustrate why this method can be problematic with an example, in the hope that someone might be able to help me with a solution.

Imagine two pairs of signals y1, x1, and y2,x2. All four signals contain the same frequency component (e.g. a 440Hz sine wave). In one pair, the sines waves are in phase w.r.t. each other. In the other pair, the sine waves are in inverse polarity.

Separate MSCs of both of these pairs of signals will yield approx. '1' at 440Hz in both results, identically, telling us that 440Hz is common to both signals in both pairs (i.e. the relative phase of each of the signals within a pair is ignored). If we now wish to calculate the average MSC, we must first compute the average PSDs and then the average CPSDs. Ideally, this average MSC for the two pairs would also yield a '1' at 440Hz. BUT the CPSD output is complex, so averaging the two CPSDs of the two pairs of signals according to the method above will result in cancellation, because the averaging is taken BEFORE the magnitude of the CPSDs is calculated.

I have offered one solution below - is this correct?

I believe I might have a solution... Editing A Donda's work above, I think a simple rearrangement of the last two lines is necessary:

sxya = (abs(sxy1) + abs(sxy2)) / 2;
ca = sxya.^ 2 ./ sxa ./ sya;

So here, the average of the magnitudes is taken, rather than the average of the complex CPSDs. Then, this average can be squared.


What is a $BT$ (Bandwidth-Time) product with reference to modulation? I would like to know what does Bandwidth-Time product mean. I understand that Bandwidth ($B$) = 1/Symbol Time ($T$), hence $BT = 1$.

But how can it vary?
What is its significance?
e.g. When we say GFSK is GMSK when $BT = 0.5$, what does that mean exactly?

The $BT$ product is the bandwidth-symbol time product where $B$ is the $-3\textrm{ dB}$(half-power) bandwidth of the pulse/filter and $T$ is the symbol duration. For different applications you will find varying recommended values. In GSM telephony for instance, a $BT=0.3$ is recommended. In satellite communications with GMSK, for near-earth missions the CCSDS recommends a $BT = 0.25$ whilst for deep-space/interplanetary missions, the use $BT=0.5$ is recommended. You find more details in this CCSDS report on bandwidth-efficient modulations. See page 2-2 and page 2-3 for the mentioned recommended values.

What does that mean ? Let's say we have $1$ bit per symbol ($T$ then correspond to bit time). For a $BT = 1$, the pulse shaping the symbol spreads over one bit period duration. For $BT = 0.25$, the spread is over $4$ bit periods, for $BT = 0.3$ the spread is over approximately $3$ bit periods, and for $BT = 0.5$ the spread is over $2$ bit periods.

This means that a smaller $BT$ product results into higher ISI and a compact spectrum. Measures need to be taken for the introduced ISI in this case much more than in the case of a higher $BT$ product where less ISI is introduced and we have much less compact spectrum.

In GMSK, one of its properties is it maintains a constant envelope and that's because of the Gaussian pulse applied prior to modulation. The GMSK pulse can be defined as in equation $(1)$ below$^1$:
$$
g(t) = \frac{1}{2T}\left[Q\left(2\pi B\cdot\frac{t-\frac T2}{\sqrt{\ln(2)}}\right)-Q\left(2\pi B\cdot\frac{t+\frac T2}{\sqrt{\ln(2)}}\right)\right]\tag{1}
$$

Where $Q(t)$ is the complementary cumulative distribution function defined as:
$$
Q(t)=\int_{t}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac 12 x^2\right)dx\tag{2}
$$
The figure below shows the pulse shapes for different values of the $BT$ product:

The spread of this pulse is inversely proportional to the $BT$ product and its peak amplitude directly proportional to the product. Again, this means that lower $BT$ results into a wider spread (over bit symbol period) and with lower peak amplitude and a high $BT$ results into a narrower spread with a higher peak amplitude. In conclusion, the pulse duration increases as the bandwidth of the pulse decreases.

Here they show a link between the $BT$ product, the filter’s $-3\textrm{ dB}$ cutoff frequency and the bit rate $R_b$ as:
$$
BT = \frac{f_{-3\rm dB_{\rm cutoff}}}{R_b}\tag{3}
$$
You can say that the $BT$ somehow determines the degree of filtering. More on fundamentals and properties of GMSK can be found in this paper $^2$ and this paper $^3$. In both papers discussions in relation to equation $(3)$ and the variations in the eye pattern as a result of $BT$ product values are given.

You can find extra stuff here, here, and here.

$[1]$: John G. Proakis, Digital Communications, 4th Edition, McGraw-Hill, 2000.

$[2]$: A. Linz and A. Hendrickson, "Efficient implementation of an I-Q GMSK modulator," in IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, vol. 43, no. 1, pp. 14-23, Jan 1996.

$[3]$: K. Murota and K. Hirade, "GMSK Modulation for Digital Mobile Radio Telephony," in IEEE Transactions on Communications, vol. 29, no. 7, pp. 1044-1050, Jul 1981.


Determining false positives in feature matching I extracted features from an image using SIFT and SURF. Then I applied few transforms on the image and again calculated the features. While matching the features, I noticed that there are some erroneous matches. For example, if using image of a face, some points on face are matched with points not on face in the transformed image.

Is there a way to detect those false positives ?

A simple way to reduce the number of erroneous matches is to tweak matching parameters, including trying different matching strategies. If you are using a simple distance threshold to determine what is a match, you can tighten that. Alternatively, you can use the ratio test (used in David Lowe's demo SIFT implementation) to eliminate ambiguous matches. There you can tightened both the distance threshold and the ratio threshold.

Unfortunately, the world is not perfect, so even if you tweak every parameter, you will still have some false matches. Then you need to use additional information and/or make additional assumptions about your problem. For example, you can assume that the images are related by a rigid transformation (e. g. affine). You can estimate the parameters of this transformation from your matches using an algorithm like RANSAC. Then you can reject the matches that do not fit this transformation as erroneous.

Here is an example of removing erroneous matches by assuming a similarity transformation, using the Computer Vision System Toolbox for MATLAB.


Dynamic time warping in percentage I'd like to compare 2 speech similarity and show it in terms of 0-100%.
Here's my approach:

1: Convert 2 audio waveform into audio spectrum by Fast Fourier Transform (FFT)
2: Use Dynamic Time Warping to compare 2 spectrums

Everything is fine when converting them into audio spectrum
But how to show the similarity in terms of percentage?

I've read this site: https://stackoverflow.com/questions/2168027/dynamic-time-warping-to-compare-two-audio-recordings

I'm working on exactly the same things with him.
However I have no idea about the answer.

Could anyone please help me?
Thanks!

There are no doubt plenty of ways you could define one of the measures in the linked question, but the question remains, after doing so, is the measure in any way meaningful?

This all depends on what you intend to do with the metric after computing it: If it's a perceptual estimate (quality of a Codec,  subjective similarity and so on), you will need a set of measured perceptual data to which you then attempt to fit (mathematically) your measure.


Using Neural Network for locating shifting resonant frequencies I have multiple FFT's taken from a sample at different pressures, through different analysis I can see that the resonant frequencies are shifting in the spectrum for each FFT at a different pressure.
Using conventional peak tracking has been difficult as the peaks increase/decrease in magnitude within the FFT as well as shifting in the spectrum.
Is it possible for a neural network to 'detect'/'pick out' these frequency values?
Any help or guidance is appreciated :)
Thanks!

Is it possible for a neural network to 'detect'/'pick out' these frequency values?

Yes, that sounds generally possible.
There's the Universal Approximation Theorem that says that a sufficiently large Neural network¹ can approximate any continuous function on anything isomorphic to $\mathbb R^N$ (and your FFT output is that), including things like

map this FFT output to a vector of resonant frequencies.

Now, what that says that if such a continuous function exists, then a neural network can reproduce the same.
It's then "just" a problem of training the weights of that neural network. That training usually happens through putting in data and judging the output (e.g. by how "correct" the output is, which requires you to already know the correct output, or by making sure the output probability distribution is as close as possible to what the distribution of the true phenomenon is).
However, whether or not that function exists is more of a physical / stochastic problem: For example, when I want a function

map the FFT of the acellerometer signal to the weight of the fattest cat of the oldest cat-owning person on this planet

no neural network could approximate that – simply because there's no relationship between the signal and the cat.
So, a Neural Network can't do magic: if things are stochastically independent, or ambiguous, then no estimation can be implemented, by a Neural Network or anything else. That's kind of an important thing to keep in mind, especially when reading the hype literature about Neural Networks that still occasionally show up:

If there's no mutual information between in- and output of the proposed mapping, then there's nothing any algorithm could do – be it a Neural Network, or anything else.

But in your case, I think you have good reason to assume there's a solid amount of relation between your input data and the output you want, so you should be fine.

¹ "sufficient" can be very large; and you need activation functions that aren't polynomials; what we call "neural networks" usually doesn't use polynomials.


For a standard second order transfer function, what is the equivalent time domain significance of $\zeta>0.707$? In frequency domain $\zeta>0.707$ implies no resonant peak in frequency response. But I am unable to correlate what is the exact time domain effect of this? If  $\zeta=0.707$ acts like an 'important damping value' in the frequency domain, shouldn't its importance be reflected in time domain as well. But in time domain $\zeta=1$ is where system behavior changes.

You are absolutely right that the damping value $\zeta=1/\sqrt{2}$ is important for the frequency domain behavior of a second-order system. However, as you've also noted, there is no one-to-one correspondence in the time domain. Generally we distinguish over-damped ($\zeta>1$), under-damped ($\zeta<1$), and critically damped ($\zeta=1$) systems, and these categories have distinct time-domain behaviors. Over-damped and critically damped systems only have real poles, whereas under-damped systems have complex poles, which results in oscillations in the impulse and step responses.

Note that only under-damped systems are useful for implementing frequency-selective filters, because the pole angles are related to the cut-off frequency, i.e. you need complex conjugate poles. So for frequency-selective filters you only consider under-damped systems, and then you can find different categories of under-damped systems according to their behavior in the frequency domain. A second-order Butterworth low pass filter has $\zeta=1/\sqrt{2}$, i.e. it is maximally flat at $\omega=0$, a Bessel low pass is even smoother and has larger damping ($\zeta=\sqrt{3}/2$), and a Chebyshev low pass filter is less damped than both of these filters. All of these filters satisfy $\zeta<1$, i.e. they are under-damped.


how to apply Short-time fourier transform to code in matlab Please Help, I have a white noise signal that I created and need help to apply the Short-time fourier transform to my code so that it can do the bandpass filtering. I am trying to put it in the code so that I don't have to use the FDA tool. Also be able to plot the graphs from the output of the STFT being applied.

L = 5000; %Sample length for the random signal
Pause = 10000; %Sample Pause Gap
mu = 0;
sigma = 2;

%Need to see left signal is not displaying
Left_signal = sigma*randn(L,1) + mu;
Right_signal = sigma*randn(L,1) + mu;

Long_signal = [Left_signal zeros(L,1); zeros(Pause,2); zeros(L,1) Right_signal];

%Plots subplots in graph
figure
subplot(211);
plot(Left_signal, 'b'); grid on;
subplot(212);
plot(Right_signal, 'r'); grid on;

this is very simplistic, but it

inputs a wav file and uses only one channel of it
breaks it up into windowed frames,
transforms it to the frequency domain,
does nothing to it in the frequency domain,
transforms back to the time domain,
overlap-adds that frame back into a time-domain array, and
eventually writes it back out to a different .wave file.

try it out.  can't guarantee there ain't typos.

if ~exist('inFile')
    inFile = 'input.wav';
end

if ~exist('outFile')
    outFile = 'output.wav';
end

if ~exist('frameWidth')
    frameWidth = 4096;          % size of FFT frame, better be a power of 2
end
frameHop = frameWidth/2;

analWindow = hanning(frameWidth);

[inBuffer, Fs] = wavread(inFile);

x = [inBuffer(:,1); linspace(0, 0, frameWidth)'];                   % use left channel only, zeropad one frame at the end

clear inBuffer;

numSamples = length(x);

numFrames = floor((numSamples-frameWidth)/frameHop);

y = linspace(0, 0, numSamples)';

n = 0;                              % init sample pointer.  unlike MATLAB, i like counting from 0

for frameIndex = 1:numFrames

     xWindowed = x(n+1:n+frameWidth) .* analWindow;     % get and window the input audio frame

     X = fft(fftshift(xWindowed));              % do the FFT

     Y = X;                         % copy the input spectrum to output

                                % do whatever processing to Y that you like

     yWindowed = fftshift(real(ifft(Y)));           % convert back to time domain, toss the imaginary part

     y(n+1:n+frameWidth) = y(n+1:n+frameWidth) + yWindowed;

     n = n + frameHop;
end

wavwrite(y, Fs, outFile);


Steps for successful image recognition I have asked the same question in stackoverflow ,https://stackoverflow.com/questions/17794772/steps-for-successful-image-recognition/17840975?noredirect=1#17840975,  but since I am not getting considerable responses , I'm going to ask it again here .

In order to do successful image recognition what are the steps one has to follow ,for example denoising , use of LDA , PCA , then use of neural network . I am not looking for any algorithm , just a brief outline of each of the steps from capturing an image to test an input image for similarity.

First of all, your question is board and it is hard to be answered thoughtfully. In general, image recognition itself is a wide topic. It might refer to classify a given image into a topic, or to recognize faces, objects, or text information in an image.

Depending on the objective of image recognition, you may use completely different processing steps. However, you may write the following general steps:

Training

get training data set
preprocessing (e.g. denoising, normalization, scaling, contrast
equalization, binarization)
feature extraction (e.g. sift, hog, lbp, raw pixels, contours, edge
map, DCT, DFT)
machine learning (e.g. logistic regression, neural network, svm)
parameter tuning (need a dev. set)
check results on how good you are, and see whether you need more data. If so, go to step 1 and repeat the process again.

Recognition:

the same preprocessing in training
the same feature extraction in training
feed your features into the trained classifier
output classification results.

Note: depending on your objective, preprocessing and selected features can be very different.

Asking someone to mark data manually is expensive. So people often use free dataset, each of which is often created for a specific purpose. Thus, these dataset may or may not satisfy your needs. Meanwhile, images from different dataset are often inconsistent, e.g. one set maybe of small grayscale images, while another maybe of large color images.


How to extract vocal part from stereo audio signal? I'm now processing MP3 file and encounter this problem. My MP3 is stereo encoded. What I want to do is extract vocal part for further processing(whatever mode of output signals, mono or stereo are both OK).

As far as I know, audio is encoded into different dis-joint sub frequency bands in MP3. I think I can limit the signals to the vocal range through high-pass/low-pass filter with cutting-off frequency proper set. However, result must contain parts of pure music signal in this case.
Or after googling, I think I may calculate the background signals first(by inverting one channel adding with signals from the other channel assuming vocal part is centered in the stereo audio called phase cancellation). After this transformation, the signal is mono. Then I should merge the original stereo into mono from which extracting the background signal.

Given the effectiveness, which one is preferred(or any other solutions:)? If the 2nd one, let two channels A and B, will (B-A) or (A-B) used when compute the background? As with merging two channels, does the arithmetic mean accurate enough? Or I can downsample each channel by a factor of two and interleave the downsampled signals as mono result?

Thanks and best regards.

First of all, how the data is encoded in a mp3 file is irrelevant to the question unless you aim at doing compressed-domain processing (which would be quite foolish). So you can assume your algorithm will work with decompressed time-domain data.

The sum / difference is a very, very basic trick for vocal suppression (not extraction). It is based on the assumption that the vocals are mixed at the center of the stereo field, while other instruments are panned laterally. This is rarely true. L-R and R-L will sound the same (the human ear is insensitive to a global phase shift) and will give you a mono mix without the instruments mixed at the center. The problem is, once you have recovered the background, what will you do with it? Try to suppress it from the center (average) signal? This won't work, you will be doing (L + R) / 2 - (L - R), this is not very interesting... You can try any linear combinations of those (averaged and "center removed"), nothing will come out of it!

Regarding filtering approaches: the f0 of the voice rarely exceeds 1000 Hz but its harmonics can go over that. Removing the highest frequency will make consonants (especially sss, chhh) unpleasant. Some male voices go below 100 Hz. You can safely cut whatever is below 50 or 60 Hz (bass, kick), though

Some recent developments in voice separation worth exploring:

Jean Louis Durrieu's background NMF + harmonic comb > filter model. Python code here.
Rafii's background extraction approach. Straightforward to code and works well on computer-produced music with very repetitive patterns like Electro, Hip-hop...
Hsu's approached based on f0 detection, tracking and masking. "A Tandem Algorithm for Singing Pitch Extraction and Voice Separation from Music Accompaniment" (can't find accessible PDF).


Is convolution distributive over multiplication? Is there any formula or expansion for

$$a(t)*[b(t) \cdot c(t)]$$

$$a(t) \cdot[b(t)*c(t)]$$

where $*$ denotes the convolution?

By expansion I mean something like $a(t)\cdot[b(t)+c(t)]=a(t)b(t)+a(t)c(t)$.

Addition is (by nature) additive. Binary operations like multiply and convolution are both somehow of multiplicative structure. Thus, I suspect there is no generic distributivity (apart for Boole or Boolean algebras, where the restriction on allowed values shrinks the problem, see below).

However, there is a sort of commutativity of multiply and convolution, as follows. It is interesting to first look at such operator  properties from an algebraic point of view. Properties are for instance, for a generic operator binary $\bigcirc$:

commutativity: $a \bigcirc b = b \bigcirc a $
associativity: $a \bigcirc( b \bigcirc c) =  (a \bigcirc b) \bigcirc c  $
alternativity: $a \bigcirc( a \bigcirc b) =  (a \bigcirc a) \bigcirc b  $ or $(a \bigcirc  b) \bigcirc b =   a \bigcirc (b \bigcirc b)$
flexibility: $(a \bigcirc  b) \bigcirc a =   a \bigcirc (b \bigcirc a)$

Distributivity of a binary   operator $\diamond$  over operator $\bigcirc$  is more involved: it can be left-distributed: $$a \diamond (b \bigcirc c) = (a \diamond b) \bigcirc (a \diamond c)$$ or right-distributed: $$(a \bigcirc b) \diamond c = (a \diamond c)  \bigcirc (b\diamond c)$$

Division is right-distributive over add, not left-distributive. Distributivity of multiply over add is at play in rings and fields. But the converse does not apply. There are mathematical structures where
two binary operators distribute over each other, like Boolean or switching algebras. Apart from that, I cannot remember of standard cases of multiply-like $\ast$  or $\cdot$ that would distribute. However, I will details some specific cases, useful to signal processing, ending with the Hilbert transform and Bedrosian theorem.

In the case of multiplication and convolution, or $a$ scalar, this is associative, as :

$$a.[b(t)\cdot c(t)] = [a.b(t)]\cdot c(t)= b(t)\cdot [a.c(t)]$$
and
$$a.[b(t)\ast c(t)] = [a.b(t)]\ast c(t)= b(t)\ast [a.c(t)]$$

You see that you cannot directly"distribute it", as this would become "bilinear".
Convolution and point-wise product are linear. So, if you consider that $c(t)$ (or  $b(t)$ by symmetry), there exists a linear operator $\Lambda$ such that:

$$\Lambda(c(t)) = a(t)\ast[b(t)\cdot c(t)] $$

Whether $\Lambda$ has a simple closed form expression is a difficult matter. As a side note, where $\mathcal{H} $ denotes the Hilbert transform, the Bedrosian theorem (A product theorem for Hilbert transforms states that (under some conditions):

$$\mathcal{H}(a(t)e^{i\theta t}) = a(t)(\mathcal{H}(e^{i\theta t}))\,.$$

where $a(t)$ is an (low-pass) envelope, and $e^{i\theta t}$ a modulation. More generally, if a low-pass signal $x_\flat$ and a high-pass  $x_\sharp$ have non-overlapping spectra, then under Bedrosian theorem:

$$\mathcal{H}(x_\flat(t)\cdot  x_\sharp(t))=x_\flat(t)\cdot \mathcal{H}(x_\sharp(t))$$

Now, let us remind that the Hilbert transform can be seen as the convolution with the distribution ($\operatorname {p.v.}$ is the Cauchy principal value):

$$h(t) =\operatorname {p.v.} {\frac {1}{\pi t}}$$

Thus:

$$h(t) \ast (x_\flat(t)\cdot x_\sharp(t))  =  x_\flat(t) \cdot(h(t)\ast x_\sharp(t))$$

This is the only practical case I can remember of a switching of operators in signal processing, more a form of limited commutativity.


Can a Kalman Filter Be Applied Using Measurement-Space Dependent Sensors? I'm currently attempting to apply a Kalman filter to track the angular position, velocity, and acceleration of a bike wheel, and I'm having a lot of trouble, so I want to check if I'm even applying the right tool for the job.  As you will see, this is not a "classical" Kalman implementation, and I want to check if I can still use it.

I have a wheel with n reed switches evenly spaced around the circumference of the wheel.  The way that I have the filter set up is that every time one of the reed switches is triggered by a magnet sitting on the bike's fork, a new measurement is fed to the filter that has the position and velocity (z = [x_prev + wheel_circumference/n ; wheel_circumference/(n*dt)] where x_prev is the previous position, wheel_circumference is the circumference of the bike wheel, and dt is the time since the last measurement).  (I have tried this system using only position measurements, and I am still having trouble.)

Unlike the classic application of the Kalman filter, the measurements in this system occur periodically in position, not in time, i.e. the sensor(s) sends a measurement to the filter iff the wheel's position reaches a certain position.  Compared to a traditional sensor, this implies a heavy correlation between position, velocity, and inter-measurement period.  Does this correlation negatively affect the Kalman filter?

Yes, sorta kinda.  Basically, if you can model your system accurately and the resulting model fits the Kalman filter paradigm, then yes.

However, I think it's a mistake to give the filter position and velocity -- were I doing this, I'd just give it the wheel position and let the Kalman filter determine the velocity -- which it will do well to the extent that your model matches the actual plant.

If it'll work at all, it'll be because you take the typical Kalman system model, $\mathbf{x}_k = \mathbf{F}_k \mathbf{x}_{k-1} + \mathbf{B u}_k + \mathbf{w}_k$, and you (probably, unless you have some prior idea of the acceleration being applied) let $\mathbf{B} = \mathbf{0}$, let $\mathbf{w}_k$ reflect your maximum expected acceleration, and let $k$ be the triggering instance.

The first fly in this ointment is that you have to calculate $\mathbf{F}_k$ and $\mathbf{w}_k$ according to the amount of time elapsed since the last trigger instant.  So for $\mathbf{x} = \left [ \begin{matrix} \theta & \omega \end{matrix} \right ]^T$, where $\theta$ is angular position and $\omega$ is angular velocity, you'd get something like $$\mathbf{F}_k = \left [
\begin{matrix} 1 & t \\ 0 & 1
\end{matrix}\right ]$$ and $$\mathbf{w}_k = \left [
\begin{matrix} t^2 & t \\ t & 1
\end{matrix}\right ] \sigma_\alpha^2$$ where $\sigma_\alpha$ is your expected variation in acceleration.

(Note that I'm not giving the derivation of $\mathbf{w}_k$, for which I only sorta-kinda-apologize -- it's tedious for me to translate my hen-scratches to properly formatted math, but it'll be good for you to derive it yourself.  Start with the definition of the Kalman-Bucy filter and grind through the math to find $\mathbf{P}(t)$ for a constant $\mathbf{w}(t)$ -- the above expression for the effective sampled-time $\mathbf{w}_k$ should drop out).

Note that the second fly in the ointment is that as the speed of the wheel drops the effect of the acceleration uncertainty ($\sigma_\alpha$, above) will get bigger and bigger compared to the position change from sample to sample.  At the point where the likelihood of overrunning the sample gets significant, you can't really use a "pure" Kalman filter anymore -- logically, this is because if you haven't triggered a pulse, then you know darned well that your position and velocity excursion are limited, but there's no way to effectively communicate that to a pure Kalman application.

If you need to track a really slow wheel, then you need to come up with some other, nonlinear, state estimator.

PART II

On reflection, I thought of a possibly-better way.  It's still not the best way, because you're jamming a nonlinear system with non-Gaussian uncertainty into a model of a linear system with Gaussian noise -- but it has the advantage both of being conceptually simpler, and of working down to zero speed.

The system model that you use to build a Kalman filter doesn't have to be time-invariant.  So model the system as $$\begin{align}
\mathbf{x}_k = \mathbf{F} \mathbf{x}_{k-1} + \mathbf{B u}_k + \mathbf{w} \\
\mathbf{y}_k = \mathbf{H} \mathbf{x}_k + \mathbf{v}_k
\end{align}$$

Note that $\mathbf{v}_k$ is modeled as time-varying.  In truth, the system is nonlinear (you get a pulse at specific values of $\theta$), but we'll just pretend.

Run the Kalman at a constant rate in time, with $\mathbf{F}$, $\mathbf{B}$, and $\mathbf{w}$ $\mathbf{H} = \left [ \begin{matrix} 1 & 0 \end{matrix} \right ]$ all constant.

But for each sample where a pulse is triggered, set $\mathbf{y}$ equal to the known position of the wheel, and $\mathbf{v}_k$ small, to reflect that you darned well know where the wheel is at that moment (size $\mathbf{v}$ to reflect any uncertainties in position due to mounting or irregular triggering or whatever).

For each sample where a pulse is not triggered, set $\mathbf{v}_k$ large, to reflect the fact that you don't know where between the triggers the wheel may be.  I'd start with $\mathbf{v}_k = \left( 3 \theta_0 \right)^2$, with $\theta_0$ equal to the spacing between trigger points -- and note that the multiply by three is coming directly out of my ear, so feel free to experiment with it.  Set $\mathbf{y}$ equal exactly halfway between the last known position of the wheel and the expected next position of the wheel (I'm assuming you know the wheel direction).

The net effect of this should be that if the wheel is turning fairly quickly the filter will have a darned good idea of the position and velocity, with "coasting" between readings.  If the wheel is turning slowly, it should coast OK if it's going at a constant rate, it should update at least fairly well if it speeds up, and if it slows down or stops, the fact that you're telling the Kalman filter that it's between positions "somewhere" should cause the Kalman's velocity estimate to slowly decay to zero (which is your best guess, intuitively), while the position estimate slowly settles on halfway between triggers (which again, is intuitively your best guess).


auto-correlation of differentiator How to do this entirely in the time-domain, without using frequency domain?
Let's say I have a continuous-time system that is a differentiator, and X(t) is WSS:
$$x(t) \rightarrow \boxed{\frac{d}{dt}} \rightarrow y(t)$$
I can change the differentiator operator into convolution as follows:
$$X(t) \rightarrow \boxed{\frac{d}{dt} \delta(t)} \rightarrow Y(t)$$
$$h(t) = \frac{d}{dt}\delta(t)$$
$$Y(t) = h(t) * X(t)=\frac{d}{dt}x(t)$$
$$\begin{matrix}R_{XX}(\tau)\end{matrix}\longrightarrow \boxed{\begin{matrix}h(\tau)\end{matrix}}\longrightarrow \begin{matrix}R_{XY}(\tau)\end{matrix}\longrightarrow \boxed{\begin{matrix}h(-\tau)\end{matrix}}\longrightarrow \begin{matrix}R_{YY}(\tau)\end{matrix}$$

cross-correlation of X(t) and Y(t):
$$R_{XY}(\tau) = h(\tau)*R_{XX}(\tau)$$
$$R_{XY}(\tau) = \Big(\frac{d}{d\tau}\delta(\tau)\Big) * R_{XX}(\tau)$$
$$R_{XY}(\tau) = \frac{d}{d\tau} R_{XX}(\tau)$$

autocorrelation of Y(t):
$$R_{YY}(\tau) = h(-\tau)*R_{XY}(\tau)$$
$$R_{YY}(\tau) = h(-\tau)* \frac{d}{d\tau} R_{XX}(\tau)$$
$$h(\tau) = \frac{d}{dt}\delta(\tau)$$
$$h(-\tau) = \frac{d}{dt}\delta(-\tau)$$
$$\delta(\tau) = \delta(-\tau)$$
$$h(-\tau) = \frac{d}{dt}\delta(\tau)$$
$$R_{YY}(\tau) = \Big(\frac{d}{d\tau}\delta(\tau)\Big) * \Big(\frac{d}{d\tau} R_{XX}(\tau)\Big)$$
$$R_{YY}(\tau) = \frac{d^2}{d\tau^2} R_{XX}(\tau)$$
However, the book says the answer should be:
$$R_{YY}(\tau) = -\frac{d^2}{d\tau^2} R_{XX}(\tau)$$
Somewhere I lost a sign...
and the reason why this sign is important because if I take the Fourier of $R_{YY}(\tau)$, it needs this sign to make $S_{YY}(\omega)$ positive as is required for PSD.

The derivative of the Dirac delta impulse is written as $\delta'(\tau)$. This helps with notation because the mistake you made is to write $h(-\tau)=\frac{d}{d\tau}\delta(-\tau)$, which is not the case because $\delta(\tau)$ is an even (generalized) function, whereas the derivative operator $\delta'(\tau)$ is an odd (generalized) function:
$$\delta'(\tau)=-\delta'(-\tau)\tag{1}$$
Consequently, the auto-correlation function of the output of the differentiator is
$$\begin{align}R_{YY}(\tau)&=R_{XX}(\tau)\star\delta'(\tau)\star\delta'(-\tau)\\&=-R_{XX}(\tau)\star\delta'(\tau)\star\delta'(\tau)\\&=-\frac{d^2}{d\tau^2}R_{XX}(\tau)\tag{2}\end{align}$$


data rate in digital communications I am a beginner in digital communications and I have question for the data rate with CDMA. Let, the length of the spreading sequence is $L_c$, The modulation is BPSK and the symbol time is $T_s$.

To calculate the data rate, I make the following reasoning: I have $1$ bit for $1$ symbol. After, the spreading, the length is $L_c$.

Here my questions:

What is the duration of the spreading sequence?
On the internet, I found that the data rate for CDMA is $D=\frac{1}{L_cT_s}$. Why $L_cT_s$?
Many people calculate the data rate with the length of a frame, can you explain to me how to do that?

Spreading refers to spreading of the frequency spectrum of the signal, and is achieved by dividing up the symbol interval $T_s$ into many shorter intervals called chip intervals ($L_c$ chip intervals in your notation) and

transmitting successive bits of the spreading sequence ($L_c$ bits long, remember?) in these $L_c$ chip intervals if the data bit you want to transmit in the symbol interval under consideration happens to be a $0$
transmitting the complements of the successive bits of the spreading sequence in these $L_c$ chip intervals if the data bit you want to transmit in the symbol interval under consideration happens to be a $1$

Symbolically, if $L_c = 3$ and the spreading sequence is $011$, then in a $T_s$-second symbol interval, you get to transmit either $011$ or $100$ in that $T_s$-second interval which we have mentally divided into $3$ subintervals of length $T_s/3$ each. So, the bits transmitted during each of these chip intervals are of shorter duration (by a factor of $3$) and thus require more bandwidth (by a factor of 3 too!) to transmit than we would need if we were transmitting just one data bit every $T_s$ seconds.

"But, but, but,..." you splutter, "I don't have data bits, but rather these beautifully shaped RRC pulses that are of duration $T_s$ seconds. How do I create a CDMA signal from these?"  Well, you have to go back to the purely digital domain where bits are bits, use a frequency multiplier on your clock to get clock pulses at a rate of $L_c/T_s$ Hz so that you can replace each data bit at the $1/T_s$ Hz clock rate with $L_c$ chip bits at the $L_c/T_s$ Hz   clock rate, and tell the analog guys to design their circuits to produce RRC pulses of duration $T_s/L_c$ to accommodate your higher rate digital input signal that is going to send them bits at rate $L_c/T_s$ Hz.


How to calculate SNR of dead-band quantizer? Suppose we have a dead-band quantizer with $\ N = 5 \ quantizing \ levels $  which uses $\ b = 3 \   bits $ for each produced symbol. The quantizer gets as input a uniformly distributed signal $\ X $ at the interval $\ [-10,10] $. How can I calculate the signal to noise ratio $\ SNR $ of the quantizer ?

Time align the quantized signal with your input signal and subtract to create an error signal (match the amplitudes such as to minimize the error signal). The SNR can then be determined as the ratio of the standard deviation of the input signal (after amplitude matching) to the standard deviation of the error signal (20Log10 of this ratio to be in dB).

The top diagram that I have shows this for a complex signal example, but it is also equally applicable to real signals.


Tutorial for a stereo Phaser in C++? I was searching for a tutorial on how to implement a stereo Phaser in C++. All i´ve found was sourcecode with little or no explanations which did´nt help me much. Does anyone know a Tutorial with which i can learn and understand how to implement a stereo Phaser in C++? By the way im using the Juce Plugin.

Have you seen this MusicDSP Phaser code and discussion by Ross Bencina? It is commented, and says that its implemented to be "clear, not efficient". Cheers.


Recursive Implementation of the Gaussian Filter I'm trying to implement an IIR form to approximate the Gaussian Blur Filter.
I'm working with the article "Recursive Implementation of the Gaussian Filter" by Ian T. Young and Lucas J. van Vliet.

They suggest a form and way to calculate the coefficients as given by:

I'm trying to reproduce their example for $ q = 5.0 $.

There is the MATLAB code I wrote:

qFactor = 5;

b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;

normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);

vDenCoeff = [b0Coeff, b1Coeff, b2Coeff, b3Coeff] / b0Coeff;

vXSignal = zeros(61, 1);
vXSignal(31) = 10;

vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));

figure();
plot(vYSignal);

Now, the result I get is this:

Namely, I get a filter which isn't stable.
Yet it seems I get the same coefficients as they get in their example.

What am I missing?
Has anyone ever implemented this method?

Thank You.

The answer was simple, the article uses the coefficients value on one hand where the MATLAB implementation on the other.
Namely, a minus sign should be added.

Here's the correct code:

qFactor = 5;

b0Coeff = 1.57825 + (2.44413 * qFactor) + (1.4281 * qFactor * qFactor) + (0.422205 * qFactor * qFactor * qFactor);
b1Coeff = (2.44413 * qFactor) + (2.85619 * qFactor * qFactor) + (1.26661 * qFactor * qFactor * qFactor);
b2Coeff = (-1.4281 * qFactor * qFactor) + (-1.26661 * qFactor * qFactor * qFactor);
b3Coeff = 0.422205 * qFactor * qFactor * qFactor;

normalizationCoeff = 1 - ((b1Coeff + b2Coeff + b3Coeff) / b0Coeff);

vDenCoeff = [b0Coeff, -b1Coeff, -b2Coeff, -b3Coeff] / b0Coeff;

vXSignal = zeros(61, 1);
vXSignal(31) = 10;

vYSignal = filter(normalizationCoeff, vDenCoeff, vXSignal);
vYSignal = filter(normalizationCoeff, vDenCoeff, vYSignal(end:-1:1));

figure();
plot(vYSignal);


FFT averaging with different block sizes - how to scale the amplitudes I need to analyze a real-valued time signal with a length of 300000 in the frequency domain. The sampling frequency is 100000Hz. In order to increase the signal-to-noise ratio  reduce the "fine-grain" fluctuation when your number of bins is too large I decided to split the whole signal in smaller blocks, calculate the FFT of the smaller blocks and average the spectra. I used Gnu Octave to implement this:

function A = fft_average(x, win_length, samp_freq)

%overlapping samples
step = win_length / 2;

% calculate frequency range
freq = 0:1/win_length*samp_freq:(win_length-1)/win_length*samp_freq;

%fft_win = ones(win_length,1);
%fft_win = hanning(win_length);
%fft_win = hamming(win_length);
fft_win = blackman(win_length);

y = zeros(win_length,1);
start = 1;
stop = win_length;
nsteps = 0;
while (stop <= length(x))
    nsteps = nsteps + 1;
    ytemp = abs(fft(x(start:stop).*fft_win));
    start = start + step;
    stop = stop + step;
    % add (and scale)
    y = y + ytemp;
end
% average
y = y./(nsteps);

%copy to result array, omit negative frequencies
A(:,1) = freq(1:win_length/2);
A(:,2) = y(1:win_length/2);

endfunction

Here are the results with 3 different block sizes of 1024, 4096 and 16384.
Additionally, I tried 4 different window functions. I only post results for rectangle and Hann Window due to lack of reputation.

I did not apply any scaling to the amplitudes yet. What I would expect from a correct scaling is consistency of the results for different block sizes. Hence no offset like in the second image for the Hann window.
If I scaled the amplitudes by the square root of the block size, the results for the Hann-windowed FFTs match. But I have 3 problems with this approach:

What would be the mathematical foundation to scale by the square root of the block size? Or is it the square root of the block size because I made an obvious mistake?
The amplitudes for the rectangle window do not match if I apply this scaling. Is this rather caused by spectral leakage?
Do I have to apply an additional scaling to account for the shape of the window function, i.e. so the integral remains the same?

I do not immediately see any SNR advantage to splitting and averaging assuming you are dealing with additive white noise. I will explain why in case there are other possible approaches- but before detailing in simple terms, the SNR for the case of additive white noise is directly related to the bin width in frequency, and the bin width in frequency is inversely related to the length. For the rectangular window, the equivalent noise bandwidth is 1 bin wide. If we halve the number of bins, the bin width will double and therefore noise will double (+3 dB). If we average the result of two bins, the noise will go down by 1/2 relative to the signal (-3dB). So we ended up doing a lot of work to get back to where we started.
Before going into further details on that, to answer your scaling question the FFT scaling in magnitude is the length of the sequence (So $1/\sqrt{N}$ in power).
And yes you must apply additional scaling to account for the window shape. See one of my favorite papers which details this for all the common windows; harris on Windows, and I draw attention to specifically the equivalent noise bandwidth and coherent gain metrics for each window.  Regarding coherent gain, the mathematical basis for this is simply the sum of the window compared to the rectangular window. For Hann, this comes to 0.5 as the length of your window approaches infinity, and 0.49902 for a length of 512, so basically 0.5 overall. So to normalize in this case you would scale in magnitude by 2/N, compared to the rectangular window where you would scale by 1/N.  (This window adjustment factor changes with an overlap-add approach based on the overlap, which is further detailed in fred's paper).
More Details:
The FFT is an algorithm that efficiently calculates the DFT with no difference in the result, so reviewing the DFT formula can be very insightful to the scaling as well as the SNR effects:
$$X[k]=\sum_0^{N-1}x[n]e^{-j (k\omega_o) n}$$
Notice what is going on here: In doing this computation, we are summing over N samples of $x[n]$ for each value of k as $e^{-j (k \omega_o)n}$ counts from 0 to N-1. This is a correlation of x[n] to $e^{+j (k \omega_o)n}$, what I view as a "spinning phasor"  ($e^{j\phi n}$ is a vector with magnitude 1 and phase $\phi n$, and n is increasing with time). Note that any "signal" components that are correlated to the spinning phasor (meaning spinning in the equal but opposite direction), will grow in magnitude at rate N. Consider our "DC" bin 0 as this is the simplest example: If x[n] is constant as M (and $e^{-j0} = 1$), the output from the summation will be NM. For any other frequency component in x[n] the same result occurs, as our spinning phasor essentially moves that signal component to DC by spinning in the opposite direction prior to summing.
Thus we see all correlated signals grow in magnitude (not power!) at rate N; since it is a magnitude quantity they grow at rate $20Log_{10}N$.  Noise on the other hand, if uncorrelated from sample to sample, will sum in power, so the noise component grows at rate $10Log_{10}N$. By doing this "correlation" we get a net processing gain (SNR increase) of $10Log_{10}N$, as the "signal" component went up by 20dB per N while the noise component went up by 10dB per N. In fact, we see from this exactly why "correlation" works in general, which is the process of multiplying and accumulating (or multiply and integrate in the analog domain as is done with the Fourier and Laplace transforms)! By correlating, to the extent that the noise from sample to sample is not correlated, we get a $10Log_{10}N$ advantage.
NOTE: Additional updates to this response are required to be complete in order to show two different cases: Averaging the complex FFT values directly (resulting in scalloping through the spectrum from the delay+add) vs averaging the absolute magnitudes as the OP has later indicated he has actually done.


simulating directional interference to evaluate model performance I have a model for audio processing. The model receives as input, a source, convolved with a room impulse response (RIR). I would like to test the performance of the model against inputs containing directional interference. The sources and interferences are simply audio recordings (for example human speech, a dig bark, car engine, etc.). both source and interference are convolved with two different RIRs (with respect to angle) and then a gain $A$ is applied to the noise in order to control the desired SIR.

The signals power: $P_s(t)=mean(x^2(t))$

The interference (noise) power: $P_n(t)=mean(n^2(t))$

The amplified interference power: $P_{An}(t)=mean(A^2n^2(t))=A^2mean(n^2(t))=A^2P_n(t)$

And if we want to isolate $A$, we have:
$$SIR=10 \log\left(\frac{P_s(t)}{A^2P_N(t)}\right)$$
$$A=\sqrt{\frac{P_s(t)}{10^{\frac{SNR}{10}}P_N(t)}}$$

Now, considering I am working from a dataset of audio files with variable length, will this not cause a problem due to the averaging process? For example, if the interference is of length $T_1$ and the source is of length $T_2>T_1$. The interference will be of worst than the desired SIR in the interval $0<t<T_1$ and better in the interval $T_1<t<T_2$.

Is my thought process correct?
Is there a known method of dealing with such a problem?
How do I set the correct SNR for such an experiment?

You do need to be careful how exactly you define the SNR. I think you are actually talking about SIR (signal to interference ratio). The way you describe it, it sounds like you plan to just compute the power (time averaged energy) emitted by the signal and interference and use that to get after the SIR. This would be fine but, as you mentioned, sometimes the interference is not present. I might think about using a different definition such as instantaneous SIR, since it sounds like sometimes you have signal + interference and sometimes it is just signal. For this reason, I'm not convinced that computing the power as you described is the right thing to do.

I would suggest doing this: take the source audio file and interference audio file, figure out which is longer, zero pad the shorter one so they are of equal length. Now compute the metric, "sample-by-sample SIR", by doing $\frac{s[n]^2}{i[n]^2}$. Just be careful when $i[n]=0$. This is where knowledge of the noise (real noise, not interference) would be helpful so that you could ideally report SINR (signal to interference plus noise ratio) and avoid the denominator being $0$ problem.


i have an audio signal. i have to find whether my audio signal is from given avi scene I have audio signal. I want to find whether my audio signal is from a given avi file . I tried using correlation and mfcc. But i'm not getting accurate output. So please help me in this.

The information about your problem is quite incomplete. If you want to find where is small signal is in large signal, then correlation give a best results usually.
If your problem is more complicated you can see on very popular application Shazam for music search by audio signal. Googling with string "shazam algorithm" give you a lot of information.  See for example article of author - An Industrial-Strength Audio Search Algorithm.

In brief - you create "fingerprint" of signal by using frequency spectrum of shorts chunks of signal. Fingerprints are used for searching in catalog of songs. You can use such method for your problem.


FFT on non-periodic signal and signal power at each harmonic - is there a more robust approach? I have a non-periodic signal collected from a force plate, it represents the foot/ground contact of a human when running.

My aim is to calculate the signal power of each harmonic up to 50Hz. I will break down what I have done so far.

Firstly, I calculate the fundamental frequency by dividing the sampling frequency by the length of the signal. I then find the other harmonics by multiplying the fundamental frequency by 2,3...etc (not shown here, this example just includes the fundamental frequency (FF).

Fs = 1000;
FF = Fs/(length(Fz_data));

I run the FFT on my signal.

dt = 1/Fs;
t=(0:1/Fs:(length(Fz_data)-1)/Fs)';
NFFT = length(Fz_data);
Y = fft(Fz_data,NFFT);

For each harmonic (identified in 1 - FF in this instance),I zero out frequencies above and below the harmonic and then run ifft to convert the data back to the time domain which leaves me with the data for the first harmonic only(?).

F = ((0:1/NFFT:1-1/NFFT)*Fs).';
Y(F<FF-0.1) = 0;
Y(F>FF+0.1) = 0;
y = ifft(Y,NFFT,'symmetric');

I then calculate the power of this signal using pwelch.

[P,~] = pwelch(y,ones(NFFT,1),0,NFFT,Fs,'power');
pwr = sum(P)

Questions

From a programming perspective, this works, but I have reservations that this approach is not based on sound theory. I will, therefore, present some questions below, and I hope the community may be able to offer some solutions to improve the robustness of this approach.

Is it suitable to calculate the fundamental frequency from the length of the signal?
I have read that it is not suitable to use an FFT on a non-periodic signal. Is there a way to adapt the FFT method to make it more suitable e.g. some sort of windowing function?
Is there an easy way to remove drift/bias from the signal?
I would also welcome comments on any other aspects of this process if you feel it would make it more robust.

Thanks in advance.

I'll try to address your questions one at a time.

  Is it suitable to calculate the fundamental frequency from the length of the signal?

Your signal is not periodic, so speaking of "fundamental frequency" is not correct. The fundamental frequency corresponds to the lowest frequency of a periodic waveform, which is not the case here. So, regarding your question, no: you cannot define a fundamental frequency from the length of a signal - two things not related at all.

  I have read that it is not suitable to use an FFT on a non-periodic signal. Is there a way to adapt the FFT method to make it more suitable e.g. some sort of windowing function?

FFT can indeed be used for non-periodic signals. In case you are not too familiar with DFT, I'll try to explain how it works in a nutshell.

You have a discrete signal. You window it (you only take $N$ of its values). You perform the DFT on those values. What does "perform the DFT" mean? The signal whose DFT will be calculated is assumed to be $N$-periodic, and you get a discrete frequency representation of that signal (not the original non-windowed one). I believe that this is the reason why you state that FFT is not suitable on non-periodic signals. Nevertheless, you can still use it, we do it all the time. Maybe some answers in this question can help.

If you are interested in which frequencies are present at each instant of time, one can use several windows an overlap them, and take the FFTs of those signals. See STFT.

  Is there an easy way to remove drift/bias from the signal?

Do you mean DC value? In that case, you can just substract the mean value of the signal, or apply a high-pass filter with sufficiently low cutoff frequency.

  I would also welcome comments on any other aspects of this process if you feel it would make it more robust.

At some point in your script, you just set some FFT coefficients to zero. Note that this approach is generally not a good idea.


Prove the dirac delta contains all frequencies I'm looking for a mathematical proof that the dirac delta contains all frequencies.

I just read in a text book that the frequency spectrum of a dirac is just a horizontal line of amplitude 1, whereas some things I've read online show the Fourier transform should look more like a complex exponential.

Either way, would someone be able to clear up for me the way in which the dirac has a broad band frequency spectrum such as white noise?

some things I've read online show the Fourier transform should look more like a complex exponential.

Don't believe everything you read online!

A signal "contains all frequencies" is a vague description; it can apply to any finite-energy signal $x(t)$ whose Fourier transform $X(f)$ is nonzero for all values of $f$, $-\infty < f < \infty$.  On the other hand, a finite-energy signal necessarily has the property that $|x(t)| \to 0$ as $t \to \infty$ as well as when $t \to -\infty$. Absent this property, $\int_{-\infty}^\infty |x(t)|^2 \, \mathrm dt$ will not converge and so $x(t)$ could not be called a finite-energy signal.

Dirac deltas are not finite-energy signals, but it is common engineering practice to treat them as if they are, and apply all the results in Fourier theory them. That this is legitimate requires more high-power math than I understand and so I, like most other engineers, simply accept the more delicate analyses as being understood. With that, consider that the Fourier transform of $x(t) = \delta(t)$ is $X(f) = 1, -\infty < f < \infty$ and so the Fourier transform is nonzero at all $f, -\infty < f < \infty$, that is,
the Dirac delta "contains all frequencies" as you would like it to. If you would like more generality, then the Fourier transform of $x(t) = \delta(t-\tau)$ is
$$X(f) = \int_{-\infty}^\infty x(t)\exp(-j2\pi ft) \, \mathrm dt = \exp(-j2\pi \tau f)$$
which is a function of $f$ (not $\tau$, that is just a constant), is a complex exponential function of $f$ as you have read online, and is nonzero for all $f, -\infty < f < \infty$, that is, $\delta(t-\tau)$ also "contains all frequencies" just like $\delta(t)$ does.  Exercise: What is the complex exponential in the case $\tau = 0$?

  Either way, would someone be able to clear up for me the way in which the dirac has a broad band frequency spectrum such as white noise?

The definition of white noise is a zero-mean process $\{X(t)\colon -\infty < t < \infty\}$ whose power spectral density is
$$S_X(f) = \frac{N_0}{2}, -\infty < f < \infty$$
which is not just nonzero but also constant for all $f$. That is, not only does the noise "contain all frequencies" but it contains all frequencies in "equal quantities". If the white noise is filtered through an LTI system with transfer function $H(f)$, then the power spectral density of the output is $\frac{N_0}{2}|H(f)|^2$ and the output noise power is $\frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$. In particular, if the filter is an ideal bandpass filter with bandwidth $B$, then the output noise power is $N_0B$ regardless of where the band is in the frequency spectrum.

In this sense,
the Dirac delta is like white noise: its spectrum $X(f)$ has constant value $1$ for all frequencies (when $\tau = 0$) and more generally, $|X(f)|$ has constant value $1$. The energy spectral density $|X(f)|^2$ of the Dirac delta also has value $1$ for frequencies $f$ (which is similar to the result for the power spectral density of white noise). Also, if the Dirac delta is filtered through an LTI system with transfer function $H(f)$, then the energy of the output signal is $\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$ which is similar to result for output noise power of the LTI filter when the input is white noise. Finally, if the LTI system is an ideal bandpass filter of bandwidth $B$, then the output is a finite-energy signal of energy $2B$, regardless of what the center frequency of the filter is, which is similar to the comparable result for output noise power from the LTI filter driven by white noise.


Spectral density of signal written as quadratures Suppose I have a real-valued random signal $f(t)$ written as
$$f(t) = I(t)\cos(\Omega t) - Q(t)\sin(\Omega t)\, .$$
What is the relationship between the spectral densities of $I$ and $Q$ and the spectral density of $f$?

One approach is to use the fact that
$$\int_0^\infty \frac{d\Omega}{2\pi} S_f(\Omega) \cos(\Omega \tau)=
\langle f(t)f(t+\tau)\rangle_t \, .$$
But so far I haven't understood how to make progress.

This answer is related to my answer to this question of yours. The short answer is that in general the random process $f(t)$ has no power spectral density because it is not wide sense stationary (WSS), even if $I(t)$ and $Q(t)$ are.

I will first introduce some notation. I use the following definition of the autocorrelation function of a (possibly complex) random process $X(t)$:

$$R_X(t,\tau)=E\{X(t+\tau)X^*(t)\}$$

where $E\{\cdot\}$ denotes the expectation operator, and $^*$ denotes complex conjugation. Note that this definition with the complex conjugate is consistent with the usual definition of the inner product of complex functions.

The cross-correlation of two (possible complex) random processes $X(t)$ and $Y(t)$ is defined by

$$R_{XY}(t,\tau)=E\{X(t+\tau)Y^*(t)\}$$

Consequently, the cross-correlation of $X(t)$ and its complex conjugate $X^*(t)$ is given by

$$R_{XX^*}=E\{X(t+\tau)X(t)\}$$

If $X(t)$ is wide-sense stationary (WSS), then its autocorrelation function only depends on the time difference $\tau$ and is written as $R_X(\tau)$. Similarly, if two processes $X(t)$ and $Y(t)$ are jointly WSS, their cross-correlation function only depends on the time difference $\tau$ and is written as $R_{XY}(\tau)$.

For WSS processes, the power spectrum is defined as the Fourier transform of their autocorrelation function:

$$S_X(\omega)=\int_{-\infty}^{\infty}R_X(\tau)e^{-j\omega\tau}d\tau$$

Similarly, for two jointly WSS processes, their cross-spectral density is given by the Fourier transform of their cross-correlation function:

$$S_{XY}(\omega)=\int_{-\infty}^{\infty}R_{XY}(\tau)e^{-j\omega\tau}d\tau$$

Now for the answer of the question. To see that $f(t)$ is generally not WSS, define a complex-valued WSS process $X(t)$ by

$$X(t)=I(t)+jQ(t)\tag{1}$$

with real-valued and jointly WSS processes $I(t)$ and $Q(t)$.
The process $f(t)$ is then given by

$$f(t)=\Re\{X(t)e^{j\Omega t}\}=\frac12(X(t)e^{j\Omega t}+X^*(t)e^{-j\Omega t})\tag{2}$$

The auto-correlation function of the real-valued process $f(t)$ is

$$R_f(t,\tau)=E\{f(t+\tau)f(t)\}\tag{3}$$

Plugging (2) into (3) gives (after some straightforward manipulations)

$$R_f(t,\tau)=\frac12\Re\{R_X(\tau)e^{j\Omega\tau}\}+
\frac12\Re\{R_{XX^*}(\tau)e^{j\Omega(\tau+2t)}\}\tag{4}$$

From (4) it is clear that the auto-correlation function of $f(t)$ depends on $t$, and, consequently, the power spectral density of $f(t)$ is not defined. However, if $R_{XX^*}(\tau)=0$ and if $E\{X(t)\}=0$ then $f(t)$ is WSS and has a power spectral density.

The cross-correlation of $X(t)$ and $X^*(t)$ is given by

$$\begin{align}R_{XX^*}(\tau)&=E\{X(t+\tau)X(t)\}\\&=E\{(I(t+\tau)+jQ(t+\tau))(I(t)+jQ(t))\}\\&=\ldots\\&=R_I(\tau)-R_Q(\tau)+j(R_{IQ}(\tau)+R_{IQ}(-\tau))\end{align}\tag{5}$$

where $R_I(\tau)$ and $R_Q(\tau)$ are the auto-correlation functions of $I(t)$ and $Q(t)$, respectively, and $R_{IQ}(\tau)$ is the cross-correlation of $I(t)$ and $Q(t)$. From (5) it is obvious that the condition $R_{XX^*}(\tau)=0$ is satisfied if

$$\begin{align}R_I(\tau)&=R_Q(\tau)\\
R_{IQ}(\tau)&=-R_{IQ}(-\tau)\end{align}\tag{6}$$

is satisfied, i.e. if $I(t)$ and $Q(t)$ have the same auto-correlation function (and thus the same power spectrum), and if $R_{IQ}(\tau)$ is an odd function, the latter implying $R_{IQ}(0)=0$, i.e. $I(t)$ and $Q(t)$ are uncorrelated when sampled at the same instant.

The auto-correlation function $R_X(\tau)$ can be expressed in terms of $R_I(\tau)$, $R_Q(\tau)$, and $R_{IQ}(\tau)$:

$$\begin{align}R_X(\tau)&=E\{X(t+\tau)X^*(t)\}\\&=\ldots\\&=R_I(\tau)+R_Q(\tau)-j(R_{IQ}(\tau)-R_{IQ}(-\tau))\end{align}\tag{7}$$

If the conditions (6) are satisfied, this becomes

$$R_X(\tau)=2(R_I(\tau)-jR_{IQ}(\tau))\tag{8}$$

Pluggin (8) into (4) gives for the auto-correlation function of $f(t)$

$$R_f(\tau)=R_I(\tau)\cos(\Omega\tau)+R_{IQ}\sin(\Omega\tau)\tag{9}$$

(assuming the conditions (6) are satisfied). From (9) the power spectral density of $f(t)$ can be expressed by the power spectral densities of $I(t)$ and $Q(t)$, and by the cross-power spectral density of $I(t)$ and $Q(t)$:

$$S_f(\omega)=\frac12[S_I(\omega-\Omega)+S_I(\omega+\Omega)]+\frac{1}{2j}[S_{IQ}(\omega-\Omega)-S_{IQ}(\omega+\Omega)]\tag{10}$$

which is the final result. Note that $S_f(\omega)$ is of course real-valued because due to the second condition in (6), the cross-power spectral density $S_{IQ}(\omega)$ is purely imaginary and, consequently, the second term on the right-hand side of (10) is real-valued.


What is "feature indexing" in image processing? What does "feature indexing" mean in image processing? does it mean that we assign a number to each row/column of our feature vectors? any precise explain would appreciate.

The idea is that each image generates a series of features, and you want to be able to search all images based on these features.  That way, features in one image can be quickly associated with close features in another image.

For example, this slide deck has decided to use four "visual words" to say that each block of images in the figure are "the same" in some sense.

I'm not quite sure what you mean by:

  we assign a number to each row/column of our feature vectors

Each figure vector is a series of numbers.  What that slide deck does is assigns an identifier (visual word) to each block of similar images (features).

What indexing then does is finds the "middle" of each visual word (group of similar figures) so that other images can be indexed against these words. The green dots in the left-hand of the diagram are the middles for three of the words.