Implementing inverse FFT using forward FFT I am trying to implement an inverse FFT using the forward FFT. For clarity:

Let S[t] be a signal in time, and S[w] the transformed signal. As per this site, it seems one can reverse S[w], use the forward FFT routine, then reverse the resulting signal again and this should give S[t]. I won't go into why should it work but it's all in the link provided.

I attempted to try this method, and it seems to recreate the signal when I use a high sample rate. However, it seems that even if I do not go over the niquist frequency, I encounter some weird effect. Here is my code in numpy:

# Generate a signal of s cosine with 200 [Hz]
f = 200
Fs = 10
t = np.linspace(-10,10,Fs*20)
s = cos(2*np.pi*f*t)

# Take fft
u = np.fft.fft(s)

# Reverse in time
u = u[::-1]

# Transform again
u_t = np.fft.fft(u)

# Reverse and normalize

s_new = np.divide(u_t[::-1],s.shape[-1])

# Finally slice for easier viewing
plt.plot(t[1:100],s_new[1:100])
plt.plot(t[1:100],s[1:100])

This code yields the following graph:

I am a little confused, the theory looked sound and I can't think what I did wrong here.

Note: I tried calling np.fft.fft and then np.fft.ifft and the reconstruction goes as planned. Therefore I believe this problem is not due to aliasing.

Edit: I made way by simply taking the complex conjugation instead of reversing (using np.conj() where I reverse). This solves the problem, but I still do not understand why reversing in time does not accomplish a conjugation, so I'd love someone to explain it to me. I leave the modified code in case it helps someone else:

# Generate a signal
f = 200
Fs = 10
t = np.linspace(-10,10,Fs*20)
s = cos(2*np.pi*f*t)

# Take fft
u = np.fft.fft(s)

# Conjugate
u = np.conj(u)

# Transform again
u_t = np.fft.fft(u)

# conjugate and normalize

s_new = np.divide(np.conj(u_t),s.shape[-1])

# Finally slice for easier viewing
plt.plot(t[1:100],s_new[1:100])
plt.plot(t[1:100],s[1:100])

In your link I cannot find a reference to the time-reversal. However, here are some additions:

The time reversal $y(t)$ of a continuous signal $x(t)$ is given by $y(t)=x(-t)$. Similarly, for a discrete signal we have $y[n]=x[-n]$. However, here the delicate issue occurs that $n=0...N-1$, where $N$ is the signal length. Now, noting that in the finite discrete domain every signal is considered to be periodic with $N$, you get
$$y[n] = x[-n] = x[(-n)_N]$$

where $(x)_N$ is the modulo operation of $x$ by $N$. Hence, time-reversal in the discrete domain is not just reading your vector from back to front. Instead, the first element remains the same ($n=0$) and the others are read from back to front.

Now coming to the time-reversal property wtih the Fourier Transform:

$FFFF=I$ where $F$ is the (unitary) Fourier transform operator and $I$ is the identity. I.e. applying 4 times the (forward) Fourier transform yields the original signal.
Hence, logically $FFF=F^H$, where $F^H$ is the inverse Fourier Transform (yes, it is also the Hermitian of the forward transform). However, calculating inverse FFT by 3 times forward FFT is not very efficient...
Further, we know $FF=T$ where $T$ is the inflection operator $(Tx)[n]=x[-n]$, again read with the modulo interpretation.

Hence, we find $F^H=TF$, i.e. in order to caculate inverse FFT, you can do forward FFT and time-reverse the result:

f = 0.1
Fs = 10
t = np.linspace(-10,10,Fs*20)
s = np.cos(2*np.pi*f*t*(t+1))

# Take fft
U = np.fft.fft(s)

# Transform again
u_t = np.fft.fft(U)

# Reverse and normalize

s_new = np.divide(u_t[-np.arange(u_t.shape[0])],s.shape[-1])

# Finally slice for easier viewing
plt.figure(figsize=(12,4))
plt.plot(t,s_new)
plt.plot(t,s)


Are the components of a Bandpass signal namely "inphase" and "quadrature" real functions in time domain? A band pass signal representation goes by the generalization as
$X(t)=XX(t)*e^{j \cdot2\pi \cdot ft}$
where $f$ be the carrier freq. and $XX(t)$ be the complex envelope.
On further decomposition it boils down to:
$X(t)=XXi(t)\cos(2\pi ft)-XXq(t)\sin(2\pi ft)$

Here is it implicit that $XXi(t)$ and $XXq(t)$ are real functions and their spectra be even symmetric?

$X(t) = Re\{XX(t)e^{j2\pi ft}\}$. So $XX(t)$ can be any generic complex baseband signal with real I and Q components - $XX(t) = XX_i(t)+j XX_q(t)$ and $e^{j2\pi ft} = \cos(2\pi ft) + j\sin(2\pi ft)$

After the complex multiplication, and taking real part, you get $X(t) = XX_i(t)\cos(2\pi ft) - XX_q(t)\sin(2\pi ft)$.

So $XX_i(t)$ and $XX_q(t)$ are indeed real functions and their spectra is even-symmetric.

If $\hat{X}(t)$ is the Hilbert transform of $X(t)$, then $X^+(t) = X(t)+j\hat{X}(t)$. $XX_i(t)=X(t)\cos(2\pi ft)+\hat{X}(t)\sin(2\pi ft); XX_q(t)=\hat{X}(t)\cos(2\pi ft)-X(t)\sin(2\pi ft)$


Question about MRI signal construction As I understand it, in MRI, a line in k-space is usually acquired while the readout gradient is applied. So during the application of this gradient the protons spin at different frequencies along the gradient.

Now, I sample the MRI signal during the application of this gradient. My confusion is that I can sample the signal at arbitrary resolution and the inverse FFT will produce the image with as many pixels as the sampling resolution.

To make the question a bit concrete. let's say I have simulated 128 spins at a distance of 1 mm apart. This will give me an object with the FOV of 0.128 meters.

Now say during the readout I sample 10000 equidistant points. So, effectively I have collected 10000 complex k-space values. An inverse FFT would also give me an image with 10000 pixels. Now, how do I resolve this 10000 pixel image into a 128 pixel image where each pixel corresponds to the corresponding spin.

I would expect that you simply cut the outer parts of the image away. Certainly, there is always a hassle to figure out the proper pixel in the center, but that is a detail of the FT-algorithm employed.

Remeber: Sampling denser than you need increases the FOV. Hence, the image you simulate should simply be larger with the original image (that you get by simulating 128 sampling points in the setting you described above) sitting at its center. Cutting off the outer void parts is therefore the technique you need, and not an interpolation technique.


In Matlab: how do I quickly create smoothing matrices? Say I have 5 measurements of a signal as a column vector: x.
I can smooth the signal by multiplying it with a smoothing matrix:

$\ \mathbf{xs=S*x} $

For instance I might use [1 2 1] as a smoothing kernel in which case I get:

S =
     2     1     0     0     0
     1     2     1     0     0
     0     1     2     1     0
     0     0     1     2     1
     0     0     0     1     2

So far I have been using the function sparse() to build the matrix diagonal by diagonal,
but this quickly becomes tedious.

How can I quickly create a mxm smoothing matrix from a 1xn convolution kernel in Matlab?

You could do the following-

kernel = [1 2 1];
s = conv2(eye(numMeasurements), kernel, 'same')


Fun problem: Dividing a dataset into parts depending on mean value Ok, so I have some recorded data over several hours that looks like this (values on y, time on x).
See this image http://imgur.com/fdakUlk
I need a way to use matlab to divide the data into sections depending on their mean value. The sections needs to be atleas a couple of hours. Just by looking at the data they can easily be spotted, see this image: http://imgur.com/svrP79h
Is there a good method of doing this? What do I need to research?

Thanks

Try total variation denoising. Here I found a piecewise constant toolbox: http://www.maxlittle.net/software/pwctools.zip


Asymptotic bode plots I am really struggling to see how the lecturer took this transfer function and produced the bode plot that I have in my notebook. It is not the plotting so much that is confusing, I just don't understand how he manipulated the transfer function into this form that he did and why it needed to be in that form in the first place when it looked solvable in its original form.

  $\textrm{(b)}\quad T(s) = 180\displaystyle\frac{s(1+0.01s)}{(1+0.05s)(1+0.001s)}$

Looking at the transfer function shown I would have thought that all one would simply need to do is evaluate the gain at $20\log(180)$ and then plot the bode plot using the two zeros and two poles.

Supposedly to plot the asymptotic bode plot for this function correctly one needs to get it into the following form:

$$3600\times\frac{\left(\frac s{20}\right)\left(1+\frac s{100}\right)}{\left(1+\frac s{20}\right)\left(1+\frac s{1000}\right)}$$

Sorry for the bad format but basically the answer provided has it rising to its pole at $20$ on the $71\textrm{ dB}$ line and going flat until it hits $100 \textrm{ rad/s}$ and then rising up again at $20\textrm{ dB}$ per decade until it meets its pole at $1000$. I don't understand what is going on and everything I have written hitherto is just my observation of how the lecturer demonstrated.

If anyone with any knowledge on the mechanics of plotting bode plots by hand could shed some light on this for me I'd greatly appreciate it.
Simon.

Quick Answer cause I'm at work and don't have time to do the plotting but basically, in equation b you already have the expression i what i like to call "The Bode Form", to only a gain and terms of the form $(1+s/A)$ with $A$ some value. I'm sure the professor explained this, that when you take logarithm you add and subtract all this different terms because the log turns multiplication into addition and division into subtraction, etc. This individual terms, When $s < A$ are approximately equal to 1 (a constant line at 0db Gain in the plot), and when$ s > A$ are approximately equal to $s/A$ ( a 20db/decade drop in the gain).

So the plot would be an increasing line of 20db/decade from  up to your first pole at 20 rad/s, where the gain would be $20\log_{10}(3600)$, and then it would remain at that until 100 rad/s where the next pole appears which again makes the gain drop at 20db/decade, and then at the next pole at 1000 rad/s it drops at 40db/decade and that's it, that's your plot.

Could you maybe put up the image you have in your notebook so everyone could see what the professor did.

Attempt at illustrating below. This video may be instructive.


Upsampling PCM audio: from 6 kHz to 8kHz What would a straight-forward way to convert an audio payload (PCM 16bit) from 6 kHz to 8kHz? I understand that this is an interpolation problem.

I fear that a linear interpolation would introduce too much noise. If so, what would be a better interpolation? I was looking at Lanczos resampling but I've seen it used in image scaling only and not sure it's worth it for audio.

I'm using Java.  I've found the JSSRC library butI would have prefer something more lightweight. Beside, I'm not sure this is really necessary for low-frequency audio (which happens to be decoded from ADPCM).

BTW, I'm not interested in creating a WAV having a 6 kHz frequency header, I really want to resample.

If processing cost is not a problem, you may use the following algorithm for sample rate converting your data from 6 khz to 8 khz.

1- expand your signal x[n] by 8, producing xe[n] (insert 7 zeros after every sample of x[n])
2- design a low pass filter hlpf[n] with a gain of 8 and cutoff frequency pi/8
3- process the expanded signal with this lowpass filter = xef[n]
4- compress output of the filter by 6 via selecting every 6th sample: y[n]=xef[6n]
(NOTE: You can also do this with 4:3 upsample and down sample, with a new filter,  which would require less processing cost)

Lowpass filter must be sharp enough to ensure proper anti-aliazing. If this is too much overburden for  simple rate conversion then try high order polynomial resampling. That could also be costly however.


Generating digitized white noise: uniform vs normal sampling Consider the following two ways of generating noise in the time domain for audio applications:

Generate samples from a uniform distribution [-amplitude, +amplitude], where amplitude is in the clipping range.
Generate samples from a Gaussian/normal distribution with μ = 0 and some σ > 0, followed by clipping the samples to the clipping range.

I always thought that the important property of white noise is that it is an uncorrelated process, and expected (1) to produce white noise in the sense of having a flat frequency spectrum. I was never sure why some implementations prefer to use a Gaussian distribution and just came across this comment:

White noise is a continuous process from any uncorrelated random process, like uniform or normal. However, if you digitize it, you must apply a bandpass filter at the Nyquist frequency, otherwise your approximation of the continuous process contains aliasing. It turns out that bandpassing white noise results in a discrete random process where each sample is picked from a Gaussian/normal distribution.

I do not fully understand that comment, and it triggers a bunch of questions:

Does it mean that generating noise via (1) suffers from aliasing, and thus, its frequency spectrum isn't flat, i.e., it should not be considered white noise?
If so, what is the frequency spectrum of (1), i.e., how does aliasing manifests itself here?
If proper white noise generation requires (2), what about the clipping issue of sampling from a Gaussian? Isn't it an issue that for large σ the tail of the distribution gets chopped off?

I'd appreciate any hint that clarifies above It turns out that... statement.

Your question is an interesting project for you to research it on your own. Well, maybe with a little help from your friend and the SE community. And, as your question goes, start with generating samples from (1) a uniform distribution, and (2) a Gaussian distribution. I will accompany you in your first steps.
First, I generate waveforms of 1,048,576 samples for each type of noise. Here in the pictures the ranges of 512 samples long are shown.

Uniform noise

Gaussian noise
The pictures look quite similar.
Then I compute FFTs of waveforms for each of type of noise and plot the distributions of values of the noise waveform (blue graph) and of the corresponding FFT (orange graph).

Uniform noise. Distribution of sample values in blue, distribution of FFT values in orange

Gaussian noise. Distribution of sample values in blue, distribution of FFT values in orange
I do not quite understand what you mean talking about "a flat frequency spectrum". Both signals are realizations of stochastic processes. But here they are, and the trends in the graphs of sample value distributions clearly betray their origin. At the same time, distributions of FFT values (real parts of FFT values) for both signals have a distinct look of Gaussian distribution. To prepare yourself to follow the explanation of this behavior, I recommend you to study one of the proofs of the central limit theorem of probability theory, the proof, which uses Fourier transform. With this knowledge you may even find the explanation on your own. If not, still it is a useful thing to know when following the explanations in the literature.
Notice also, in the graphs of the distribution of FFT values, the central peaks, which are two times higher then the interpolated distribution at zero.
UPDATE
While modeling white noise in computations, under certain conditions you are free to select any method that will generate random values for samples in your simulation. The generated values can follow uniform, Gaussian, any other distributions with a zero mean value. Even a coin toss random generator, randomly producing +ampl/-ampl with a 50/50 chance, will fit. Of these "certain conditions", one is just mentioned -- a zero mean value. Another requirement is a high sampling frequency, the higher the better. Within any sample interval between adjacent signal samples you must generate a great many noise samples, densely filling the time interval. This, admittedly nonobvious, procedure is justified by the essence of "white noise signal": it is a stochastic process, not a classical function. In an ordinary language, it means that we cannot know a noise value at any given moment, but noise values averaged over an arbitrarily small time intervals can be "measured" in experiment, or "simulated" in computations following a known formula, this "known formula" following from the central limit theorem of probability theory.
With this approach, the signal sampling (and, consequently, Nyquist) frequency is much lower than the noise sampling frequency. To avoid aliases in the noise frequency domain, you must bandlimit the noise. In this process, the noise is averaged over the time interval between the adjacent signal samples in a manner dependent of the filter used. By force of the central limit theorem, the averaged noise values tend to approach a Gaussian distribution -- the more samples used in averaging process, the closer the distribution to Gaussian. But this fact suggest the idea to directly generate samples from a Gaussian distribution and avoid oversampling at frequency much higher than the signal sampling frequency.
The method generating a random variable with a Gaussian distribution is readily available in matlab, numpy and scipy (I believe). In the general-purpose computer languages you can easily implement the Marsaglia polar method or the Box–Muller transform.
Summing up, the "It turns out that..." statement from you question is the corollary of the central limit theorem. You can see it either following the central theorem proof or even better generating a number of i.i.d. random variables and computing the distribution of their sum, repeating the process with an increased number of random variables and examining the results.

The generation of noise via (1) never suffers from aliasing in principle, but the computed spectra can be aliased if an inadequate Nyquist frequency is used for computation. With the Nyquist frequency fixed by the signal sampling frequency, you must adjust the noise sampling frequency.
Aliasing manifests itself via a phenomenon called Spectral Leakage.
The standard deviations resulting in the signal values exceeding the maximum values of the variable types used in the program do chop the tail of the distribution. Using float's rather than UINT16's is a safe option to fight this chopping.

NUMERICAL EXPERIMENT: NOISE GENERATION
Start with a silent (zero-valued) signal sampled into a data array of 65536 samples. First, add a uniform-distributed noise, one noise sample per one signal sample. The first 256 samples of realization (pure noise) are shown here:

Compare the distribution of sample values against a Gaussian distribution of equal power:

The FFT of the realization:

The PSD of the realization:

The zoomed in PSD of the realization 256 samples long:

Back to the original zero-valued pure signal, and add a uniform-distributed noise again, but this time insert two noise samples per signal sample. To produce a realization of the noisy signal sampled at the original pure signal Nyquist frequency, we sum the adjacent pairs of samples. To keep the noise power constant, we scale down the noise amplitude by the square root of oversampling (=2). The first 256 samples of the new realization are shown here:

Compare the distribution of sample values against a Gaussian distribution of equal power:

With just two-fold oversampling, the distribution is already close to the Gaussian distribution, as it is seen in the difference plot of the two distributions:

The zoomed in PSD of the realization 256 samples long:

Not much different for the plots of purely uniform noise.
Finally, calculate with the oversampling of 4096 noise samples per signal sample:
The realization after averaging (first 256 samples):

The difference of the two distributions, uniform-folded and Gaussian:

The PSD of the realization:

UPDATE on AWGN:
The phrase "but the computed spectra can be aliased if an inadequate Nyquist frequency is used for computation" may lead to confusion. I'll try to explain it without explicitly referring aliasing (the aliasing thing still implicitly leaks through the references, but these are presumably more trustworthy than my explanations).
See it the other way round: when, in your attempt to generate noise, you generate one random uniform-distributed variable per sound sample (the sampling frequency for generated noise is 44.1KHz), the PSD is also random and very uneven. To approximate the white noise's PSD behavior, you can resort to N-oversampling, increasing the sampling and Nyquist frequency to N·44.1KHz: you spread each value of sound sampled at 44.1KHz among adjacent N samples of a new (N·sample_count)-sized vector (array) of noised samples and add uniform noise.
std::random_device rd;
std::mt19937 gen(rd());
std::uniform_real_distribution<double> unirand(-1.0, 1.0);
for(int isample=0; isample < vecSig.size()/*44.1-samples_count*/; ++isample)
  for(int ix=0; ix < N; ++ix)
    vecOversampled[isample*N + ix] = vecSig.at(isample) +  unirand(gen);

The PSD of this signal averaged over N adjacent frequency bins is a way more even. Seeing that the frequencies higher than 44.1KHz are sort of computation crutches, we filter out these frequencies from the signal and arrive at the equivalent of the Gaussian-distributed random variables added to the original 44.1KHz samples, because the central limit theorem holds for the averaged sample values. Alas, the PSD of the filtered signal is ragged anew, for the PSD computation is a non-linear operation, and the PSD plot can be smoothed only with enlarged frequency bin sizes. But this time we know that this raggedness is due to windowing function and we have a leaking in process rather than leaking out.
This PSD raggedness is the result of sampling of the white noise "signal".
Summing up: to spare memory and CPU time, you can do without oversampling, generate a Gaussian-distributed noise from the beginning
std::random_device rd;
std::mt19937 gen(rd());
std::normal_distribution<double> normalrand(-1.0/std::sqrt(M_PI), 1.0/std::sqrt(M_PI));
for(int isample=0; isample < vecSig.size()/*44.1-samples_count*/; ++isample)
    vecSigPlusNoise[isample] = vecSig.at(isample) +  normalrand(gen);

and arrive at the same results as with the oversampled uniform-distributed noise. You cannot do without oversampling when modeling the noise with uniform-distributed random variables -- without delving into the depths of probability theory, I refer you to the Additive White Noise definition (https://en.wikipedia.org/wiki/White_noise)

In discrete time, white noise is a discrete signal whose samples are regarded as a sequence of serially uncorrelated random variables with zero mean and finite variance; a single realization of white noise is a random shock. Depending on the context, one may also require that the samples be independent and have identical probability distribution (in other words independent and identically distributed random variables are the simplest representation of white noise). In particular, if each sample has a normal distribution with zero mean, the signal is said to be additive white Gaussian noise. [my emphasis]

whereas each sample in a non-oversampled uniform-distributed noise model has a uniform distribution.
I presume that at this stage you need an AWGN simulation for your project.
It may be just an intuition in the requirement that the "ideal" noise model must exhibit similar patterns at increasingly small scales, a so-called self-similarity feature. You see that if you start with a uniform-distributed or any other non-Gaussian-distributed noise at some fine frequency resolution and "zoom out" this noise to a larger, coarser scale, you arrive at the noise tending to approach a Gaussian-distributed noise -- whereas the Gaussian-distributed noise reproduces the Gaussian distribution of sample values at any zoomed resolution (with a properly adjusted standard deviation parameter). This self-similarity feature is often used for the analysis and computation in noise-related development projects. As soon the white noise model is applicable, you can safely assume it is AWGN in discrete time and use the self-similarity feature for analysis.
In addition to the self-similarity of noise value distributions among zoomed frequency resolutions, another argument for Gaussian distribution versus any other distribution can be derived from a process of direct construction of Fourier spectrum to satisfy a "flat PSD" requirement.
The straightforward function of frequency which gives the most flat, indeed constant, PSD, is a constant, i.e., independence of frequency. Unfortunately, the signal, of which the Fourier spectrum is independent of frequency, is a delta function, and the delta function is everything but the noise.
To endow the constant Fourier spectrum with noise-like traits, one can follow your recipe for producing the noise in time domain: each sample of frequency spectrum is made a random uniform-distributed variable. But, as we have seen, the Fourier spectrum constructed in this way gives the time-domain realization with Gaussian-distributed sample values!
We have discussed two possibilities to simulate white noise. In one approach, we start with the uniform-distributed samples in a time domain and have the Gaussian-distributed samples in a frequency domain. In another approach, we start with the uniform-distributed samples in a frequency domain and have the Gaussian-distributed samples in a time domain. Interestingly, this means that there exist realizations with Gaussian-distributed samples in the frequency domain for which the procedure of inverse Fourier transform gives uniform-distributed samples in the time domain. But the share of these realization in the signal space is negligible (a set of measure zero), and for an overwhelming majority of signal realizations, starting with the Gaussian-distributed sample values in the time (frequency) domain, one obtains the Gaussian distributed values in the frequency (time) domain.
The Gaussian distribution appears to be the most natural distribution for generating white noise realizations. In the signal space, the subspace of non-Gaussian-distributed white noise realizations has measure zero.


Proper rounding when converting between digitized values and floating point I have some very large files that are unsigned 8 bit (0 to 255), real valued.  I need to do some simple tuning and decimation.  In the past, when I've dealt with low bit rate samples, I've run into issues with artifacts.  I want to make sure I get the rounding correct.

The data will be processed in Matlab using all native double precision.  The first question is should I subtract 127.5 or 128 when converting, or does it depend?

To keep things simple, I'm not doing anything fancy with the filter gain or multiple passes to make sure I don't clip.  Instead I'll notify the user if there was a problem.  The output also needs to be the same format, so I can't use any additional effective bits gained by decimation.  The next thing I'm unsure about is do I round 255.3 down to 255 or do I consider this an overflow condition?

EDIT

I thought I could just process the data with the offset in there doing something like this at the output.

y = y * scale + (1-scale)/2*255

The trouble is the tuning operation may filter out the offset before I correct for it.  So instead, I decided the simplest thing to do is subtract 127.5 and add it back in at the end.  Then round as I described above.

I'm not sure how the different methods of adjusting for the offset and scaling effect the resultant artifacts, but as arnfinn said there is probably nothing I can do about in any case (unless I'm willing to use dither).

EDIT 2

This is actually more along the lines of what I was thinking would be an issue.  In the distant past I have experienced this exact problem with requantization.  Certain combinations of quantization parameters can really mess up quality (worse artifacts).  My situation is a little different because I'm processing the data before requantization.  But I believe it can still be an issue.

A requantization algorithm for the transcoding of JPEG images

Signal Processing: Image Communication.
Volume 21, Issue 1, January 2006, Pages 13–21.
Jae Won Moona, Jong Seok Leeb, Nam Ik Cho.

I think you don't have to worry too much about so-called 'transparency' (ability to represent accurately the integer value as a floating point value) unless you do normalization. See the libsndfile FAQ.

As quantization is a lossy operation (many-to-one, or surjective, mapping) you will always have artifacts when re-quantizing from floating point to integer values. The standard method for de-correlating the quantization error is to use dither. For small word lengths, such as 8 bit, the classic solution of adding noise with a triangular distribution might cause the dither noise to be unaccptably high, so you might have to look around for some semi-optimal (in terms of quantization error moments) dithering method that suits your application better. Or just re-quantize and accept the artifacts, they might be small, depending on your data and what you are planning to do with it.


Determine poles and zeros of a specific filter design I have a design question and need to determine poles-zeros of a IIR bandpass filter.
Requirements are:

1-order should be 3

2-causal

3-real valued impulse response

4-stable

5-cutoff freq: pi/3-pi/2

I have tried butterworth design on matlab to determine poles and zeros but I do want to learn how are they placed analytically or by inspection. Also a problem of matlab is it designed a filter with order=4. Is it possible to solve this with n=3.

Thanks

There's a good reason why band pass filters usually have an even filter order. As you may know, real-valued filters have either real-valued poles and zeros, or their poles and zeros occur in complex-conjugate pairs. If the order is odd, there must always be at least one real-valued pole and zero. In the case of a band pass filter, a real-valued pole is not very useful, because poles create peaks in the frequency response, but a real-valued pole can only create a peak at DC (frequency zero), or at Nyquist (the maximum frequency given the sampling rate). However, a band pass must attenuate the signal at DC and at Nyquist, so peaks at those frequencies are undesirable.

So a reasonable choice would be to place a complex conjugate pole pair at (or close to) the center frequency of the band pass filter, and leave the remaining real-valued pole at $z=0$ (where it doesn't affect the magnitude response). The three zeros should all be placed on the unit circle to create frequency response zeros. A reasonable choice would be to have one zero at DC and a complex pair of zeros close to Nyquist.

The figure below shows the magnitude response and the pole zero diagram of a possible choice of poles and zeros for a third order band pass filter with the given pass band.


Root raised cosine pulse shaping filter What are disadvantages of root raised cosine pulse shaping filter in digital communications and why does it need to be improved?

Links:

Square Root Raised Cosine Fractionally Delaying Nyquist Filter
Square Root Raised Cosine Filter
Generalized raised-cosine filters

The primary reason that root raised cosine filters are used is because they can be expressed in closed form and they are simple to understand. There are far better root Nyquist filters out there.

Usually they are optimized through least squares or remez exchange to improve stopband attenuation or some other channel specific criterion.

The ISI-free response (obtained through the cascade of a root Nyquist, the channel, and another root Nyquist) is only accurate in a flat channel (although it is also a good approximation in an almost flat channel). The shortcomings of root Nyquist filters have largely been addressed by DMT/OFDM, SC-FDM, and other proposed (and still experimental) technologies. These technologies result in ISI-free responses even for channels that the root Nyquist filters do not.

However, all of this is very dependent on channel models. If you really want to study the advantages and disadvantages then you need to provide a channel model, because the limitations and drawbacks of any method will be dependent on the channel model you are considering.


Output of an FIR filter has a peak at the beginning of the signal I have a signal  $x[n]$ and I filter the signal using an FIR bandpass filter to receive $y[n]$. (implementation in MATLAB using built-in function firpm).
The filtered signal $y[n]$ is compensated for its delay of half the filter order.

The problem is that the filtered signal $y[n]$ has a major peak at the beginning that does not exist in the original signal.

What is the cause of that peak?
Is it possible to attenuate the peak?

The peak (or response) is shown in the second plot of the filtered signal between 0-0.5 seconds (Note that I only show the first few seconds of a several minute signal)

The signal is sampled at 250 samples/sec. The desired filter passes frequencies between 10-20 Hz. The code for the implementation of the filter:

fs=250; % sampling rate;
n= 120; % filter order;
freqs = [0 5 10 20 25 fs/2]/(fs/2); % normalized frequencies in Hz;
amps = [0 0 1 1 0 0 ]; % amplitudes;
b = firpm (n,freqs,amps); % define filter coefficients;
filtered_signal = filter (b,1,original_signal); % apply FIR filter;

It's the initial transient at the filter output due to the implied zeros of the input signal for $n$ less than $0$.

Let me simulate the effect here, since you do not post the actual data I will model it as close as necessary, with the following Matlab/Octave code:

Fs = 250;                    % Sampling frequency in Hz.
Ts=1/Fs;                     % Sampling period in seconds.
t = [0:Ts:4];                % Simulation interval for the signals.
L = length(t);               % Simulation signal lengths in samples.

x = 8E4 + 1E4*exp(-t)  + 1E3*filter(fir1(43,0.03),1,randn(1,L));
% Note: x[n] is roughly modeling your input signal sampled at 250 Hz.
% I have not included those spikes, as they're irrelevant for this analysis

b = fir1(65,[10/125 20/125],'bandpass');   % a BPF with passband [10,20] Hz.
y = filter(b,1,x);           % Obtain the filtered output signal.

% Now, observe the outputs...
figure,subplot(2,1,1);
plot(t,x);title('Input signal x[n] of duration 4 seconds')
subplot(2,1,2)
plot(t,y);title('Output signal y[n] obtained from BPFing of the input');

As it can be seen below:

So the peak at the beginnig of the output is clearly observed.

It's intuitive to look at the bandpass filter's impulse response :

Looking at the BPF impulse response (especially the second subplot) reveals the fact that the transient is inline with the impulse response; their durations are the same and their shapes are very similar.

In fact, looking at the transient's duration (which is about 0.25 seconds) it matches with the number of samples the impulse response has:
$ N \times T_s = 0.25 s \rightarrow N = 0.25 \times 250 = 62.5 \text{ samples}$ This coincides with the BPF impulse response which was designed as a $65^{th}$ order filter having $66$ coeffcients.

So we have the clue that the peaking at the begining is indeed the impulse response (a transformed version of it) at the output.

The reason of this can be understood by looking at the convolution sum:
$$y[n] = \sum_{k=0}^{L-1} x[k]h[n-k] $$ which shall be evaluated for $n \geq 0$

now in the sum it can be seen that signal values for $k < 0$ are implied to be zero. Lets show this effect on the flip and drag graphical evaluation below:

Now first obervse that in its initial part, the signal $x[n]$ has an almost constant value. Then it can also be seen that until $n=65$ the overlap between the signal and impulse response is partial and therefore there won't be enough cancellation from the positive and negative samples of the impulse response to produce (apprx) zero output. But when $n \geq 65$ there is full overlapp between the impulse response and the signal (which is almost constant) therefore in the summation the positive and negative samples of the impulse response will cancel each other yielding a small, (zero-dc), output which represents the expected output after $n > 65$.


What exactly is a complex envelope? I have seen this be mentioned a couple of times in some books I read, so I want to make sure. Is the complex envelope simply the summation of the real and quadrature components of a signal, whereby the absolute value is the (real) envelope? I have read this wiki page but I am not sure I completely understand. Is the complex envelope simply the real and imaginary parts of a passband signal combined? Thank you.

Depends on context.

For strictly real-valued signals it sometimes refers to creating an analytic signal (see: http://en.m.wikipedia.org/wiki/Analytic_signal ) from a signal with no imaginary component, and taking the magnitude (or complex with both magnitude and phase) envelope of that hypothetical now complex creation.


How i can a 1 Hz resolution fft, from a limited data (between 0 to 100 msec)? How i can have a 1 Hz resolution fft, from a limited data (between 0 to 100 msec)?
sampling frequency is about 20 kHz or step time is 0.05 msec.
i try it by using zero padding but the resolution not satisfy me.

Your DFT's resolution is always $\frac{f_\text{sample}}{L_\text{FFT}}$, by definition of what the DFT is.

So, for any given frequency resolution, you'd need to transform the inverse of that frequency worth in samples. So, for 1Hz resolution, you need 1/(1Hz) = 1s worth of samples; you're giving it less than 100ms, i.e. 1/10 of that – you can zero-pad, but for an FFT-based frequency estimate, that few samples simply won't be enough to give you a low-variance estimate.

So, garbage in, garbage out.


Calculating phase of DFT sinewave? I have been attempting to make a basic, slow, DFT in Matlab and have noticed peculiar behavior that I don't understand.

I have been trying to plot the phase of a 100Hz sinewave captured at 250kHz sample rate with a record length of 10k points. I find that the real part of the DFT is negative at 100 Hz, which I think makes sense due to a sinewave having negative correlation with a cosine wave. However, the imaginary part at 100Hz is also negative, which again makes sense because a sinewave is shifted -90 degrees from a cosine wave.

When I take the inverse tangent of the real/imaginary, obviously this is a positive number. Resulting in a +90 phase shift instead of a -90. Am I misunderstanding something, or do I just have incorrect results?

edit: I've realized the issue is that due to spectral leakage and the fact that this operation there are non-zero magnitudes that need to be filtered out before calculating the phase.

Is there a good way to choose a threshold for when to zero out DFT results? It seems a hard value to choose as if I throw in an arbitrary waveform then I have no information beforehand to choose a proper threshold value. My initial idea is to calculate the bin size based on the time length of the input, and then use a lookup table to change the threshold to be 10%, 1%, 0.1%, etc, of the maximum magnitude in the DFT.

4 full cycles of a sine wave (starting at a phase of 0) is an odd function in aperture, so the real part of DFT result bin 4 should be zero (or nearly so with numerical noise).

atan2(-1,0) should result in -pi/2


Spatial receive diversity: a different view It is well known that spatial diversity can combat multipath small scale Rayleigh  fading in a SIMO system by offering the receiver multiple options of independent receptions using more than one antenna ( this is called selection combining, which is one of several spatial diversity schemes). The basic idea is that the probability of several antennas all simultaneously seeing destructive fading is much lower than the same probability in case of only one antenna
Since by adding more and more antennas we get closer and closer to an AWGN channel, this means that the channel becomes more and more independent from multipath fading, and this is also called channel hardening (channel converges more and more to a deterministic value)
The question is the following: Is there a mathematical derivation that shows how can a random Rayleigh PDF (probability distribution function)  in case of one antenna converge more and more to a deterministic value in case of a large number of receiver antennas

multiple options

So, we're talking about Selection Combining (SC); there's other combining methods, too!

independent receptions

Good that you're stating that, this will be handy in a minute.

Is there a mathematical derivation that tackles this subject from this point of view?

The very definition of a CDF tackles that, but only if the channel realizations have a distribution that have a finite magnitude maximum, which is usually not the case.
Let $|h_i|$ be the magnitude of the $i$th of the $N$ individual channels. Then, the cumulative density function of the "best of all $N$ channels" is simply defined as $\eqref{def}$
\begin{align}
F_{\text{SC},N}(x) &= P\left(\max_{i=1,\ldots,N}|h_i| < x\right) \label{def}\tag{1}\\
&=P\left(|h_1|<x \,\wedge\, |h_2|<x\,\wedge\,\ldots\,\wedge\,|h_N|<x\right)\tag2{}\\
&=\prod_{i=1}^N P\left(|h_i|<x\right)\tag3&\text{due to independence}\\
&=\prod_{i=1}^N F_{|h_i|}\left(x\right)\tag3 & \text{definition of "CDF"}\\
&=\prod_{i=1}^N F_{|h_1|}(x)\tag4 &\text{all identically distr.}\\
&= \left(F_{|h_1|}(x)\right)^N\tag5\label{ord}
\end{align}
I'm almost certain you've seen $\eqref{ord}$ already – you're looking for the concept of diversity order.
Now, let's see about your infinite number of antennas. Two ways to say the same thing:

You know that this doesn't make sense. If I'm sampling any continuous distribution, such as an absolute normal distribution, an infinite amount of times, I will surely hit any non-zero-mass interval. So, even from basic definition of probability, the answer is already clear: you can assume that you get an arbitrarily large channel coefficient! Congratulations!
$\lim_{N\to\infty} k^N = 0\quad\forall |k|<1$, no matter what $k$ actually is. But if $k=F(x)$ is a cumulative probability function (CDF), it's always $<1$ (unless you hit the maximum value), by definition of CDF. In other words: you surely never get anything but 0; that means, the probability of getting a value out of your diversity combining that is smaller than any value that's not the maximum of your possible channel realizations is 0, which means you always get the maximum possible channel.

This makes no sense physically, because it especially means that with finite transmit energy, one of your infinitely many receive antennas gets infinitely much receive energy, if you assume any fat-tailed distribution of channel coefficient magnitudes (e.g. normal). That's (one of the many reasons) why your considerations don't matter for communication engineering.


Benefit to know Fourier series for image processing? I know there's a benefit of knowing the Fourier Transform for image processing, but is there a benefit to know Fourier series, or could you just skip them? Would you recommend skipping Fourier series or would that make it too abstract? What's your opinion on this?

In image processing, you are always dealing with discrete data, so the transform isn't really a Fourier transform, it is a discrete Fourier transform (DFT with a sum rather than an integral).  The DFT can be understood as providing the coefficients of a Fourier series from which a periodic version of you're original image can be reconstructed.  The "checker board" image is always implied in the inverse of the DFT which would construct an infinite checker board of squares made up of your original image.

Normally the area outside the boundaries of you're original image is just ignored (treated as undefined), but if for some reason that area must be considered when reconstructing from a DFT (you transformed a sub-section of an image for example) that area will be filled with a periodic replication of of the transformed portion of you're image.

I would say it is worth the little bit of time it takes to understand Fourier series before moving on to Fourier transforms so you understand the underlying assumptions and mechanics of the transform.


Shifting and phase difference Assume we have a series v[k] of length M=16. It's shifted to the right in a circular manner by k0=2 samples.

Which phase difference between the DFT spectrum of the original signal and that of the shifted version is to be expected?

The phase shift will be a factor of $e^{j2\pi k_0 m/M}$, where $m$ is the DFT bin number


Smallest FFT buffer size given zero-padding I am currently working on a project that involves processing an audio signal by dividing it into chunks of size B. The buffer is then zero padded to have a size of 44100 points, after which FFT is applied so that the peak frequency of interest is detected. The reason why I'm zero-padding is to have a 1-Hz FFT resolution.

I would appreciate your input on two questions, please:

Given this process, what's the smallest FFT buffer B (before zero-padding) that I could have?
Is there a relationship between the size of the buffer B (before zero-padding) and the accuracy of detecting the peak?

Zero padding does not increase frequency resolution, but interpolates the spectrum you have by adding more samples in between as an interpolation. The frequency resolution is well known to be set by the size of your buffer before zero padding with the relationship:

$$f_{res} = 1/T$$

Where T is the length of your buffer in time (notice that this relationship is independent of the sampling rate!). However this classic relationship is the width of each bin by their 3.0 dB bandwidths or half power point (viewed as a filter-- it is helpful and accurate to view the FFT as a bank of filters), while to completely resolve two frequencies of equal magnitude, due to their coherent addition in the FFT, the frequencies need to be separated by greater than the half magnitude width of each bin (or the 6.0 dB bandwidth).  The 6.0 dB bandwidth for a rectangular window (no window) is 1.21, so to resolve two tones of equal magnitude and any arbitrary phase the length N must be greater than:

Minimum number of samples, two tones same magnitude (no Window):

$$N > 1.21 f_s/f_{\Delta}$$

Where $N$ is the number of samples, $f_s$ is the sampling rate and $f{\Delta}$ is the frequency separation both in Hz, using the relationship:

$$T = \frac{N}{f_s}$$

Note that using any window will decrease this resolution further (the best resolution is obtained with the rectangular window, which is no window) at the needed benefit of dynamic range. Thus choosing a window is always a trade between frequency resolution and dynamic range, and therefore the window is chosen based on the anticipated maximum power difference between closest frequencies and number of samples needed.  For example, the 6.0 dB BW of the Hamming window is 1.81 which modifies the above formula accordingly:

Minimum number of samples, two tones same magnitude (Hamming Window):

$$N > 1.81 f_s/f_{\Delta}$$

Please refer to this paper by fred harris On the Use of Windows for Harmonic Analysis that provides detailed considerations for window selection and resolution bandwidth for resolving tones, particularly the section G on p 178 "Minimum Resolution Bandwidth" with reference to his figure copied below, along with Table 1 on page 176 that lists detailed metrics for all common windows.

In the case of multiple tones at different power levels and in the presence of background noise, it is helpful to look at the complete magnitude response of the kernel (frequency response of the "filter" for each bin). Note again, since the tones in different bins add coherently they will sum in magnitude based on their relative phases (and subtract!), while noise components will sum in power. This is shown below by comparing the normalized magnitude plots of the rectangular and Hamming kernel, with a horizontal line showing the -6 dB half magnitude point to relate to the relationships derived above. The width of the main lobe of the rectangular window is 1.21 and for the Hamming window it is 1.81, and with reference to the figure above we see how this is the minimum separation for two tones of equal magnitude, being the boundary condition where a single maxima can exist.  But we also see the details related to dynamic range considerations.

Specifically consider two tones separated by 20 dB, the rectangular window will require greater than 5.36 bins ($f_{\Delta} > 5.36 f_s/N$) while the Hamming window will be able to resolve with separations greater than 3.32 bins in this condition:

Minimum number of samples, two tones 20 dB difference (no Window):

$$N > 5.36 f_s/f_{\Delta}$$

Minimum number of samples, two tones 20 dB difference (Hamming Window):

$$N > 3.32 f_s/f_{\Delta}$$

The paper referenced above details this further with consideration to multiple tones and in the presence of noise, where due to the effect of the side-lobe roll-off the window chosen has a significant impact.

For consideration of noise, margin is added to your maximum separation in power between tones; theoretically for AWGN you can statistically always encounter a condition where the noise added to one tone in magnitude compared to the noise subtracted from the next is beyond your designed separation, therefore it becomes a probability of false detection problem, where you design a margin based on the probability of not separating two tones in your designed resolution. I suspect this would quickly be insignificant if you are working with high SNR's in a bin (and remember the FFT has a processing gain of 20LogN decreasing the noise in one bin compared to your overall noise). This is best demonstrated with an example: Consider the case of adjacent tones of equal magnitude with a noise limit we bound to 1.5 (which could also then reduce the adjacent tone by 0.5), this would result in a worst case separation in the tones of 20Log(1.5/5)= 9.5 dB.  If we therefore added 9.5 dB to our total maximum designed separation in magnitude between adjacent tones (and then use the process we used above to compute number of FFT samples needed) we will ensure that we will always be able to detect the adjacent frequencies for all cases that the noise itself does not cause a 3x magnitude separation between bins. (Placeholder here to show awesome graph of SNR vs false alarm rate in this case: This would be an interesting joint probability problem to solve at some point if I or someone else has the time and energy to work through it assuming not readily available somewhere).

So to answer your question, your smallest buffer size is set by the actual frequency resolution you need, with consideration to total number of tones in close proximity and dynamic range needed (max distance in magnitude between tones). Zero padding is done if you want to interpolate more samples (without adding any more resolution; it is just a smoothing), which you may do for convenience (to get the FFT out to the closest power of 2 for example, or to fill in more samples on a plot.

In the case of detecting the frequency of a single tone, the zero padded FFT will more readily point out the precise frequency location (but this is not different from what you can do more efficiently through interpolation using the information derived from the closest bins if you did not zero pad), but if you had multiple frequencies together within the span of the frequency resolution defined above, you would not be able to resolve them (they would appear as one frequency).

As @OlliNiemitalo has correctly pointed out that zero padding WILL allow you to visibly observe the obtained frequency resolution described above (which is what I have done to plot the responses in between bins in the plot above).

In this post is an example plot showing the difference between zero padding to increase the number of samples versus adding more actual samples:

What happens when N increases in N-point DFT


IIR filter with circular buffer - Arduino I am trying to implement a butterworth lowpass filter with a circular buffer. The code below does not filter as intended and I am not really able to understand what is wrong. Any suggestion?

#define PERIOD_MICROSECS 1000

static uint32_t lastRead = 0;

int analog_pin = 0;
int32_t analog_input0 = 0;
double analog_input0_lp_filtered = 0;

// The coefficients are calculated offline. I set the sample frequency to 1000 Hz, the cutoff to 50 Hz.
// The filter has order = 4.
double c[] = {0.000416599204407, 0.001666396817626, 0.002499595226440, 0.001666396817626, 0.000416599204407};
double d[] = {1.000000000000, -3.180638548875, 3.861194348994, -2.112155355111, 0.438265142262};

int m = 10;
double x[10];
double y[10];

int n = 0;

double lp_butterworth(double input_sample, double *x, double *y, int m, int n)
{
    x[n] = input_sample;
    y[n] = d[0] * x[n] + d[1] * x[(n-1+m)%m] + d[2] * x[(n-2+m)%m] + d[3] * x[(n-3+m)%m] + d[4] * x[(n-4+m)%m]
                     - c[1] * y[(n-1+m)%m] - c[2] * y[(n-2+m)%m] - c[3] * y[(n-3+m)%m] - c[4] * y[(n-4+m)%m];

    return y[n];
}

// ************************* //

void setup() {

  Serial.begin(115200);
}

void loop() {

  if (micros() - lastRead >= PERIOD_MICROSECS) {
        lastRead += PERIOD_MICROSECS;

        analog_input0 = analogRead(analog_pin);
        analog_input0_lp_filtered = lp_butterworth(analog_input0, x, y, m, n);
        n = (n + 1) % m;

        //Check the original and filtered signals with the serial plotter
        Serial.print(analog_input0);
        Serial.print(" ");
        Serial.println(analog_input0_lp_filtered);
  }
}

Below you can see the result of the code in the serial plotter using as analog input an accelerometer at steady state. You can notice a huge offset and not really a filtering effect.

First of all :

have you initialized all the previous data to 0 ? The x and y buffer should be initialized to 0. Otherwise, you could get some unexplained transients at the beginning. Maybe this is the offset that you see at the beginning. Maybe the compiler will set the initial data in the data at 0, but you should not rely on the compiler, it should be explicitely reset to 0.

Secondly : I think you have mixed up the c and d coefficients. I think the d coefficients should be used with the output (y) and the c coefficients for the input (x). I say that because d[0] = 1 which makes me think that this the a0 coefficient.  Usually the a0 coefficient is 1 in an IIR. It can be different than 1, but usually we try to normalize the coefficients by setting a0 to 1. https://en.wikipedia.org/wiki/Infinite_impulse_response

Third : Have you tried the coefficients in Octave/Matlab/Python ? Use freqz or Bode in Octave/Matlab. Do you get the expected frequency response?

Finally : This is not related to the matter at hand, but I notice that you use the modulo function. The modulo function is slow unless you're using powers of 2. You should consider having a circular buffer with a size of 4 instead of 10. You don't need more since you only go 4 samples in the past.


Sampling Frequency in a OFDM Technique I am a new learner of OFDM. One source I refer explained OFDM in this way:

  Consider 4 symbols $X_0$, $X_1$, $X_2$, $X_3$ $\in \mathbb C$. They
  are modulated using Multi-Carrier Modulation. The time domain equation
  is given as

  $$c(t)=\sum_{n=0}^{3}X_ne^{j2\pi f_nt}$$

for the duration $0\leq t
> \leq 4T$, and $f_n=\dfrac{n}{4T}$.

Now, he samples the signal at a sampling rate, $F_s=\dfrac{1}{T}$ so that

$$c(mT)=\sum_{n=0}^{3}X_ne^{j2\pi \dfrac{n}{4T}mT}=\sum_{n=0}^{3}X_ne^{j \dfrac{2\pi nm}{4}}$$

This means that the sampled version of $c(t)$ is the IDFT of the symbols $X_i$, $i=0,1,2,3$.

What I am thinking is, since the spectrum of $c(t)$ is the shifted version of sinc pulses, considering only the main lobe of each sinc pulse, the spectrum is band limited to $\dfrac{-1}{4T}\leq f \leq \dfrac{4}{4T}$. So the sampling frequency should be twice the maximum which is $2*\dfrac{4}{4T}=\dfrac{2}{T}$. This is not taking me to the IDFT equation. Can you please correct where I am doing wrong?

What you're doing wrong is considering only the main lobe. Don't do that!

The sinc as spectral shape is elegant, because it has its zeros in all the other carriers' exact frequencies. So don't truncate; also, truncating a sinc after the main lobe ignores a lot of its energy.


Online filter for EEG My goal is to build a real-time filter for EEG data. The data are being recorded in 512Hz and I wish to filter it (in a real-time, from 0.5Hz to 30Hz in the most optimal way).

Unfortunately I'm not expert in DSP stuff, so I'm not sure which type of filter I should use, an IIR or FIR? What is the recommended filter type butter, cheb1, or cheb2?

All I understand is there are some properties to set when designing a filter. Among of them are the transition band, the order of the filter, a group delay, Gibbs-effect, linear phase, etc.

Currently, I'm really puzzled about these lists of properties. I also believe that some of the properties mentioned above are dependent on each other. I couldn't find any recipe that is "easy to use " on the web of how to build such a filter and what considerations should be made when building such filters.

I'll try to guide your research a bit rather than giving you a direct answer. One-click recipes aren't too good for obtaining optimal results, and this requires some careful insight. After all, you are dealing with biological signals.

To develop a real-time filter you need to process each sample/window before the next one comes up. You need to take into account the CPU time available for your calculations (depending on the platform) and the number of signals you must filter at once (with the subsequent impact in used memory, of course). After all, a digital filter is no more than a bunch of sums and products.

Now you should study the kind of signal you need to filter. You seem to be handling frequencies from delta to beta waves. You need an insight in the characteristics of each signal so that you keep them after the filtering (e.g. if studying brain activity during sleep, you should make sure that you don't miss out k-complexes and spindles during Phase II). This way you can look into some values: the maximum band pass ripple you are allowed, the width of the transition band, etc.

Having this information, after some research (e.g. Digital filter design by Burrus & Parks; Discrete-time signal processing by Oppenheim & Schafer) you can choose the kind of filter you need, and go into further details about implementation.

Fortunately, with a simple google query you can find out lots of things. There are Matlab toolboxes available to aid your filter design:

FDATool, Matlab filter design toolbox.
you can also take a look on ERPLab, which seems a good starting point.

And lots of free papers: Google query for "EEG filtering".

The higher the order of the filter, the more accurate it is, but the more operations per second it needs. So the main constraint here is that you need to keep it real-time for a bunch of signals.


Perfect reconstruction ( Single Channel ) I was trying to solve the following problem .

  Given an input x[n], consider upsampling by 2, followed
  by interpolation with a filter having z-transform H(z) for magnification of the signal.
  Then, to recover the original signal size, apply filtering by a decimation filter G(z) followed
  by downsampling by 2, in order to obtain a reconstruction ˆx[n].

using the results for upsampling as $Y(z) = X(z^2)$

and downsampling as $Y(z) = X(z^{1/2})+X(-z^{1/2})$

I obtained

$\hat{X}(z) = 1/2\bigg(G(z^{1/2})H(z^{1/2})X(z) + G(-z^{1/2})H(-z^{1/2})X(-z)\bigg)$

however, the solution is listed as ,

$\hat{X}(z) = 1/2X(z)\bigg(G(z^{1/2})H(z^{1/2}) + G(-z^{1/2})H(-z^{1/2})\bigg)$

Where have i gone wrong ?

I solved in the following process ..

$A(z) = X(z^2)$

$B(z) = H(z)A(z)$

$C(z) = G(z)B(z))$

$\hat{X}(z) = 1/2[C(z^{1/2})+C(-z^{1/2})]$

and substituted them in a chain.

HINT: Note that if

$$Y(z)=X(z^2)$$

then obviously

$$Y(z^{1/2})=X(z)$$

but also

$$Y(-z^{1/2})=X((-z^{1/2})^2)=X(z)$$

from which the desired result follows in a straightforward way.


What frequencies are present in the Fourier transform of the Dirac impulse? When I do the Fourier transform of the Dirac impulse I get a pure sinusoid (or complex exponential, however you wanna call it) but I read in several places that all frequencies are present in the dirac impulse and all of them with the same amplitude. How is this possible? Am I wrong when I perform the transform?

A Dirac impulse $x(t)=\delta(t-d)$ has the continuous-time Fourier transform $X(\Omega)$ of
$$\mathcal{F}\{\delta(t-d) \} = 1 e^{-j\Omega d} $$

whose magnitude is $$|X(\Omega)| = 1 ~~~, \text{ for all } \Omega $$ and a phase of $$<X(\Omega) = -\Omega \cdot d $$

Note that it's incomplete to think of the real or imaginary parts of the Fourier transform alone. Rather the magnitude and phase point of views are more reflective of the nature of the result.

So in this case the magnitude is $1$ and hence it's said to contain all frequencies of magnitude $1$.
Note that these are differential amplitude components of continuum frequency range as opposed to a finite amplitude of discrete set of frequency components, aka line components.


Why does shortening window widen its Fourier Transform? On page 76 of the book Discrete Time Speech Signal Processing by Thomas F.Quatieri, the author states "shortening the window widens its Fourier Transform".

Can anyone explain why?

As per my knowledge, you can use the context of Duality property of Fourier to understand this. So if you reduce the time window size then the spectrum (Fourier transform)should widen(Duality). Hope this helps!


Aliasing Issues in Sub Band Coding The general structure of Sub-Band Coding (SBC) is shown below:

The diagram shows that after each of the analysis filters, decimation is performed by a factor M which corresponds to the number of channels. Would this not result in aliasing when this decimation by N is performed for the higher frequency sub-bands and why is DDC not required in each path to translate each sub-band to baseband before performing the decimation and thus avoiding aliasing?

I think you're just misjudging what the [CHANNEL] block does: this is just saying "all the bits that fall out of the different quantizers/encoders get transmitted".
It's especially not an adder! Most books actually show this just as a "frame packer" or "interleaver".


Fourier coefficients of product of two periodic signals question:
If $x(t)$ and $y(t)$ are two periodic signals(both with period T) with Fourier coefficients $c_{n}$ and $d_{n}$ respectively then, Fourier coefficient of $z(t)=x(t)\cdot y(t)$ is:

(a) $\dfrac{c_{n}d_{n}}{T}$
(b) $T c_{n}d_{n}$
(c) $c_{n}d_{n}$
(d) $c_{n}*d_{n}$

my attempt :

let, $x(t)=\sum_{-\infty}^{\infty} c_{n}\cdot e^{jn\omega_{0}t}$ and $y(t)=\sum_{-\infty}^{\infty} d_{n}\cdot e^{jn\omega_{0}t}$.

Then, fourier cofficient of $z(t)$ will be equal to

$$\begin{align}
\dfrac{1}{T}\displaystyle\int_{0}^{T} z(t) e^{-jn\omega_{0} t}dt
&=\dfrac{1}{T}\displaystyle\int_{t=0}^{T} x(t)\cdot y(t) e^{-jn\omega_{0} t}dt \\
&=\dfrac{1}{T}\displaystyle\int_{t=0}^{T} \left(\sum_{l=-\infty}^{\infty} c_{l}\cdot e^{jl\omega_{0}t}\sum_{m=-\infty}^{\infty} d_{m} \cdot e^{jm\omega_{0}t} \right)e^{-jn\omega_{0} t}dt\\
&= \dfrac{1}{T}\displaystyle\int_{t=0}^{T}\left(\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}c_{l}d_{m}\cdot e^{j(l+m)\omega_{0}t}\right)e^{-jn\omega_{0}t}dt
\end{align}$$

[in above $l$, $m$ are dummy variables]

after changing order of integration and summation i got this:

$$\begin{align}
\displaystyle\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}c_{l}d_{m}\left[\dfrac{1}{T}\int_{t=0}^{T}e^{-j(n-l-m)\omega_{0} t}dt\right]=\displaystyle\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}c_{l}d_{m}\cdot\delta[n-l-m]\end{align}$$

where

$$\delta[n] \triangleq \begin{cases}
1 \qquad & \text{for } n = 0 \\
0 \qquad & \text{for } n \in \mathbb{Z}, n \ne 0 \\
\text{undefined} \quad & \text{otherwise} \\
\end{cases}$$
is the Kronecker delta.

But from here on I'm unable to proceed,

i know answer is option (d)i.e, $c_{n}*d_{n}$ but i am interested in proving the result formally.

You did everything right, and you're almost there. Your last equation becomes

$$\sum_lc_l\sum_md_m\delta[n-l-m]\tag{1}$$

Due to the Kronecker delta the sum over $m$ reduces to a single element:

$$\sum_md_m\delta[n-l-m]=d_{n-l}\tag{2}$$

because $\delta[n-l-m]$ is non-zero for $m=n-l$ and zero otherwise. Plugging $(2)$ into $(1)$ gives the desired result:

$$\sum_lc_ld_{n-l}$$


Making sense of Haar transform and inverse Haar transform Yesterday, my teacher taught us about the Haar transform (for pictures). We had to implement it in Python and when I finished I was wondering: what was the point?
I can't explain the purpose of what I'm studying: what a shame !

Can someone explain me what is the main purpose of this algorithm for images ?

I would like to tell myself that I can explain everything I study to a random person.

I fully agree with your option: I consider that "to learn about something, you should teach it first". The Haar transform is important at several levels, here are a few ones:

As an orthogonal basis, it was thought to provide useful (in some strict mathematical sense) expansions of functions
As a 2D dimensional discrete operator, it emulates an oriented singularity detector, allowing to reveal horizontal and vertical orientations
As a regularity absorber, it deals well with locally linear variations
As a multiscale device, it provides an invertible way of looking at data in a multiresolution power-of-two scheme
As an algorithmic tool, it mostly requires adds and subtracts, hence it is very fast and cheap to compute.

On only one level, the result looks as below

The original image is reversibly cast to a lower resolution version, and three panels exhibiting singularities. They are quite grayish, showing that locally, pixels don't vary that much, except around edges. And you can iterate this again and again.

For locally smooth images with edges, (opinion) it is the most natural idea to represent, detect and compress information at different levels, with very low cost.


Overlapping power-of-cosine windows with constant weight I would like to refer you to a relevant question I posed on the Math StackExchange: A generalisation of the Pythagorean trigonometric identity.

As can be seen from the question and the answer, what this solution means is that the $(2n)$th-power-of-cosine window will exhibit a constant weight of $2n\binom{2n}{n}4^{-n}$ with an overlap of $1 - \frac{1}{2n}$, i.e.

$2$nd-power-of-cosine window exhibits constant weight of $1$ with overlap $50\%$,
$4$th-power-of-cosine window exhibits constant weight of $\frac{3}{2}$ with overlap $75\%$,
$6$th-power-of-cosine window exhibits constant weight of $\frac{15}{8}$ with overlap $83.\overline{3}\%$, etc.

I trust it is clear without explanation why constant-weight overlapped window schemes are great for application in signal analysis-synthesis. As such, it seems to me this would be an important result in signal processing and, since I am sort of a newbie to signal processing, I find it unlikely that I have discovered something important - it is rather more likely that I am just "reinventing the wheel".

However, I couldn't find anything relevant online - for example, the most relevant result in a Google search for "overlap power of cosine window constant weight" seems to be a page which, in passing, mentions that:

  [w]hen four-times-overlapping Blackman windows are plotted (...) a neatly constant sum appears again(...)

So, where can I find material which discusses issues similar to above in detail and, specifically, is there a paper which provides the above result?

your result is interesting but you should take a look at :

F. J. Harris, "On the use of windows for harmonic analysis with the discrete Fourier transform," in Proceedings of the IEEE, vol. 66, no. 1, pp. 51-83, Jan. 1978.

doi: 10.1109/PROC.1978.10837

Abstract: This paper makes available a concise review of data windows and their affect on the detection of harmonic signals in the presence of broad-band noise, and in the presence of nearby strong harmonic interference. We also call attention to a number of common errors in the application of windows when used with the fast Fourier transform. This paper includes a comprehensive catalog of data windows along with their significant performance parameters from which the different windows can be compared. Finally, an example demonstrates the use and value of windows to resolve closely spaced harmonic signals characterized by large differences in amplitude.
keywords: {Discrete Fourier transforms;Fourier transforms;Frequency;Harmonic analysis;Oceans;Parameter estimation;Signal processing;Signal resolution;Signal sampling;Smoothing methods},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1455106&isnumber=31261

There a number of Figures of Merits (FOMs) that are discussed for windows including $cos^k(),\;
 k=\{1,2,3,4\}$ but not at the generalized overlaps you show.  Looking at the main table in Harris's paper, as $k$ increases, many of the FOMs like the ENBW get worse.  Also while the constant energy is preserved for your result, the coherent gain of the DFT decreases.  The only FOM that touches on overlap is Overlap Correlation.

If you may permit me being blunt, but $cos^{20}$ as windows go, doesn't seem to offer much advantage in general terms. I think that your result is interesting and someone might publish it, but you might hold off going down to the Porsche dealership.


Air-coupled microphone Does anyone know what the heck is a coupled-air microphone? I know that they were used for recording internal body sounds in papers like this and this. I googled the term, but I couldn't find any circuitry schema or place where to buy one. Are they some kind of diaphragm, condenser or contact microphone?

Any help would be greatly appreciated.

Problem you are aiming at is connected with Phonocardiography (recording bodily sounds). Traditionally contact sensors (such as stethoscopes) were used for this task, but problem is that you must raise and lower the device between two measurements. In case of air-coupled microphone (meaning that it is not touching the body) you can move it constantly, performing a continuous measurement.

Obviously there is some disadvantage, such that it is prone to environmental noise, but shielding (tubes) are used, as well as microphone arrays for noise cancellation. Problem of coupling the membrane with patient body is very hard to tackle, as big mismatch exists between air characteristic impedance and body impedance, but using a conical shape is helping to transform these impedances.

What kind of microphones were being used? These are almost the same used widely in recordings:

Condenser
Electret
Piezoelectric

Some of these were adapted for recording of infra-sounds, i.e. Brüel & Kjaer, model 4117 or other great, 4135. Although people are also using quite cheap electret microphones for this task.

This article can shed some light on that:

  Algernon D., et al. - Signal Processing for Air-Coupled Impact-Echo
  using Microphone
  Arrays.


Does constant Q transform have linearity property in the transformed domain? With FFT, I sometimes take advantage of the fact that I can pre-calculate a signals Fourier Transform, and then you can add the noise in the spectral domain:

$$ \mathscr{F} (x[n] + h[n]) = \mathscr{F}(x[n]) + \mathscr{F}(h[n]) $$

This is useful when we want to discard of the original signal, but we still want to noise it with a spectrum. If we know the spectrum of the noise beforehand, this is also faster computationally (just adding together the two spectra).

Does the same relation hold for constant Q transforms?

It’s equivalent to a bank of linear filters, so yes it is linearly additive because each component filter is linearly additive with respect to its inputs/outputs.


Elegant way of calculating minimal sampling frequency without aliasing for annoying signals Given a symmetric Fourier transform with finite support as below, there's an algorithmic procedure to determine the minimal frequency for sampling without aliasing:

Find an integer $L>1$ satisfying $f_2=L(f_2-f_1)$
The desired frequency $f_s$ is given by $\frac{2f_2}L$.
If no such $L$ exists, expand to the right or left and repeat.

The problem I was given is to find this minimal frequency for the signal below.

However, I can only manage to find the answer (which is 4) manually. Is there some algorithmic procedure again to find the answer?

Let me start by pointing out that for a band pass signal with a single band, there is no need for an iterative procedure. Taking your first figure as an example, in order to avoid aliasing, the sampling frequency $f_s$ must satisfy

$$\frac{2f_2}{n+1}\le f_s\le \frac{2f_1}{n}\tag{1}$$

for some integer $n$. The largest integer $n$ for which the range defined by $(1)$ is non-empty results in the smallest possible sampling frequency. This largest value of $n$ is given by

$$n_{max}=\left\lfloor\frac{f_1}{f_2-f_1}\right\rfloor\tag{2}$$

where $\lfloor\cdot\rfloor$ denotes the floor function. So you just need to determine $n_{max}$ from $(2)$ and then evaluate $(1)$ with $n=n_{max}$ to determine the possible range for $f_s$.

For multi-band signals, you usually need an iterative procedure. Let's take your second figure as an example (assuming that the spectrum is only shown for positive frequencies, and that it is symmetrical around $f=0$). Let's define $f_0=1$, $f_1=5$, and $f_2=6$. First consider the band pass component in the range $[f_1,f_2]$. From $(2)$ we get $n_{max}=5$. From $(1)$ we get the following "range" for $f_s$:

$$2\le f_s\le 2$$

Let's also evaluate the other ranges for values $n<n_{max}$:

$$\begin{align}n=4:&\quad 2.4\le f_s\le 2.5\\
n=3:&\quad 3.0\le f_s\le 3.3\\
n=2:&\quad 4.0\le f_s\le 5.0\\
n=1:&\quad 6.0\le f_s\le 10.0\end{align}\tag{3}$$

Let's now consider the low pass component. In any case, due to $f_0=1$, the minimum sampling frequency must satisfy $f_s\ge 2$. However, the aliases of the low pass component must not interfere with the band pass component. In a manner completely analogous to the single band band pass case, this leads to the following inequalities:

$$kf_s+f_0\le f_1,\quad (k+1)f_s-f_0\ge f_2$$

resulting in

$$\frac{f_0+f_2}{k+1}\le f_s\le \frac{f_1-f_0}{k}\tag{4}$$

The maximum value of $k$ for which the range defined by $(4)$ is non-empty is

$$k_{max}=\left\lfloor\frac{f_1-f_0}{2f_0+f_2-f_1}\right\rfloor\tag{5}$$

For the given values this gives $k_{max}=1$, resulting in the range

$$3.5\le f_s\le 4.0\tag{6}$$

Comparing $(3)$ and $(6)$, the only overlap is $f_s=4.0$, as you've already figured out.

Since this "overlap" is only a single frequency at the edge of two ranges, this solution is a purely theoretical one. Any imperfection in the implementation will cause aliasing. In practice, you would choose a larger sampling frequency. Considering $k=0$ in $(4)$ (meaning that the upper limit becomes infinite), we get $f_s\ge 7$. Comparing this with the ranges in $(3)$, we can conclude that a practical sampling frequency should satisfy $7+\delta<f_s<10-\delta$ with some reasonably chosen $\delta>0$.

In sum, for a signal composed of a low pass and a band pass signal (as in your example), the procedure is as follows:

determine the maximum values $n_{max}$ and $k_{max}$ from $(2)$ and $(5)$
determine the possible ranges for $f_s$ from $(1)$ and $(4)$ for decreasing values of $n$ and $k$, until you find overlapping ranges.

A general iterative procedure for finding the minimum sampling frequency which guarantees zero aliasing for arbitrary multi-band signals can be found in this paper.


What is the name of this color transformation? This is a question regarding color transformation in images. I have this color transformation matrix I am using to convert an RGB image to a color space whose name I do not know:

T = [(1/3) (1/3) (1/3); (1/2) 0 (-1/2); (-1/2) 1 (-1/2)]

What is the name of this color transform? I think I once saw somewhere that it is called KL transform because it was close to a KLT computed over a large collection of images, but I am not sure...

I do not know the name of your transformation $T$, but it can't be the Karhunen-Loeve transform: The Karhunen-Loeve transform is image specific and its transformation matrix is a rotation matrix (see for example the KLT in Mathematica).

Given the values in your matrix, it seems to attempt to separate between intensity (first row) and color (the other two rows, spanning the remaining hyperplane) while being  optimized for speed.


Does the Nyquist frequency of the Cochlear nerve impose the fundamental limit on human hearing? The bandwidth of human hearing by empirical data is $20 \; Hz$ to $20 \; kHz$. A cochlear implant stimulates the auditory or acoustic or Cochlear  nerve directly so that the hearing can be improved in the case of stimulation mechanism upstream of the Cochlear nerve has degraded.
Let us assume that the ear mechanism has not degraded (such as in a young and healthy adult). The Cochlear implant can likely improve hearing, even in this case, by increasing the bandwidth by amplifying the effect of the ear drum vibration (sensor actuation). However the neurons connecting the Cochlear nerve to the hearing region of the brain have an upper limit on the sampling rate on the order of $1 \; kHz$.

Does the the Nyquist sampling theorem limit the superhuman hearing and sound localization capability made possible by a Cochlear implant?

Does the Nyquist frequency of the Cochlear nerve impose the fundamental limit on human hearing?

No.
A quick run-through the human auditory system:

The outer ear (pinnae, ear canal), spatially "encodes" the sound direction of incidence and funnel the sound pressure towards the
ear drum, which converts sound into physical motions, i.e. mechanical energy
The middle ear (ossicles) is a mechanical transformer (with some protective limiting built-in) that impedance matches the air-loaded ear drum to the liquid-loaded oval window of the
Cochlea (inner ear). The vibration excites a bending wave on the basilar membrane. The membrane is highly resonant and transcodes frequency into location: for any given frequency the location of the resonance peak is in a different spot. High frequencies wiggle very close to the oval window, low frequencies towards the end of it. This motion is picked up by the
Cochlea neurons, which transmit the intensity of the excitation at their location to the brain. About 20% of the neurons are efferent (come out of the brain) and are used to actively tune the resonance with a feedback loop (which causes tinnitus if misadjusted)

So in essence the Basilar Membrane performs sort of a mechanical Fourier transform. The frequency selectivity of the Neurons is NOT determined by the firing pattern but simply by their location. A neuron at the beginning of the basilar membrane is sensitive to high frequencies and a neuron at the end detects low frequencies. But they are more or less the same type of Neurons.
The Nyquist criteria doesn't come into play at all since no neuron is trying to pick up the original time domain waveform. The couldn't anyway: human neurons have a maximum firing rate of less than 1000 Hz and average firing rates are way below that. The firing rate of a cochlea neuron represents "Intensity at a certain frequency" where that frequency is determine by the location of that specific neuron.
So you can think of it as a short term Fourier Transform. Instead of a single time domain signal you get a parallel stream of frequency domain signals where each individual signal has a much lower bandwidth.
A cochlea implant basically does the short term Fourier transform internally and then connects the output for each frequency range to the "matching" neurons in the cochlea nerve. Theoretically you can create ">20 kHz" hearing with an implant that can actually receive and process higher frequencies and simply routes them to existing neurons, i.e. you could feed 40 kHz activity to the 10 kHz Neuron. The human would have a sensation when exposed to 40 kHz but it's unclear what they could do with that: they would have "relearn" how to hear. Aside from the highly questionable practical and ethical issues, it probably wouldn't be useful. In order to get to 40 kHz you'd have to give some other frequencies, and chances are that evolution has chosen the current "normal" range for humans pretty carefully.


Develop a simple image classification system by using Laplacian of Gaussian (LoG) I am tring to write a simple CAD system to classify some images within two groups by using Laplacian of Gaussian (LoG). I am using scikit-image for this tasks and I want to use DNN in keras to train a classifer.

The problem is that I am not sure how to create feature list by using LOG.

My images could have different sizes so LOG could have different feature set.

I am using below codes

from math import sqrt
from skimage import data
from skimage.feature import blob_dog, blob_log, blob_doh
from skimage.color import rgb2gray

img1T = Image.open( "C:\Users\Public\h1T.jpg")
img2T=Image.open( "C:\Users\Public\h2T.jpg")
data_img1T = np.asarray( img1T , dtype='uint8' )
data_img2T= np.asarray( img2T, dtype='uint8' )
data_gray_img1T = rgb2gray(data_img1T)
data_gray_img2T = rgb2gray(data_img2T)

blobs_log_img1T = blob_log(data_gray_img1T, max_sigma=30, num_sigma=10, threshold=.1)
blobs_log_img2T = blob_log(data_gray_img2T, max_sigma=30, num_sigma=10, threshold=.1)

when getting the LOG feature for tow images:

>>> blobs_log_img1T.shape
(73L, 3L)
>>> blobs_log_img2T.shape
(51L, 3L)

now these two images have different shapes, but I need to find a way to use identical number of features for all images. Do I have to resize images to be of identical sizes or I have to take another approach?

Thanks,

The two simplest ideas are probably to either pad the smaller images with a constant pixel value to a common shape, or crop the bigger images to a common shape.

In the first case, you will however likely end up with a network biased towards padded images/samples with constant pixel value areas. In the second case, you loose possibly prominent information with respect to your classification task.

You could also use the cropping approach with several crops per image in order to avoid loss of information in your data set, leading to an increased proportion of the classes associated with bigger images in your training set. This also biases your neural network (this could be more or loss controlled thanks to appropriate weights attached to the classes' samples during training, if you can predict the imbalance before training, cf. this towardsdatascience article to read about training with imbalanced classes, with Keras code).

Those questions were raised and discussed in this SO post.


what is the frequency selective channel and doubly-selective channel Could anyone please help me to explain the frequency-selective channel and doubly-selective channel and the difference between them  in brief and easy way. I have tried to check on internet,  I couldn't get a clear explanation about them.

Example. For a channel that can be modeled by LTI system, we send an impulse-like signal and receive channel impulse reponse $h(\tau)$ where $\tau$ is delay.

Frequency-flat channel means $h(\tau) \sim \delta(0)$ or for discrete-time version, $h[n] \sim \delta_0$.

Frequency-selective channel means $h(\tau) \neq \delta(0)$ or for $\tau > 0,h(\tau) \neq 0$.

Now if we do the measure at time $t=t_1$ and $t=t_2$, $h(\tau,t_1) = h(\tau,t_2)$, the channel is called non time-selective (anyone remember the name?).

If we do the measure at time $t=t_1$ and $t=t_2$, $h(\tau,t_1) \neq h(\tau,t_2)$, the channel is time-selective.

Doubly selective is simply both time-selective and frequency-selective.


How to outline a defect in an image? The images given below are a result of fruit defect detection algorithm i have made myself...

Just wanted to know is there a way of outlining these defects?

Original image 1:

1st Image:

Original image 2:

2nd Image:

What i want to achieve is this :

        =>      2.

As far as I understand what you call defects are detcted by changes in color of the fruit.

I would proceed as follow:

edge detection algorithm, (as you did), to detect the fruit region in the image.
In this region, compute the median color, that is the color that matches the most pixels (for example, it would be some kind of yellow for a banana)
For each  pixel in the fruit region, check the distant between the median computed at step 2 and the color of the pixel. If the distance is high, mark the pixel as defect
At the end, apply some morphological operator to remove outliers

Hope this helps.


Can you predict how long a discontinuity will take to reach its target if it is low pass filtered at a certain frequency? Let's say you have a signal that relatively slowly rises from 0 to 1, then suddenly from one sample to the next drops back to 0, creating a sharp discontinuity. Let's say you run it through a one pole low pass filter at a given frequency, say 200 Hz to smooth it.

Is there any way to predict or graph or calculate how long it will take after the discontinuity for the signal to reach 0 or near 0?

Infinitely long. It's an exponential decay, so it never goes back to 0.

The time constant for a 200 Hz filter is $\frac{1}{2 \pi f} \approx 0.8ms $. About every 0.8ms it drops by about $2/3$ or every 1.84 ms it drops by a factor of 10.

If you define "almost zero" as one millionths you need to wait $6*1.8ms$ or about 11ms.


what does "fixed multipath channel model" means? I want to check the effect of this channel model on my signal transmitted with "impulse response-Ultra Wide Band" transmission (use very short pulses ~0.5 ns each), but i don't know where to find this channel impulse response or how can I generate it?
if there is only one line of sight (LOS) path and one non-line of sight (NLOS) path, if I assume the channel impulse response something like h=[1 0 0 0 .4 0 0] and then filter(h,tx signal), it would be correct?
best regards

Your "h" is a valid multipath channel impulse response.  It is a special case, though, in that it represents a multipath delay that is an exact multiple of your sample period.  Usually the delay will not be an exact multiple of the sample period.

I am not an expert on multipath channels, so hopefully someone else can chime in, but I think that the easiest way to generalize your multipath channel is to make the gain of the multipath complex rather than real.  The phase of the gain would alter the phase of your signal, simulating a non-integer delay.

If you really wanted to be hardcore about it you could use this answer on how to make fractional delay filters to create a filter that delays the signal by as much or as little as you want, and then add it to [1 0 0 0 0...] (the LOS channel) to get the complete multipath channel.


OFDM with beamforming I'm trying to compute the the interference that one OFDM user with beamforming causes to another user. The two systems are fully overlapped in frequency and time. I know how to generate an OFDM signal in MATLAB and compute its BER in AWGN and fading channels, I just don't know how to implement beamforming. The beamform is a circular sinc, I believe I need multiply my signal with the gain of the main lobe of the sinc function. I don't know how to calculate the interference caused on the other user. I'm not asking for anyone to provide me with a code, I'm just hoping that someone here can provide some directions on how to approach this or if there are any other resources that I can check.

The interference you cause to the other user is determined by:

The total transmit power of your transmitter
The gain of your beampattern in the direction of the other user
The path loss to the other user
If the other user uses receive beamforming, then also by the gain of its beampattern in the direction of the transmitter

So if you transmit with a power of 10dBm, have a gain in the users direction of -20dB and a pathloss of 60dB, the other user receives -70dBm of interference power from the transmitter. Whether or not this is bad then depends also on how powerful the useful signal of the other user is.

More information is required here for a more detailled answer.


How to choose the "best" measurment (from a given set) as input for a kalman filter? Problem:
For an object tracking scenario with multiple objects and multiple tracks, I want to choose the "best" assignment object<->track. Therefore, I need a parameter indicating how suitable an object is for a track.
(Each object refers to a measurement and each track is a Kalman filter)

Idea:
Given the classical Kalman filter with the matrices state $x$ and measurement transition $H$
We can calculate the residual $y_t$ between the predicted measurement $x_{pred}$ and the measurement $z$ by:

$y_t = z-H x_{pred}$

However, this will be a vector. Is there a general approach to reduce this to one parameter? (e.g distance when having just position information). I assume that some kind of normalisation is needed to get all the residuals for each state in the same range? Building the l2-norm would be my first approach.

Question:
Which parameter is suitable to indicate how "good" the measurement fits to the Kalman filter?

Is the residual suitable? If so, how to handle different value ranges, and how to "compress" the vector to just one value?

Question: Which parameter is suitable to indicate how "good" the
  measurement fits to the Kalman filter?

To estimate a quality of association you can use likelihood function. The likelihood considers not only residual but also uncertainty and represented as scalar value:

$$\mathcal{L} = \frac{1}{\sqrt{2\pi S}}\exp [-\frac{1}{2}\mathbf y^\mathsf T\mathbf S^{-1}\mathbf y]$$

  Problem: For an object tracking scenario with multiple objects and
  multiple tracks, I want to choose the "best" assignment
  object<->track.

The question which you are interested is called data association and as was said there are a lot of solution and this question is too broad. Here are some of methods (you can google it and find more information):
 - Global nearest neighbour (in Euclidean space);
 - Global strongest neighbour;
 - Joint probability data association;
 - Multiple hypothesis tracking (the best choice).


Choice of Gaussian kernel parameters when lowpass filtering before image resampling? I need to decimate a signal by a factor of q.

More specifically my signal is a 3D "image": $\ I(x_i,y_j,z_k)$,
which I need to downsample by a factor of two in the z direction.

I want to do lowpass filtering before decimation by convolving with a Gaussian kernel of size n.

I create my Gaussian kernels 2 standard deviations below and above 0 since this accounts for 95 % of the distribution.

I am looking for a rule of thumb that tells me how large n should be.

Am I right in thinking that the Gaussian should filter out all frequencies above fN/q,
where fN: Nyquist frequency of original signal?

I know that the Fourier of a Gaussian with standard deviation $\ \sigma $ is another Gaussian with standard deviation  $\ \sigma^* =1/\sigma $. I am not sure how strict my lowpass filter should be. Should the cut frequency, fc, be at 2 or 3 standard deviations?

What is then the equation for the cut frequency, fc, of a Gaussian kernel with size n: fc(n)=?

Below is the frequency response of some Gaussian kernels calculated in Matlab:

My actual problem involves q = 2, and from this figure I see that n=5 should work nicely.
It would have been good to have a rule of thumb do, so I did not have to do this for each q I encounter.

I guess what you want to do is similar to image resize. If so, then multiple build-in functions in MATLAB can be used for this purpose, for example imresize.


Conjugate Pole Pairs I am currently studying for an exam for a Speech Processing course I am taking and stumbled across this question in one of the sample exams:

I began to plot the pole-zero configuration placing my poles on the unit circle at points +j and -j as well as placing my zero in the 0.5 location.

When I compared my plot to his, I noticed that he placed the poles instead in the halfway point of the Imaginary axis and was wondering why he decided to place them there?

Here is my professors answer:

Thanks

If the dotted-line circle is the unity circle and both the real and imaginary axis are on the same scale, then your professor's answer is wrong.

Your solution is the correct one. It would be something like this:

where the dotted-line circle has a radius $r=1$.


Is it meaningful to find linear convolution of just two random sequences? Could it bring any meaning? Is this kind of convolution be useful is solving anything?

No, if the two sequences are random then convolving them is not useful.


Difference between causality and memorylessness I have found the particular definitions.

Causality means that the output of the system does not depend on future inputs, but only on past input.

Memory-less - does not depend on previous values of either input or output.

I am looking for example of systems which are either non-casual or memoryless but not both.

Also i want to check whether I understand it well.

  $ T(x[n]) = x[n − n_{0}] $

  $T(x[n]) = ax[n] + bx[n − 3]$

Are those an examples of non-casual and memory-less systems?

A causal system does not need to know the future in order to compute its output. A memoryless system computes the output only from the current input. A memoryless system is always causal (as it doesn't depend on future input values), but a causal system doesn't need to be memoryless (because it may depend on past input or output values).

The system $$y[n]=x[n]+2x[n+1]$$ is non-causal because it needs to look into the future (by $1$ sample) to compute its output. The system $$y[n]=3\big(x[n]\big)^2$$ is memoryless (and necessarily causal) because it only needs the current input sample $x[n]$ to compute the output.

The systems in your question are both causal and have memory (if $n_0>0$).


Understanding LTI systems graphically I'm currently working on this problem from Oppenheimer's book:

Given $x_{1}(t)$ and $y_{1}(t)$ I should figure out $y_{2}(t)$ and $y_{3}(t)$ given that the system is LTI.
My progress so far:

I'm specifically asking about the graph for $y_{3}(t)$ because it looks like my solution doesn't agree with what's in the solution manual as for $0<t<1$ it's $y_{3}(t)=2$ instead of $y_{3}(t)=4t$ so am I missing anything?

Suppose that you could write $x_.$ as $x_.(t) = a_1x_1(t-n_1)+a_2x_1(t-n_2)$, then by LTI hypothesis, you can derive $y_.(t) = a_1y_1(t-n_1)+a_2y_1(t-n_2)$ with very little computation.
You can get a first insight to that method  by looking at the shape of $x_2$ which is self-evident, and a little more imagination gives you the structure for  $x_3$.
I find that this exercise is very pedagogic in the progression of difficulties. Some other exercises of the book might have the same gradation, don't hesitate to look back before jumping into computations at first.


How are 12 outputs being reconciled so that a Homomorphic filter could be applied? Take a look at this paper. It writes in the section 4.1:

  We selected 12 not overlapping filters, to analyze 12 different directions, rotated with respect to 15° each other. A homomorphic filter is then applied to enhance lines and to produce a dark uniform background. Finally a threshold is applied to obtain a binary mask. The threshold value is computed as mean plus standard deviation of the image intensity (see Figure 3(b)).

I could not understand this.

They began the process with a single image. Then, they applied 12 filters to that input image so that they can produce 12 outputs from one image.

How are they selecting one output image from those 12 outputs so that they can apply the homomorphic filter?  I mean, out of those 12 outputs, to which output image they would apply the homomorphic filter?

If they apply homomorphic filter and thresholding to 12 images, when did those 12 images become one image?

From my reading, the authors don't combine the 12 images in a single image: each generates a directional scratch mask, each of them being processed independently. If I read too fast, tell me. Indeed, they first remove a low-pass version, and apply a logarithmic homomorphic filter. But you could combine the binary masks with an OR to get a single image of potential scratches.

However, it would be possible with filters alone to recombine them by addition.

Suppose that you only have two filters $f_1$ and $f_2$, that are applied to an image $i$. $\mathcal{F}$ is the Fourier transform: $I=\mathcal{F}(i)$ and $F_k=\mathcal{F}(f_k)$.

The hypothesis that they do not overlap is related to their frequency behavior, with spectra $F_1$ and $F_2$. Their supports $S_{F_1}$ and $S_{F_2}$ (imagine 2D patches in the full 2D frequency domain $S_F$) do not overlap if their intersection is void: $S_{F_1} \cap S_{F_2}= \emptyset$. What is not explicit is that, often, the filters are also complementary, their union covers the whole frequency domain: $S_{F_1} \cup S_{F_2}= S_F$. Thus, the spectrum $I$ of $i$ is cut into two complementary parts, and $I = F_1\cdot I+F_2\cdot I$. Then, you could recover $i$ with
$$i = \mathcal{F}^{-1}(F_1\cdot I+F_2\cdot I) = \mathcal{F}^{-1}(F_1\cdot I)+\mathcal{F}^{-1}(F_2\cdot I)\,.$$

Equivalently, this amount to directly write, in the original domain:

$$i = f_1\star i + f_2\star i\,.$$

Visually, the frequency domain is a pie, and each filter is a slice. When you have the the slices, you can reproduce the original pie. In practice, people don't use exactly non overlapping filters in frequency, but orthogonal filters in space, that sum to unity to equally cover all frequencies.  This is a little more general as non overlapping filters, as $F_1\cdot F_2 = 0$, thus $f_1\star f_2 = \mathcal{F}^{-1}(F_1\cdot F_2) = 0$,  are a subset of orthogonal filters.


What Is the Difference Between the Terms Color Intensity and Color Saturation? In Digital Image Processing, what is the difference between the intensity and the saturation of a color?

How would you explain the difference to a newbie?

When you take an RGB Image matrix and convert the color into HSV Color Model the color is represented on Cylinder.
Now, the intensity (Lightness / Value) is the height on this Cylinder which is going from black to white and basically sets the Gray Color of the neutral color (One which blends RGB in the same intensity).

Saturation is the Radius and intuitively is how far the color is from being neutral (Pure Gray).
Namely, how dominant is one color component versus the others.