Calculating power and energy of a signal in matlab I am a little new to matlab and would like to get some help on how to write a program or a function that can calculate the energy/power of a signal depending on whether the signal is an energy or a power signal.
Here is the program the I tried:

function [ ] = power_or_energy( s,t)
%UNTITLED Summary of this function goes here
%   s is the input signal
%   t is the independent variable

s2=s.^2;
energy=int(s2,t,-inf,inf);

if energy~=inf && energy~=-inf
    disp ('the signal is an energy signal with energy equal to ' );
    disp(energy);
end

power=limit((int(s2,t,-t/2,t/2))/t,t,inf);

if power~=inf && power~=-inf
    disp('the signal is a power signal with power equal to ');
    disp(power);
end

But the problem here is that when i want to check it for a function like

t=0:0.01:10;
s=sin(t);
power_or_energy(s,t);

then i get the following error:
"Undefined function 'int' for input arguments of
type 'double'."

So how to i actually calculate the energy for such signals?

Also the above function works fine with the below input:

syms t;
s=sin(t);
power_or_energy(s,t);

Please help me in writing a function that can be used to find the energy of any input signal.

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

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

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

Refer this link for calculation of Energy and Power of signal: http://www.mathworks.com/matlabcentral/fileexchange/11202-automatic-gain-control


Fourier Series of Aperiodic convolution of periodic functions we were given the following classic exercise:

  Given two periodic signals $x(t), y(t)$ with fundamental period $T$
  with Fourier series coefficients $c_m^x, c_m^y$ respectively, find the
  Fourier coefficients of the signal $z(t) = x(t) * y(t)$ with relation
  to $T, c_m^x, c_m^y$.

Now, this can easily be solved when the aforementioned convolution is the circular convolution (integral over a period only).

However, in class our professor noted that it can be solved even when we have an aperiodic convolution (that is, convolution as an integral from $-\infty$ to $+\infty$). We argued that, in this case, that infinite integral doesn't converge, and he responded that, even though the convolution integral doesn't converge (i.e. might be infinite), the Fourier Series coefficients are still finite and can be calculated!!

Is this true? If yes, then is the relation the usual one: $c_m^z=Tc_m^xc_m^y$ or another and how do you prove that? If not, why? Austere mathematical proofs would be appreciated.

okay, now that the popcorn is eaten and my fingers are clean...

$$ x(t) = \sum\limits_{m=-\infty}^{+\infty} c_m^x e^{j 2 \pi (m/T) t} $$

$$ y(t) = \sum\limits_{m=-\infty}^{+\infty} c_m^y e^{j 2 \pi (m/T) t} $$

then

$$\begin{align}
 z(t) &= \sum\limits_{m=-\infty}^{+\infty} c_m^z e^{j 2 \pi (m/T) t} \\
\\
  &= x(t) \circledast y(t) \\
\\
  &= \int\limits_{-\infty}^{+\infty} x(\tau) y(t-\tau) \, d\tau \\
\\
  &= \int\limits_{-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x e^{j 2 \pi (m/T) \tau} \sum\limits_{n=-\infty}^{+\infty} c_n^y e^{j 2 \pi (n/T) (t-\tau)}  \, d\tau \\
\\
  &= \sum\limits_{n=-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x c_n^y  \int\limits_{-\infty}^{+\infty} e^{j 2 \pi (m/T) \tau} e^{j 2 \pi (n/T) (t-\tau)}  \, d\tau \\
\\
  &= \sum\limits_{n=-\infty}^{+\infty} \sum\limits_{m=-\infty}^{+\infty} c_m^x c_n^y \int\limits_{-\infty}^{+\infty} e^{j 2 \pi ((m-n)/T) \tau} \, d\tau \ e^{j 2 \pi (n/T) t} \\
\end{align} $$

this integral: $\int\limits_{-\infty}^{+\infty} e^{j 2 \pi ((m-n)/T) \tau} \, d\tau$ ain't converging when $n \ne m$ and certainly not when $n=m$.

you will not be getting finite values for $c_m^z$.


How can I identify Arabic letters in an image using Java? I want to identify Egyptian car number plates that look like this:

How can I do this in Java ?

I plan to do this in Android, as the user will catch the number plate with a camera and the results will be submitted to server.

You have to do following methods:

First step: Number plate extraction

You can use OpenCV for this purpose. It has Java api. Also it is ported to Android.

Second step : OCR

OCR is used to recognise the letters on notice board.

You can train your own OCR for arabic letters and digits using OpenCV.

Or you can use tesseract-OCR. I think it supports arabic also.

Following paper will be useful for you: Paper 1

Also, this site includes complete code for number plate extraction, but in C++. But it will be a good start.

The link provides an algorithm to detect text from an image in Java... it works for English..

You can try for something similar for Arabic.
http://code.google.com/p/textdetection/source/browse/trunk/src/com/jpprade/text/detection/MainApp.java?r=7


Decompose video frames into slices I am working on a video motion feature implementation. I am implementing a work in a research paper. What I want to implement is this:

  The frames contained in a video clip can be thought of as a volume obtained by considering all the frames in time. This volume can be decomposed into a set of two 2-D temporal slices I(x,t)and I(y,t).

What's the methodology of implementing this i.e slicing the volume of frames in time?
Relevant link of paper: http://www.cs.cmu.edu/~yaser/RasheedSheikhShah_CVST_2005.pdf

Each frame of your movie is a 2D plane indexed by $x$ and $y$.
When you consider several consecutive frames, you can stack them in order to have a 3D cube, where each slice of the cube is a frame.
This space-time volume is now indexed by $x,y$ and $t$.

You can manipulate this volume as any 3D volume, just replace i nyour mind $t$ by $z$ if you find it confusing.

Once you have this volume, if you change the normal axis used to define the 2D slices (remember it was the $t$ axis initially), you obtain the different slices that you want:

choosing the $x$ axis as normal, you obtain 2D $(y,t)$ slices
choosing the $y$ axis as normal, you obtain 2D $(x,t)$ slices

By the way, the classical term for these kind of slices is epipolar-plane image analysis if you are loking for more references.

As the different type of videos have different size, you should first specify your video format. These videos formats use different types of compression which makes you decode them before stating the processing. So, researchers usually use uncompressed (raw) format of videos named YUV. If you are using this format I can give you the m files which extracts each frame of the video sequence as a matrix in a cell. You should create a cell which its elements are the video frames.


Criteria to distinguish FIR and IIR filters from difference equation? What criteria should be used to safely decide if it is an IIR filter or FIR filter from a difference equation? FIR filter are always stable, meaning all poles are inside the unit circle AND have pole-zero cancellations i.e have equal number of poles and zeros at the same point in the z-plane. IIR filters are conditionally stable, i.e always have poles on/outside the unit circle. Is this good enough to distinguish from the difference equation? what other criteria should be considered?

PS: Please don't say IIR filters have recursive difference equations, and FIR doesn't. It is not true.

If the Z-transform of the feedforward section is divisible by the Z-transform of the feedback section, the filter is FIR.

Consider your example: $y[n] = y[n-1] + x[n] - x[n-3]$.

The Z-transform is $\mathrm Y(z)- z^{-1}\mathrm Y(z) = \mathrm X(z) - z^{-3}\mathrm X(z)$, and the Z-transform of the response is $\mathrm H(z) = \mathrm Y(z)/\mathrm X(z) = (1 - z^{-3})/(1 - z^{-1}) = 1 + z^{-1} + z^{-2}$. This means that the filter can be realized as $\mathrm y[n] = \mathrm x[n] + \mathrm x[n-1] + \mathrm x[n-2]$, which is clearly FIR.

As far as i know, unless you want a deep philosophical analysis fo filters, there is absolutely no working distinction between a FIR/IIR filter or a diference equation, in the sense that a filter has "ALWAYS" an associated diference equation that you can write.

Having said that though, if you randomly pick up a piece of papaer that has a difference equation written on it, you can't automatically say "Oh they are working with filters and signal processing here", because difference equations are used in lots of other applications where you aren't even processing a signal. And the solution to the equation (giveng initial conditions or an input $X[n]$) might mean something completely different to the output of a filter.


Why does FFT output need to conjugate processing? I came across some matlab code as the below,

fftout   = fft(adc_in, 1024);
fftout_p = fftout.*conj(fftout);

I think fftout is formatted as complex form. from here I'm ok.

Next stage that signal is multiplicated by itself which is conjugated.

I think it means that it want to get the real number not complex number.

From here, I'm not sure and want to know that practically what does Imaginary part affect to real signal processing?

Can we ignore like above that?

In real practically engineer field, What do conjugating and Imaginary part mean and affect to the real engineer field?

The real part represents cosine correlation, and the imaginary coefficient represent the sine correlation.  Cosines are even or symmetric functions in the DFT aperture, and sines are odd functions (or anti-symmetric) in the DFT aperture.  So you can only ignore the imaginary coefficients for exactly symmetric input.

Conjugation and multiplication is “just” a mathematical trick (identity) to compute sin^2 + cos^2 which is the hypotenuse squared (Pythagoras).

Any Fourier Transform (and in particular the DFT computed by the FFT) contains information regarding the amplitude and phase of signals (or amplitude and phase response of systems). The output of the Fourier Transform has a complex number associated to every frequency considered. This complex number has a magnitude and phase. The magnitude corresponds to the amount a (complex) sinusoidal component at that frequency will be amplified (or attenuated). The phase angle of the complex number corresponds to the phase shift that will be observed by passing a (complex) sinusoidal component at that frequency.

For example, if I have a frequency response at frequency $f_0$ of $1 + j$ and I pass a signal $x(t) = \exp(j 2 \pi f_0 t)$ through the system, then the output observed will be $y(t) = \sqrt{2} \exp(j 2 \pi f_0 t + j\frac{\pi}{4})$. This is a direct consequece of the fact that $1 + j = \sqrt{2}\exp(j\frac{\pi}{4})$. This is the information that you receive from the FFT.

By multiplying the FFT output by its conjugate, you are losing the phase information and obtaining the square of the magnitude information. Sometimes, this is all you are after. Presumably, this is the case in the algorithm code that you are reading.

Note that the magnitude and phase are real numbers. I believe this answers your question. Please let me know if you still have confusion.


What is the state of the art for vocal elimination/reduction? I have no background in signal processing but I would like to know what is the actual state of the art for vocal elimination/reduction. Obviously, I am not talking about the infamous "stereo elimination", but other, more sophisticated techniques that not everybody is aware of.

I understand that for many projects there might be no software or online application available, but I'd find very useful to hear a demonstration.

Some recent ideas (2-3 years old, I haven't followed recently):

Design an elaborate generative model of the spectrogram of music signals, combining a non-negative mixture of a small set of background frames (models the background music) + a harmonic comb filtered by a non-negative mixture of smooth basis functions (the source filter model of the lead instrument), with constraints on the pitch transitions. Estimate model parameters from the audio clip. Use Wiener filtering to resynthesize the background music/vocals from the estimated spectra. See A musically motivated mid-level representation for pitch estimation and musical audio source separation by Durrieu et al. Such Bayesian approaches can be made more complex ad nauseam, by integrating prior knowledge about note transitions/durations - at the cost of increasingly hairy inference procedures. Demo time!
Techniques inspired by background subtraction in video analysis. The idea is that the accompaniment (especially in "bland" pop music) is invariant by a translation of a few bars (the duration of the "loop" or "riff") - while the vocals are not. Once the structure of the piece of music is known and can be described in terms of repeated groups of bars (a well-studied topic in itself), one can take the median of translated versions of the spectrogram to suppress the vocals. See REpeating Pattern Extraction Technique (REPET): A Simple Method for Music/Voice Separation by Rafil & Pardo. Demo time!

In the commercial space, you can download a trial for ADX Trax.

I looked at this about 5 years ago when I was first hosting karaoke. I didn't find anything that did anywhere near a decent job. The generally had a knob that would allow you to remove more or less vocal, but they all removed/garbled the instrumentals along with it.


Continuous-time mathematical formula for deconvolution filters I have an elementary function $p:\mathbb{R}^2\rightarrow\mathbb{R}$ which (locally) represents an image. It's a polynomial, and its the result of the following 2D convolution:

$$p=f\star G\star \mathbf{1}_{A}$$

where $G(x,y)$ is a 2D Gaussian, and $\mathbf{1}_A$ 'box' function:

$$G(x,y)=e^{-\frac{(x^2+y^2)}{2\sigma^2}}, \quad\mathbf{1}_A(x,y)=1 \text{ for }x \text{ and } y < a, 0\text{ otherwise.}$$

Is there a closed form for a function with which I can convolve $p$, which results in $f$?

Currently I am 'avoiding' the use of such a filter by performing Fourier deconvolution, but mathematically this opens a huge can of worms.

The FT of the convolutional inverse is the point-wise inverse in the FD. So the convolutional inverse of a gaussian might have numerical problems. The FT of a gaussian is a gaussian. Its point-wise inverse thus can be large in many places.

For a box function, the situation doesn't look all that good either. Its FT is a variant on $\operatorname{sinc}$, which has many zeroes and loci of low magnitude.

You can certainly determine the convolutional inverse of a gaussian, but the inverse of a box function is not even a function. Convolving with a box will cancel any function that is periodic with a period multiple of the edge length.

The problem with deconvolution is that you amplify quantization noise and other stuff and the inverse might not even have a stable representation.  Look up "Wiener filtering": it is least-square error inverse filtering with a signal and noise model and can usually be done easily enough in the frequency domain.


Why does an image get blurry beyond a limit of zooming? As we have came cross a phenomenon where a picture get blurry and more blurry as we zoom in a image beyond a limit.

Why does it happen?

Digital zooming is commonly done by using some form of interpolation, which commonly represents making up some form of smooth change between actual data samples.  Assuming the data was low pass filtered (or roughly near bandlimited) before sampling, the information about where any sharp (non-smooth) changes were located between samples has been lost.  Adding back any non-smoothly-changing details would be arbitrary (guessing from context, etc.), and thus possibly incorrect (non-corresponding to the optical field before sampling).

Smooth changes or interpolations across a large area of pixels looks blurry to humans, since this effect does not occur in real life vision often unless one's eye or optical instrument is out-of-focus or not focused on the feature of interest.

Lens structure can no longer focus the image to the image sensor. Every lens system has a depth of field and focus distance. Once you are outside the focus distance and depth of field, photons from the target area are no longer hitting the same pixels and start to hit adjacent pixels which you perceive as blurry. As you come closer impact will be more severe.

That's the optical explanation. If you enlarge an image more than its resolution the application you use will have to interpolate the missing information. It creates new pixels using some interpolation such as bicubic. If it doesn't interpolate, you will see pixelation. Either way the info is estimated therefore contrast is lost.


Sampling $x(t)=\cos(4\pi t)+\cos(2\pi t)$ Imagine that we sample the signal $x(t)=\cos(4\pi t)+\cos(2\pi t)$ with a certain sample frequency $f_s$ and we obtain $x[n]$. Now, by ideal interpolation, we get $y(t)$ from $x[n]$.

How can we know the used sample frequency $f_s$ by looking at $y(t)?$

As a particular case, which $f_s$ verifies that $y(t)=2\cos(2\pi t)?$

How can we know the used sample frequency fs by looking at y(t)?

You can't.

The sampling theorem states that any sampling higher than twice the highest signal frequency allows for perfect reconstruction. In you case any sampling frequency higher than 4 would result in the same $y(t)$.

I'm struggling with this a bit.  Is this a posed problem somewhere?  Maybe you didn't quote it quite right.  For instance the relationship between $n$ and $t$ is not given but assumed to be $n=t$ at integer values.

The only sampling rate (assuming sampling starts at zero) at which $ \cos(4\pi t)+\cos(2\pi t) $ appears the same as $ 2 \cos(2\pi t) $ would be once per second, or whole number multiples of seconds.  But that would also appear to be a DC signal with a value of 2.  So this fails "ideal interpolation" unless that term has some kind of trick meaning.

There is no other $f_s$ which would "verify" your assertion.  A sampling rate of four per second would allow reconstruction of the original function, but only because you have pure cosine functions with zero phase.  The Nyquist criteria is you must sample faster than four per second.

P.S. Personal pet peeve: Sampling is a rate with a typical unit of samples per second, not a frequency with units of Hz.  Yes, most textbooks are sloppy on this.


How does averaging increase the Signal to Noise ratio? When attempting to spot activities on an EEG waveform given a stimuli, it is recommended that several trials be taken and their average is done to increase the SNR.

Intuitively, since noise is not stationary or localized, averaging different trials will have a smearing effect on the noise.

But how exactly does it increase the SNR since the noise is still present (but averaged out)? Can we more rigorously define the effect of averaging on the SNR?

Intuitively this is true, because averaging a zero mean noise processes approximates its expectation value - which is zero. More rigorously:

If the signal $x$ that you want to observe (estimate, actually) is constant for all observations $y$ we can write the $n\text{th}$ observation as
$$
y_n = x + r_n
$$
where $r_n$ is the noise which is different for every observation, of course. The SNR $\gamma_1$ of $y_n$ is
$$
\gamma_1=\frac{\operatorname{E}[|x|^2]}{\operatorname{E}[|r_n|^2]}
$$
where $\operatorname{E}$ is the expectation operator. By averaging $N$ observations of $y_n$ we get
$$
\hat x = \frac{1}{N}\sum_{n=1}^N( x+r_n) = x + \frac{1}{N}\sum_{n=1}^N r_n
$$
The SNR $\gamma_2$ of the averaged observations is
$$
\gamma_2=\frac{\operatorname{E}[|x|^2]}{\frac{1}{N^2}\operatorname{E}[|\sum_{n=1}^N r_n|^2]}
$$
To simplify that expression, we assume that the noise sequences $r_n$ are uncorrelated (1). In that case, the expectation of the squared sum of $r_n$ is equal to the sum of the expectations of the squared $r_n$, i.e.
$$
\operatorname{E}\left[|\sum_{n=1}^N r_n|^2\right] = \sum_{n=1}^N \operatorname{E}\left[|r_n|^2\right]=N\operatorname{E}\left[|r_n|^2\right]
$$
Inserting this in the expression for $\gamma_2$ yields
$$
\gamma_2=N \frac{\operatorname{E}[|x|^2]}{\operatorname{E}[|r_n|^2]}
$$
Comparing the SNR $\gamma_1$ of the single observation and the SNR $\gamma_2$ of the averaged observation we see that averaging increases the SNR by a factor of $N$.

(1) This holds for a Gaussian noise process, for example. Right now I'm not sure for what classes of stochastic processes it holds. Comments on this are welcome.

The above analysis's are correct if the signal is constant/DC; if the signal is varying but you know the shape you have to weight the samples.  In this case you are trying for the amplitude of the waveform with slow edges but you know the shape; and you are trying to improve your signal to noise by looking at the edges as well as the "top"; actually in the case I have in mind the "top" was never reached.
In this case the correct form of combining Gaussian contaminated  samples is:

$\sum_{i=0}^n(S_i\cdot \frac{W_i}{N_i^2})$

Where $W_i$ is the pre-calculated relative value of the waveform at point j; typically you would make the sum of the weights one.  The ratio $\frac{W_i}{N_i^2}$ is pre-calculated and $S_j$   is the sample value.
As usual in signal in signal processing scaling is required but the formula is robust in calculating $(\frac{S}{N})^2 $
As a matter of fact you can see that the formula is a weighted sum of the squares of the S/N across the samples.
I deduced it via a symbolic algebra program a long time ago; if somebody wants I could reconstruct the proof.  As I recall the proof wasn't enlightening and the result looked accidental.  Sometime later I convinced myself that geometric vector space argument would work in an obvious manner but never wrote that down.


Impulse response of a time scaling system Assume a bandlimited signal $X(t)$. Given that the output for this signal is $X(t/2)$, what will be the impulse response $h(t)$ of such a system?

\begin{array}{l}
X( \omega ) \ =\ \int ^{\infty }_{-\infty } x( t) e ^{-j\omega t} dt  \\
Y( \omega ) \ =2X( 2\omega )\\
\end{array}

As the system won't be LTI the output won't be the convolution of the input and the system response. How do I proceed with finding the response of such a system? Any guidance in the direction of the solution will be very helpful.

Edit: Thanks for all the answers, i had reframed this question which was on a similar line. The question goes like this ->
Given an input pulse of duration T, the system scales it to a pulse of 2T, design such a system? It would be of great help if someone could point out the right direction to proceed.

The input-output relation of such a system can be written as

$$y(t)=\int_{-\infty}^{\infty}\delta(t/2-\tau)x(\tau)d\tau\tag{1}$$

Note that the system is linear but time-varying, and such systems can generally be described by a two-dimensional impulse response $h(t,\tau)$:

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

In the given example we simply have

$$h(t,\tau)=\delta(t/2-\tau)\tag{3}$$

For linear time-invariant (LTI) systems, $h(t,\tau)$ only depends on the difference $t-\tau$.

Because you have this information about the fourier transform of the output that $Y(\omega) = 2X(2\omega)$, so I will suggest you to use the time-frequency scaling property of Fourier Transform itself to compute response of the system to any input $x(t)$.

Time-Frequency scaling of the Fourier Transform says that:
$$x(at) \leftrightarrow \frac{1}{|a|}X(\frac{\omega}{a})$$

Notice that the arrow is double sided and so we know that the system's output for any input $x(t)$ in time domain will be inverse fourier transform of $Y(\omega)$. Here $a=\frac{1}{2}$ and so the time-domain output $y(t)$ will be given by:
$$y(t) = \mathcal F^{-1}\{Y(\omega)\} = x(\frac{t}{2})$$.

You can use this relation between input and output to find the response of this system.

You can not call it impulse response of the system. It is just input-output relationship of the given system. But if you could give an impulse as input to the system then you would have gotten $y(t) = \delta(\frac{t}{2}) = 2\delta(t)$ as output.


Doppler Shift and the Pulse Repetition Interval Suppose that I have a signal $x(t)$ consisting of $N$ pulses at a given frequency $\omega_{c}$. I space them $T$ seconds apart, ie the Pulse Repetition Interval (PRI) for the waveform is $T$.

Now this signal bounces off a target and experiences a large Doppler shift, such that the compression of the pulses are evident. My question is -- does the PRI change as well?

In other words, are each of the waveforms still $T$ seconds apart, or does the PRI scale with the Doppler shift?

I guess the question really is, should the model of the second pulse be:

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t -T - \tau}{v}) $$

or

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t -T}{v} - \tau) $$

or

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t}{v} - T - \tau) $$

Where $\tau$ is the target delay and $\frac{1}{v}$ is the Doppler shift.

Note: A similar question is asked in the DSP Stackexchange question
Doppler Shift of Phase Coded Pulse Compression Waveform

Actually the PRI also changes. There is a complete explanation about this topic in "Cobbold - Foundations of Biomedical Ultrasound" chapter 10 (Pulsed Methods for Flow Velocity Estimation and Imaging).

Consider a case where PRI of emitted pulse is 1 s and the wave speed is 9 m/s and the object is 10 m away and moving with speed of 1 m/s toward the antenna. The first echo will be received 2 seconds after the emission (after 1 seconds, the object moved 1 meter toward the antenna and placed 9 meter away from the antenna and the transmitted wave traveled 9 meter away from the antenna). For the next pulse the echo will be received 1.8 seconds after the emission (considering the object is moved closer to the antenna during the PRI between first and second pulse). So the first echo received at t=2 and the second echo received at t=2.8 and the PRI would be 0.8 s (the same is true for the next pulses until the object passes over the antenna).

Also the transformation you have here is something like x(t - (x0-vt)/c ).

The Pulse Repetition Interval (PRI) is chosen based on the maximum unambiguous range you want to detect. The waveform of a pulsed radar is usually represented by the pulse duration $\tau$, the pulse repetition interval (PRI) $T$, or the pulse repetition frequency (PRF) $f_p$. PRI is the time interval between two adjacent pulses.

Target range is the distance between the target and the radar. It can be determined by measuring the time required by the pulse to travel to the target and return. The target range is given by:

$$R = c \tau_d/2$$
where $c$ is the speed of light, $\tau_d$ is the time delay, and the factor 2 accounts for the round-trip delay.

If the time delay is smaller than the PRI, the range determined by the above equation is unambiguous. If the PRF is too high or the PRI is too short, the echo pulses from the target might arrive after the transmission of the next pulse, and then the range measured might be ambiguous. The maximum unambiguous range is given by:
$$R_{\rm unamb} = c (T-\tau) /2.$$


Multiple cheap IMUs vs one expensive IMUs In my project, there is a possibility of either buying one expensive IMU or 5 cheaper IMUs. They have to be cheaper due to budget limitations, and the idea is to use them in order to be able to do "voting" - i.e. vote out outlier measurements of one IMU if 4 other IMUs show a different reading - and also for redundancy. Question: does anybody know if 5 cheaper IMUs (on the order to 300 dollars each) can create a signal at least as smooth, via sensor fusion, as one expensive IMU (on the order of 1500 dollars)?

That really will depend on the variance (standard deviation) of the measurements you get from the five vs the one.

Suppose your expensive IMU gives a measurement distributed as:
$$
x_\mbox{expensive}  \sim  N\left(x_\mbox{truth}, \sigma^2_\mbox{expensive}\right)
$$
and suppose your cheaper ones give measurements distributed as:
$$
x_\mbox{cheap}^i  \sim  N\left(x_\mbox{truth}, \sigma^2_\mbox{cheap}\right)
$$
for $i=1\ldots 5$ and that you form the average of these cheaper measurements to get the "real" measurement:
$$
\hat{x} = \frac{1}{5} \sum_{i=1}^5 x_\mbox{cheap}^i
$$
then $\hat{x}$ is distributed as:
$$
\hat{x} \sim  N\left(x_\mbox{truth}, \frac{\sigma^2_\mbox{cheap}}{5}\right)
$$
So if
$$
\frac{\sigma^2_\mbox{cheap}}{5} \le \sigma^2_\mbox{expensive}
$$
you should be OK with the cheaper sensors.

Of course, there may be other problems with your cheaper sensors: they may be biassed (so that the mean of their distributions will all be different from $x_\mbox{truth}$ ).

And I'm not quite sure what you mean by "voting"? Do you vote off the one furthest from the mean of the 5?  What if the mean changes for the remaining four?  I'll do a bit more digging to see if there's a viable algorithm there, but the direct mean of all five might be the easiest.

If you are going for the 5 IMU option, you should definitely use Kalman filtering with them. It is extremely easy to implement with as many IMUs you like.
If you do not understand the whole concept yet, I suggest watching the following playlist.
Special Topics - The Kalman Filter

Very down-to-earth lectures and if you write everything down from his whiteboard, you can just "C&P" the mathematical formulas to your code and fine-tune from there.


What is a good way to distinguish aircraft noise from other sounds? The FAA's NextGen routes, which were modeled, not measured, have created a bunch of aircraft noise for tens of thousands of people. I want to make it easy for those people (including me) to buy a cheap USB microphone, plug it into a laptop or PC, detect aircraft noise, correlate that with publicly available flight-path data, and file complaints with the FAA, airport(s), and legislators. With enough actual (as opposed to theoretical) data, we might be able to get them to roll back the route changes. Or not, but it's worth a try.

To do that, I need some help with the "detect aircraft noise" part of the problem. It's not enough to trigger off noise above a certain dB threshold, because lawnmowers, leaf blowers, trucks, etc. also make noise. I can weed out some over-threshold noise that doesn't last long enough or lasts too long (aircraft take 5-10 seconds to pass overhead, vs. minutes for lawnmowers and a couple of seconds for trucks), but that's a pretty crude approach. And I still need to be able to factor multiple noise sources into individual components.

Correlation would probably be easier with a pair or array of microphones, to establish noise direction that could be tied to flight paths, but that adds a lot more hardware and complexity. Ideally we'd just need one microphone.

Aircraft include airplanes (jets, turboprops, and standard propellor-driven), as well as helicopters.

I'm a good coder (Java expert, former C expert), but have zero experience with DSP or training algorithms. Can you help me get started?

You need to train a classifier. The way I would approach it is to allow users to submit audio recordings through your app, and indicate which sections of the recordings contain the plane. This would enable you to quickly crowd source a lot of data with which to build a model. Then you would find the STFT and volume of the signal over time. Finally, train the classifier (such as an SVM or neural network) using the amplitude spectrum and volume features. Once you have trained the classifier, you can use it to detect planes by extracting the same features on live audio (STFT isn't computationally expensive, so it should be feasible to do on the phone).

Here's a backgrounder on Audio signal classification.

Your objective is to measure Sound Pressure Level (SPL) at some position on the ground and compare it against safe levels.

There already exist instruments that are doing exactly that and are called Sound Level Meters. These instruments are not cheap because of all the hard work that is required to produce a valid measurement. Their functionality cannot be simply replaced by a "cheap USB microphone". For example, to derive SPL from the microphone's output, you need to have the microphone calibrated (not a cheap process) or use a measurement microphone (not a cheap microphone). To further ensure the validity of the measurement, you need to make sure that amplification and conversion to digital form has been linear and remains linear throughout the frequency spectrum of interest irrespectively of changing factors such as temperature, voltage, etc (not cheap hardware).

What you could do to get some indication of the problem in a relatively simple way is this:

Use a Sound Level Meter with data logging capability.
Set the meter up in the desired area
Examine the recording of the sound level meter and look for regions that are over safe levels. Hopefully, these are not going to be many.
Cross reference the timings that safe levels are exceeded with aircraft arrivals / departures (obviously, adjusted for the location of the mic)

In this way, you have the opportunity to estimate a background SPL (in between aircraft noise) that residents are constantly exposed to and the additional contribution to it by aircraft. By "cross correlating" aircraft timings with SPL measurements, you are doing away with all the problems of trying to distinguish the type of aircraft engine that contributes to the sound level. SPL peaks should correlate more with the aircraft schedule and less with the occasional leaf blower.

To further rule out false positives, that is occasions where the SPL exceeds safety levels because of any reason other than an aircraft, you could set up a cheap sound recorder too and automate the whole process by using the SPL curve to automatically preview the sound recording for what caused the peak at that time.

However, be prepared that the position of the sound level meter could be challenged and the "simplest" way out of this is to obtain many measurements at different positions to reduce its effect.

Hope this helps.

Potentially Useful Links:

Noise and Emissions

Aircraft Noise Issues


What is the advantage of multiple antenna in receiver in energy harvesting ?harvested more power? What is the advantage of multiple antenna in receiver in energy harvesting ?

Assume there are $3$ transmitted antennas in the A device,and $1$ received antenna in the B device,now assume A device use beamforming with $3$ antenna to transmit some signal to B device,and B device harvest $0.1W$ power from the signal.

Now if i add one more received antenna to device B,that is, there are $2$ antennas in B device now,and assume these two antennas are very closed,how much power will device B harvest now? still $0.1W$? or $0.1W\times2$ ?

I am confused about this question,because if device B still harvest $0.1W$,then what is the advantage  of multiple antenna in receiver in energy harvesting ?

If device B harvest $0.1W\times2$,and this mean if i want to create power,i just use $1$ antenna to transmit 1W RF power to a receiver with $\infty$ antennas receiver,then i can get $0.1W\times \infty $ power ,but it is impossible .

After thinking,i think when the antenna number of device B become 2,the harvested power of B will not become $0.1W \times 2$,WHY?

Because when there is just one antenna in B device,the beamforming direction just aim that antenna,because just that antenna has to receive or harvest the signal,so that antenna can harvest all the RF energy

However,if the antenna number become two,device A has to let both of two antenna of B receive or harvest the signal,so the beamforming direction has to change,maybe the middle of these two antenna,so now first antenna will harvest about $0.05W$,and the other will receive $0.05W$ as well.Or perhaps $0.06W$ and $0.04W$,but the total power will not be larger than 0.1W


What is the opposite of sampling called? We call the conversion from a continuous time signal $f(t)$ to a discrete time signal $f_s[k]$ "sampling". Is there a name for the reverse operation, i.e. creating a continuous time signal from a discrete time signal?

In mathematics it is called interpolation.  If you fill data in between data points (in signal processing we call them samples) it is called interpolation.  If you try to fill missing data outside the range of the collected data/samples, it is called extrapolation.  In signal processing, the operation of filling the data in between samples is called interpolation and the blocks that do this operation are called interpolation filters.


Audio signal processing I am trying to code audio signal processing software in java.

Purpose: It will use audio signal processing library to filter out noise from the raw PCM based input (.wav) of heartbeat recording.

Steps involved

Select input .wav file.
Apply Low Pass filter at 33 Hz.
Apply High Pass Filter at 2 Hz.
Apply Amplification with Max Gain.
Store back it as .wav file

There are many C++ library which will cover most of the requirement. But looking for Java for this particular project. Need input from Audio Signal Processing engineers for their input if these steps will be sufficient to able to hear the heartbeats. Also please let me know which order of filtering and which window will be most effective in this scenario.

Input is taken from

http://www.lehman.edu/faculty/hoffmann/itc/techteach/audio/sound2.html

Your low pass should be higher - probably 40 or so.
Your high pass can be higher - around thirty
Your gain should be variable else you will amplify it into distortion - what you want is called normalisation.  You will need to find the maximum volume of your signal and calculate a scaling factor that will bring this up to full scale.  (Say full scale is 1.0 and your maximum value is 0.01 then multiply all samples by 100.)

When I experiment with heart sounds, I find the beat itself to be around 34Hz - I usually do a band pass from 30 to 40 in order to isolate just the beat.

If you do this, you will have a clean recording of ONLY the heartbeat.  If all you want is to count heartbeats or listen to them you will be fine.  If you need to listen for heart defects or other things, then you will need to change the parameters.  You will need higher frequencies and a wider bandwidth for those types of things.

One other thing you need to consider:  30 to 40 Hz doesn't really come through headphones and ears all that well.  Even when you isolate the hearbeat and amplify it to full scale, it is still hard to hear.  Earphones will attenuate it some, and your ears even more so.

As to actually doing the filtering:

You can use IIR filters at higher orders to isolate the heart sounds.  Avoid using chebyshev filters, though.  The "ring" badly and can make one beat "smear" into the next.
If you use an FIR, you may want to downsample first.  It depends on your sampling rate.  If your sampling rate is high (say the typical 44100) then an FIR than can bandpass something around 30Hz will have to be VERY long.  If you downsample to 1000 Hz, then shorter filters can be used.  You could get away with maybe 100 taps at 1000Hz.  Downsampling consists of doing a low pass filter to remove everything above half of your new sampling rate, then "throwing away" a lot of your samples - downsampling by two would throw away every second sample, downsampling by 4 keeps one out of four samples.

When I do audio in Java, I usually import libpd, and use pd to develop the actually processing stuff.  Then the same pd "patch" can be used with libpd in the java programm.
pd  libpd

The order in which you apply the filters (for your  bandpass) shouldn't matter.

I think you're mixing up things here.

Acoustically, 2 Hz is inaudible, and 33 Hz is barely audible. The audible range is 20 Hz to 20.0000 Hz. Heartbeats are somewhere between 40 BPM and 180 BPM, so 0.7 Hz to 3 Hz. So your filter seems to match neither well.

Now the main challenge here is exactly what "noise" you have to filter. For medical purposes, some sounds are clinically relevant which to laymen sound like noise. But if you just need to count the beats, you don't care about noise at all - just determine the envelope, resample it to 100 Hz or so, pass several seconds to an FFT, done.


Alternative Methods to Zero Padding When using a sequential neural network to classify signals with their frequency spectrum, we would need to normalize the length of the signals to perform vectorization. The only method I'm aware at this point is zero padding. Does zero padding have any disadvantages? If so, are there any alternative methods to that?

The other method would be to capture more data in time out to the length desired. This approach would improve the actual frequency resolution achieved, where here frequency resolution means the ability to distinguish in frequency the existence of two closely spaced tones. Zero padding decreasing the bin spacing, meaning how many Hz each bin represents but does NOT increase the resolution: two closely spaced tones would still appear as one tone (spread over multiple bins) if the frequency resolution was not sufficient and the length was only increased by zero-padding.

In general, without further windowing, the frequency resolution in Hz is related to the time duration of the signal according to:

$$f_{BW} \approx 1/T$$

Where $f_{BW}$ is the equivalent noise bandwidth in Hz (which is the same power that would result from a brickwall filter of that bandwidth if white noise was at the input) and $T$ is the duration of the data in seconds.

This relationship applies universally to continuous and discrete systems given the relationship between a rectangular window and its Fourier Transform in both the CTFT and DFT (A Sinc for the CTFT and essentially an aliased Sinc function in the DFT, approaching the Sinc as the number of samples increases).

I don’t know much about sequential neural networks, but depending on the features that you are looking for, could it make sense to resample the waveforms to a common length?


Does negative frequency actually exist or it is just theoretical? This below questions gives good explanation regarding negative frequency

What is the physical significance of negative frequencies?

But I have one query

Does negative frequency only exist theoretically or does it actually exist in nature and can be measured?

Consider a wheel rotating counter-clockwise at one revolution per second. Its frequency of rotation is 1 Hz. If it rotates clockwise, its frequency of rotation is -1 Hz. It's as simple as that.

With respect to a known given fixed point in time (and space if you assume a time-space coordinate system), a negative frequency sine wave is simply the negative of a sine wave or sine() function basis vector that starts with a phase of exactly zero at that same exact point in time.

The fixed point in absolute time can also be defined to be periodic with reference to a given fixed frequency, give any one fixed point.  Or you can use a given phase of another given reference sinusoid at the identical frequency to create you train of periodic fixed phase reference points in absolute time and space.

Added: For phase starts other than zero at the reference point: Only the odd component (in an even/odd signal/function decomposition) of a sinusoid with respect to the absolute reference time/position indicates positive or negative frequency. The even component is symmetric around the reference point and thus makes to difference with respect to the direction of time.

The convention if you don’t have or know some exact point in time for your phase zero reference is to call the frequency of some sine wave or basis vector positive.

As with special relativity, two observers using two different reference frames (or clocks) can end up with two different measurements of phase, thus two different ratios of positive to negative frequency.


FIR filters: direct form / transposed FIR I am learning about FIR filters and I'm confused.
I am trying to find out different types of FIR filters.

Is direct form and n-tap FIR filter the same?
What does transposed FIR filter do?

A finite impulse response (FIR) digital filter implements the following convolution sum

$$y(n)=\sum_{n=0}^{N-1}h(k)x(n-k)\tag{1}$$

for each output sample $y(n)$, where $x(n)$ is the discrete-time input signal, $h(n)$ is the filter's impulse response, and $N$ is the filter length. The values $h(n)$ are also called filter taps, and $N$ is then referred to as the number of taps. The filter described by Equation (1) is also called an $N$-tap filter.

Direct-form and transposed direct-form are just different implementations, i.e. different ways to compute the sum in (1). In theory they are identical, but when computed with finite precision, there can be differences between the different implementations. The direct-form FIR structure is also called tapped delay line or transversal filter.

The two realizations below are the direct-form structure (transversal filter, tapped delay-line) and the transposed structure (from Oppenheim and Schafer, Discrete-time Signal Processing):

Recall the various types of linear phase FIR filters,in which the coefficients h(n) of the transfer function

are assumed to be symmetric or antisymmetric. Since the order of the polynomial in each of these two types can be either odd or even, we have four types of filters with different properties.


Store a FFT with a minimal amount of data I have an array x of length 1024 (stored as 16 bits integers, named for example np.int16 in numpy/python), i.e. the size of x is 1024*2 = 2048 bytes.

(Remark : x comes from an audio .wav file, stored as 16 bits integers, as it is very common. But it is also very common to interpret it as a float array, with values in $[-1, +1]$ by doing: x = x * 1.0 / 2^16)

When I take fft(x), as the input was real, there is some symmetry that makes that I only need to store half of the array fft(x), that's often also called rfft(x) : real fft.

This means that, by taking fft, I translated 1024 real numbers into 512 complex numbers (i.e. can be viewed as 1024 real numbers again) : in a mathematical point of view, we have the same amount of data :

1024 real coefficients -- rfft --> 1024 real coefficients

But in a programmaing point of view, is it possible to store, losslessly* and without compression, the fft of an array of 1024 elements of type int16 (using 2048 bytes) with 2048 bytes maximum ?

If not, what is the minimum number of bytes required to store the fft of such an array?

remark (*) : by losslessly I mean that the original x can be recovered later

In general, quantizing an FFT result is a lossy process.  Given that the FFT twiddle factors are transcendental functions, any finite storage size will result in adding some quantization noise; and that noise increases as the resulting format decreases in bit size.  For extremely sparse FFT results, it may be possible to compress those results, but those are statistically very rare cases.

This is not exactly an answer (I hope it will become it soon) but a result based on experience :

import numpy as np
from scipy.io.wavfile import read, write
from scipy.fftpack import rfft, irfft

(fs, x) = read('nowomannocry_16bit_mono.wav')

y = rfft(x)
M = max(abs(y))
y *= 1.0 * (2**23-1) / M               # mapped into [-2^23, 2^23 -1]
y = np.round(y)
y = np.int32(y)   # here we store the fft by using the full width of 24 bits

z = irfft(y) * M / (2**23-1)
z = np.round(z)
z = np.int16(z)
print z - x, max(abs(z-x))
# [0 0 0 ..., 0 0 0] 0

According to this small test (done on Bob Marley's No woman no cry, 16 bit, mono ;)), 24 bits seem to be enough to store the FFT of the original signal, in such a way that the original signal can be recovered by RFFT.

I tried with 16, 20, 22 bits, and it wasn't enough.

These are only experiences, but can something be proven in this direction ?


Recommendation for book - Writing DSP code in C I am looking for some good book, that simply show how you actually write a code in C, to do all the main DSP methods .

FFT.
Low-pass and high-pass filters.
Auto-correlation.
Noise processing.

And all the basics of DSP, from theory into a real code in C.

For example, I have got 1000 samples, now i would like to compute FFT of it, remove noises, then bring back to the time axis.

Is there something good that cover all of these?

I can recommend you two books about DSP for C language.

  Embree P. M. - C Language Algorithms for Digital Signal Processing

It is old and you can easily get it second-hand for a decent price. It covers pretty much all 4 topics that you described.

The other one I recommend is:

  Malepati H. - Digital Media Processing: DSP Algorithms Using C

It covers way more topics (including error correction algorithms and image/video processing algorithms). The nice thing about that one is that test signals are provided to test your implementation.

I just recalled that there is another book on my desk which has some C code for DSP in it.

  Press W. H. - Numerical Recipes: The Art of Scientific Computing

There are two chapters about the Fourier Transform and its applications.

One last suggestion - if you can, please use an existing C library instead of writing your own tools. There is no point in re-inventing the wheel.

You can look to the source code of openSAL and Octave sources of Matlab signal processing functions (not C, but useful). This is what I did.


Gaussian noise generator algorithm trouble I stumbled upon the following algorithm for Gaussian noise generation and I just can't figure out how this is supposed to work at all.
The algorithm is as follows:

First of all: $j = \arg\min_i (q_1 - p[i])$ does not depend on $q_1$ and therefore yields the same result every time - I assume $\arg\min_i (|q_1 - p[i]|)$ is what's really meant here.
But then there's more: I can safely assume that $q_1$ is distributed linearly within $[0,1]$, otherwise I could just scale that, round afterwards and wouldn't have so much hassle.
In this case I got another problem: As $q_1$ is linearly distributed, choosing $j$ near to zero (which is what $\mu=0$ would implicate) is not particularly likely, so I suspect the resulting distribution is not a normal distribution.

This algorithm reminds me of inverse transform sampling ( http://en.wikipedia.org/wiki/Inverse_transform_sampling); however it does not appear to be a variant of it. I'm a bit confused here because that's somewhat easy to compute and if p was a cumulative distribution here, I'd understand this one right away.

So basically I'm stuck with understanding this one and looking forward to any hints as to where I'm wrong.

Thank you folks in advance!

I agree with you, this looks like a flawed description of an inverse transform generator.

The $p[i]$ shouldn't be like described, but should be the prefix sum of these, corresponding to the cumulative distribution.
The $\arg\min$ should indeed be taken on the absolute value (or allow positive values only).
As the pdf is computed on the positive $i$ only, the range of the uniform drawing should be $[0,\frac12]$ (or the $p[i]$ scaled by $2$.)

(When the $p[i]$ are not cumulated, they are in $[0,\frac1{\sigma\sqrt2\pi}]$, and the range $[0,1]$ doesn't make sense.)

For better efficiency, you can precompute a lookup-table of the $i$ indexes for sufficiently dense $q_1$.

Yes, the algorithm is distinctly poor.  I've made a quick attempt to implement it (as stated) in R, and this is the histogram of the resulting values.  Please let me know if there's a problem with it (as the algorithm is stated, not trying to fix the algorithm).

It looks no different from a uniformly distributed random variable.

R Code Below

#21675

q_21675 <- function()
{
  #1
  sigma_i <- 3
  #2
  G <- 256
  p <- rep(0,G)
  for (i in 1:G)
  {
    p[i] <- 1/(sigma_i*sqrt(2*pi))*exp(-i*i/(2*sigma_i*sigma_i))
  }
  #3
  q1 <- runif(1,0,1)
  diff <- rep(0,G)
  for (j in 1:G)
  {
    diff[j] <- q1-p[j]
  }
  j <- which.min(diff)
  #4
  q2 <- runif(1,-1,1)
  addition <- q2*j
  return(addition)
}

x <- rep(0,1000)
for (t in 1:1000)
{
  x[t] <- q_21675()
}

hist(x)


Power and Energy spectral densities of a linear system Here is a schematic of my desired circuit.

I want to calculate Power and Energy that the load dissipates during 0< t<8 You can see my attempts:

Therefore

How can find load's energy in Frequency domain for desired time
domain interval (0 < t < 8) ?
Did I obtain power correctly in
    frequency domain?

The relationships, with $x$ the input voltage to the network, $y$ the output voltage from the network, $R$ the load resistance, $p(t)$ the instantaneous power dissipated at the load, and $E$ the energy dissipated at the load and $t_0=8$

$$Y(f)=H(f)X(f)\\
y(t)=\mathscr{F}^{-1}\{Y(f)\}\\
p(t)=\frac{y^2(t)}{R}\\E\bigg|_{t=0}^{t_0}=\int_0^{t_0}p(t)dt=\int_0^{t_0}\frac{y^2(t)}{R} dt$$

IF the cable or network is a lossless cable (!), the example turns trivial:
$$y(t)=x(t)\\
p(t)=\frac{x^2(t)}{R}\\
E\bigg|_{t=0}^{t_0}=\int_0^{t_0}\frac{x^2(t)}{R} dt$$


[2D fourier transform]: Most people can't explain this I am confused with 2d magnitude plot of frequency spectra. So we have 2 images, the first one is shown at the top, and the dilated or enlarged version of the white box is shown at the second row.

For the dilated version, in the spatial domain when travelling in the $x$-direction, you see a fast change in intensity, and when you travel in the $y$ direction, the change is slower, so that the $v$ direction frequency component is smaller. In the frequency domain, $x$ maps to $u$, and $y$ maps to $v$.

I do understand that when you enlarge in one domain, you must do the opposite to the other domain.
The thing I don't understand is that in the dilated version, when you travel in $x$, and you still see the same rate of change as in the original image. Since $u$ maps to $x$ in the frequency domain, we expect that the frequency spectrum looks the same in the $u$ direction. However, the result is that, in the $u$ (or the horizontal) direction of frequency domain, its frequency component shrinks; the $v$ direction remains intact.  But we just change its $y$ direction!

Update:
To be more precise, please refer to the area of the red box. In the second figure, it gets shrink.

The scaling property of the Fourier transform tells you that if  $\mathcal{F}(x(t)) = X(f)$, then  $\mathcal{F}(x(at)) = \frac{1}{|a|}X(\frac{f}{a})$, with $a\neq 0$. So when you dilate in time, you contract in frequency, and vice-versa.

Now, for continuous 1D signals, the Fourier transform of a uniform box will be a cardinal sine. If you dilate the box, the cardinal sine will be contracted. This could be coarsely observed from the location of its zeros.

Here, you deal with a 2D discrete signal. So the result should be a discrete (or periodic) version of a cardinal sine. Its spread in one direction can be visually estimated by the horizontal and vertical black lines, since your 2D discrete Fourier transform is separable.

Now look at one horizontal line in your red boxes. The blackest dots that separate the horizontal gray/white bumps have the same location on both images. Now look at several lines. The black dots align vertically. So you see that the vertical lines are spaced in the same manner on the top and the bottom, which is exactly what you expect.

The result in the vertical direction is the controverse. Horizontal black lines are closer to each other in the bottom row.

What is sometimes misleading in Fourier plots is that some features in the frequency domain look orthogonal to those in the spatial domain. You can see more of this effect in An Intuitive Explanation of Fourier Theory.

If you want to gain additional conviction, crop a line from each spatial image, and an horizontal line from the spectrum, and superimpose the plots.

You increased the length of the white box. Now when you traverse in x direction , the rate of change remains the same. But now you have more number of pixels where intensity changes from black(0) to white (255). This has resulted in more energy along that direction. check your image. The u component is more brighter in expanded image.

This energy can be represented using few harmonic coefficients. Check your image. You would see a denser grid along u direction.


Difference Between Convolution and Multiplication I read that multiplication is convolution in frequency domain. I also understand that convolution is just polynomial multiplication. Can somebody explain what are the advantages of doing convolution in frequency domain? Also, i would like to know how to multiply two signals in time domain and what is the significance of it.

Multiplication in your first sentence is term-by-term multiplication: $z[n] = x[n]y[n]$ for all $n$. Convolution, for discrete-time sequences, is equivalent
to polynomial multiplication which is not the same as the term-by-term
multiplication. Convolution also requires a lot more calculation: typically $N^2$
multiplications for sequences of length $N$ instead of the $N$ multiplications
of the term-by-term multiplication.

The key point of Fourier analysis is that term-by-term multiplication in one
domain is the same as convolution in the other domain.  So, in order
to calculate the results of a convolution, you can either do it directly,
using $N^2$ multiplications, or transform to the other domain, do a term-by-term
multiplication, and transform back. This requires two transformations
to go from one domain to the other, $N$ multiplications in the other domain,
and one inverse transformation to come back to the domain
where the convolution result is needed. This more complicated process
can, in fact, require less computation because the transformations can
be done very efficiently via the Fast Fourier Transform (FFT) algorithm
which requires about $N \log N$ multiplications. So, compare the computational
effort in calculating the three transforms and doing the term-by-term
multiplication
to $N^2$ to see if the FFT gives you a more efficient method of computing
a convolution.

The convolution operation can be imagined as swiping a credit card ,each instance of both the signals are multiplied with each other.
In which case a bar code scanner is more like multiplication ,only two instances of each signals are multiplied  .


The Fourier Series Of This Triangle Wave I am using matlab to study digital signalling and have come across a problem which i was wondering if anyone with more experience could help me with.

I need to work derive the Fourier series of a triangle wave that i have generated, I just do not know how to actually go about this problem in Matlab.

I am generating a 100hz Triangle signal using the following code:

t = 0:1/10000:1;
f=100;
x1 = sawtooth(2*pi*f*t, 0.5);
plot(t,x1);
axis([0 0.10 -1 1]);

Now how should i go about deriving the Fourier series of this signal, i am completely lost.

Any help would be appreciated.

Check out Wikipedia's Fourier Series page.  Their "Example 1" shows how to derive the Fourier series of a sawtooth wave.  All you have to do is normalize the results for your particular time and amplitude values.

I assume you want to calculate the FFT of this signal in MATLAB? This would be like this:

N_fft = 2^15; %Just an FFT size
x1_fft = fft(x1, N_fft); %Take the FFT.
plot(abs(x1_fft));


Text in image is low or high frequency component? Text present in captured images is high frequency or low frequency component?Why?

The edges of the words will probably contribute to a high frequency, the background of the text and the portion within the boundaries of the letters of the text will be otherwise.

Actually this is very general. It can also depends on the clarity of the text.


How is it possible to Filter out a person's voice out of 100 of other voices? So, I've just learned that the human voice is not a single sine wave its a bunch of unlimited Sine waves each having a different frequencies,

According to Wikipedia,

  The voice consists of sound made by a human being using the vocal
  folds for talking, singing, laughing, crying, screaming, etc. Its
  frequency ranges from about 60 to 7000 Hz.

So if the human voice is a composite signals, it contains various frequencies ranges between 60Hz to 7KHz.

Suppose if there is a group of people singing the same song all together, each person has its own ranges of voice frequencies,

For example,

If a person A has the following frequencies,
100Hz , 250HZ 6KHz, 10Hz, 87Hz, 52Hz, 2KHz.......

and a Person B has the following,
217Hz1, 11Hz, 12Hz, 2323Hz, 839Hz, 4KHz, 100Hz, 10Hz.....

there must be so many frequencies which are similar in both the person A & B, like in above example the frequencies 100Hz and 10Hz are common between two persons.

I was watching a TV Show name "Fringe" where they filter out the particular Man's voice from an audio file while there were other people voice present there too.

So how does exactly they filter someone's voice out of the voice of 100s of people if there are so many frequencies common among all of them does it have to do something with the amplitudes of person's frequencies ?

If the signal is recorded using just one microphone, you can use methods such as spectral subtraction. This method is more suitable for "constant" noise, like the noise from a fan or an idle engine. Other methods rely on statistics and perceptual models of speech. If the signal is recorded with several microphones, you can use blind source separation for separating the (speech) signals. As it stands today, you won't get perfect results. The typical end-result is always a trade off between "noise" and clarity of the speech signal of interest. More "noise" suppression --> more degradation of the signal of interest.

okay, lets say we have a sound file of two people talking.. if they were not talking in unison, it would be possible to separate the tonal element of their speech. it might be trickier to separate the noise elements of their speech (ssss or fffff sounds), but again if they were not speaking in unison it might be possible. this would be much easier to do if we could perform spectral analysis on sound without losing the phase information.

basically, spectral analysis takes a waveform (sound as single wave or line) and separates all the individual tones so you can see them from low to high and left to right in time. when doing this, most of the information showing the rise and fall of the single waveform is lost. if we could preserve the phase information, this would be much easier because all individual component frequencies of one voice would be harmonically related to each other, thus their phase would line up.

at present, i don't know of any algorithm that achieves this, but i believe it to be theoretically possible. i remember reading some article about melodyne's creater peter neubäcker working on a way to do this, so that two singers singing together could be separated, but i don't know how i would find this article again


If the frequency response of a system is time-independent, is this system LTI? Intuitively, I feel like a time-variant system would necessarily have a time-dependent frequency response, and vice-versa. So, is the time-independence of the frequency response necessary and sufficient to say that a sytem is LTI ?

LTI means Linear Time-Invariant systems. The system has to satisfy two conditions.
(1) Linear and (2) Time-Invariant.
(1) Linear means, if the response of the system due to load Px and Py is Rx and Ry respectively, then for the load (Px+Py), the response of the system will be (Rx+Ry).
(2) Time-Invariant means, the parameters of the system does not vary with time.

If the frequency response of the system is Time-Independent, then it suggests that the parameters of the systems are not changing over time, i.e., the system is time-invariant.

But it does not say anything about nonlinearity. Hence, in a nutshell, if the frequency response of a system is time-independent, this system is NOT an LTI?

I think the answer to your question is "yes" and this is why:

A Linear Time-Variant (or possibly time-variant) system is fully described (from the POV of output $y(t)$ in terms of input $x(t)$) by this convolution integral:

$$ y(t) = \int\limits_{-\infty}^{\infty} h(t,u) \, x(u) \, \mathrm{d}u $$

$h(t,u)$ is the (possibly time-variant) impulse response function, evaluated at time $t$ of an impulse that was applied to the input at time $u$. That is what $y(t)$ is when $x(t)=\delta(t-u)$.  In general, impulse responses of impulses applied at different times are different.

However, if the Linear system is also Time-Invariant (LTI) then

$$ h(t,u) = h(t-u) $$

which means that the only difference between impulse responses of impulses applied to the input at different times, $u$, is that they are delayed by the same time that the input impulse is delayed, otherwise exactly the same.

If the system is LTI, this impulse response that was dependent on the time that the impulse is applied, and is a function of two variables, $h(t,u)$, becomes a function of one variable $h(t)$.  It doesn't matter what time, $u$, the impulse is applied other than the impulse response happens the same time later: $h(t-u)$.  Knowing the impulse response for an impulse applied at time $0$ which is $h(t,0)$ tells you exactly what the impulse response will be for an impulse applied at any other time.

The LTI frequency response, $H(f)$ is the Fourier Transform of the impulse response of an impulse applied at time $0$:

$$ \begin{align}
 H(f) &\triangleq \mathscr{F} \big\{ h(t) \big\} \\
      &= \int\limits_{-\infty}^{\infty} h(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t \\
      &= \int\limits_{-\infty}^{\infty} h(t,0) \, e^{-j 2 \pi f t} \, \mathrm{d}t \\
      &= e^{j 2 \pi f 0} \int\limits_{-\infty}^{\infty} h(t,0) \, e^{-j 2 \pi f t} \, \mathrm{d}t \\
\end{align}$$

That factor $e^{j 2 \pi f 0} = e^0 = 1$ is in there for a reason.  This is the Fourier Transform of an impulse response of an input impulse that is not delayed by any time.  But if the impulse response is of a delayed impulse, everything is delayed by that time and you get a linear phase shift factor because of that, which I want to remove, for the purpose of discussing time independence of the frequency response.

$$ \begin{align}
 H(f,u) &\triangleq \mathscr{F} \big\{ h(t+u,u) \big\} \\
      &= \int\limits_{-\infty}^{\infty} h(t+u,u) \, e^{-j 2 \pi f t} \, \mathrm{d}t \\
      &= e^{j 2 \pi f u} \int\limits_{-\infty}^{\infty} h(t,u) \, e^{-j 2 \pi f t} \, \mathrm{d}t \\
\end{align}$$

$H(f,u)$ is the time-dependent frequency response of the linear system around time $u$.  You can tell that if the linear system is time-invariant (LTI) then $h(t,u)=h(t-u)$ and consequentially $H(f,u)$ does not vary with time $u$.  So, at least, the converse is true: If the system is LTI, the frequency response is independent of time.  But I cannot see how any other case of a Linear system where $H(f,u)=H(f)$ that can allow $h(t,u)$ to be dependent on $u$.


Tool for Harmonic Analysis & Synthesis I've been using a program for identifying harmonics for a wave after doing a fourier analysis.

http://www.phon.ucl.ac.uk/resource/sfs/esynth.htm

I would like to be able to tweak some of the harmonics involed in a wave and also to add harmonics.

The program mentioned above permits this kind of operations but it isn't precise.

Does anyone know of a good tool for harmonic analysis and synthesis like the one mentioned above?

ps:It's important that the program exports the harmonics to a file.

thanks.

Try NCH Software

From Wikipedia:

  NCH Software is a software development company founded in 1993 in
  Canberra, Australia. The Colorado office was started in April 2008 due
  to the large U.S. customer base. NCH Software primarily sells to
  individuals via their website.

And further

  Browser Toolbar

  During the installation process NCH Software also by default (but you
  can uncheck the option during install) installs Google Chrome. This
  can be unselected during install or removed after install. In older
  versions of their software a browser toolbar powered by Conduit,[11]
  was installed as well as making changes to the default homepage and
  search engine.


Drawing a Spectrogram off of FFT - what are the decibels relative to? I am trying to draw a spectrogram from FFT data I get from audio files.
The problem is that I can't figure out is what the yielded power (decibels) are relative to. I'm getting some arbitrary -100db-ish numbers, but I don't understand how I can have them represent a finite intensity-scale.

If I have a gradient which will represent the intensity of these numbers, where silence is black and the maximum is yellow I need to know what the maximum is. I'm really not good in the DSP-area so help would be appreciated.

For audio files you want to use dBFS, dB relative to full scale:

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

dB are in fact meaningless units on their own. They are a ratio unit, so for dB to be meaningful, it needs to be computed within a context. For example, dB SPL (sound pressure level) is computed with respect to the pressure of 20 micropascals, which is the threshold of audibility in terms of sound pressure.

http://artsites.ucsc.edu/ems/music/tech_background/te-06/teces_06.html


Power Spectrum in dBFs I'm trying to get the spectrum of my signal in dBFS units. I'm not sure what is the refernce that I should devide by. I computed the power spectrum as follows:

1) FFT normalized to the number of samples

p = abs(fft(signal))/(No_of_samps/2);

2) and then I calculated the power spectrum for half of the samples (one-sided)

p = p(1:Nsamps/2).^2;

3) I  tried to represent it in dBFS as follows:

p_dBFs = 10*log(abs(p)/(vref/2)^2);

Where vref is my reference p-p voltage in the ADC.

I'm just not sure whether I should devide by the reference voltage to get values in dBFS, like when I apply a FS signal I would like to see the highest peak in the spectrum at exactly 0 dBFs.

Note: The time-domain signal is also multiplied by a hamming window before computing the FFT.

Any help is appreciated.

Sorry for the long post. Thank you guys very much.

Hello again guys after a lot of thought, I finally broke this down and I would like to share it with you hoping it'll help.

So what I did was:

1- I used matlab's pwelch to get the power spectrum like this

[AutoPowerSpectrum ,F] = pwelch(signal, Nsamps, 0, Nsamps, fs, 'power');

you can read about the different parameters for pwelch in matlab's documentation right HERE.

2- I calculated Vpeak squared, because as far as I know pwelch produces Vrms values.

AutoPowerSpectrum = 2.*AutoPowerSpectrum;

3- got the Vpeak by squaring the above equation

AutoPowerSpectrum = sqrt(AutoPowerSpectrum);

4- converted it to dBFs

p_dBFs = 20*(log10(AutoPowerSpectrum./(vPerLevel*2^(N-1))));

where: vPerLevel is the voltage corresponding to each level in my ADC. It is calculated like this Vref/2^N.

and N is the resolution of the ADC.

Finally used plot() to display it. This solution worked with me and gave me the signal in dBFS, i'm opened to edits and suggestions.

Hope this helped.

Cheers.


Image processing coding I was wondering which is the best language for image processing? I know MATLAB has a good library and user community for it. I havent checked it out, its there for python as well.

What is the most efficient language for this?

It's not only about programming language but library you are using. I can think of the following:

MATLAB - image processing capabilities are quite ok, but for more advanced and real time processing you would need some low-level stuff. Additionally, it does not offer very good portability.

Mathematica - good for prototyping and quick visualization, but that is all in my humble opinion.

OpenCV - I think this is the most popular library in IP community. Great capabilities (GPU computing, Machine Learning module, GUI - what else do you need?), fast, and still under development (so minor bugs are being removed very quickly). Regarding community - it's big! Mainly for C/C++ programming, but also Python (probably suits you).

JAI - Java Advanced Imaging - only if you like Java. Personally, I do not like it.

ImageMagick - you can use it with many programming languages, check the API.

CxImage - good if you want to create something better than MS Paint and with some functionality of the Photoshop.

CImg - obviously to be used with C++, but OpenCV is still better.

PIL - Python specific library with lots of functionality. If you like it then you can take a closer look.

SimpleCV - basically it is OpenCV python binding with some tweaks. Very easy to use and quite efficient.

scikit-image - also Python library, but in my opinion worse than SimpleCV (although my friend would not agree to that). One advantage is that they included DAISY features extractor - quite useful if you need dense features.

GIL - part of Boost, but less functional that OpenCV. Although if you like and use Boost then for some basic functionality it should be ok.

ResIL - continuation of DevIL project. A major advantage is operating on many file formats, including some game files.

PINK - originally written in C++. I recently tried it out and it was very pleasant to work with. Great integration with Python and fast in embedded C. Apart from that they do have plenty of nice algorithms.

Obviously, there is plenty more Image Processing libraries, but these are ones I had contact with. So if you want to start with something, then choose OpenCV (preferably with C++ framework) - you will not regret! On the other hand, if your programming skills are not strong, then you might want to think about using Python-based libraries - really easy to learn and set-up.

If you want to roll your own algorithms or code, then include OpenCL, as many image processing kernels are "embarrassingly parallel".

Or, on Big Iron, Fortran still beats most other languages at numeric vector performance.


Algorithm(s) to mix audio signals without clipping I'd like to mix two or more PCM audio channels (eg recorded samples) digitally in an acoustically-faithful manner, preferably in near-real-time (meaning little or no peek-ahead).

The physically "correct" way to do this is summing the samples. However when you add two arbitrary samples, the resulting value could be up to twice the maximum value.

For example, if your samples are 16-bit values, the result will be up to 65536*2. This results in clipping.

The naive solution here is to divide by N, where N is the number of channels being mixed. However, this results in each sample being 1/Nth as loud, which is completely unrealistic. In the real world, when two instruments play simultaneously, each instrument does not become half as loud.

From reading around, a common method of mixing is: result = A + B - AB, where A and B are the two normalized samples being mixed, and AB is a term to ensure louder sounds are increasingly "soft-clipped".

However, this introduces a distortion of the signal. Is this level of distortion acceptable in high-quality audio synthesis?

What other methods are there to solve this problem? I'm interested in efficient lesser-quality algorithms as well as less-efficient high-quality algorithms.

I'm asking my question in the context of digital music synthesis, for the purpose of mixing multiple instrument tracks together. The tracks could be synthesised audio, pre-recorded samples, or real-time microphone input.

The physically "correct" way to do this is summing the samples. However when you add two arbitrary samples, the resulting value could be up to twice the maximum value. ... The naive solution here is to divide by N, where N is the number of channels being mixed.

That's not the "naive" solution, its the only solution.  That's what every analog and digital mixer does, because it's what the air does, and it's what your brain does.

Unfortunately, this appears to be a common misconception, as demonstrated by these other incorrect non-linear "mixing" (distortion) algorithms:

Mixing digital audio (the wrong way)
A quick-and-dirty audio sample mixing technique to avoid clipping (don't do this)

The "dividing by N" is called headroom; the extra room for peaks that's allocated above the RMS level of the waveform.  The amount of headroom required for a signal is determined by the signal's crest factor.  (Misunderstanding of digital signal levels and headroom is probably partially to blame for the Loudness war and Elephunk.)

In analog hardware, the headroom is maybe 20 dB.  In a hardware DSP, fixed-point is often used, with a fixed headroom; AD's SigmaDSP, for instance, has 24 dB of headroom.  In computer software, the audio processing is usually performed in 32 bit floating point, so the headroom is enormous.

Ideally, you wouldn't need to divide by N at all, you'd just sum the signals together, because your signals wouldn't be generated at 0 dBFS in the first place.

Note that most signals are not correlated to each other, anyway, so it's uncommon for all the channels of a mixer to constructively interfere at the same moment.  Yes, mixing 10 identical, in-phase sine waves would increase the peak level by 10 times (20 dB), but mixing 10 non-coherent noise sources will only increase the peak level by 3.2 times (10 dB).  For real signals, the value will be between these extremes.

In order to get the mixed signal out of a DAC without clipping, you simply reduce the gain of the mix.  If you want to keep the RMS level of the mix high without hard clipping, you will need to apply some type of compression to limit the peaks of the waveform, but this is not part of mixing, it's a separate step.  You mix first, with plenty of headroom, and then put it through dynamic range compression later, if desired.

A + B + {
    (|A| = A) = (|B| = B) = true: -AB;
    (|A| = A) = (|B| = B) = false: AB;
    else: 0
}

That is, if both A and B share a sign, apply a limiting offset.  The magnitude of the offset is the product of A and B.  The direction of the offset is opposite to that of A and B.

If A and B do not share a sign, no limit is applied, as there is no way to overflow.


What are True Positive and False Positive with reference to image processing Say for instance, face detection algorithms - I see researchers testing their algorithms and table their results under TP and FP. What are these?

A ,,true'' positive would be an image that has a property (in the context of face detection maybe: ,,contains a face'') and that is recognized by a program as such. A ,,false'' positive then is one that does not have the property but is recognized anyway.

These terms are validation metrics used for verifying quality of a segmented image. In a scenario where you want to compare a segmented image with ground truth, then taking the ground truth image as base of comparison you can make assumption of taking foreground as "white" pixels and background as "black" pixels in ground-truth. The terms that you referred then would mean:

True positive (TP)  : pixels correctly segmented as foreground
False positive (FP) : pixels falsely segmented as foreground
True negative (TN)  : pixels correctly detected as background
False negative (FN) : pixels falsely detected as background

These metrics are then used to calculate sensitivity, specificity and accuracy as:

sensitivity :
The sensitivity tells us how likely the test is come back positive in someone who has the characteristic. This is calculated as TP/(TP+FN).

specificity :
The specificity tells us how likely the test is to come back negative in someone who does not have the characteristic. This is calculated as TN/(TN+FP).

accuracy : (TP+TN)/(TP+FP+TN+FN)

Hope this helps !!


Dimension of Matrix in MIMO with multi-path In a MIMO system, the received signal is supposed to be $y = Hx + n$, where $H$ is the channel of dimension $N_t \times N_r$, $x$ is the transmitted signal with dimension of $N_t\times 1$ , and $n$ is the noise, finally the resulted signal $y$ represents the received signal of dimension of $N_r\times1$.

Assuming that I need to use channel $H$ as multi-path channel whose number of taps is $e.g.$ 5 .. The question is:

Will the dimension of channel $H$ will be $5N_t\times N_r$?
If so, how can we perform the matrix multiplication of $Hx$ in  $y = Hx + n$, so the received signal $y$ will be of dimension $5N_r\times1$?

The short answer is: if you transmit one block of $N$ symbols (no data before it for at least $L$ symbol times) over a frequency-selective channel of length $L$, then the channel matrix will be of dimension $NN_TN_R\times N$.

The long answer: start with $N_T=N_R=1$. Then the $n^{\text{th}}$ received sample can be written as

$$y_n=\sum_{l=0}^Lh_lx_{n-l}+z_n$$

If you arrange these samples in a vector-matrix form you get

$$\mathbf{y}_{N\times 1}=\mathbf{H}_{N\times N}\mathbf{x}_{N\times 1} + \mathbf{z}_{N\times 1}$$

Now let $N_T=1$ while the number of receive antennas $N_R$, then the $n^{\text{th}}$ received sample over $m^{\text{th}}$ receive antenna is

$$y_n^{(m)}=\sum_{l=0}^Lh_l^{(m)}x_{n-l}+z^{(m)}_n$$

for $m=1,\,2,\,\ldots,\,N_R$. If you arrange these sample in a vector-Matrix form, you get

$$\mathbf{y}_{NN_R\times 1}=\mathbf{H}_{NN_R\times N}\mathbf{x}_{N\times 1} + \mathbf{z}_{NN_R\times 1}$$

Similarly, you can do the general case for $N_T$ transmit and $N_R$ receive antennas.

As you can see, to make this work on block-by-block basis, you need to padding $L$ zeros before each block to eliminate the inter-block interference. Alternatively, you can use cyclic-prefixed (CP)-OFDM, which is more efficient at the receiver side.

Will the dimension of channel $H$ will be $5N_t$x$N_r$ ?? If so, How can we perform the matrix multiplication of $Hx$ in  $y = Hx + n$ ? ... So the received signal $y$ will be of dimension $5N_r$x1, right?

No!

$\mathbf H$ is defined to work as following: With $\mathbf x$ being the transmit vector, and $\mathbf y$ being the receive vector, and in the absence of noise, it is

$$\mathbf y = \mathbf H\mathbf x\text.$$

Thus, since the dimensions of $\mathbf x$ and $\mathbf y$ are fixed (at $N_r$ and $N_t$, respectively), the only choice for the dimension of $\mathbf H$ is $N_t\times N_r$.

A multipath channel can't be represented with a simple construct such as 2D Matrix. So, channel Matrix and multipath will not work together unless you redefine $\mathbf x$ and probably also $\mathbf y$.

I like your idea of "inflating" the transmit vector i.e. you could argue that individual delay slots are active orthogonally to each other and could restrict $\mathbf H$ to be causal. However, I'm personally not aware of anyone who's done that – so, finding a good provable notation that works out like the actual multipath MIMO channel would be worth giving it a shot. I'm pretty sure that the result won't be a simple $5N_t \times N_r$ Matrix, but probably be a tensor of higher dimensionality ($N_t \times N_r \times 5 \times 5$) to describe the convolution that happens between every path, but it'll be certainly interesting.

In short: You will need a more complicated system model.

In fact, describing and correcting multipath MIMO channels is so complicated that in practice we don't. Instead, modern communication systems that use MIMO (that being WiFi and LTE) use techniques like OFDM to split the wideband multipath MIMO channel in many narrowband single-path MIMO channels and then deal with these individually, because these single-path MIMO channels are the MIMO channels we know how do deal with using the well-established "boring" math.


Power of summed signals Pretty much a newbie here, but would like to understand a pretty basic concept.

We all know that a quick and dirty way of calculating signal power of a part of a discrete sampled system (such as audio signals) is the age-old power formula for RMS power, where we simply break the signal into chunks of N samples and sum the squares of the amplitudes, then divide by the number of samples

$$
P_\mathrm{x} = \dfrac{1}{N}\sum_{n=1}^{N}|x[n]|^2
$$

So let's say I have 10ms samples of two completely uncorrelated signals at 48khz (480 samples), and I have calculated the RMS power of these samples using the above formula and have come up with P1 and P2, the power of my small segment of each signal.

Then if I sum the two signals (simple software mixer function), what is the resultant power of the two summed signals?  Is it simply P1 + P2?

If so, do I simply keep going?  So if I had four signals with P1, P2, P3 and P4, if I sum all four of them together, is the resultant RMS power P1+P2+P3+P4?

The formula you gave is an estimator for the "statistical power", i.e. the expectation of the square of a continuous  underlying random process:
$$ P_X = E(X^2(t)) $$
Computing the latter theoretically would require to know the probabilistic behavior of $x(t)$, which you don't have access to in real life, especially since:

you get sampled data (not continuous),
you generally have a limited number of realizations.

This is even more complicated when you add variables, because you should know their paired probability behavior. A classical result is that, for variables $X$ and $Y$,
$$\mathrm{Var}(aX+bY) = a^2\mathrm{Var}(X)+ b^2\mathrm{Var}(Y)+ 2.a.b.\mathrm{Cov}(X,Y)$$
This quantity can take different values, depending on the last term: $\mathrm{Cov}(X,Y)$. For uncorrelated variables, the covariance vanishes, and the variables add as:
$$\mathrm{Var}(aX+bY) = a^2\mathrm{Var}(X)+ b^2\mathrm{Var}(Y)$$
So, for uncorrelated signals, powers add. And this generalizes to several signals.
Note that when  you sum over a time frame to estimate such quantities, you somehow assumes that you can estimate quantities by averaging in time, with assume that the signals are correctly sampled for a stationary and ergodic process.
For $K$ signals of  identical power $P$, the combined power can take value:

minimum: $0$ when signals sum to zero,
$KP$ when uncorrelated,
maximum: $K^2 P$ when signals are identical.

Additional pointers:

Addition of signal variance and noise variance
Variance of Uncorrelated Variables
Signal averaging
Superposing Signals

For two signals in general:

$ {1 \over N} \sum_{n=0}^{N-1}  (x_1[n] + x_2[n])^2 = {1 \over N}(\sum_{n=0}^{N-1} x_1[n]^2 + \sum_{n=0}^{N-1} x_2[n]^2 + \sum_{n=0}^{N-1} 2x_1[n] x_2[n])$

For uncorrelated signals the last sum must be zero. This leads to:

$ {1 \over N} \sum_{n=0}^{N-1}  ={1 \over N} ( \sum_{n=0}^{N-1} x_1[n]^2 + \sum_{n=0}^{N-1} x_2[n]^2)$

The same would work for more than two uncorrelated signals.


Is inverse matched filter useful? I am aware of matched filter and its application. Now, wondering if there is any application of inverse matched filter? What I mean my inverse matched filter is that convolution of matched filter and the inverse matched filter would lead to close to delta function.

In your question, you postulate the existence of two filters, let's say $p(t)$ and $g(t)$, such that $p(t) \star g(t)=\delta(t)$. This assumption is problematic because it implies that $G(f)=1/P(f)$, which requires $P(f) \neq 0$ for all $f$. Furthermore, for small values of $P(f)$, $G(f)$ will take arbitrarily large values (see Dilip Sarwate's comments above and below this answer).

Assuming (for the sake of discussion) the existance of such a pair of filters, there may not be any advantage to using them. I'll describe a possible application from digital communications. Let's say you want to transmit a number $A$ over an analog channel. You choose an appropriate analog pulse $p(t)$ and transmit $$s(t) = A\delta(t)\star p(t).$$ (The pulse $p(t)$ could be chosen to fit the channel bandwidth, and/or to have a certain energy). The receiver could perform the "inverse filter" operation on the signal $s(t)$: $$s(t) \star g(t)=A\delta(t),$$ where $g(t)$ is your "inverse matched filter". However, in practice it turns out that detecting $A$ in this way is not optimal, because it does not have the best signal-to-noise ratio. In other words, when you add noise to the transmitted signal, the received signal becomes $$r(t)=s(t)+n(t).$$ In this case, filtering $r(t)$ with the "inverse matched filter" $g(t)$ results in a worse signal-to-noise ratio than filtering with the filter matched to $p(t)$.

If you have access to it, I highly recommend "An Introduction to Matched Filters", by G. Turin, IRE Transactions on Information Theory, June 1960. Your question is answered in page 318.


Counter design for audio signals I am designing a schematic in which a Counter is to be increased whenever a particular band of frequencies is received at the Audio input.
We are receiving the audio input as electrical Signal and I have already designed the bandpass filter with the required passband. How can I design the Counter for the same purpose?

def doppler(x):
    """
    Parameters
    ----------
    x : array-like
        Domain of x is in (0,1]

    """
    if not np.all((x >= 0) & (x <= 1)):
        raise ValueError("Domain of doppler is x in (0,1]")
    return np.sqrt(x*(1-x))*np.sin((2.1*np.pi)/(x+.05))

def blocks(x):
    """
    Piecewise constant function with jumps at t.

    Constant scaler is not present in Donoho and Johnstone.
    """
    K = lambda x : (1 + np.sign(x))/2.
    t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    h = np.array([[4, -5, 3, -4, 5, -4.2, 2.1, 4.3, -3.1, 2.1, -4.2]]).T
    return 3.655606 * np.sum(h*K(x-t), axis=0)

def bumps(x):
    """
    A sum of bumps with locations t at the same places as jumps in blocks.
    The heights h and widths s vary and the individual bumps are of the
    form K(t) = 1/(1+|x|)**4
    """
    K = lambda x : (1. + np.abs(x)) ** -4.
    t = np.array([[.1, .13, .15, .23, .25, .4, .44, .65, .76, .78, .81]]).T
    h = np.array([[4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 2.1, 4.2]]).T
    w = np.array([[.005, .005, .006, .01, .01, .03, .01, .01, .005, .008, .005]]).T
    return np.sum(h*K((x-t)/w), axis=0)

def heavisine(x):
    """
    Sinusoid of period 1 with two jumps at t = .3 and .72
    """
    return 4 * np.sin(4*np.pi*x) - np.sign(x - .3) - np.sign(.72 - x)
In [5]:
x = np.linspace(0,1,2**11)
dop = doppler(x)
blk = blocks(x)
bmp = bumps(x)
hsin = heavisine(x)
In [6]:
fig, axes = plt.subplots(2, 2, figsize=(10,10))
ax1 = axes[0,0]
ax2 = axes[0,1]
ax3 = axes[1,0]
ax4 = axes[1,1]

ax1.plot(x,dop)
ax1.set_title("Doppler")

ax2.plot(x,blk)
ax2.set_title("Blocks")

ax3.plot(x,bmp)
ax3.set_title("Bumps")

ax4.set_title("HeaviSine")
ax4.plot(x,hsin)

for ax in fig.axes:
    ax.tick_params(labelbottom=False, labelleft=False, bottom=False,
                   top=False, left=False, right=False)

fig.tight_layout();


discrete sine transform using FFT i have made a subroutine to calculate 1d fft of a series of size 2^n. I need to calculate discrete sine transform of a series using the fft subroutine. I am not sure what to do.

well, you can do a little Googling and Wikipeding and find the definitions of the DST and pick one of them.  The get the definition of the DFT (the FFT is an efficient method to compute the DFT), use Euler's formula on the DST definition and see where a DFT pops out.

i might get back to this tonight and spell it out, but if you're not sure where to go next, the above (hit the links) tells you where to approach the problem.


Kiss fft does not work after giving it more than 32 samples I am trying to take data from an accelerometer and apply Kiss FFT to the samples. I'm using a Freescale Kinetis FRDM-K22F board. I want to use 64 samples, but when I run the program I get an error saying "kiss fft usage error: improper alloc" I started turning down the sample size and saw that the FFT does work with 32 samples, but giving it 33 samples the program just stops and returns no errors. Giving it any more samples gives similar results.

I played around with how I set up the FFT and followed a few websites and forum posts:

KissFFT output of kiss_fftr

http://digiphd.com/programming-reconstruction-fast-fourier-transform-real-signal-kiss-fft-libraries/

Kiss FFT on a dsPIC33

From what I can see, I haven't done anything different from what the above websites and forums have done. I've included my code below. Any help or advice is greatly appreciated.

void Sample_RUN()
{
    int size = 33;
    kiss_fft_scalar zero;
    memset(&zero,0,sizeof(zero));
    kiss_fft_cpx fft_in[size];
    kiss_fft_cpx fft_out[size];

    kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,NULL,NULL);

    signed short samples[size];

    for (int i = 0; i < size; i++) {
        fft_in[i].r = zero;
        fft_in[i].i = zero;
        fft_out[i].r = zero;
        fft_out[i].i = zero;
    }

    printf("Data Collection Begins \r\n");

for(int j = 0; j < size; j++)
{
    for(;;)
    {
        dr_status = My_I2C_ReadByte(STATUS_REG);
        dr_status = (dr_status & 0x04);
        if (dr_status == 0x04)
        {
            //READING FROM ACCEL OUTPUT DATA REGISTERS
            AccelData[0] = My_I2C_ReadByte(OUT_X_MSB_REG);
            AccelData[1] = My_I2C_ReadByte(OUT_X_LSB_REG);
            AccelData[2] = My_I2C_ReadByte(OUT_Y_MSB_REG);
            AccelData[3] = My_I2C_ReadByte(OUT_Y_LSB_REG);
            AccelData[4] = My_I2C_ReadByte(OUT_Z_MSB_REG);
            AccelData[5] = My_I2C_ReadByte(OUT_Z_LSB_REG);

            // 14-bit accelerometer data
            Xout_Accel_14_bit = ((signed short) (AccelData[0]<<8 | AccelData[1])) >> 2;     // Compute 16-bit X-axis acceleration output value
            Yout_Accel_14_bit = ((signed short) (AccelData[2]<<8 | AccelData[3])) >> 2;     // Compute 16-bit Y-axis acceleration output value
            Zout_Accel_14_bit = ((signed short) (AccelData[4]<<8 | AccelData[5])) >> 2;     // Compute 16-bit Z-axis acceleration output value

            mag_accel = sqrt(pow(Xout_Accel_14_bit, 2) + pow(Yout_Accel_14_bit, 2) + pow(Zout_Accel_14_bit, 2) );

            printf("%d \r\n", mag_accel);

            samples[j] = mag_accel;

            break;

        } // end if
    } // end infinite for
} // end for

for (int j = 0; j < size; j++)
{
     fft_in[j].r = samples[j];
     fft_in[j].i = zero;
     fft_out[j].r = zero;
     fft_out[j].i = zero;
}

printf("Executing FFT\r\n");
kiss_fftr(fft, (kiss_fft_scalar*) fft_in,  fft_out);

printf("Printing FFT Outputs\r\n");
for(int j = 0; j < size; j++)
{
    printf("%d \r\n", fft_out[j].r);
}

kiss_fft_cleanup();
free(fft);
} // end Sample_RUN

It is possible that the program expects a power of 2 samples. Does it work with 64 samples? The FFT by definition usually refers to power of two series.


Why do we choose 44.1 kHz as recording sampling rate? Peoples' ears can hear sound whose frequencies range from 20 Hz to 20 kHz. Based on the Nyquist theorem, the recording rate should be at least 40 kHz. Is it the reason for choosing 44.1 kHz?

It is true that, like any convention, the choice of 44.1 kHz is sort of a historical accident. There are a few other historical reasons.

Of course, the sampling rate must exceed 40 kHz if you want high quality audio with a bandwidth of 20 kHz.

There was discussion of making it 48.0 kHz (it was nicely congruent with 24 frame/second films and the ostensible 30 frames/second in North American TV), but given the physical size of 120 mm, there was a limit to how much data the CD could hold, and given that an error detection and correction scheme was needed and that requires some redundancy in data, the amount of logical data the CD could store (about 700 MB) is about half of the amount of physical data. Given all of that, at the rate of 48 kHz, we were told that it could not hold all of Beethoven's 9th, but that it could hold the entire 9th on one disc at a slightly slower rate.  So 48 kHz is out.

Still, why 44.1 and not 44.0 or 45.0 kHz or some nice round number?

Then at the time, there existed a product in the late 1970s called the Sony F1 that was designed to record digital audio onto readily-available video tape (Betamax, not VHS). That was at 44.1 kHz (or more precisely 44.056 kHz).  So this would make it easy to transfer recordings, without resampling and interpolation, from the F1 to CD or in the other direction.

My understanding of how it gets there is that the horizontal scan rate of NTSC TV was 15.750 kHz and 44.1 kHz is exactly 2.8 times that. I'm not entirely sure, but I believe what that means is that you can have three stereo sample pairs per horizontal line, and for every 5 lines, where you would normally have 15 samples, there are 14 samples plus one additional sample for some parity check or redundancy in the F1.  14 samples for 5 lines is the same as 2.8 samples per horizontal line and with 15,750 lines per second, that comes out to be 44,100 samples per second.

Now, since color TV was introduced, they had to bump down slightly the horizontal line rate to 15734 lines per second. That adjustment leads to the 44,056 samples per second in the Sony F1.

Look [http://batmobile.blogs.ilrt.org/audio-analysis-on-an-iphone for description
.
 .A theorem called the Nyquist sampling theorem states that in order to sample a signal of X Hz without significant loss of quality, you need to sample at 2X the frequency.  The limit of human hearing is approximately 20kHz, which hence requires a sample rate of approximately 40Khz.  This is why CDs are sampled at 44Khz.  i.e. each second of recording in a CD contains 44,000 measurements of the highest possible frequency contained in the recording.


Detecting blobs that connects to any other blob, maybe with OpenCV In the image, I connects to universe, but II and III not.

I need to detect both II and III, also I if possible.

Is it possible with current computer vision libraries?

Flood fill  the background, then hough transform to pick out rough circles II and III