Increasing SNR and Dynamic Range using Oversampling How much gain in dynamic range and SNR can be expected if we are to oversample a signal with fixed analog input bandwidth. For Example if I have a analog filter at the input which limits the bandwidth to 100 MHz. The noise power will also be limited to the 100 MHz bandwidth ($kTB$ formula for noise).
Now if I were to oversample the data by a factor of 2 how much gain should I expect and why? I know we get 3 dB gain due to spreading of quantization noise, will a similar gain be achieved also due to limiting the noise bandwidth at the input?

OP clarified that the question in the comments as follows:

If we ignore any modulation for now and assume that we are receiving
pure tones plus the band limited noise and we try to improve the SNR
in post processing how much improvement can we expect by oversampling
and is there a limit to it? My original question was about this
aspect.

First consider the case of sampled white noise (no band-limiting of the sampled signal other than anti-alias filtering prior to the A/D converter), the improvement in SNR due to oversampling a pure tone would be the processing gain in a matched filter such as a correlator to $e^{j\omega n}$ given as $\sum x[n]e^{-j\omega n}$:
$$PG = 10\log_{10}(N) \text{ dB}$$
This is because as we sum the samples in the correlation process, the signal component of each sample which is correlated from sample to sample will increase at $20\log(N)$ dB while the noise components will increase at $10\log(N)$ dB: The total signal magnitude would be $x_1+x_2+x_3 +\ldots$ while the total noise rms magnitude would be $\sqrt{x_1^2+x_2^2 + x_3^2 + \ldots}$
But with band-limiting, as we increase the sampling rate beyond the Nyquist frequency for the bandwidth of the signal, there will be no further increase in SNR: since the signal is bandlimited (no longer white), adjacent samples will be correlated both in signal and in noise, and as we sum the samples in an attempt to increase SNR, both the signal and noise components of the signal will increase at the same rate ($20\log_{10}(N)$).
Assuming a constant noise out to $f_b$ and then rolling off after that, we would increase SNR after correlation up to the point where $f_s/2 = f_b$ (for a real signal) Increasing the sampling rate beyond that will not increase SNR. Further we have a practical limit for long averaging situations from ADC spurious free dynamic range (SFDR), phase noise from clocks and local oscillators, and similar $1/f$ noise sources where the signal is no longer stationary and increasing the averaging duration further will begin to degrade SNR.  The Allan Deviation (ADEV) is an excellent statistical tool for determining the optimum duration to average given we have no other constraints on the time interval. (For more on ADEV for this application see What determines the accuracy of the phase result in a DFT bin? )

If the input signal has a fixed analog bandwidth, and you are properly not limited by the quantization noise, and front-end filtering has been sufficient to not have aliasing to occur, then oversampling will not have any effect on the SNR in band.
The idea of oversampling is to spread the quantization noise out so that it is no longer swamping out the noise floor of your signal which should be the limiting factor to SNR, and to allow for appropriate filter design to eliminate noise out of band from folding in at a lower sampling rate. These are the two most prominent reasons to oversample: increasing the effective resolution by spreading out the same total quantization noise over the full sampling bandwidth, and increasing the total frequency space in which we can do digital radio processing of signals outside our primary band of interest, for filter rejection and multi-carrier operations. For more details on the first and it's limitations please see:  What are advantages of having higher sampling rate of a signal?
It is assumed that the proper receiver design will filter out all noise that is outside of your signal bandwidth (the matched filter), so given you are sampling sufficiently high according to Nyquist (2x the bandwidth, plus some margin for realistic filter realization) then there is no reason to sample any higher once the signal is properly filtered (and this only increases power dissipation, increased resources, etc).
The strategy should be to sample sufficiently high at the front-end of the digital receiver for the filter designs involved and in consideration with the analog filtering that may be done for anti-aliasing, then once out of band noise, interference and other channels are filtered out, reduce the sampling rate to be as small as possible (typically 2 samples/symbol).
This may add a further intuitive explanation: Consider a band-limited random noise process with the sampled time domain signal shown below along with a histogram of those samples along the right axis, and then immediately below that is the spectrum. As long as we sample that spectrum more than twice the bandwidth shown, no matter how much faster we sample we will still get the same histogram, from which we get the standard deviation of that noise.
If we had additional noise energy outside the signal's bandwidth, then oversampling can of course help us to eliminate that interference through filtering where it could otherwise fold into band with a lower sampling rate-- but this isn't improving SNR by oversampling, this would be implementing proper receiver design by filtering out interference (which that task alone may require a higher sampling rate).


LPF and HPF Confusion - Rectangular and Circular Cut-Offs Interpreting some homework, I'm a little confused on what it's asking for. I'm asked to filter in the frequency domain (easy enough, just a fft2 and ifft2 when I'm done) a 256x256 8bit grayscale tiff with the following filters:

"Rectangular low pass filter using cutoff frequency of
(uc=N/4,vc=M/4) where N and M are the row and column sizes of the
image."
"Circular low pass filter using cutoff frequency of Circular
filter the cutoff frequency is given by Rc=N/4= sqrt[(uc)^2+(vc)^2]
where N and M are the row and column sizes of the image"
"Rectangular high pass filter that has the same cutoff as part a."
"Circular high pass filter that has the same cutoff as part b."

Now, I understand some of this well enough, but the "Rectangular/Circular" throw me.  We didn't really cover anything like that.  What we did cover were the ideal low/high pass filters and the butterworth low/high pass filters.  Is this perhaps a synonym for the two, the ideal with its immediate cutoff and the buttworth with its tapering? My book does not mention anything in this regard (though perhaps I have not found it...the book is verbose and I am an impatient student at times), only that the general form it gives for the Ideal filters is a circular cutoff, which is what prompts my second supposition: that the circular is the method the book uses, and the rectangluar means to filter anything inside a submatrix one quarter the size of my image, centered, in the frequency domain.

Or are we talking two different kinds of ideal LPFs and HPFs,

Yes, this is homework for a course, but I'm not looking for easy answers.  I'm looking for understanding.

You're probably best off asking your professor for clarification. As you reasoned, it sounds like you're being asked to multiply the image in the frequency domain by either a rectangular or circular mask that eliminates all frequencies outside the mask. This is not typically how you would apply a filter, however, as "ideal" frequency-domain filter masks have unsatisfactory time-domain properties (e.g. a truly ideal rectangular filter has an infinitely long impulse response).

As you noted, filters with a tapered response (of which the Butterworth filter is an example) are more often used. However, the problem specification doesn't seem to contain the information that you would need to design such a filter, such as the type of filter and its order.


Discrete Harmonics - Why multiplying digital frequency by k does not get next harmonic For continuous time $ e^{jk\Omega_0t} $ gives a complete set of orthogonal harmonics for fourier decomposition but for discrete  $ e^{jk\omega_0n} $ does not form a complete set orthogonal basis set to decompose a signal because... ?

No Answer yet on 21 Nov 2018,

If you know please help :(

First the continuous-time:

Consider the complex exponential: $ e^{j \Omega_0 t} $ where $\Omega_0$ is the continuous-time radian frequency, $t$ is the continuous time itself. Now from elementary calculus, the period of this complex exponential can be found from $x(t) = x(t+T)$ as :

$$ e^{j \Omega_0 t} = e^{j \Omega_0 (t+T)} = e^{j \Omega_0 t} e^{j \Omega_0 T} $$

Now in order for this equality to hold, the last exponential term must be $1$. From the algebra of complex-exponentials it's known that it will be $1$ whenever $\Omega_0 T = 2\pi m$ for some integer $m$ and real $T,\Omega_0$ which indicates that

$$e^{j \Omega_0 T} = 1 \implies \Omega_0 T = 2\pi m $$ from which we find the period of the continuous-time complex exponential as:

$$ T = \frac{2 \pi m}{\Omega_0} = \frac{2 \pi}{\Omega_0} = T_0$$

where in the left we used the fact that the fundamental period $T_0$ is the minimum (nonzero) real number $T$ that satisfied the equality for $m=1$.

Conclusion: any continous-time complex exponential $e^{j \Omega_0 t}$ is periodic for any real value of $\Omega_0$... This is so because $t$ is a continuous variable which admits the period $T_0$ to be continuous as well. Since the period is allowed to be continous, it can always be found to make $ \Omega_0 T$ an integer multiple of $2\pi$.

Furthermore, in the continous-time case the harmonic family of the complex exponential $e^{j \Omega_0 t}$ is defined to be:
$$\phi_k(t) = e^{j k \Omega_0 t}$$ for integer $k = 1,2,...,\infty $.

The particular period associated with the k-th harmonic $\phi_k(t) = e^{j k \Omega_0 t}$ is $T_k = \frac{T_0}{k}$, nevertheless its fundamental period is $T_0$.

Now, since the period in continous-time is a real variable it can take any value possible, as you can see, as the harmonic member index $k$ increases the particular period decreases like $T_k = \frac{T_0}{k}$. As the member index $k$ goes to infinity, the member period goes to zero but is a valid value. Hence we see that there are an infinite number of such members; i.e., the harmonic family in continuous-time has infinite members. for each $k=1,2,... , \infty$

Part-II: The discrete-time case :

The complex exponential is $e^{j \omega_0 n}$, where $n$ is an integer and $\omega_0$ is the discrete-time radian frequency (real variable).The headache comes beacuse the index $n$ is not continous but an integer, and therefore, admits only integer periods $N=1,2,...$ . The smallest allowed period is $1$ (unlike in the CT case where the smallest period goes to zero taking any real value)

Now for a discrete-time complex expoeential to be periodic, you should have $x[n] = x[n+N]$ for some integer $N$ as shown:

$$ e^{j \omega_0 n} = e^{j \omega_0 (n+N)} = e^{j \omega_0 n} e^{j \omega_0 N} $$

Again from the complex exponential algebra it's seen that :
$$ e^{j \omega_0 N} = 1 \implies \omega_0 N = 2\pi m $$ for some integer $N$ and $m$ . This implies that :

$$ \omega_0 N = 2\pi m  \implies N = \frac{ 2 \pi m }{\omega_0} $$

Now, in order for a discrete-time complex exponential $e^{j \omega_0 n}$ to be periodic with period $N$, it must be true that its frequency $w_0$ must be a rational multiple of $\pi$ as indicated by :

$$ \omega_0 = \frac{2 \pi m}{N} $$

Hence if $\omega_0$ does not satisfy the above condition that that complex exponential cannot be periodic at all. This is never the case with continuous-time complex exponentials which are always periodic for any $\Omega_0$.

After finding the condition on the periodicity of the discrete-time complex exponential, let's also observe their harmonic family:

In the discrete-time case, the harmonic family of the periodic complex exponential $e^{j \omega n} $ with a period of $N$ is defined to be:
$$\phi_k[n] = e^{j k \omega_0 n}$$ for integer $k = 0,1,...,N-1$.

Another difference between the CT and DT harmonic families occur on the particular period $N_k$ associated with the k-th harmonic: it is not $N_k = \frac{N}{k}$, as that won't be an integer for any $N$ and $k$. The correct period for the k-th member is $N_k = \frac{N m}{k}$ where integer $m$ is chosen to make $N_k$ the minimum integer for given integer $N$ and $k$

Now the last point, why are there just a finite number of harmonics in DT unlike the CT where there are infinite harmonics? This can most easily explained by the following:  Let integer $k = mN + r$ for some integers $r,m$ and let $r < N$ where $N$ is the period of the fundamental member $e^{j \omega_0 n}$. Then we have the k-the member as:
$$ \phi_{k}[n] = \phi_{m N+r}[n] = e^{j \omega_0 (mN + r ) n} = e^{j \omega_0 m N n } e^{j \omega_0 rn } $$

Now since $N$ is the period, we have $e^{j \omega_0 m N n } = e^{j 2\pi l m N n} = 1$ and therefore we see that

$$ \phi_{k}[n] = \phi_{m N+r}[n] = e^{j \omega_0 rn } = \phi_{r}[n] $$

The $k = mN + r$ th member is identical to the $0\leq r \leq N-1$ st member. Hence we conlude that in discrete-time, for a periodic complex exponential $e^{j\omega_0 n}$ of period $N$ there are only $N$ distinguishable members of the harmonic family.

I hope this clarified it.


Spectrogram Interpretation I collected vibration data (amplitude and time) on a component that was attached to an engine. Using MATLAB I plotted a spectrogram to see what frequencies it vibrated at.

The code:

windsize = 1000;
beta = 25;
overlap = 900;
thresh = -140;
spectrogram(a,                      ...
            kaiser(windsize, beta), ...
            overlap,                ...
            windsize,               ...
            'MinTreshold', thresh,  ...
            'yaxis', fs             ...
);

I am not sure how to interpret the results of the plot. Order analysis of the engine expected a maximum output of 100 Hz, where as here the strongest signals are at the higher frequencies and at high harmonics.

Have I done something incorrectly?

Abstracting this a bit:

This doesn't look too bad for a signal!

So, first of all, you might be driving your sensor and digitizer too hard – artifacts due to nonlinearity would be at multiples of the fundamental tone! Idea is the following: If the transfer function of my system isn't $y = \alpha x$ but (also contains) a function $y=\alpha x^2$, then a sine in $x$ would lead to two new signals: one at 0 Hz, and one at twice the original sine's frequency. That's intermodulation as understandable by the classical trigonometric addition / multiplication theorems!

So, first thing: make sure your system is nice and linear over the whole range; for systems that take vibrations, exciting the thing with a large acceleration of known amplitude might be hard, so a first start would be plotting a histogram of the amplitude of your signal. Does it look nice and somewhat smooth with few values at the extreme values and most at low amplitudes?

If that's not the case, you might simply need to use a sensor / sensing system that can work with a larger range of values. That might just mean your sensor is working fine and linearly, but your ADC needs to get the values attenuated to map them to the whole range (you'd see that as sharp clipping if you plot your signal over time!); but based on your relatively benign "ghost signals", my guess is that your acceleration sensor simply isn't that linear. Maybe it's not meant for that range of acceleration, maybe you're not supplying it with a supply current or voltage high enough to represent these values? Maybe it's not even specified to do anything linear in frequency ranges above a certain cutoff frequency?

Another thing: Your plot goes up to what I presume is half the sampling rate, and I see no cutoff of any filter that was applied to the signal.

Are you properly anti-aliasing filtering your signal prior to sampling it? If your signal contains e.g. the oscillation of a MEMS fork at 5 kHz (wild guess!), but you're sampling at 4.8 kHz, then that oscillation would appear at 0.2 kHz, and a lot of multiples of 0.2 kHz. Do not mess with Nyquist unless you want aliases! Quite possible, the accelerometer actually generates a signal that contains the intermodulation of MEMS oscillation and signal to be observed:

$$\sin(2\pi f_{machine} t+\varphi_{machine}) \cdot \sin(2\pi f_{oscillator} t)$$

Which in turn would be

$$\frac12 \left(\cos(2\pi(f_{machine}-f_{oscillator})+\varphi_{machine})-\cos(\cos(2\pi(f_{machine}-f_{oscillator})\right)\text.$$

Staying with the number from my example, $f_{oscillator}=5\,\text{kHz}$ and your signal at 100 Hz, we'd thus should see components at 4.9 kHz and 5.1 kHz – except we don't, because we're not sampling fast enough, so at our 4.8 kHz sampling rate, these would end up at 100 and 300 Hz; and then, combine that with the non-linearity intermodulation, you get the funkiest frequency content.

So, first of all, if you know you're not expecting oscillations above 100 Hz, use an analog filter to get rid of anything above. If you still get signal components above that frequency, your filter isn't sufficient or your ADC is driven into nonlinearity.

On the other hand, without any education on the mechanical side of things: doesn't seem so unlikely to me that materials aren't perfectly rigid, and that something oscillating at 100 Hz might be "fluttering" at higher frequencies, especially if things like spring elements aren't the perfect linear springs I was told to assume in school physics.

In conclusion:

You're looking at something non-linear in nature. Whether it's your measurement system or the system you're observing is something we can't tell you, but systematic elimination of error sources will probably yield results pretty quickly.

Start with reducing sensitivity of your measurement system.

Also, make double sure you're properly anti-alias-filtering your signal prior to digitization.


Matlab: plotting magnitude and phase response I've got the frequency response $H(e^{j\omega}) = \frac{2 + 3.125e^{-2j\omega}}{1-0.9e^{-j\omega}+0.81e^{-2j\omega}}$, and I'm trying to figure out how to plot the magnitude response and the phase response in Matlab. Is there a simple way to do this?

Thank you in advance. This is my first post here, so please tell me if there's anything I forgot to add.

Hal

@Hal: Unlike MBaz, I realize you weren't asking for information on the theory of signal processing. But, rather, you were asking a Matlab implementation question. I only have an old version of Matlab so I couldn't try Anthony's code. (Anthony was smart to insert that zero-valued coefficient.) But to plot the frequency mag. and phase responses of a discrete version of your filter, try these:

figure(1), freqz([2,0,3.125],[1,-0.9,0.81]);

or

freqz([2,0,3.125],[1,-0.9,0.81], 512, 'whole');


Mutual Information, is this calculation correct or should I contact the author about fixing it? Background information

According to Cover's text book on information theory  Mutual Information is calculated as:

$$(1)  I(W;C) = \Sigma_k\Sigma_i P(C_k,W_i)log(\frac{P(C_k,W_i)}{P(C_k)P(W_i)})
$$

Now if we were interested in a single Word W_i we would just calculate this with one sum:

$$(2)   I(W_i;C) = \Sigma_k P(C_k,W_i)log(\frac{P(C_k,W_i)}{P(C_k)P(W_i)})
$$

OR we can rewrite this as

$$ (3)   I(W_i;C) = \Sigma_k P(C_k,W_i)log(\frac{P(W_i|C_k)}{P(W_i)})
$$

because
$$ (4)  P(W_i|C_k) = \frac{P(W_i,C_k)}{P(C_k)}$$ (unless i am incorrect)

Question
While reading a paper on feature selection using information theory measurements I came across A segment of the text:

Should the equation be more like (3), and I should contact the Author, or have I mis-understood/interpreted mutual information?

No there is nothing wrong with the formulas in the paper.

MI(S) is the average mutual information of set S over all classes c in C. It is NOT a reformulation of the formula for MI(W,C).
MI(w) is identical to your formula for I(Wi,C)  (except it uses slightly different notation).


How does a matched filter know about input signal beforehand? I know that a matched filter at time t0 is just the input signal reversed in time about that time. But how does the filter know about input signal. Received signal is transmitted signal plus the noise and if filter already has knowledge about input signal then why is there any need of filtering? You already know the input signal!

A matched filter is used when the shape of the "desired" part of the received signal is known. Under the additive white Gaussian noise (AWGN) channel, a matched filter is the optimal detection scheme in the maximum-likelihood sense. This is a commonly-modeled case, so matched filtering is a typical signal-processing operation.

Based on your question, I think you're misunderstanding what a matched filter is: it doesn't require full knowledge of the exact transmitted signal. Instead, it helps to sift out the transmitted signal from the noise that is received with it, allowing the receiver to make more accurate decisions on what was actually transmitted.

Consider a baseband communication system that transmits a rectangular pulse with amplitude $1$ to communicate a zero and $-1$ to communicate a one. The appropriate matched filter in this case would be a rectangular pulse that matches the shape that is used by the receiver.
When the train of transmitted rectangular pulses is passed through the receiver's rectangular matched filter, the result is a train of triangular pulses (the autocorrelation of a rectangular pulse is a triangle pulse).
The receiver then makes symbol decisions based on the amplitudes of the triangular peaks that it observes. At this point, the SNR is maximized, as the matched filter tends to pass the signal of interest while averaging out any additive noise.

This is a simple, but representative example. The matched filter doesn't carry any of the transmitted information in itself (as you noted, that would be nonsensical); instead, it is used as an optimal means of extracting the transmitted information from the observed waveform.

I would also call your attention to this very good answer to a previous question about matched filters. It contains many of the diagrams that would make the above discussion easier to understand but I didn't have time to produce. It's a great resource for additional details.


Finding the input bit rate of an analogue input Suppose I have 4 of 300bps synchronous digital input and 1 of analogue input with a bandwidth of 500Hz. The analogue input is sampled at 4bit using PCM. All the channels are multiplexed with TDM.

As one of the input is of analogue, I need to convert that to digital input first before TDM can be carried out. Since the analogue signal will be sampled at 4bits, my initial idea was to use the Nyquist sampling theorem to have 500Hz × 2 = 1000 samplings per second. But on a second thought, this is merely the bandwidth and I cannot do this unless it is said that the highest frequency is 500Hz.

However, then how should I find out the input bit rate of the analogue input channel for the TDM multiplexing?

Theoretically speaking you need 2 samples per second for 1 Hz of bandwidth (regardless whether the signal is centered around DC or not), so in this case a sample rate of 1 kHz would be enough.
Practically speaking, you always need more since the theoretical requirements is based on the assumption of an ideal low pass (or bandpass) filter.
The correct choice of sample rate depends on how the band limiting is achieved and what bandwidth you are really interested in and what noise and sampling accuracy you need. Most audio systems use 44.1 kHz or 48 kHz to achieve 20 kHz usable bandwidth.
Another consideration is that you already have some digital inputs. These are derived from some master clock anyway. If you can use the same clock to sample your input, you can make the subsequent processing a lot easier.


Transfer function estimation from frequency response Let's assume that we know that we are dealing with a SISO second order system for which we have the frequency response (magnitude and phase for a known frequency range ω). What methods would people use to fit the frequency response to a transfer function (i.e. transfer-function estimation)? How does this process look like?
Also, this looks like a one-liner in Matlab (see tfest). Is there a Python equivalent?

Here you find an extensive discussion about transfer function estimation and even source code.
Your problem can be expressed as $F(2j \pi f_i) \approx g_i$, where $f_i$, $g_i$ are your measurements, and you also you  can expresss $F(s) = N(s) / D(s)$ as a parametric function, then the parameters may be adjusted with curve fit.
Python provides curve_fit, that can be directly applied to this problem as in this answer, with the difference that you want to apply for x in the imaginary axis. This can be improved by also passing the gradient in terms of the parameters as well.
Maybe you prefer to express your function in terms of poles and zeros, this may be specially convenient if you want to ensure stability (you can add to your cost function the Lagrangian Multiplier for $real(p) < 0$)
def pztf(omega, h, p, z):
  num = np.prod(1j*omega[None,:] - p[:, None], axis=0);
  den = np.prod(1j*omega[None,:] - z[: None], axis=0);
  return h * num ./ den

Maybe you want to fit the logarithm of this, this may produce a better fit in the low gain frequencies.
There are also more analytic approaches. The so called Levy method is interesting,
$F(s) = \frac{N(s)}{D(s)}$ is equivalent to $F(s) D(s) - N(s) = 0$ this can be expressed as least square fitting.
This is biased in the sense that the errors close to the poles receives less weight, this can be mitigated by solving using $D_0(s) = 1$, and iteratively refining it by iteratively fitting
$$\frac{1}{D_{i-1}(s)}\left( F(s) D_i(s) - N_i(s) \right) = 0$$
when it converges after some iterations we have $D_{i-1}(s) \approx D_i(s)$ we have
$$\frac{1}{D_{i-1}(s)}\left( F(s) - N_i(s) \right) \approx F(s) - \frac{N_i(s)}{D_{i-1}(s)}$$
That corresponds to the original problem (without the $D(s)$ factor).
If you want better accuracy in different frequencies. You achieve this by multiplying the whole expression by a given $W(s)$, then te fitting iteration becomes
$$\frac{W(s)}{D_{i-1}(s)}\left( F(s) D_i(s) - N_i(s) \right) = 0$$


Color and circle detection in image I try to detect and count circles in image (for example smarties)

I use HSL color space. But I am not able to distiguish colors od the same color if they are in touch. I try to erode and dilate picture. But result is the same, I have only some blobs(connected components) od the same color. Do you have some general algorithm for this problem?

(I try to do that with EMGUCV library)

Thanks in advance.

There is no need to Dilate. Here is my segmentation results:

And here is my code using OpenCV:

#include <iostream>
using namespace std;
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
using namespace cv;

int main()
{
        int ht,wt;
    IplImage *im = cvLoadImage("2Dpg6.jpg");
        ht=im->height,wt=im->width;
    cvShowImage("RGB",im);

    IplImage *img = cvCreateImage(cvGetSize(im),8,1);
    cvCvtColor(im,img,CV_RGB2GRAY);

    CvRect rec;
    CvSeq *contours;
    CvMemStorage *storage=cvCreateMemStorage(0);
    IplImage *img2 = cvCreateImage(cvGetSize(im),8,1);

    cvThreshold(img,img2,30,255,CV_THRESH_BINARY);  // image segmentation
    cvErode(img2,img2,NULL,7);  // erode
    cvShowImage("BW",img2);  // show black-wihte image
    cvFindContours(img2,storage,&contours,sizeof(CvContour),CV_RETR_CCOMP);  // find contours
    while(contours){
        rec = cvBoundingRect(contours,1);
        cvRectangleR(img,rec,cvScalar(255));
        contours=contours->h_next;
    }
    cvShowImage("Gray",img);

    cvWaitKey(0);
    return 0;
}


Software to register a 2D image to a 3D volume I have been searching for software that can handle the task of registering a 2D image to a 3D volume.  The 3D volume is a microCT scan, while the 2D image is taken from a plane passing through the sample.  I am looking for something that can locate the plane of the 2D image within this 3D volume, i.e. get the 3D translation, rotation and scale transform betwee the image and volume.

First prize is for pre-compiled free Windows software, but if that isn't possible then any libraries that can handle this task.  I see that ITK from Kitware provides registration tools, but I don't see any discussion of my specific task, which doesn't seem to a common one.

The software Amira from Visage Imaging has some excellent built-in registration features. I never tried that exact problem, though.

It is unfortunately pretty pricy, but you can download a free trial version of it from the website to give it a try. http://www.vsg3d.com/amira/overview

I also found some well implemented registration algorithms in the software MedInria, an open source and free software specially developed for medical image analysis. http://med.inria.fr


Why do we say that "zero-padding doesn't really increase frequency resolution" Here is a sinusoid of frequency f = 236.4 Hz (it is 10 milliseconds long; it has N=441 points at sampling rate fs=44100Hz) and its DFT, without zero-padding :

The only conclusion we can give by looking at the DFT is: "The frequency is approximatively 200Hz".

Here is the signal and its DFT, with a large zero-padding :

Now we can give a much more precise conclusion : "By looking carefully at the maximum of the spectrum, I can estimate the frequency 236Hz" (I zoomed and found the maximum is near 236).

My question is : why do we say that "zero-padding doesn't increase resolution" ? (I have seen this sentence very often, then they say "it only adds interpolation")

=> With my example, zero-padding helped me to find the right frequency with a more precise resolution !

Resolution has a very specific definition in this context. It refers to your ability to resolve two separate tones at nearby frequencies. You have increased the sample rate of your spectrum estimate, but you haven't gained any ability to discriminate between two tones that might be at, for instance, 236 Hz and 237 Hz. Instead, they will "melt together" into a single blob, no matter how much zero-padding you apply.

The solution to increasing resolution is to observe the signal for a longer time period, then use a larger DFT. This will result in main lobes whose width are inversely proportional to the DFT size, so if you observe for long enough, you can actually resolve the frequencies of multiple tones that are nearby one another.

--

To see how this plays out, here's a plot of the zoomed-in FFT of the addition of two signals: your original sinusoid, and one that differs in frequency from it by 0 to 100 Hz.

It's only towards the 100Hz difference end of the plot (left-hand side here) that you can distinguish (resolve) the two.

Scilab code for generating the plot below.

f = 236.4;
d = 10;
N=441;
fs=44100;
extra_padding = 10000;

t=[0:1/fs:(d/1000-1/fs)]
ff = [0:(N+extra_padding-1)]*fs/(N+extra_padding);

x = sin(2*%pi*f*t);

XX = [];

for delta_f = [0:100];
    y = sin(2*%pi*(f+delta_f)*t);
    FFTX = abs(fft([x+y zeros(1,extra_padding)]));
    XX = [XX; FFTX];
end

mtlb_axis([0 1300 0 500])

figure(1);
clf
[XXX,YYY] = meshgrid(ff,0:100);
mesh(XXX(1:100,[50:90]),YYY(1:100,[50:90]),XX(1:100,[50:90]))


What is the meaning of a Nyquist diagram in a vibration signal I need to analyze some vibration data from accelerometers (g-Hz) but I don't quite understand the meaning of each component that the registered signal contains. Basically it is Magnitude and Phase. I'm able to get the Real/Imaginary parts from those two as well.

Today I plotted Real vs Imaginary parts (Nyquist plot/diagram, correct me if I'm wrong). These are the results:

The Nyquist shows like a circular pattern, sometimes there are spirals, sometimes just random points. What could this mean?

The second picture displays the Magnitude of this particular accelerometer.

In the third the raw Real/Imaginary parts are included. The Magnitude is the RSS of Re/Im parts. I've realized that I need to compare (divide by) this channel to the reference (the input or control channel).

The last one shows the chaotic Phase, it gets better if I do the same as at Re/Im parts: compare to the reference, in this case I subtract in order to get the gap between the input and my channel (output).

The pink and green dashed lines are the Re/Im parts divided by the Re/Im parts of the reference (the input). The blue is the same but with Magnitude

The blue line is the subtraction.

It would be great if someone could help me to have some better understanding and interpretation of this kind of data.

By the way, is Argand diagram the same in essence as Nyquist?

EDIT: The vibration signal comes from a UFF (Universal File Format) dataset 58, the test consists in measuring the response of 50-100 accelerometers attached all around the structure to a sinusoidal motion (varying in frequency). There is an accelerometer situated in the base, so it is measuring the real input to the system (in g-Hz).

The hardware of the test transforms internally the response from g-s (time) to g-Hz with some kind of Fourier Transform algorithm I believe.

As the vibration may not be in phase with the input motion the gap between them is also measured. For example: a gap of 180 degrees between an accelerometer and the base means that if the shaker is going up in a particular instant, this accelerometer is going down.

The "Nyquist plot" you have is really just the Fourier transform of the vibration data.  That the plot is "circular" doesn't mean much: it just means that the phase of the signal is continually increasing and that the magnitude is (roughly) constant.

I wouldn't call such a plot a Nyquist plot.  Nyquist plots are plots of the transfer function of a system. You may be able to use your measured signal and the reference to attempt to generate a Nyquist plot, but it's probably not going to be just a matter of multiplying / dividing Fourier transforms of the signals.

Generally, plotting the Fourier transform of a signal on the complex plane doesn't show much useful information.  More useful is to plot the magnitude against frequency, as this shows the dominant frequencies in the signal.

The reason the phase plots are "chaotic" (sic) is that the signal is noisy.  Taking the argument of a noisy complex signal will nearly always generate discontinuous phase plots.

I believe an Argand diagram is the same as your "Nyquist plot".


Is sum of period samples in the DSP domain? Following program outputs sum of samples of signal's period. Values in program are:

Sample rate 96 ksamples/sec.
First frequency (fundamental, 1st harmonic) 96 Hz.
Second frequency (2nd harmonic) 192 Hz.
n = 1 to 1000, representing one (the first) period of the signal.

Output (sum of period's samples) should be zero. Actually is 18.4062, apparently due to round off errors. In general, when samples are multiple of fundamental's (in this example: 2000, 3000, ... samples), sum is least. Other sums are much greater. So first harmonic is represented by the first least sum.
I have some question about, but first of all, I have to ask whether this process belongs to DSP domain.

// Outputs sum of input values.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
int n = 0; // sample.
float xn = 0; //pure number.
float xn_cos = 0; //pure number.
float xn_sin = 0; //pure number.
float sum = 0;

//Sampling frequency and samples of segment.
int samplFreq = 96000; // sample/second.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 96; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 192; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

for (n = 1; n <= 1000; ++n)
{
xn = ampl_1 * sin(n * unitAngle_1 + phase_1)
   + ampl_2 * sin(n * unitAngle_2 + phase_2);
sum += xn;
}
 std::cout << sum;
return 0;
}

Ladies, Gentlemen,

Please let me say that my answer is wrong on that lower (first) frequency is represented by first least sum. By trial and error I found that lowest frequency is represented by the fourth (absolutely) least sum, indeed given that existing frequencies are two. I apologize. However, I have solved my problem and do not need any other question. My answer to my question is that this process is DSP, but I ask Messrs Administrators view so that if I have any other question about this method, ask it here.

Following C++ program outputs sums of input values (xn) and you can see that 4th least sum represents first (fundamental) existing frequency. You can change values on the condition second existing frequency be double the first. Example values are: sample rate 96000 ksamples/sec, 1st existing freq = 100 Hz, represented by 960 samples, 2nd = 200 Hz, represented by 480 samples, n = 1 to 1970 with step 11. Fourth least sum (138131) represents 1O1.373 Hz. Next frequency, 100.209 Hz, with sum -357858, is the nearest to 100 Hz. Obviously there is accuracy error. However for step n = 1, fourth absolutely least sum is -7.40625 representing 100 Hz. Eighth least sum represents half the 1st frequency, 50 Hz.  I point out that in detecting fundamental frequency of periodic signal by this method, precise is high when ratio sample rate/fundFreq is high, in opposite to method by DFT where results are bad when ratio samplFreq/existFreq is too high.

With regards and friendsip.

// Outputs sum of input values.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
float n = 0; // sample.
float xn = 0; //pure number.
float xn_cos = 0; //pure number.
float xn_sin = 0; //pure number.
float sum = 0;

//Sampling frequency and samples of segment.
int samplFreq = 96000; // sample/second.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 100; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 200; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

for (n = 1; n <= 1970; n)
{
xn = ampl_1 * sin(n * unitAngle_1 + phase_1)
   + ampl_2 * sin(n * unitAngle_2 + phase_2);
sum += xn;
std::cout << samplFreq/n << "       " << sum << "\n";
n += 11;
}
return 0;
}


Feature Detection on a Circular Metal Surface for Character Recognition Hello I'm trying to recognize characters on the circular metal surface.

Without any preprocessing I could match cropped features with  original image.
I've actually cropped 2 characters together to increase interested points on a template.

however on different locations of the conveyor belt and with different illumination angle and perspective SURF matching fails.

I will definately crop the ROI from conveyor belt but still I am not sure whether I have to create a database of images to increase accuracy or some kind of preprocessing (I've read that for feature detection there is no require preprocessing on images)

I am open to new suggestions for this problem. Maybe a different algorithm or method that can help me out.

thanks in advance.

Do not struggle forming a database of images to match via descriptors. This would be too computationally cumbersome and would require immerse amount of training. Such a scalable solution doesn't exist out of the box yet. I would rather rely on Neural Networks or SVMs to train the possible appearances of characters.

Of course using a classifier relies on the assumption that the images appear more or less in a controlled fashion and you can segment the characters. To meet these criteria, you need to stabilize your illumination conditions. You can do this with the appropriate choice of LED lighting. After the illumination is more or less uniform and robust, the character segmentation is straightforward. You can simply threshold the image through local techniques (which posses certain illumination treatment) and apply blob analysis. You can also put all the constraints you have in this stage. Don't forget to unwarp your images prior to segmentation as described in this post.

To train your classifier do not use plain intensity values. This would be very noisy. Instead exploit the gradient information (or normalized gradients). This is proven to be robust in many works such as this one. Also, don't even use the gradients directly. It is better to quantize them (Remember canny?). Robustness of LINEMOD is also a good example of this.

I would suggest further reading on machine vision as this is a practically involved topic. Many techniques on machine vision differ substantially from the methods in computer vision. You can think of this as a practice-theory duality. Refer to this book for instance.


Aliasing in Doppler Radar my book about RADAR says that:

If a Continuous wave radar sends a sine wave at frequency $f_T$ to a moving object (at speed V), a frequency $f_R$ is received. Their difference is called Doppler frequency and it is equal to $f_D = f_R - f_T = - 2 \frac{V}{\lambda}$.

It is quite clear. Then, it says that the situation is slightly different if a radar which sends rectangular pulses modulated with a sine wave at frequency $f_T$ is used. In such a case, the continuous wave is sampled by a train of strict rectangular pulses spaced by period PRT (Pulse Repetition Time):

Then, it says that:

If $PRF = \frac{1}{PRT}$ (Pulse Repetition Frequency) is less than $2f_D$, spectral aliasing will occur (Nyquist Theorem).

First Question: why is the Nyquist theorem applied to the doppler frequency $f_D$? It should be applied to the sampled sine wave, which is at frequency $f_T$.
Then, it shows the aliased spectrum:

and says that:

The spectral aliasing causes spectrum folding at $f_D > \frac{PRF}{2}$ in the main interval ranging from $-\frac{PRF}{2}$ to $+\frac{PRF}{2}$

Second Question: why does aliasing cause folding? As far as I'm concerned, spectral aliasing means that the original (correct) spectrum is replicated and those spectra overlap (and hence it is not possible to separate them correctly).

I don't understand the link between my last two pictures.

So I took a look at the slides and it seems to be giving a cursory overview of how to measure Doppler with multiple pulses, so it's very generic. Hopefully this is a pretty straight forward thing to discuss given that your questions are suspicious of their statements, that's good!

First Question: why is the Nyquist theorem applied to the doppler frequency $f_D$? It should be applied to the sampled sine wave, which
is at frequency $f_T$.

The received frequency $f_R$ is the transmitted frequency $f_T$ with the added Doppler shift $f_D$ as $f_R = f_T + f_D$. So in order to measure Doppler unambiguously, you need to be able to sample $f_R$ above Nyquist. But the question still stands as to why the text states that we must appropriately sample $f_D$ instead of $f_R$?
As is done in modern radar systems, and as implied in the text, the received pulses are mixed-down to baseband. That is, the carrier is centered around zero frequency. This means that  $f_T = 0$ and so
$$f_R=f_D$$
Where the Nyquist criterion for the PRF on sampling Doppler is now more apparent.
It is possible to undersample around carrier frequencies that are not zero and still extract the desired signal, but I don't think that's necessarily important for this discussion. It's a slightly more advanced technique that simply achieves the same thing. You can look up "bandpass sampling" for more information.

Second Question: why does aliasing cause folding? As far as I'm concerned, spectral aliasing means that the original (correct) spectrum is replicated and those spectra overlap (and hence it is not possible to separate them correctly).

As you already stated, sampling produces spectral copies that are separated apart by your sampling frequency, in this case which is the $PRF$. Aliasing causes "folding" because as the frequency being sampled goes past Nyquist, in this case $\frac{PRF}{2}$, there is a copy incoming from the left into your un-aliased window that spans $[-\frac{PRF}{2}, \frac{PRF}{2})$.
Let's do a simple example where we have the following properties:

$PRF = 50 \space kHz$

Unambiguous spectrum span: $[-\frac{PRF}{2}, \frac{PRF}{2}) = [-25 \space kHz, 25 \space kHz)$

Doppler frequency $f_D = 10 \space kHz$

Everything is fine here, we're well below the Nyquist rate and we can recover the Doppler frequency nicely:

Due to sampling, the adjacent copies to the left and right this sinusoid exist at $f_D \pm PRF$. Meaning we have two tones that exist at the frequencies
$$f_{left} = 10 \space kHz - 50 \space kHz = -40 \space kHz$$
$$f_{right} = 10 \space kHz + 50 \space kHz = 60 \space kHz$$
As expected, the copies are outside our Nyquist zone achieved by properly sampling.
Now, what if the Doppler frequency was beyond Nyquist? This is where aliasing causes "folding". Let's make the Doppler frequency be $f_D = 35 \space kHz$. The two adjacent frequencies are now
$$f_{left} =  35 \space kHz - 50 \space kHz = -15 \space kHz$$
$$f_{right} = 35 \space kHz + 50 \space kHz = 85 \space kHz$$
We can see here that the left copy has now moved, or folded, into our Nyquist zone! Then due to this aliasing, we measure the wrong Doppler:

A note on your last image:
This is a generic representation of signals that have some finite instantaneous bandwidth. In our case, we're dealing with sinusoids which have some peculiarities when it comes to aliasing since they are ideally delta functions in the frequency domain that don't have a "width" per se.


Detection of circles (ellipses) in 2D point cloud Given a set of points (2D) i.e., point cloud (PC), the question is about a robust, accurate and computing-friendly method to find circles (or ellipses in advanced version).

The intuitive idea is to use Brute-Force Search on all possible points (as center){infinite!} and radii (again infinite!). This is ultra-extremely slow and inefficient.

As demonstrated bellow each fitted circle would be ranked based on the number of points (nn) positioned on the circle-circumference in a distance shorter than a threshold (t). So there is derr to present an average distance.

In advanced form ellipses are of interest to be fitted.

Any ideas, brain-storming, experiences, comments?

The best ideas that exactly tries to solve this problem is Hough Transform .

Basically, the signal in hough space will be r, x, y co-ordinates. Here r
 stands for radius and x,y stands for center. Every points may belong to one or many circles. So in the Hough plane go through all possible circles where this point could belong to and just do a +1. This is not a search, just collection.

Now, if a real circle exists, so many points will add and the score of such a r, x, y will be much higher than all others. Selecting such a point will allow you pick the right circles.

Here is a classical paper way back in 1971 (before i was born!) that invented this concept.

USE OF THE HOUGH TRASFORMTION  TO DETECT LINES AND CURVES IN PICTURES
By:  Richard O. Duda, Peter E. Hart
Tech report Artificial Intelligence Cente April 1971.

For Tutorial i would suggest references below:

HIPR2 -Link
Amos storkey
IDL Reference

Specifically for circle detection, you can refer to this below:

AI Shack
Tech Report from Chicago Univ.
Lecture Notes from Rochester Institute

These methods is very efficient, and very computer friendly.


Unit Impulse funciton FT How is FT of $ \delta $(t) equal to 1 ? Normal FT gives the result $\infty$. Can someone please explain? I did the normal integration and substituted the limits.

Is it because  $ \delta $(t) is a unit impulse function so as it's height is large it's width is very small so no matter what the FT will always be equal to 1 ? (I'm just trying to figure out the logic)

If follows directly from the definition of the Dirac delta distribution. It is defined so that

$$\int_\mathbb{R} \delta(x) f(x) dx := f(0)$$

for any test function $f(x)$. In other words, the Dirac distribution is the generator of the linear functional that extracts a single function value.

With this definition the Fourier transform of the Dirac distribution is simply:
$$\int_\mathbb{R} \delta(t) \exp(-2\pi i \omega t) dt=\exp(-2\pi i \omega\cdot 0)=\exp(0)=1$$

It doesn't make much sense to say the Dirac distribution is infinitely high or infinitely narrow. Just use the definition given above and apply it. If you really need to, you can understand the distribution as the limit of a sequence of certain functions, but it's not a function itself. And specifically it doesn't have a graph.


salt noise in an image Suppose an image has noise, but it is exclusively salt noise. What effect will dilation have?
I am looking for a reference so I can answer this question. anyone has any idea?

To answer your 1st question, dilation will enlargen the white spots in your image over the darker spots, which may be useful if you only have tiny dark salting you would like to remove in your image.

However, you are referring to salt noise, so if the spots were also white crumbs in a dark image, then all you are doing is only making the crumbs bigger, which is what you don't want.

What I would recommend is to use either Opening or Closing, (or both) depending on your application. These are a combination of dilations and erosions so that the minor parts of the image are removed while keeping the original portions of the image intact.

For example, here is an image for a penguin (with salting):

Notice how the image has both black and white salting.

Running both an opening and closing algorithm (OpenCV) on this penguin, we remove the black and white salts, and get a nice, smoother image without the spots:

Also, here is my OpenCV code that does this:

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>

using namespace cv;
using namespace std;

int main( )
{
     Mat image,dst;
     image = imread("saltnoise.png", CV_LOAD_IMAGE_GRAYSCALE);

     // Create a structuring element
     int erosion_size;
     cout << "What element size?" << '\n';
     cin >> erosion_size;

     imshow("Salted", image);

     Mat element = getStructuringElement(cv::MORPH_ELLIPSE,
            cv::Size(2 * erosion_size + 1, 2 * erosion_size + 1),
            cv::Point(erosion_size, erosion_size) );

    //Opening and Closing
    erode(image,image,element);
    dilate(image,image,element);

    dilate(image,image,element);
    erode(image,image,element);

     imshow("Unsalted", image);

     waitKey(0);
     return 0;
}


What is the most lucid, intuitive explanation for the various FTs - CFT, DFT, DTFT and the Fourier Series? Even after having studied these for quite sometime, I tend to forget (if I'm out of touch for a while) how they are related to each other and what each stands for (since they have such similar sounding names). I'm hoping you'd come up with an explanation that is so intuitive and mathematically beautiful that they'll get embedded into my memory forever and this thread will serve as a super quick refresher whenever I (or anyone else) needs it.

I wrote this handout as a complement to Oppenheim and Willsky. Please take a look at Table 4.1 on page 14, reproduced below. (Click for larger image.) I wrote that table specifically to answer questions such as yours.

Note the similarities and differences among the four operations:

"Series": periodic in time, discrete in frequency
"Transform": aperiodic in time, continuous in frequency
"Continuous Time": continuous in time, aperiodic in frequency
"Discrete Time": discrete in time, periodic in frequency

I hope you find these notes helpful! Please feel free to distribute as you wish.


Algorithm to change an already sampled guitarstring to sound like it was a bit damped by the hand I'm trying to simulate the damping of the hand on a guitarstring on an already recorded/sampled open guitarstring sound. I've been trying to use low pass filter and had a moving frequency range but that didn't make it sound like a palm muted string, just the loss of higher frequencies.

Could someone help me find good material on this, that a human could atleast grasp a bit? It's going to be implemented in C++ and I have been searching and found almoust everything about the karplus-strong string algorithm, but that's not what I want.. I do want the damping part implemented on a sample of an already recorded real played string.

My question was first asked here.. I hope that will not be considered as crosslinking now?? https://stackoverflow.com/questions/19667049/algorithm-to-change-an-already-sampled-guitarstring-to-sound-like-it-was-a-bit-d

Here is one thing to try: Guitar wave form typically has a sharp on set (the "pluck"), followed by some ringing. You could write an onset or transient detector that finds these onsets. Once you have one, apply a fast exponential decay to the waveform to make the ringing a lot shorter. Find the next onset and repeat. You probably also want a mild low pass filter, just something with negative frequency slope to simulate that the palm damps higher harmonics more than lower harmonics.


How to Calculate the Ambiguity Function (Maybe this should be on Stackoverflow instead ...)

What is the best way to calculate the ambiguity function? There are different implementations out there for Python (like PyTFTB or from TooringAnalytics) and Matlab but there are many difficulties:

1. Python - PyTFTB

I am not sure whether the implementation of pytftb is what I am looking for because the plots for 200 Samples of my data always look similar. Here are two examples:

I'd expect more distinct plots. The plots were generated like this:

from tftb.processing import ambiguity
waf, tau, theta = ambiguity.wide_band(signal)
cntr = plt.contourf(tau, theta, np.abs(waf) ** 2, cmap="Greys", antialiased=True)
// even less difference when I use np.abs(waf) ** 2
plt.imsave("filename.jpg", np.abs(waf), cmap="jet")

2. Python - TooringAnalytics

Unfortunately I can't make this one work. Maybe this is content for another question.

3. Matlab

Is not an option. I tried my best to translate this to R but couldn't make it work.

4. R

Unfortunately I couldn't find an implementation yet in R.

I'd really like to ask Dreamcooled how he computed his ambiguity function.

Thanks for your help

Octave is a free mostly MATLAB compatible application. You shouldn’t rule it out. Also, your signals would look more alike if you removed the low frequency trend obvious in both plots. You may have a valid ambiguity function although I think you may be calculating an ambiguity on signal plus noise not the proper calculation of just signal.

There is more than one ambiguity function, this is the wide band single transmitter single receiver version. Most people think Doppler is just a frequency shift but it is actually a time dilation, although the frequency shift model is very often sufficient, particularly when looking at steady lines like molecular emission lines.

In acoustics the motion of the propagation media can also cause Doppler.

This code should be ok as long as the interpolation is accurate, which assumes that the signal has sample more than about 5 points per cycle.

 % wide band ambiguity function
    clc
    clear all
    close all
    %
    tb=-1;
    tend=1;
    dialation=.8;
    %
    %signal=[zeros(1,100) randn(1,200) zeros(1,100)];
    %signal=[zeros(1,100) exp(1j*2*pi*(4*[0:199]/200)) zeros(1,100)];
    signal=[zeros(1,100) exp(1j*2*pi*(4*[0:199]+.01*[0:199].^2)/200) zeros(1,100)];
    conj_rev_signal=conj(signal(end:-1:1));
    time=linspace(tb,tend,length(signal));
    delay_time=linspace(2*tb,2*tend,2*length(signal)-1);
    tau=linspace(dialation,1/dialation,512);
    out=zeros(length(tau),2*length(signal)-1);
    for i=1:length(tau)
    out(i,:)=conv(signal,interp1(time,conj_rev_signal,time*tau(i),'spline'));
    end
    figure(1)
    imagesc(delay_time,tau,abs(out))


Circular vs Linear Convolution When deriving DFT from DTFT,we sample the frequency domain with uniformly spaced samples,hence adding periodicity to time domain.

But DFT requires a limited support: we take only 1 period.

Does that make circular convolution the same as a linear convolution with the signal's periodic extension?

Convolution in DFT is still circular.

Think of the DFT as taking the 1st period (in time and in frequency) of the DFS (discrete Fourier series). In DFS, both the time sequence and the frequency sequence are N-periodic, and the circular convolution applies beautifully.

I personally think all properties in terms of DFS, and then consider the 1st period when speaking DFT.


Convert normalized frequency to real frequency in AR model Let us suppose that we have modeled signal using AR model, and suppose we have following model:

I used spectral estimation function from MATLAB

pyulear

Now frequencies are given in normalized frequencies and I would like to know how to convert them back to real frequencies? From there it is clear that we have four deterministic model and also plus some white noise, actually I want to know approximate frequencies in each deterministic model. I can of course determine this frequencies using  FFT, periodogram and so on, but I am studying application of the AR/ARMA model. So in my case, I have such frequencies on following figures, how can I determine actual frequencies? Thanks in advance.

Matlab generally uses normalized frequencies in radians/sample. If the signal is real-valued then these frequencies are in the interval $[0,\pi]$ (as you can see in the plot). The Matlab frequency variable $w$ is normalized as

$$w=\frac{2\pi f}{f_s}$$

where $f_s$ is the sampling frequency. So $w=\pi$ corresponds to $f=f_s/2$.


LoG filter creating additional maxima in scale space To create a scale space, I applied a Laplacian of Gaussian filter on the following image:

After the scale space was created, I plotted circles around local maxima in scale space. However, instead of $9$ circles, I got many more as shown below(image below not to scale).

The following video shows the LoG filter responses as scale is varied from $\sigma=1$ to $16$ in steps of $0.01$ -- http://youtu.be/m4LyXBI_6w0

It shows additional maxima being created at higher scales. But this should not be the case since LoG was especially chosen so as not to create additional maxima. What could be the explanation here ?

NOTE: Thresholding the magnitude of scale space responses eliminates the additional maxima as shown below. Could it be that the additional maxima being absent holds only up to a certain degree ?

Without exactly seeing your code it is hard to know for certain, but I suspect this is because you are not tracking your maxima well through scale space.

The position of the maxima will move in scale space. At the smallest scales you get the 4 small blobs around each circle. These are due to digitization and possibly noise/irregularities in the boundary edge.

These maxima will expand as the scale increases and merge together to form the "correct" blobs around the circles.

However, the other side of the peaks will also merge together resulting in the 12 medium size blobs and finally these will merge to form the four largest blobs.

Note at each stage here the number of maxima decreases 36->21 (9+12) -> 13 (9+4)

The problem is as you only look for local maxima in scale-space you lose  a lot of the relationship between the maxima.

One option may be to look for all local maxima at each scale and track these through scales. However, I'm not sure how straightforward this would be or how you determine which are the significant peaks using this approach.

A possibly better option may be similar to what you suggested by thresholding the responses. This article does basically the same thing as you. It applies a threshold on the mean blob volume over all scales. Which it claims it better than thresholding the peak response.


Detecting edges oriented at certain angle How would I detect an edge oriented at 22.5° using a spatial filter?

You can use steerable filters.  To detect an edge at a specific orientation you need to compute the derivative in the direction normal to the edge. To do that, you first compute the x-derivative and the y-derivative. Then the oriented derivative $I_{\theta} = I_x cos(\theta) + I_y sin(\theta)$

Doing this with a simple Sobel filter  generates the following images. The top left is a contour plot of the test image. The top right is the contour plot of the result of the steerable filter.  The bottom two show the contour plots of the original horizontal and vertical filters applied to the test image.

scilab code below.

Ix =  [1 2 1; 0 0 0; -1 -2 -1];//Sobel horizontal edge kernel;
Iy =  [1 0 -1; 2 0 -2; 1 0 -1];//Sobel vertical edge kernel;

theta = 22.5/180*%pi;

I_theta = Ix*cos(theta) + Iy*sin(theta);

test_image = zeros(100,100);
test_image(1:100,1:100) = (([1:100]-50.5)'*ones(1,100)).^2 + (ones(100,1)*([1:100]-50.5)).^2 < 100;

output_image = conv2(test_image, I_theta,"same");
xset("fpf"," ")

levels = [0:0.5:1];
axes = [0 100 0 100];

clf;
subplot(221);
contour(1:100,1:100,abs(test_image),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
title('Test image');

subplot(222);
contour(1:100,1:100,abs(output_image),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
xtitle('${\rm Steered\ in\ direction\ }\theta$');

subplot(223);
contour(1:100,1:100,abs(conv2(test_image,Ix,"same")),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
title('Horizontal filter')

subplot(224);
contour(1:100,1:100,abs(conv2(test_image,Iy,"same")),levels)
a = gca();
a.isoview = "on";
a.title.font_size = 5;
mtlb_axis(axes);
title('Vertical filter')


Extracting a DC value from a sum of sinusoids I guess it might be a basic questions but here it goes anyways.

How could I extract a DC value from a sum of sinusoids, i.e.:

$$v(t) =\widetilde{v_{dc}} + \sum_{n=1}^{\infty}\sin(\omega_nt+\phi_n)$$

I tried using low pass filters in Matlab or just mean of the signal but it's not quite what I'm looking for. I'd also like to have a mathematical way of expressing this.

You can go with linear filtering as others have suggested or do something nonlinear.

The sinusoids are symmetric around zero so if the exact dc value is subtracted prior to a hard clipper like Matlab’s sign() function, the, leaky integral of the hard clipped  time series would be zero, which means that it can act as the error of a feedback loop.

This sort of thing was done a lot by circuitry in older systems to compensate for the dc offsets that occurred in A/D converters that operated in a wide temperature range.

The integrator can be a linear leaky integrator, or a counter depending on how you implement it


Digital filtering based on frequency component amplitude... why not? I've got a very general question.
I've just implemented a filter working in frequency domain that simply cuts every frequency component whose modulus amplitude falls below a certain threshold value. As to say, "just keep the main frequencies in the signal". See example below:

Here is the effect in time domain on a noisy square wave:

It could seem to have poor performance in time domain, but this 'frequency picking' worked perfectly for my particular problem.
However, it seems to me that this is not so usual in DSP, and, since I'm not a DSP expert (I'm a physicist and not an electronic engineer), I wonder if this kind of approach (i.e. select the frequencies to keep based on their amplitude in the transformed domain and discard all the rest) is a known method in literature and, if it's so, what's its name. Or, on the contrary, I would be curious to know if there are some general reasons to avoid such a criterion.
I found no reference on the web...

Thanks!

Federico

I can tell you why I avoid such methods. When you pick some frequencies above threshold and remove others, you are using a kind of hard limit. In DSP this is same as Rectangular Windowing in frequency domain. It's effect in time domain will cause ringing. Because when you multiple in frequency domain with rectangular window, in time domain the signal gets convolved with a sinc function. There is a whole wikipedia page dedicated to this (Gibbs Phenomenon) - https://en.wikipedia.org/wiki/Gibbs_phenomenon

Again it depends on our requirements though. For example, if you see your filtered signal, the rise and fall time are not steep compared to original signal because you have eliminated higher frequencies which contribute to this. If my requirement was to apply a filter without compromising on fast rise time, I would not choose this method.


How to solve this equation (about complex envelope)? I am a beginner of digital signal processing. In my research project, I need to calculate this equation which is below.

$$ e(t) = \Re \Big\{z(t) \, e^{j \omega_0 t} \Big\} $$

I know the $e(t)$ and I want to know how to calculate $z(t)$ and $\omega_0$. And I also want to know what's the meaning of $\omega_0$, if it's the same as the center frequency of $e(t)$.

I am sorry about that this question is very basic. Thanks!

If $e(t)$ is band-limited and if the frequency $\omega_0$ is chosen such that the complex signal

$$y(t)=z(t)e^{j\omega_0 t}\tag{1}$$

has no negative frequencies, then $y(t)$ is called the analytic signal, and it can be determined from $e(t)$ using the Hilbert transform:

$$y(t)=e(t)+j\mathcal{H}\{e(t)\}\tag{2}$$

This analytic signal $y(t)$ is unique. However, the complex envelope $z(t)$ is not, unless we define the frequency $\omega_0$. Even though $\omega_0$ can be called center frequency, there is no unique definition of such a center frequency. That's why $\omega_0$ is often called reference frequency, and the complex envelope $z(t)$ depends on the choice of that reference frequency.

If $\omega_0$ is not given, you cannot uniquely solve for $z(t)$. It's up to you to define a convenient value for $\omega_0$, from which the complex envelope is easily obtained:

$$z(t)=y(t)e^{-j\omega_0 t}\tag{3}$$

Usually you would choose $\omega_0$ such that $z(t)$ is a low pass signal.

Note that since

$$e(t)=\Re\{z(t)e^{j\omega_0 t}\}=\frac12\left\{z(t)e^{j\omega_0 t}+z^*(t)e^{-j\omega_0 t}\right\}\tag{4}$$

(where $*$ denotes complex conjugation), the complex envelope can also be determined form $e(t)$ by down conversion and low pass filtering:

$$2e(t)e^{-j\omega_0 t}=z(t)+z^*(t)e^{-2j\omega_0 t}\tag{5}$$

Assuming that $z(t)$ is a low pass signal and that $\omega_0$ is chosen such that it is greater than the bandwidth of $z(t)$, low pass filtering $(5)$ will yield $z(t)$.


What is causal signal? I'm using Digital Signal Processing Principles, Algorithm and Applications 4th Edition written by Proakis and Manolakis. In Chapter 3, subtopic 3.3.2 it mentions the term "causal signal" which I can't find any definition. I know what a causal system means: a system that output depend entirely on past and present input, not future input but here it's causal signal, not system

So what does causal signal definition?

I tried to look at the Index and it says that the term "Causal Signal" is mentioned in page 85 but in fact, I find nothing in page 85 mentioning it

A system is causal if its output depends only on the current input and past inputs (and not on future inputs).  As a consequence, if the system is a linear, time-invariant (LTI) system (which input/output relationship can be completely characterized by an impulse response), that impulse response, $h(t)$ must be zero for all time $t \lt 0$.

Some people define a causal signal, $x(t)$, to be one that can be the impulse response of a causal system: it is zero for all time $t \lt 0$.


frequency characterritic of activ high pass filter How can i see  the characteristic(frequency response) of  an active high pass filter?I attached  the circuit and the value of resitence .Can i see it on the computer without an osciloscope ?How can i calculate it on the paper?
Thank you very much!

There are many online calculators with plot etc. options you could use. An example: http://sim.okawa-denshi.jp/en/Fkeisan.htm


Why is it called continuous-time frequency? I'm just wondering about the CTFT. My lecturer refers to capital Omega $\Omega$ in the following as being the continuous-time frequency:

Why is it called continuous-time frequency here but in the DTFT the small omega is called discrete-time frequency? CTFT and DTFT both give continuous signals out.

$\Omega$ is the usual angular frequency in radians per second, and is equal to $2 \pi f$. It is the way to measure frequency for continuous-time signals.

In discrete-time, frequency is measured in radians per sample, and is denoted as $\omega$. Here, a frequency component with $\omega = 2\pi$ is indistinguishable from $\omega = 0$.

$\Omega$ and $\omega$ is the notation in texts such as Oppenheim-Schafer. Other authors use $\theta$ for discrete-time, and the classic $\omega$ for continuous-time signals.


how to add two images in one, next to each other I have two splited images from diference matlab files. Now i need to merge them, near to each other. I tried supblot, imadd,cat methods, but nothing works. Please help me.

rezult after sublot experiment.

combinedImage = [leftImage, rightImage];


Inverse $\mathcal Z$-transform of system with an 8th order pole Can I find the inverse $\mathcal Z$-transform of this transfer function: $$H(z)=\frac{1}{1-\alpha z^{-8}}$$
in a way other than contour integration and finding the residues of the 8 poles?
If so, how?

Applying the geometric series directly

$$
H(z) = \frac{1}{1-\alpha z^{-8}} = \sum_{n=0}^{\infty}(\alpha z^{-8})^n = \sum_{n=0}^{\infty}\alpha^n z^{-8n}
$$
If $h_n = \mathcal Z^{-1}\{H(z)\}$, then
$$
h_n = \begin{cases}
a^{n/8} \; , \;  \text{n is multiple of 8} \\
0 \; , \; \text{otherwise}
\end{cases}
$$


What are the characteristics of a "good" smoothing convolution kernel? At work we were smoothing a signal by convolving with either

f1=[0.2000 0.2000 0.2000 0.2000 0.2000]

or

f2=[0.1111 0.2222 0.3333 0.2222 0.1111].

A colleague noticed that f1 smoothed "stronger" than f2 and suggested a new filter:

f3=[0.2727 0.1818 0.0909 0.1818 0.2727]

His reasoning was that this should smooth even "stronger" than f1, and therefore be better.

My answer to this was that

if he wanted "stronger" smoothing he should instead increase the
size of the filter
a good smoothing filter should always approach zero at the start and beginning

My motivation for 2. was Gaussian filters which are often used for noise reduction in image analysis.

Is 2. correct?

Also because of the periodic nature of DFT does not f3 have the same frequency response as

f4=[0.0909 0.1818 0.2727 0.2727 0.1818],

which in turn should have a very similar frequency response as f2?

Edit: I followed Jim Clay's suggestion and zero padded and the result shows that f3 has terrible stop band ripple while f2 is quite OK in  this respect:

Thanks in advance for any answers!

Yes, in general, your #2 is correct.  That being said, both of the filters stink (with your triangle filter being a little better).

No, f3 does not have the same frequency response as f4.  To get an idea of why that is so, you generally have to zero-pad the impulse response before DFT'ing it to get a reasonable idea of what its frequency response looks like.


the ROC of a Z-transform for shifted signal I have got two different answers for the ROC of the signal.
In that PIC, I have solved it in 2 methods, but I'm getting different answer.
Which one is correct?
Also please explain how to find the ROC of shifted signal?

The $\mathcal{Z}$-transform of your signal is

$$X(z)=1+\frac12 z^{-1}\tag{1}$$

(as you figured out by yourself).

The ROC is obviously the whole $z$-plane except for $z=0$, because there's a pole at $z=0$. Note that the ROC of the $\mathcal{Z}$-transform of any finite length sequence must be the whole $z$-plane, possibly except for $z=0$ (for sequences with non-zero values for $n>0$), and $z=\infty$ (for sequences with non-zero values for $n<0$). So by merely looking at the sequence you could determine its ROC.

The reason why your second method of determining the ROC doesn't work is because when you add the two $\mathcal{Z}$-transforms you get a pole/zero cancellation:

$$X(z)=\frac{1}{1-\frac12 z^{-1}}-\frac{\frac14 z^{-2}}{1-\frac12 z^{-1}}=\frac{(1-\frac12 z^{-1})(1+\frac12 z^{-1})}{1-\frac12 z^{-1}}=1+\frac12 z^{-1}$$

So the restricted ROCs of both functions caused by the pole don't show up in the final result due to this pole/zero cancellation.

Note that you could express the unit impulse by $\delta[n]=a^nu[n]-a^nu[n-1]$ and obtain arbitrary ROCs of the separate $\mathcal{Z}$-transforms of the two terms on the right-hand side, depending on the choice of $a$. However, the ROC of the $\mathcal{Z}$-transform of $\delta[n]$ (which is $1$) is obviously the whole $z$-plane, including $z=0$ and $z=\infty$.

In general, if you have a $\mathcal{Z}$-transform $X(z)$ and you shift the corresponding sequence to the right by $k$ samples, $X(z)$ gets multiplied by $z^{-k}$. So all poles and zeros of $X(z)$ remain unchanged, and a pole at $z=0$ and a zero at $z=\infty$ are added, both with multiplicity $k$. Of course they can cancel with already existing poles or zeros. Analogously, if you shift the sequence to the left by $k$ samples, $X(z)$ gets multiplied by $z^k$, which adds a $k$-fold pole at $z=\infty$ and a $k$-fold zero at $z=0$. As an example, take the sequence $a^nu[n]$ with $|a|<0$. Its $\mathcal{Z}$-transform is

$$X(z)=\frac{z}{z-a}$$

with a pole at $z=a$ and a zero at $z=0$. If you shift the corresponding sequence to the right by one sample, the resulting $\mathcal{Z}$-transform is $z^{-1}X(z)$. The original zero at $z=0$ gets cancelled by the new pole at $z=0$, and you're left with the original pole at $z=a$ and a new zero at $z=\infty$. If you shift the sequence to the left by one sample you get $zX(z)$, which gives a double zero at $z=0$, the original pole at $z=a$, and an additional pole at $z=\infty$.


Cyclic prefix vs. guard interval. Why use one over the other? I'm currently studying some OFDM techniques, and have found that in order to avoid Inter-Symbol Interference (ISI), there needs to be some space in between symbols. Lets assume this length to be $\mu$. This can either be done by adding a cyclic prefix to the beginning of the OFDM symbol, or by adding a simple guard interval (no transmission).

It seems to me like a cyclic prefix is typically used, but I'm not sure why this is used over a guard interval. using a guard interval would save power, and since the first $\mu$ points of the received signal will be ignored anyways, wouldn't it be better to save power? I'm not sure what the advantages of using a cyclic prefix are. Thanks for any insight!

A silent  guard interval would likely increase PAPR (peak to average power), which requires a sender capable of more dynamic range (without clipping or generating power supply noise, etc.)

A long enough cyclic prefix prevents multipath echoes from one symbol interfering with the portion of the next symbol after the prefix.  And a symbol interfering with itself from only its own multipath echo can be modeled as a flat fading channel, which may be easier to equalize or to error correct over a far wider range of multipath channel responses.


why type 2 FIR cannot be used to design high pass filter Kindly explain why type 2 FIR filter cannot be used to design high pass filter.

Type 2 filters are even and symmetric.  A signal with a frequency of $f_s/2$ is just a scaled version of [+1,-1,+1,-1,+1,...].

That means, at $f_s/2$, the response of a type 2 filter will always be zero as the coefficients of the same magnitude are an even number of samples apart so they cancel each other out.

A filter with a magnitude frequency response of zero at $f_s/2$ makes a very poor high pass filter.


Is possible to get overshoot on bessel filter? I tried to use bessel filter from here,
but in the case of a large dataset, I have overshoot.  Since I know a bessel filter shouldn't have overshoot, is this code wrong or am I doing something else wrong?

// Sample rate : 1e9
// Cut off Frequency : 50e6

int s=50000;
double[] sample = new double[s];

for(int i=s/2; i<s; i++) {
  sample[i] = 1.0;
}

for (int i = 0; i<sample.length; i++) {
    xv[0] = xv[1];
    xv[1] = xv[2];
    xv[2] = xv[3];
    xv[3] = sample[i] / GAIN;
    yv[0] = yv[1];
    yv[1] = yv[2];
    yv[2] = yv[3];
    yv[3] = (xv[0] + xv[3]) + 3 * (xv[1] + xv[2])
          + (  0.3402988787 * yv[0])
          + ( -1.4152725242 * yv[1])
          + (  2.0223552191 * yv[2]);
    filteredValue[i] = yv[3];
}

All coefficients have been calculated based on source on given website, it might been messy because its used for chebyshev and butterworth filters as well.

It's not exactly clear what you mean by overshoot. Let's define it as "when filtering a square wave, the max of of the the output should not exceed the steady state of max of the square wave".

In this case your statement "as i know bessel filter shouldnt have overshoot" is wrong. A bessel filter will indeed overshoot. A few first order filters will not overshoot (butterworth, chebychev, ...) but any filter of order > 1 will overshoot.


How do we distinguish between Image Recognition,Detection and Segmentation? I am reading image processing book.while reading i am getting terms like image detection,image recognition and image segmentation. it seems that all above terms are one and the same.i am totally confused amongst them.so can anybody help me to understand or distinguish between the above?

It can be a little confusing at times, and the terms are not completely independent.

Detection:

In detection, you are simply detection the presence of something. For example, you might design an algorithm to detect beach-balls from pictures of a beach. You would feed an image into your algorithm, and it will spit out an answer 'yes', if there was a beach ball, and 'no' if there was no beach ball. This is very simply what it means to 'detect' something.

Segmentation:

Segmentation is a method - a technique - an algorithm. So for example, if someone says "How did you design your beach-ball detection algorithm?", you might say "I used a segmentation algorithm". Or you might say, "I used a template-matching algorithm". Or you might say "I used a scale and shift-invariant matching score algorithm", etc. Those are all methods of doing something, and segmentation is just a method.

In particular, segmentation simply means a thresholding of your image into various areas. In a simple case, say you want to separate background from foreground. Then, you segment the image into two regions: A region that corresponds to the background, and a region that corresponds to the foreground.

Image Recognition:

Lastly, there is image recognition, which is simply, a very broad goal.  For example, you might say, "I achieved image recognition for beach scenes, by detecting beach balls using the method of segmentation, detecting water using the method of k-means, and detecting bikinis using the method of template matching.

Of course, those are just examples, but broadly speaking, this is how you should parse the different words and terminologies you come across in this field. Hope that helped!


Why does a pole close to the unit circle result in an enhanced Q-factor? I have the following pole zero plot (in Z-domain) of a digital filter, where its poles lie far away from the unit  circle, with its corresponding frequency response:

On the other hand I have this one:

Could someone explain me why exactly the poles give such a change in Q-factor?
I expected that the closer a pole to the unit circle the less high the Q-factor would be.

EDIT:
a pole closer to the unit circle will on one hand mean an "stronger" presence of a certain frequency it won't be decaying exponentially. But on the other hand a pole on the circle would mean a marginally stable system. In other words this means having a frequence which is perfectly present results in a +/- unstable system. Isn't that contradictory?

Thank you

Suppose some $H(q)=\frac{A(q)}{B(q)}$ where q is some complex variable and $A,B$ are functions of $q$.

Whether in the s or z planes, to evaluate the magnitude of $H(q)$ at some $q$, you evalaute and sum all distances from $q$ to the locations of the zeroes (i.e. the magnitude of $A$) and similarly for $B$ and create the fraction above for that specific $q$.

When does $H(q)$ becomes big? When either its nominator is big...OR its denominator is very small. In fact, if the denominator was ever to become exactly zero, then $H(q)$ would "shoot to infinity". It would acquire a very large value.

When does $H(q)$ becomes small? When either its nominator is small...OR its denominator is very big. The smallest value $H(q)$ can acquire is zero.

  ...why exactly the poles give such a change in Q-factor?

Given the location of the poles in your diagrams, it is useful to see what happens as your $q$ (or a particular $z$, in the case of the z-plane), approaches the bearing of the pole. In other words, as the $z$, changes frequency and walks around the z-plane, it approaches the location of the pole. EXACTLY AT the bearing of the pole, the distance to the zero is maximum. BUT! the distance to the pole is at its minimum. Therefore, the denominator becomes small and the response "shoots to infinity".

When the pole is closer to the center of the circle, its distance to the particular $z$ value is large. The closer the pole moves to the edge of the circle the smaller its distance will be with a particular $z$ and the higher the output of $H$ will be.

To explain the "peakiness" (or "...Q-factor...") you have to look at the rate by which the distance to the pole is diminishing as $z$ approaches the bearing of the pole.

  ...a pole closer to the unit circle will on one hand mean an "stronger" presence of a certain frequency it won't be decaying exponentially. But on the other hand a pole on the circle would mean a marginally stable system. In other words this means having a frequence which is perfectly present results in a +/- unstable system. Isn't that contradictory?

When examining the response of a system via the z-transform, the "presence of a frequency" means that the system resonates at that frequency. At extreme conditions, the system is so marginally stable that if you happen to excite that frequency even by a tiny little amount, it will go into oscillations that might take a very long amount of time to die out. In other words, the system will "ring".

A prime example of this are 1-pole resonators. These are bandpass filters with very high Q-factor. If you pass audio through such a resonator, you will hear the filter ringing as it effectively resonates at its own frequency.

  ...How should I modify my pole zero plot if I'd like to render undershoot in stead of overshoot like here?

If you are asking how to generate a notch filter, instead of a resonator (i.e. a filter that "cuts" one specific frequency instead of boosting it), then you can simply subtract the output of the 1-pole resonator from the rest of the signal OR, you could try to create a "well". A region of the z-plane where the conditions are such, for your $H$, that it dips to zero at a specific bearing to a particular $z$.

Have a ponder about how that looks like and if you cannot "see it", have a look at this link.

Hope this helps.


What's wrong with this Average Magnitude Difference algorithm? I implemented this Average Magnitude Difference algorithm in Javascript

var abs = Math.abs
var floor = Math.floor
var POSITIVE_INFINITY = Number.POSITIVE_INFINITY
function autocorrelation (data) {
  var length = data.length
  var halfLength = floor(length * 0.5)
  var summation = null
  var smallest = POSITIVE_INFINITY
  var smallestOffset = null
  for (var o = 1; o < halfLength; o++) {
    summation = 0
    for (var i = 0; i < halfLength; i++) {
      summation += abs(data[i] - data[i + o])
    }
    average = summation / halfLength
    if (average < smallest) {
      smallest = average
      smallestOffset = o
    }
  }
  return smallestOffset
}

and I use it this way

analyser.getFloatTimeDomainData(data)
var rms = rootMeanSquare(data)
var offset, frequency
if (rms > RMS_MIN) {
  offset = autocorrelation(data)
  frequency = SAMPLE_RATE / offset
}

to find the frequency of the signal input to the analyser.

Signal input is the output from an oscillator. The oscillator's frequency varies from C0 to B8 and I can't always detect the correct frequency. Can you tell me why, please?

Please, find below two images showing oscillator's frequencies and related detected frequencies

okay, i took a look at your code and at your data.

first of all, the algorithm you are using is not autocorrelation but is the Average Magnitude Difference Function (AMDF).  it will give results similar to the autocorrelation but upside-down.

you appear to be implementing the AMDF correctly but you are suffering from the so-called "octave error" problem.  you see, a note that is, say, A-440, has a fundamental frequency of 440 Hz.  that means the waveform repeats every $\tfrac{1}{440}$-th second.  but that same waveform also repeats every $\tfrac{1}{220}$-th second and every $\tfrac{1}{110}$-th second or every $\tfrac{1}{55}$-th second.  so, mathematically, it can be construed to be a 220 Hz or 110 Hz or 55 Hz (or even a $\tfrac{3}{440}$-th second period or $\tfrac{440}{3}$ Hz waveform).

so you need to figure out how to turn your AMDF upside down and offset it to get peaks instead of valleys.  for the math relating AMDF to ASDF and that to autocorrelation, take a look at this and particularly this.  you have multiple pitch candidates and you have to elect the correct candidate.  don't do it like the Americans did last November.


Issue with AMR file structure For developing an application I have to identify a region of sound. But I’m not sure if I have my basic knowledge right. I’m using AMR-NB sound clips to identify a sound region. According to my understanding AMR files have 160 samples for each 20ms (milliseconds) making it 800 sample for 1 minute. I am basing my complete dependence on this concept. Please let me know whether this is correct if not please correct me. Thanks in advance.

Actually, if it's got 160 samples every 20ms, you get a sampling rate of 8000 (you're missing a zero), which is the standard for over-the-phone voice communication. from a Google search, I couldn't find any AMR-NB specification that uses a different sampling rate.

Also, I found this line here pretty useful:

  For AMR-NB, only an “8000” sample rate and one of the following
  bitrates are valid:
  4.75k, 5.15k, 5.9k, 6.7k, 7.4k, 7.95k, 10.2k or 12.2k

Wikipedia also mentions that AMR-NB uses a 8000 Hz sampling rate, you can be pretty sure that this is the case all across the board.


How to mimic/copy/fake someone's voice? Is there any existing application to sample someone's voice and use it to modulate any other voice or synthesize a text to resemble the original one?

For example, this Text-to-Speech Demo by AT&T lets you choose a voice and a language from presets that I guess are based on some human voice that have been sampled.

How do you call this process? Is it voice modulation? Voice synthesis?

A first note: Most modern text-to-speech systems, like the one from AT&T you have linked to, use concatenative speech synthesis. This technique uses a large database of recordings of one person's voice uttering a long collection of sentences - selected so that the largest number of phoneme combinations are present. Synthesizing a sentence can be done just by stringing together segments from this corpus - the challenging bit is making the stringing together seamless and expressive.

There are two big hurdles if you want to use this technique to make president Obama say embarrassing words:

You need to have access to a large collection of sentences of the target voice, preferably recorded with uniform recording conditions and good quality. AT&T has a budget to record dozens of hours of the same speaker in the same studio, but if you want to fake someone's voice from just 5 mins of recording it will be difficult.
There is a considerable amount of manual alignment and preprocessing before the raw material recorded is in the right "format" to be exploited by a concatenative speech synthesis system.

Your intuition that this is a possible solution is valid - provided you have the budget to tackle these two problems.

Fortunately, there are other techniques which can work with less supervision and less data. The field of speech synthesis interested in "faking" or "mimicking" one voice from a recording is known as voice conversion. You have a recording A1 of target speaker A saying sentence 1, and a recording B2 of source speaker B saying sentence 2, you aim at producing a recording A2 of speaker A saying sentence 2, possibly with access to a recording B1 of speaker B reproducing with his/her voice the same utterance as the target speaker.

The outline of a voice conversion system is the following:

Audio features are extracted from recording A1, and they are clustered into acoustic classes. At this stage, it is a bit like having bags will all "a" of speaker A, all "o" of speaker A, etc. Note that this is a much simpler and rough operation than true speech recognition - we are not interested in recognizing correctly formed words - and we don't even know which bag contains "o" and which bag contains "a" - we just know that we have multiple instances of the same sound in each bag.
The same process is applied on B2.
The acoustic classes from A1 and B2 are aligned. To continue with the bags analogy, this is equivalent to pairing the bags from step 1 and 2, so that all the sounds we have in this bag from speaker A should correspond to the sounds we have in that bag from speaker B. This matching is much easier to do if B1 is used at step 2.
A mapping function is estimated for each pair of bags. Since we know that this bag contains sounds from speaker A, and that bag the same sounds but said by speaker B - we can find an operation (for example a matrix multiplication on feature vectors) that make them correspond. In other words, we now know how to make speaker 2's "o" sound like speaker 1's "o".
At this stage we have all cards in hand to perform the voice conversion. From each slice of the recording of B2, we use the result of step 2. to figure out which acoustic class it corresponds to. We then use the mapping function estimated at step 4 to transform the slice.

I insist on the fact that this operates at a much lower level than performing speech recognition on B2, and then doing TTS using A1's voice as a corpus.

Various statistical techniques are used for steps 1 and 2 - GMM or VQ being the most common ones. Various alignment algorithms are used for part 2 - this is the trickiest part, and it is obviously easier to align A1 vs B1, than A1 vs B2. In the simpler case, methods like Dynamic Time Warping can be used to make the alignment. As for step 4, the most common transform are linear transforms (matrix multiplication) on feature vectors. More complex transforms make for more realistic imitations but the regression problem to find the optimal mapping is more complex to solve. Finally, as for step 5, the quality of resynthesis is limited by the features used. LPC are generally easier to deal with a simple transformation method (take signal frame -> estimate residual and LPC spectrum -> if necessary pitch-shift residual -> apply modified LPC spectrum to modified residual). Using a representation of speech that can be inverted back to the time domain, and which provide good separation between prosody and phonemes is the key here! Finally, provided you have access to aligned recordings of speaker A and B saying the same sentence, there are statistical models which simultaneously tackle steps 1, 2, 3 and 4 in one single model estimation procedure.

I might come back with a bibliography later, but a very good place to start to get a feel for the problem and the overall framework used to solve it is Stylianou, Moulines and Cappé's "A system for voice conversion based on probabilistic classification and a harmonic plus noise model".

There is to my knowledge no widely piece of software performing voice conversion - only software modifying properties of the source voice - like pitch and vocal tract length parameters (For example IRCAM TRAX transformer) - with which you have to mess in the hope of making your recording sound closer to the target voice.


FastICA: what happens when we have more source than channel? I am reading the original article proposing FastICA and I have a couple of question on point not covered in the article. Both answering me or providing a source will be very appreciated.

I want to use it in blind audio source separation context, so I am going to use that language.

Let's suppose that I have more source than channels.

$$x_1 = a_{11}s_1+a_{12}s_2 + a_{13}s_3+a_{14}s_4$$
$$x_1 = a_{21}s_1+a_{22}s_2 + a_{23}s_3+a_{24}s_4$$

Where $x_1$, and $x_2$ are our observed channels, $s_i$ the unknown sources and $A = (a_{ij})$ the unknown mixing matrix.

FastICA suppose that we have as many sources as channels, but does not check it in any way. If I am going to run the algorithm on my data, it would extract two sources and a 2x2 mixing matrix.

What should I expect from those sources? Are they going to be weighted mixing of my true sources? If yes: mixed how?

Thanks in advance, and sorry if I did some language mistake.

PS: I don't have a strong formation in Signal Processing, so if you can provide some source or small explanation when using specific language it would be appreciated.

In a way, you have answered this already with:

  FastICA suppose that we have as many sources as channels, but does not check it in any way. If I am going to run the algorithm on my data, it would extract two sources and a 2x2 mixing matrix.

ICA will indeed provide a separation, along the lines of section 3 ("What is independence?") which would also require some elements of section 4.2("Measures of non-gaussianity"), 4.3("Minimisation of mutual information") and 4.4("Maximum likelihood estimation")

ICA is trying to express some $x$ (in your notation), as a mixture ($A$) of a collection of waveforms ($s$).

A better way of looking at this is that $x$ is the combined result of a number of processes that move from state to state with some probabiliy. By expressing $x$ in this way and at the same time requiring that you want to find independent components, you arrive at the concept of statistical independence.

You want the waveforms in $s$ to be as different to each other as possible but at the same time, when you sum them up (with some mixing weights determined by $A$) you want their sum to amount to $x$.

This is what ICA will do and it will do it whether you feed it audio samples or some random signals that are completely out of context.

  What should I expect from those sources? Are they going to be weighted mixing of my true sources? If yes: mixed how?

This is a very interesting and at the same time difficult question to answer accurately.

You should expect that these "sources" will be "independent" to each other but not that they correspond to meaningful waveforms in an audio context.

A trivial example is this: I record 2 trumpets playing exactly the same score. Should I expect that they will be "separated" by ICA?

No.

ICA will try to find a minimal set of independent "components" as if these were the output of a number of independent processes. One possible decomposition might have one component sounding very closely to the score and the other being a totally random set of values. These two would be independent to each other and when mixed (at the right levels) would result to the recorded waveforms.

Don't expect necessarily "meaning" out of these independent components.

A better "result" can be obtained by attempting to discriminate between percussive sounds and the rest of the sounds in a recording. Again, this does not mean that it will give you the isolated drums as if you had a contact mic on the drum set and recorded only those sounds. This is because, percussive sounds have long(er) "silences" and their distributions tend to be non-gaussian (a key requirement for the independent components (section 4)).

Hope this helps.