How to detect corners in a binary images with OpenGL? I have binary 160x120 images such as:

I would like to detect corners of those white blobs. They are previously closed by mathematical morphology so there shouldn't have any inner corners. In this specific case, I would want 16 corners, like:

My first attempt was using some OpenCV functions like goodFeaturesToTrack or FAST but they are particularly slow (plus FAST is very unstable). My idea would be to do such a computation on the GPU, as my source image comes from it. I looked for ideas on the web on how to write such shaders (I'm using OpenGL ES 2.0), but found nothing concrete. Any idea how I could start such an algorithm?

I just happened to be implementing something like this on OpenGL ES 2.0 using Harris corner detection, and while I'm not completely finished, I thought I'd share the shader-based implementation I have so far. I've done this as part of an iOS-based open source framework, so you can check out the code if you're curious as to how some particular step works.

To do this, I use the following steps:

Reduce the image to its luminance values using a dot product of the RGB values with the vector (0.2125, 0.7154, 0.0721).
Calculate the X and Y derivatives by subtracting the red channel values from the pixels left and right and above and below the current pixel. I then store the x derivative squared in the red channel, the Y derivative squared in the green channel, and the product of X and Y derivatives in the blue channel. The fragment shader for this looks like the following:

precision highp float;

varying vec2 textureCoordinate;
varying vec2 leftTextureCoordinate;
varying vec2 rightTextureCoordinate;

varying vec2 topTextureCoordinate;
varying vec2 bottomTextureCoordinate;

uniform sampler2D inputImageTexture;

void main()
{
 float topIntensity = texture2D(inputImageTexture, topTextureCoordinate).r;
 float bottomIntensity = texture2D(inputImageTexture, bottomTextureCoordinate).r;
 float leftIntensity = texture2D(inputImageTexture, leftTextureCoordinate).r;
 float rightIntensity = texture2D(inputImageTexture, rightTextureCoordinate).r;

 float verticalDerivative = abs(-topIntensity + bottomIntensity);
 float horizontalDerivative = abs(-leftIntensity + rightIntensity);

 gl_FragColor = vec4(horizontalDerivative * horizontalDerivative, verticalDerivative * verticalDerivative, verticalDerivative * horizontalDerivative, 1.0);
}

where the varyings are just the offset texture coordinates in each direction. I precalculate these in the vertex shader to eliminate dependent texture reads, which are notoriously slow on these mobile GPUs.
Apply a Gaussian blur to this derivative image. I used a separated horizontal and vertical blur, and take advantage of hardware texture filtering to do a nine-hit blur with only five texture reads on each pass. I describe this shader in this Stack Overflow answer.
Run the actual Harris corner detection calculation using the blurred input derivative values. In this case, I'm actually using the calculation described by Alison Noble in her Ph.D. dissertation "Descriptions of Image Surfaces". The shader that handles this looks like the following:

varying highp vec2 textureCoordinate;

uniform sampler2D inputImageTexture;

const mediump float harrisConstant = 0.04;

void main()
{
 mediump vec3 derivativeElements = texture2D(inputImageTexture, textureCoordinate).rgb;

 mediump float derivativeSum = derivativeElements.x + derivativeElements.y;

 // This is the Noble variant on the Harris detector, from
 // Alison Noble, "Descriptions of Image Surfaces", PhD thesis, Department of Engineering Science, Oxford University 1989, p45.
 mediump float harrisIntensity = (derivativeElements.x * derivativeElements.y - (derivativeElements.z * derivativeElements.z)) / (derivativeSum);

 // Original Harris detector
 //     highp float harrisIntensity = derivativeElements.x * derivativeElements.y - (derivativeElements.z * derivativeElements.z) - harrisConstant * derivativeSum * derivativeSum;

 gl_FragColor = vec4(vec3(harrisIntensity * 10.0), 1.0);
}

Perform local non-maximum suppression and apply a threshold to highlight the pixels that pass. I use the following fragment shader to sample the eight pixels in the neighborhood of a central pixel and identify whether or not it is the maximum in that grouping:

uniform sampler2D inputImageTexture;

varying highp vec2 textureCoordinate;
varying highp vec2 leftTextureCoordinate;
varying highp vec2 rightTextureCoordinate;

varying highp vec2 topTextureCoordinate;
varying highp vec2 topLeftTextureCoordinate;
varying highp vec2 topRightTextureCoordinate;

varying highp vec2 bottomTextureCoordinate;
varying highp vec2 bottomLeftTextureCoordinate;
varying highp vec2 bottomRightTextureCoordinate;

void main()
{
    lowp float bottomColor = texture2D(inputImageTexture, bottomTextureCoordinate).r;
    lowp float bottomLeftColor = texture2D(inputImageTexture, bottomLeftTextureCoordinate).r;
    lowp float bottomRightColor = texture2D(inputImageTexture, bottomRightTextureCoordinate).r;
    lowp vec4 centerColor = texture2D(inputImageTexture, textureCoordinate);
    lowp float leftColor = texture2D(inputImageTexture, leftTextureCoordinate).r;
    lowp float rightColor = texture2D(inputImageTexture, rightTextureCoordinate).r;
    lowp float topColor = texture2D(inputImageTexture, topTextureCoordinate).r;
    lowp float topRightColor = texture2D(inputImageTexture, topRightTextureCoordinate).r;
    lowp float topLeftColor = texture2D(inputImageTexture, topLeftTextureCoordinate).r;

    // Use a tiebreaker for pixels to the left and immediately above this one
    lowp float multiplier = 1.0 - step(centerColor.r, topColor);
    multiplier = multiplier * 1.0 - step(centerColor.r, topLeftColor);
    multiplier = multiplier * 1.0 - step(centerColor.r, leftColor);
    multiplier = multiplier * 1.0 - step(centerColor.r, bottomLeftColor);

    lowp float maxValue = max(centerColor.r, bottomColor);
    maxValue = max(maxValue, bottomRightColor);
    maxValue = max(maxValue, rightColor);
    maxValue = max(maxValue, topRightColor);

    gl_FragColor = vec4((centerColor.rgb * step(maxValue, centerColor.r) * multiplier), 1.0);
}

This process generates a cornerness map from your objects that looks like this:

The following points are identified as corners based on the non-maximum suppression and thresholding:

With proper thresholds set for this filter, it can identify all 16 corners in this image, although it does tend to place the corners a pixel or so inside the actual edges of the object.

On an iPhone 4, this corner detection can be run at 20 FPS on 640x480 frames of video coming from the camera, and an iPhone 4S can easily process video of that size at 60+ FPS. This should be a good deal faster than CPU-bound processing for a task like this, although right now the process of reading back the points is CPU-bound and a little slower than it should be.

If you want to see this in action you can grab the code for my framework and run the FilterShowcase example that comes with it. The Harris corner detection example there runs on live video from the device camera, although as I mentioned the reading back of corner points currently occurs on the CPU, which is really slowing this down. I'm moving to a GPU-based process for this, as well.


Is there a mathematical method to determine if noise is Gaussian? Is there a mathematical method to determine if a signal's noise is Gaussian?

The only way I know so far is to analyze the histogram and layover a Gaussian distribution to visually determine if the distribution is Gaussian. I would like to know if there is a mathematical way to determine if the noise is Gaussian and how accurate the result is.

There are several statistical tests if a time series is Gaussian, although in statistics, the term "tests for normality" is usually how you search for them.

The Nist EDA site is a good place to look and the probability plot is better for shorter data sets than the sample histogram.

http://www.itl.nist.gov/div898/handbook/eda/section3/probplot.htm

Near the bottom of the page, there are references to q-q plots, KS, Chi squared, and other goodness of fit tests.  You can find ample information about them on the web and replicating here isn't going to add anything.

Matlab has qqplot and prob plot in the Statistics toolbox, and the qqplot with a single argument is specific to Gaussian distributions. SAS has all these tests. R has the tests.

I recommend this book, written by 2 Engineers, and they cover several tests including for things like independence, and stationarity.  The book is oriented towards the practical, minimum of mathematics.

  Bendat, Julius S., and Allan G. Piersol. Random data: analysis and
  measurement procedures. Vol. 729. John Wiley & Sons, 2011.

The wrinkle of these tests is that they don't conform to a Signal plus Noise scenario.  The tests generally assume that the time series is all Gaussian or not. A constant mean isn't a problem. Signals are not usually Gaussian and a simple test can't tell the difference.

Signal processing operations such as a DFT, tend to manifest central limit theorem effects on data, so you need to be aware that even linear transformations will not preserve a non Gaussian pdf.

It should be also noted that from a practical perspective, Gaussianity isn't black and white.  Algorithms that have Gaussian assumptions usually work well even if the Gaussianity assumption is not strictly valid.  Things like bi-modality and non-symmetry are more important to know about. Cauchy (heavy tails) like noise and multiplicative noise are also important to know about.


Envelop detection with low sample rate? I have a signal up to $3\textrm{ MHz}$. The ADC that sample it has a rate of $1.5\mu s$. So a full $T$ of the signal is $0.3\mu s$, and I can only sample each $1.5\mu s$. It sounds not enough, but I don't need to reconstruct it, but to create an envelop detector from it, over time, based on lets say $1000T$ (periods).

So, over $1000$ periods of $3\textrm{ MHz}$, I will have: $300 \mu s/1.5 \mu s=200 \textrm{ samples}$. From these $200$ samples, I need to create some envelop detector, or continues curve to then later check if it has some large amplitude changes.

How can I chose where to sample the input signal so that I can get the "right points" of it- means mostly its maxes, where my sample rate is much slower then the signal ?
Should I use a moving average to get this curve ?
Is there another good approach expect from taking more periods?($>1000$)

Is the ADC sample-and-hold fast enough to even capture any envelope peaks or even half cycles?  If not, all bets are off.

If your ADC does have a fast enough capture time, then randomizing your sample times might be a good bet at capturing a few near extrema.  If you know, a priori, the approximate shape of your waveform, then Monte-Carlo sampling might allow one to estimate the waveform amplitude from the statistical distribution of enough sufficiently time-randomized samples.  Assuming the envelope changes slowly enough.


Do signals with a Fourier transform with discontinuities or zero amplitude (in some frequencies) have Laplace transform? I am reading a book on Laplace transform, and in the section on the convergence of Laplace transform for various signals the following theorem is stated, without any proof :

If a signal's Fourier transform be zero in some frequencies, or have discontinuities, it will not have Laplace transform (like sinc function, or periodic signals)

I can not see the proof. For signals with a discontinuous Fourier transform I know that  $\int_{-\infty}^{\infty} |x(t)|dt$ may not converge leading to discontinuities in the Fourier spectra (again, like sinc), but I don't think it will lead to not having a Laplace transform at all.

What is the proof for the above statement? (if correct)

First of all, it is important to distinguish between the unilateral and the bilateral Laplace transform. For causal signals we can use the unilateral Laplace transform. If for such a signal the Fourier transform exists, then also its Laplace transform exists. This is simply the case because here the Fourier transform is equivalent to the Laplace transform for $s=j\omega$, and if the Fourier integral converges, the Laplace transform must exist with a region of convergence $\Re\{s\}\ge a$, for some real-valued $a\le 0$.

If a two-sided signal is considered, we have to use the bilateral Laplace transform. Here the problem is that for $\Re\{s\}>0$ the signal is damped for $t>0$ but amplified for $t<0$ (and vice versa). Consequently, for two-sided signals with a constant envelop (e.g. sinusoids), the Laplace transform does not exist, even if the Fourier transform exists, because the Laplace integral cannot converge for signals with an exponentially increasing envelope. The same is true for all periodic signals (because they are sums of sinusoids). But even for two-sided signals which decay for increasing $|t|$ the bilateral Laplace transform might not exist if they decay too slowly, which is the case for the sinc function (decaying only as $1/t$).


How to Remove Double Lines Detected Along the Edges by Edge Detector? I'm trying to take real time input for hand gestures with web cam, then processing the images to feed them to a neural network. I wrote this processing function to make the hand features look prominent:
def image_processing(image, count):
    roi = image[42:338, 2:298]
    cv2.imwrite('a/'+str(count)+'.png', roi)

    img = cv2.imread('a/'+str(count)+'.png')

    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    blur = cv2.GaussianBlur(gray,(5,5),2)

    th3 = cv2.adaptiveThreshold(blur,10,cv2.ADAPTIVE_THRESH_GAUSSIAN_C,cv2.THRESH_BINARY_INV,11,2)
    ret, res = cv2.threshold(th3, 225, 255, cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
    res = cv2.Canny(res,100,200)
    cv2.imshow("Canny", res)
    cv2.imwrite('a/'+str(count)+'.png', res)

--Edit--
The input and the output images are as follows :

It's obvious that double lines, instead of one, are detected along the edges. I want to make them single. If I apply just Canny edge detection algo, then the edges are not very prominent.

I ran the following code:
clear();
close('all');
mI = im2double(imread('bCfdb.png')); %<! Loading the image
vBlurStd = [0, 0.1, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2];

mII = sum(cat(3, 0.299, 0.587, 0.114) .* mI, 3); %<! Y (Luminosity like channel)

hFigure = figure('Position', [100, 100, 1200, 900]);
hTiledChartLayout = tiledlayout(3, 3);

kk = 0;
for ii = 1:3
    for jj = 1:3
        kk = kk + 1;
        if(kk == 1)
            mE = edge(mII, 'canny');
        else
            mE = edge(imgaussfilt(mII, vBlurStd(kk)), 'canny');
        end
        nexttile();
        imshow(mE);
    end
end

set(hTiledChartLayout, 'TileSpacing', 'tight', 'Padding', 'tight');

This is the result:

What Does the Script?

We set an array of different blur radius to evaluate the algorithm at - vBlurStd.
We load the image and convert it into Double in the range [0, 1] - mI = im2double(imread('bCfdb.png'));
We extract the Y channel from the YCbCr color model - mII = sum(cat(3, 0.299, 0.587, 0.114) .* mI, 3);.
The algorithm is basically blurring the image with Gaussian Blur and then applying Canny Edge Detector - edge(imgaussfilt(mII, vBlurStd(kk)), 'canny');.

The images are left to right, top down by their blur radius.


How is wavelet time & frequency resolution computed? Mallat gives analytic wavelet time & frequency widths/uncertainties as
$$
\begin{align}
\sigma_{ts}^2 &= \int_{-\infty}^{\infty} (t - u)^2 |\psi_{u, s}(t)|^2 dt
               = s^2 \sigma_t^2 \tag{4.51} \\
\sigma_{\omega s}^2 &= \frac{1}{2\pi}\int_{-\infty}^{\infty} \left(\omega - \frac{\eta}{s} \right)^2 |\hat\psi_{u, s}(\omega)|^2 d\omega
               = \sigma_\omega^2 / s^2 \tag{4.54}
\end{align}
$$
where
$$
\begin{align}
\psi_{u,s}(t) = \frac{1}{\sqrt{s}} \psi^* \left( \frac{t-u}{s}\right),\ \ \hat\psi_{u,s}(\omega) &= \sqrt{s} \hat\psi(s\omega) e^{-i\omega u} \tag{4.53} \\
\eta = \frac{1}{2\pi} \int_0^\infty \omega |\hat\psi (\omega)|^2 d\omega \tag{4.52} \\
\end{align}
$$
$$
\sigma_t^2 = \int_{-\infty}^{\infty}t^2 |\psi(t)|^2dt,\ \ \sigma_\omega^2 = \frac{1}{2\pi}\int_0^\infty (\omega - \eta)^2 |\hat\psi(\omega)|^2 d\omega
$$
But how would one actually code this? Continuous math can't just be ported directly; the following considerations arise:

How to treat $\psi(t)$ for wavelets defined in frequency domain? We ifft, which now mixes the discrete and discretized. Any normalization considerations?
How to define $t$ and $\omega$, particularly in light of (1)? What are their ranges? For $\omega$ one can presume $0$ to $\pi$ - but what of $t$? Even if we go to infinity, one must still define the interval over which the wavelet decays.
This answer discusses center frequency computed for wavelets sampled directly in time-domain; is the approach the same here?
How to tell whether the implementation(s) is correct? One can always code integrals to spew numbers.

The key is units, and understanding wavelet behavior in context of application (in this case CWT). Full implementations for all discussed here is available at squeezepy.
This answer assumes analytic wavelets as per the question ($\hat \psi (\omega \leq 0) = 0$), but all ideas apply to general time-frequency atoms (e.g. windows for STFT) with some code adjustments.

Time Resolution
Firstly, note we're computing a "width", which is inversely-related to "resolution". The width is the standard deviation of the time-domain wavelet. Standard deviation is well-defined for normally-distributed data, and has desirable properties if data is a probability distribution.
By squaring modulus of $\psi(t)$ and normalizing it to integrate to 1, we interpret the result as a probability distribution, and define the time variance:
$$
\sigma_t^2 = \frac{\int_{-\infty}^{\infty} t^2 |\psi(t)|^2 dt}{\int_{-\infty}^{\infty} |\psi(t)|^2 dt},  \tag{1}
$$
where $\psi(t)$ is assumed zero-mean. We thus interpret $\sigma_t$ as as the time-span of ~68% of wavelet's energy.
The purpose of the denominator is to enforce validity as probability distribution, i.e. normalizing $\psi$ if we haven't already. This is valuable for coding the integral, as the wavelet at different scales may not integrate to $1$ automatically. It also reminds us of the probabilistic context rooting  $\sigma_t$. Mind that $\sigma_t$ is half the said span, as we compute the uni-lateral ('radius') variance.
The question is, how to define $t$? We can't take infinity, and even if we could, what finite width would span the non-decayed wavelet? Answer to latter is rooted in how DFT works, discussed in "Are we right?" below, but there's an extra concern: should t be zero-mean? No; it should be defined such that it aligns with the wavelet as if it were inputs to generate the wavelet.
Most wavelets will peak at $t=0$, for which t=0 should align with wavelet's peak (i.e. same array index), negative to left, positive to right. Such t will be nearly zero-mean, but not exactly, and requires different handling for even- and odd-length t.

Frequency resolution
Much of what's said for $\sigma_t$ also holds here; we define the frequency variance as
$$
\sigma_\omega^2 = \frac{\int_0^{\infty} (\omega - \tilde\omega_\psi)^2 |\Psi(\omega)|^2 d\omega
                       }
                       {\int_0^{\infty} |\Psi(\omega)|^2 d\omega},  \tag{2}
$$
where $\tilde\omega_\psi$ is the energy center frequency (see "Center Frequency"). We thus interpret it as frequency-span of ~68% of wavelet's spectral energy. In CWT filter bank terms, a small $\sigma_\omega$ for a given scale means that most of input signal's frequencies will be band-passed about $\tilde\omega_\psi$, and the rest attenuated/rejected.

Center Frequency
The cited answer's approach won't work here - at least not directly. Either way it's suboptimal; since discretization favors us in frequency domain, integrating is best. Ref [1] defines three meaningful frequencies one can associate with a wavelet's scales, and corresponding center frequencies; two of those for center frequency are the energy center frequency and the peak center frequency, latter defined as frequency at which $\Psi(\omega)$ peaks, and former as:
$$
\tilde\omega_\psi = \frac{\int_0^{\infty}\omega |\Psi(\omega)|^2 d\omega}
                         {\int_0^{\infty}|\Psi(\omega)|^2d\omega}  \tag{3}
$$
Here we don't assume $\Psi(\omega)$ is zero-mean (in fact it can't be). If $\Psi(\omega)$ is even-symmetric, then energy center frequency == peak center frequency. For finite and coded $\Psi(\omega)$, this distinction becomes significant and crucial to any wavelet regardless of shape ("Edge Cases").

Are we right?
Three sanity checks:

Can results be interpreted meaningfully?
Is $\sigma_t \sigma_\omega \approx 0.5$?
Does $\sigma_t \sigma_\omega$ scale correctly with varying wavelet parameters ($\mu$, scale, etc)?

#2 is our main check; if the product is much greater or lesser than 0.5, we're off.
To cut it short, I had to 'cheat' to find the proper range for $t$. I started assuming $\omega_\psi$ to be correct since it made sense and correctly gave center frequency as the exact location of the mode=peak. For $\sigma_t$ I ranged $t\in [-.5, .5)$, as that's interpreted as fraction of frame's length. $\sigma_t \sigma_\omega$, however, did not multiply to $0.5$.
Working backwards, $\sigma_t = 0.5 / \sigma_\omega$, I found the correct $t$ to range $t \in [-N/2, N/2)$! And then it all clicked: units.
The units of $\sigma_\omega$ are samples per cycle. Then, to cancel in Heisenberg's relation, $\sigma_t$ must have the reciprocal units: cycles per sample. This is unified by below relation for DFT:
$$
\begin{align} f_p \left[ \frac{\text{cycles}}{\text{second}} \right] & = \left( f_{\text{DFT}} \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right) \end{align}  \tag{4}
$$
which can also be used to obtain resolutions in physical units (Hertz, seconds, etc).
Above is also a bit simplified; the full units are:
$$
\begin{align}
\left[\left[ \sigma_\omega \right]\right] &= \frac{\text{cycles}\cdot \text{radians}}{\text{samples}}  \tag{5a} \\
\left[\left[ \sigma_t \right]\right] &= \frac{\text{samples}}{\text{cycles}\cdot \text{radians}}
\tag{5b} \\
\end{align}
$$
Radians in $(5b)$ are mandated by those in $(5a)$ to cancel in product to yield unitless Heisenberg area. They're also sensibly interpreted once considering the nondimensional variants.

Edge cases
Discretizing problems don't end at finite precision and scaling conversions; consider what happens with scale pushed to extremes.
High scale
Freq-domain not doing too poorly - mainly finite precision imperfections. But time-domain is trouble: it doesn't fit the frame! Two problems here:

How do we handle $\psi (t)$? Should we extend the frame by adding samples until decay?
How does our choice in 1 affect CWT?

#2 is key to #1. Recall in CWT with freq-domain wavelets, we do a form of discretized convolution theorem, getting sort of the exact continuous-time result, then discretizing. An imperfect time-domain wavelet is characterized by not fitting in the frame - but what if the frequency-domain wavelet fits just fine?
I've yet to meet a time-domain wavelet that doesn't decay sufficiently if the frequency-domain wavelet does; problems arise when the frequency-domain wavelet doesn't. At the high scale extreme, the freq-dom wavelet is defined by a single sample, which is a time-domain pure sine with $\sigma_t=\infty$. Proper decay demands even symmetry of the freq-domain wavelet, which is bare minimally attainable by two adjacent bins at equal values. How exactly various scenarios are handled can be found in squeezepy's tests/props_test.py.
Low scale
Now the time-domain wavelet fares fine, while freq-domain is clearly troubled: half (or more) of it is entirely missing! We have the "trimmed bell" problem, contrary to "insufficiently sampled bell" for high scale case. The center frequency distinction becomes critical here; peak will always equal $\pi$, while energy will lie to its left. Investigating the behavior of the $\sigma_\omega$ integral, I found the energy variant to yield stabler and at times more sensible results, but depending on application, peak can also be used (with due re-interpreting of $\sigma_\omega$).
Nondimensional measures
Doing everything right above, one can take the Morlet wavelet and vary its mu only to notice that time and frequency resolutions don't change! What's the deal?
In short, we need a nondimensional measure; this is obtained by normalizing above as follows (thus completing formulae in ref[2]): $\sigma_t \cdot= \omega_\psi,\ \sigma_\omega /= \omega_\psi$ in relevant equations above, where $\omega_\psi$ is peak center frequency. As to how this is interpreted or why it works - it's already a lot for one question, so briefly: try making sense of $\sigma_t$ relative to center frequency's time-domain period, and $\sigma_\omega$ relative to spread around center frequency in log space (or simply ratio-wise).
Non-dimensional measures change with wavelet parameters, but not with scales, so each has its own use and interpretation.
Where's the $2\pi$?
Mallat's equations divide by $2\pi$, we don't, but there's no contradiction: this is per Parseval-Plancherel's theorem relating time-domain and frequency-domain energies. Since we explicitly divide by said energies in $(1)$ and $(2)$, the $2\pi$ is automatically accounted for.
Further reading
Not everything was explained here, nor can be in a single post. ssqueezepy contains a complete implementation, with accompanying tests that provide further explanation - particularly on handling scale extrema, which is a whole topic of its own.
References

Higher-Order Properties of Analytic Wavelets - J. M. Lilly, S. C. Olhede
Wavelet Tour, Ch1-4 - S. Mallat


FFT, padding, IFFT and plot in time domain I performed an FFT with some zero padding and would like to transform it back to time domain. When I plot the signal it looks wrong. Is it because of the padding used in the FFT? How can I make it right? here is my code for doing the IFFT:

NFFT=4301;
AvIFFT = ifft(FFTavg.dat,NFFT);
plot(real(AvIFFT))

And I get this plot:

Any help would be greatly appreciated! Thank you!

ps: I have asked this question on matworks but no success as yet.

EDIT:

As suggested in the comment and answer below, I tried fftshift:

NFFT = 4301;
FFTshifted = fftshift(FFTavg.dat);
AvIFFT  = ifftshift(FFTshifted,NFFT);
plot(real(AvIFFT))

Here is the plot:

user3406207: I'm not sure what you're trying to do because I don't know the meaning of your words "wrong" and "right."  But we can say, if you perform time-domain zero padding on an original unpadded $x1[n]$ input sequence you'll produce some longer-length $x2[n]$ sequence. The $X2[m]$ fft of your $x2[n]$ sequence will be an interpolated version of (more freq-domain samples than) the $X1[m]$ fft of your $x1[n]$ sequence.  And, of course, the inverse fft of $X2[m]$ will be $x2[n]$ and the inverse FFT of $X1[m]$ will be $x1[n]$. I hope that helps. (By the way, you were smart to plot just the real part of your inverse fft.)


What is sensor drift? From Wikipedia:

  If the output signal slowly changes independent of the measured property, this is defined as drift

So my question is, if I have a gyroscope with sensor drift does that mean that $E\left\{f(x)\right\} \neq 0$ ? I have an excercise where it is written, that the sensor has drift but also that the $E\left\{f(x)\right\} = 0$. But to me one statement seems excludatory of the other one.

I assume you set the first momentum of your true measure to $\mu=0$, so your  question actually reads

  gyroscope with sensor drift does that mean that $E[f(x)] \ne \mu$?

No. It might still be the case that the sensor value, in the long term, drifts centred around the true value.

The thing is that $f(x)$ is a random variable with non-zero variance and a time dependency. That means: If there's drift, that random variable is not a stationary process.


Step response of third-order continuous-time transfer function I have a transfer function of the form:

$$H(s) = \frac{b\omega_n s^2 + a\omega_n^2 s + \omega_n^3}{s^3 +     b\omega_n s^2 + a\omega_n^2 s + \omega_n^3}$$

If it matters, $a=b=2$.  Is anyone aware of a simple explicit form of the step response for this transfer function (i.e. something directly in terms of $\omega_n$, $a$, and $b$)?

It does matter that $a=b=2$, because this gives a relatively nice looking solution. As mentioned in a comment, you should split up the given transfer function $H(s)$. But first, let's introduce a normalized variable $p$:

$$p=\frac{s}{\omega_n}\tag{1}$$

Now we can write the given transfer function as

$$\hat{H}(p)=\frac{2p^2+2p+1}{p^3+2p^2+2p+1}\tag{2}$$

The denominator can be factored as

$$p^3+2p^2+2p+1=(p+1)(p^2+p+1)\tag{3}$$

Partial fraction expansion of $\hat{H}(p)$ gives

$$\begin{align}\hat{H}(p)&=\frac{1}{p+1}+\frac{p}{p^2+p+1}\\&=\frac{1}{p+1}+\frac{p}{(p+\frac12)^2+\frac34}\tag{4}\end{align}$$

The Laplace transform of the step response equals the transfer function divided by $p$:

$$\begin{align}\hat{A}(p)=\frac{\hat{H}}{p}&=\frac{1}{p(p+1)}+\frac{1}{(p+\frac12)^2+\frac34}\\&=\frac{1}{p}-\frac{1}{p+1}+\frac{1}{(p+\frac12)^2+\frac34}\tag{5}\end{align}$$

The inverse Laplace transform of $(5)$ is (table)

$$\begin{align}\hat{a}(t)&=u(t)-e^{-t}u(t)+\frac{2}{\sqrt{3}}e^{-t/2}\sin(\sqrt{3}t/2)u(t)\\&=\left[1-e^{-t}+\frac{2}{\sqrt{3}}e^{-t/2}\sin(\sqrt{3}t/2)\right]u(t)\tag{6}\end{align}$$

Finally, denormalization according to $(1)$ taking into account that $\hat{A}(p)=\omega_nA(s)$ (because of division by $p$) gives the desired step response:

$$\begin{align}a(t)&=\hat{a}(\omega_nt)\\&= \left[1-e^{-\omega_nt}+\frac{2}{\sqrt{3}}e^{-\omega_nt/2}\sin(\sqrt{3}\omega_nt/2)\right]u(t)\tag{7}\end{align}$$


QAM Amplitude Levels I'm working on a project to stream video between 2fpgas by using QAM. While mapping symbols with different amplitude pulses(PAM) a question popped up in my mind. In real life how does communicating wireless devices know which amplitude levels to agree on for digitally modulating-demodulating bitstreams? Is there a standard for that?

n real life how does communicating wireless devices know which amplitude levels to agree on for digitally modulating-demodulating bitstreams?

You can't "agree" on levels, since the received amplitude depends on the attenuation on the channel, and that depends on distance and a lot of other factors.
So, your receiver depends on detecting something that it knows the right amplitude (or power) of – for example, a preamble, or an average squared symbol, or something like that.
How that is done depends on the system you look at – generally, this is part of the synchronization and equalization efforts every system has to make. And it really differs! For example, some continuously transmitting system will not insert any specific symbols for that, but depend on the receiver doing statistics. Others have a "silence" period followed by known symbols; other systems are packetized and depend on preambles with specific autocorrelation properties (but not necessarily 100% defined content). Other systems depend on a pilot tone, and there's many more options.


Power Spectrum of the Current of Compact Fluorescent Lamp (CFL) Compared to the AC Input Is it possible for an AC current (current passing from CFL) to have a harmonic component greater in magnitude than its fundamental?

In my country fundamental frequency which carries power is 50 Hz. But, when I pass this power to CFL(compact fluorescent lamp) and get the FFT of the current passing through it (through MATLAB). The spectrum shows me fundamental component(50 Hz) magnitude less than the third harmonic(150 Hz) roughly twice of the fundamental.

Hope this will clarify my question.

The simple explanation is that the Compact Fluorescent Lamp (CFL) acts as Non Linear System and hence generates signals with different frequency than the input.

Another thought would be that the sampling channel as a whole has different response to those frequencies and it damps 50 [Hz] / amplify 150 [Hz] (Or both).

You take into consideration the whole chain until you get the digital signal.


Single-sided Z transform with difference equations and the system function I am working on this problem:
Given an impulse response,

find the system function,
find the difference equation representation,
find pole-zero plot,
find output $y[n]$ if the input is $x[n] = 0.25^n * u[n]$

Here is what I have so far (MATLAB code)

Find H(z)
$$H(z) = \frac{5}{1 - 0.25z^{-1}}\quad,\qquad|z| > 0.25$$

For part 2 I used the fact that
$$H(z)=\frac{Y(z)}{X(z)}$$
I cross multiplied, then took the inverse $\mathcal Z$ transform and got:
$$y[n] - 0.25y[n] = 5x[n]$$

For part 3 I used the zplane function:

 b = [5 0];
 a = [1 -0.25];
 figure; zplane(b,a);

Part 4 is where I am getting confused. I took the one-sided $\mathcal Z$-transform of the difference equation to get
$$Y^+(z)-0.25[y[-1]+z^{-1}Y^+(z)] = 5X^+(z)$$
Assuming the system is casual (due to unit-step function in impulse response and input $x[n]$, also no initial conditions given in the problem),
$$y[-1] = 0$$
So, after rearranging, $$Y^+(z)=\frac{5}{1-0.5z^{-1}+0.0625z^{-2}}\quad, \qquad |z| > 0.25$$
Using the residuez function for partial fraction decomp

 b1 = [5 0 0];
 a1 = [1 -0.5 0.0625];
 [R p C] = residuez(b1,a1)

 R =

      0
      5

 p =

     0.2500
     0.2500

 C =

      0

Then taking the inverse $\mathcal Z$-transform of the terms generated by residuez,
$$y[n] = 5(1/4)^nu[n]$$
However, when I check, this is not the case. For the check, I first create a MATLAB generated $y[n]$ sequence using filter

%MATLAB check for part 4
n = 0:50;  % compare first 50 samples
x = (1/4).^n;
y = filter(b,a,x);
Then, create my check sequence calculated from the inverse $\mathcal Z$-transform and the residuez function

ycheck =  5*(0.25).^n;
error = max(abs(y-ycheck))
The sequences are not the same, and error nowhere close to nominal. Can anyone tell me where I went wrong please?

You misinterpreted the result of residuez. Note that there is also a fourth output argument, namely the multiplicity of the poles (which is $2$ in your example). In your case, the residuez function doesn't help much, because you know the poles already, and the result is already in a partial fractions form:
$$Y(z)=X(z)H(z)=\frac{5}{\left(1-\frac14 z^{-1}\right)^2}\tag{1}$$
You can either look up the inverse $\mathcal{Z}$-transform of $(1)$ in a table, or, maybe even more straightforward, compute the result by convolution in the time domain:
$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\tag{2}$$
Solving $(2)$ for the given sequences is very straightforward.


Using FFT to Unit Test a Music Player I'd like to validate changes to our app music player in our Continuous Integration unit tests (more like integration tests, really). My thoughts were to start with an audio file that contains only a pure tone, or a series of pure tones. Our automated test would stream this file over the internet and play it through our player. We would capture the output using MTAudioProcessingTap. This output would go through and FFT and the pass/fail test would be check for the existence of only the original frequencies.

Would this work? I'm thinking this would catch pops and glitches, but not necessarily discontinuities in the output.

The phase vocoder algorithm used with overlapped FFT frames might easily catch certain kinds of discontinuities in pure sine waves.

If the audio file is compressed, note that this process itself may add certain kinds of distortion.  So you would need to find an experimental threshold.


Nyquist (Nth) digital filters I am working in something were I should use a upsampling filter. I have decided to use a Nyquist filter(Lth filter). I know that there are two constraints. The first The frequency vector values must mirror each other in pairs around $\pi/2$. The second is the amplitude vector values must mirror each other in pairs around a magnitude of 0.50 What I am looking for are references to actually design these types of filters. I cant find any references that show how to implement them.

So dose anyone have a decent reference on hoe to design these filters? Does anyone have a reference on implementing these in Verilog?

This Mathworks documentation gives a good overview about the different parameters: http://uk.mathworks.com/help/dsp/examples/fir-nyquist-l-th-band-filter-design.html?requestedDomain=uk.mathworks.com

This tutorial goes into more detail: http://www.analog.com/media/en/training-seminars/tutorials/MT-002.pdf

And this one may be useful to clarify any doubt in the concepts presented before: http://www.lumerink.com/courses/ECE697A/docs/Matlab%20Filter%20Design%20and%20Implementation.pdf

hope this helps


Non Local Means Filtering Technique Compared to Gaussian Blur and Total Variation Denoising How does the Non Local Means Filter different from other approaches such as:

Gaussian Smoothing.
Total Variation Denoising.

Gaussian smoothing and TV denoising are local filters that only consider the neighborhood around each target pixel; on the other hand, non-local means updates a pixel's value by using a weighted average of the pixels judged to be most similar.

The reason why it is called non-local is because it does not consider the spatial relationship between pixels, instead, the similarity between different patches and the patch around the target pixel is computed to determine the weight mentioned above (inversely proportional to the distance between its patch and the reference patch of the target pixel).

Yet some studies also showed that Non-Local means is a local image denoising algorithm, in which experimental results indicated that the bias of the non local means estimator is an increasing function of the radius of the similarity searching zone (patch).

The paper (A.Buade et al's On image denoising methods) provided more detailed comparisons between non local means and Gaussian, TV from the aspect of (1) the method noise, (2) the mean square error, and (3) the visual quality of the restored images by simulation. You can skip the theoretical part by reading section 6 for a brief conception of the effect of those methods, which may also be helpful in better understanding the meaning of those mathematical formulas in previous sections.


Question about symbol detection in ISI channels (use of APP decoder in Matlab) I am studying the problem of digital symbols detection in ISI channel. In the figure below, the output symbol is represented as a weighted sum of input symbols.

I am trying to use MAP decoder comm.APPDecoder to estimate the transmit symbols. To construct the trellis of such figure, I am using the Matlab function poly2trellis. The function can be used as follows:

trellis = poly2trellis(ConstraintLength,CodeGenerator)

My question is how can I enter (or consider) the real-valued coefficients $\Gamma_{-L_2}, ..., \Gamma_{L_1}$ in the second input argument CodeGenerator. If these coefficients are binary, then it is straightforward to convert this binary sequence to octal. But, I have no clue with real-valued coefficients. Any hints?

As you mentioned, applying linear equalization (e.g. filtering with the inverse of the channel) will not yield the optimal performance. Though, what you are searching is not the comm.APPDecoder (which is solely for convolutional codes, i.e. bits with binary $\Gamma$), but you are looking for maximum likelihood sequence estimation (MLSE).

Here, MATLAB offers the comm.MLSEEqualizer object, which you can e.g. use like this:

Constellation = [1+1j, 1-1j, -1+1j, -1-1j];
Channel = [1 0.7, 0.3, 0.2]';
sigma2 = 0.3;

tx = Constellation(randi(4, [100,1])).';

noise = sqrt(sigma2/2)*complex(randn(size(tx)), randn(size(tx)));
rx = filter(Channel, 1, tx) + noise;

mlse = comm.MLSEEqualizer('Channel',Channel,'Constellation',Constellation);

estimated = step(mlse, rx);

[estimated tx abs(estimated-tx)]

program output (only a few rows of it):

   1.0000 - 1.0000i   1.0000 - 1.0000i        0
  -1.0000 + 1.0000i  -1.0000 + 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0
  -1.0000 - 1.0000i  -1.0000 - 1.0000i        0
  -1.0000 - 1.0000i  -1.0000 - 1.0000i        0
   1.0000 + 1.0000i   1.0000 + 1.0000i        0

as you see, the MLSE equalizer completely recovers the transmitted sequence.


How to extend line until end of image Lets say, i have a line and i want to extend the line until the end of the image.
 The line is displayed in the first photo with the green is my basic line, i want to extend it like the second image, until the end of the image.
 I want to do this for all the lines that i have detected.
 Does anyone knows how to do this?

i think i found it

X= [x1, x2]
Y=[y1, y2]

c = [[1; 1]  X(:)]\Y(:);
slope_m = c(2)
intercept_b = c(1)

aLine = [slope_m,-1,intercept_b];
points = lineToBorderPoints(aLine,size(negfillHoles))
line(points([1 ,3]),points([2 ,4]),'LineWidth',1.5,'Color','r');
 plot(x1,y1,'*','LineWidth',2,'Color','g');


How to get a generic 3d Face model? I'm working on image processing that involves facial features.

Where can I get a generic 3D model of a human face? (any format is fine as long as I can manipulate the data using C or Matlab)

You might have to make your own from a 3d face database, which may not be easy to do. This website has a good listing of public databases. http://www.face-rec.org/databases/


How do I remove motion blur? Is there a generically accepted method of removing motion blur from an image. For a simple case, we can assume that motion happens in a straight line. I assume that it's a two part process composed of motion estimation and then deconvolution, but how is it actually done?

Yes, deconvolution.  This page describes a number of deconvolution methods and methods for estimating the point spread function:

Removing Motion Blur from Astrophotographic Images

They say the deconvolution literature is "extremely extensive".   They choose Lucy-Richardson algorithm for deconvolution and develop their own motion estimation algorithm for determining the point spread function.


What's the shortest possible code that correct all the single bit errors? Suppose the original code is a 2-bit binary coding with 3 possible alphabet.
\begin{align}
a&\rightarrow 00\\
b&\rightarrow 10\\
c&\rightarrow 11
\end{align}
Is that $\operatorname{Hamming}(7,4)$ coding? Isn't 7-bit too long?

In order to correct 1-bit errors, we need to have a Hamming distance (the number of bits that must be flipped) of 3 between all the symbols: when we flip one bit at random, then flipping another bit by random must either get us back to the original symbol (if the same bit was flipped twice) or an invalid symbol.
With only two symbols we can encode as 000 and 111 because this gives us the required Hamming distance, so clearly we can encode 3 symbols with 6 bits. Can we can do better though?
If we use 4 bits to encode one symbol as 0000 then 3 flips result in 0111, 1011, 1101 or 1110. These all differ from each other by no more than 2 bits, so clearly 4 bits is not enough.
On the other hand, if we use 5 bits to encode one symbol as 00000 then 3 flips give us 00111, 01011, 01101, 01110, 11100, 11010 and 11001. If we then encode the 2nd symbol as 00111, we're left with the following valid choices for the 3rd symbol: 11100, 11010 or 11001. These are all distance 3 from the first symbol and distance 4 from the second symbol.
So one of the many possible encodings for 3 symbols using 5 bits is:
a = 00000,
b = 00111,
c = 11100.


Non Circulant Translation Using Fourier Transform The translation property of Fourier Transform (FT) for a two dimensional image $f$ is as
$$
f(x-x_0, y - y_0) = F(u, v)e^{-j2\pi(ux_0/M+vy_0/N)}
$$

Using this equation, the following code (in Matlab) computes and shows translation of the input image.

im = im2double(imread('cameraman.tif'));
F = fft2(im);
[r, c] = size(F);
x0 = 20;
y0 = 40;

% [u, v] = meshgrid(1:c, 1:r);
[u, v] = meshgrid(-c/2:c/2-1, -r/2:r/2-1);
A = fftshift(exp(-1i * 2 * pi * (u * x0 / c + v * y0 / r )));

Y = F .* A ;
y = ifft2(Y);
imshow(real(y), [])

The output of above code is as follows

but I need output image be as following

My question:
How I can produce The second image using FT?
Thanks.

You can obtain an apparently linear translation using DFT by

first zero padding the original image into a large enough size according to the translation you want.
then applying DFT for unavoidable circulant translation
and finally cropping the original size of the output

Very inefficient using FT but if it's what you want...


Algorithm implemented in the timeseries\resample function of Matlab Does anyone know what algorithm is implemented in the Matlab function resample for timeseries?

I had a look on the manual but it just gives the instruction to use it: http://it.mathworks.com/help/matlab/ref/resampletscollection.html

Since the syntax is different from the one used for the usual resample (http://it.mathworks.com/help/ident/ref/resample.html) I assumed is not the same algorithm. Am I right?

See http://www.mathworks.com/help/matlab/ref/timeseries.setinterpmethod.html
Seems to a choice between a few fairly simple interpolation schemes. The standard resample() uses a polyphase filter interpolation which is much better for most applications.


How to determine fundamental frequencies (beats/minute) of heartbeat? (matlab) read = matfile('ECG.mat');
[cor, lags] = xcorr(read.ecg, read.ecg);
read.Fs;

So i have read.ecg (signal) and read.Fs. How can i find/count beats/minute in matlab?

Here is plot

I would suggest processing the ECG signal through a standard QRS detection algorithm (e.g. Pan-Tompkins). Once you have isolated the QRS complexes in your ECG signal, you can simply count them, divide by the duration of your ECG signal, and convert to beats/min. Here is one implementation of Pan-Tompkins: https://www.mathworks.com/matlabcentral/fileexchange/45840-complete-pan-tompkins-implementation-ecg-qrs-detector. I am sure you can find others on physionet.


Most generally-accepted audio term, for "pre" and "post" lengths of a symmetric filter's impulse response? Consider a FIR audio filter with symmetric impulse response.
In my case, the filter kernel width is 15 frames on each side (so including the center frame, 31 frames total).  When examining the audio impulse response of this filter, it includes (A) the initial interval where it approaches the center, then (B) the centered impulse output frame itself, then (C) the subsequent interval following the center.
I have seen numerous terms referring to intervals (A) and (C), including

(A) lookahead, ramp-in, processing delay, pre-echo, positive filter width, ring-in, pre-response;
(C) decay time, ramp-out, post-delay, negative filter width, ring-out, causal response.

... and many others.
In the audio world, is there a term that is considered most accurate and/or most clearly descriptive for these two intervals?

That's the causal (comes after the impulse) and anticausal part.


What is the exact meaning of the output of the Discrete Fourier Transform I'm fairly new to the subject, but so far my understanding that this would be a transform you could use to go from a discrete set of data, say [1, 0, 1, 2] to a continuous sinusoidal function in the frequency domain, if I have fundamentally got this wrong maybe this is where my confusion comes from but once you have applied the DFT to this signal, and then find its magnitude. What is the meaning of the outputs, are they coefficients of a sinusoidal function? Or is the just the representation of that data as an amplitude in accordance to its frequency?

The Discrete Fourier Transform (DFT) can be understood in different contexts.  What the meaning of the values are depend on the context.  Ideally you should understand them all.

The DFT is usually taught as a Dirac delta sampled case of the continuous FT with a rectangular window.  This can be confusing if you don't know Calculus, or haven't had a little bit of Real Analysis.  In this case, the values of the DFT can be considered coefficients of complex exponential functions that will reconstitute your original signal.

Another approach is to understand the DFT as a matrix multiplication in a Linear Algebra sense.  In this case the DFT values can be considered coordinates in a N dimensional vector space.

My favorite, in terms of a concrete understanding, particularly for real valued signals, is the DFT bin value is a weighted average of a set of points (Roots of Unity) on the unit circle on the complex plane.  You can also think of it as a center of mass.  You can find a write up of this interpretation in my blog artice "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity".

If you are not familiar with complex numbers, the Roots of Unity, or Euler's Equation, I recommend that you read my first blog as well.  Understanding these concepts is key to understanding the DFT in all contexts.

Hope this helps,

Ced


Why normalize the data set before applying Direct Linear Transform Direct Linear Transform (DLT for short) is a method of homography estimation, it solves the overdetermined linear system via SVD $$Ah=b$$to find a solution $h$ under constraint $\|h\|=1$. Actually it finds the least square solution which minimize $\|Ah - b\|$.

I understand the basic idea of this algorithm, but it is recommended to normalize the data set before applying DLT on it, and here is a intro about how to do the normalization. It is lectured that data normalization is important to DLT, without normalization the results from DLT is not stable.

I wonder why? Just because DLT involves solving the linear system using SVD and $A$ might be singular?

The normalization is basically a preconditioning to decrease condition number of the matrix $A$ (the larger the condition number, the nearer the matrix is to the singular matrix).

The normalizing transform is also represented by a matrix in the case of homography estimation, and this happens to be usable as a good preconditioner matrix. The reason why is that is more elaborate and is explained briefly in H&Z book (4.4.4, p. 107: Why is normalization essential?) or in more detail in the paper "In Defense of the Eight-point Algorithm".

Put it simply, the matrix $A$ consists of products of image coordinates which can have different scale. If the scale differs by factor of $10$, the products differ by a factor of $10^2$.

The source and target coordinate data are usually noisy. Without normalization, the data from source would can have two orders of magnitude larger variance than from target (or vice versa).

The homography estimation usually finds parameters in a least-squares sense - hence the best statistical estimate is found only if variances of the parameters are the same (or known beforehand, but it is more practical just to normalize the input).

Direct solvers do not like poorly scaled problems because numerical instabilities appear (e.g. dividing very large number by a very small number easily leads to numerical overflow).

Iterative solvers struggle with badly conditioned matrices by needing more iterations.

So normalization is essential not only for numerical stability, but also for more accurate estimation in presence of noise and faster solution (in case of iterative solver).


What is the standard way to use z-score with EEG? I can't understand how to use normalization technique for EEG.

Is it used for results of FFT or for raw EEG signal?
Are there different methods with common name "$z$-transform"?
And what z-score used in this article?

Need a bit of clarity in question. There can be numerous normalization techniques for EEG depending on how they are used. That can be before processing (i.e. raw) or even after processing (for example, with FFT or other transforms).

Z-transform is ONE method. Period. There is nothing else called z-transform. You can have a look on Wikipedia entry of z-transform here.

Z-score, as used in your article, is a statistical measure, also called the standard score. More info here, again.

To sum up, z-transform and z-score are NOT linked.


UPDATE : How to continue Computing the Convolution $$
x(t):=\begin{cases}
1&\text{if $0<t<T$}\\ \\
0&\text{if otherwise}
\end{cases}
\qquad\text{and}\qquad
h(t):=\begin{cases}
t&\text{if $0<t<2T$}\\ \\
0&\text{if otherwise}
\end{cases}
$$
Compute $(x*h)(t)$.
My approach:
We first note that :$$\displaystyle x(\tau)=\begin{cases}
1&\text{if $0<\tau<T$}\\ \\
0&\text{if otherwise}
\end{cases}\qquad\text{and}\qquad\displaystyle h(t-\tau)=\begin{cases}
t-\tau&\text{if $0<t-\tau<2T$}\\ \\
0&\text{if otherwise}
\end{cases}$$ We observe that :
$$
0<\tau<T
$$
$$
t-2T<\tau<t
$$
Thus, $(x*h)(t)\in\mathcal{R}[\max(0,t-2T),\min(t,T)]$ (i.e. Riemann integrable). Therefore :
\begin{align*}
(x*h)(t)&:=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)\;\text{d}\tau\\
&=\int_{\max(0,t-2T)}^{\min(t,T)}t-\tau\;\text{d}\tau\\
&=\left. t\tau-\frac{\tau^{2}}{2}\right|_{\max(0,t-2T)}^{\min(t,T)}\\
&=\left[t\min(t,T)-\frac{\min^{2}(t,T)}{2}\right]-\left[t\max(0,t-2T)-\frac{\max^{2}(0,t-2T)}{2}\right]\\
&=t[\min(t,T)-\max(0,t-2T)]-\frac{1}{2}[\text{min}^{2}(t,T)-\text{max}^{2}(0,t-2T)]
\end{align*}

Problem: I want to further solve this to obtain the final solution as a piecewise function (if it has such form)

I present to you the MATLAB simulation of the convolution

First one corresponds to $T=1$ and the second one corresponds to $T=2$. Here is the algorithm used:
t = linspace(0,10,1000);   % assume a time span for "t"
T=1;
m = min(t,T);              % find min of (t-1 , 5) for each "t"
M = max(0,t-2*T);            % find max of (t-3 , 3) for each "t"
vind = M < m;              % find "valid" range of "t"

y = t.*(m-M)-0.5*(m.^2-M.^2);% EVALUATE the CONVOLUTION
y = y.*vind;               % force non-valid range to zero.

figure,plot(t,y);         % DISPLAY:
title('convolution is :');

Credits to $\text{@Fat32}$ for the code in his answer in the link provided in the comment

Some general hints for this type of problem

If the function is defined piecewise than chances are your calculation and solution will also have to be done piecewise in sections.
Determine first where the result is zero, i.e. where $x(\tau)$ and $h(t-\tau)$ don't overlap. This determines the start and end of the sections that you need to consider. Hint: if both functions are causal than $t=0$ is always a boundary. The length of a convolution is always equal (or smaller) than the sum of the lengths of the signals.
Then find the section boundaries: these are where the mins() and maxes() changes from one argument to the other.
In each section replace the max()/min() with the correct value for this section. Solve  each section individually and splice the solutions together.

It's also super helpful if you can develop some graphical intuition for this type of problem: draw both functions with a few different overlaps and see what happens when they start overlapping, when they are fully overlap and when they stop overlapping.
You already came a long way. Give that a shot and post back if you are still stuck. Happy to help more.
EDIT
More help on the piecewise solution. For simplicity let's call the convolution $y(t)$

$y(t)$ must be zero for $t < 0$. The convolution cannot start earlier than any of the signals
$y(t)$ must be zero for $t > 3T$. The convolution cannot be longer than the sum of the length of the signal.
That means you only need to focus on the interval $0 < t < 3t$
$min(t,T)$ is  $t$ for $ 0 < t < T$ and $T$ for $T < t < 3T$
$max(0,t-2T)$ is  $0$ for $ 0 < t < 2T$ and $t-2T$ for $2T < t < 3T$

That means that you need to split the solution into three section: $[0,T]$, $[T,2T]$ and $[2T, 3T]$. Within each section, you can replace the min/max functions by the actual value since it doesn't change over the section. Example: in the first section $[0,T]$, we have $max(0,t-2T) = 0$ and $min(t,T)=t$ so you would be integrating over $[0,t]$.


Avoiding clicks with changing biquad coefficients I have a DSP application that filters a signal with five cascaded biquads whose coefficients are determined by the Audio EQ Cookbook formulae (Thanks for this work, Mr. Bristow-Johnson btw). Each filter's coefficients can be changed individually via RS232 from a matlab app, however doing so produces audible click noises.

To prevent clicking, my first idea was to lowpass filter the individual coefficients to smoothen the transition. A friend of mine told me that this could cause instability. He is way more experienced in the field than I am, also I still lack the proficiency to prove/-disprove his claim or finding out which circumstances can cause instability in this case, so I take his word for granted.

However, his solution (Smoothing parameters instead and then calculating the correct coefficients) doesn't work for me. Calculating intermediate values in Matlab and transmitting them is not an option because of bandwidth constraints, calculating them on the DSP is not an option because of limited CPU time.

One solution would be to employ two parallel filters and crossfade between them. But I don't particularly like it, because it doubles my computational load, also two changes in rapid succession would lead to audible clicks anyways.

The next idea I had was to employ Gold-Rader structure (see below) for the biquad denominator. Since this structure uses the pole radius, which is always kept below 1 (as long as I filter it with a filter that doesn't overshoot...), it can't become unstable... except for the fact that the radius is multiplied with $sin(\phi)$ and $cos(\phi)$ before introduced in the signal flow graph. So if I'm transmitting the calculated products $r \cdot sin(\phi)$ and $r \cdot cos(\phi)$, I'm again not sure about how filtering those coefficients would affect stability.

Each of the five biquads must be tunable across 20 Hz to 20 kHz, -18 to +18dB gain (where applicable) and a Q from 0.5 to 10 for the Peaks. The type of filter stays the same though (i.e. a highpass will always stay a bypass, a peak filter will always be a peak filter...).

So my question is: what is the most efficient way to avoid clicks while not risking instability of the filter?

Another possibility is to send the coefficients as poles and zeros and interpolate/smooth in the z-plane instead of the coefficient domain. This is always guaranteed to be stable. You can either interpolate in real/imaginary presentation or, if you want to get fancy, in magnitude and phase. The first one is a lot easier, since the phase interpolation has to handle wrap around.

Coefficient calculation from real and imaginary part is straight forward and suitable for real time. Only exception may be if you have a dual real pole/zero instead of a complex conjugated pair.


infomation about Stockwell Transform I need some information about Stockwell transform (also known as the S-transform):

How can I implement it in MATLAB?
*Does it give you the damping ratio $\zeta$ of a signal like the complex Morlet continuous wavelet? $$x(t)= e^{-\omega_{n}\zeta}\sin(\omega_{d}t+\theta )$$
What are its advantages compared to CWT?

The S-transform defined as follows:
$$s( \tau ,f)= \int_{-\infty}^\infty h(t) \frac{\ |f|}{\ \sqrt{2 \pi } } e^{ \frac{-(t- \tau )^{2} f^{2}}{2}  }    e^{-j 2 \pi f} dt$$ (This is a good tutorial about the S-Transform)

These are the first few question that I can think about at the moment. It seems no one asked about the S-transform before, so put your questions in this post by editing it, and answer them if you know them. Someone with enough reputation put a tag for S-transform please. Thanks a lot.

The Matlab code of the S-transform can be found here for example: http://en.verysource.com/code/1180181_1/st.m.html

The main advantage compared to CWT is its simplicity. In fact, it can be seen as equivalent to a Morlet wavelet, see for example paper by Ventosa, Schimmel, Simon, Dañobeitia and Mànuel "The S-transform from a wavelet point of view", IEEE Trans on Signal Processing, 2008.


Is there any processing needed to playback audio? I am new to audio processing. I wrote a simple real-time playback program (when you speak on the mic your voice is immediately played on the output), however it is not working. I am using a 8Hz ADC with 128 frame size, and I output these 128 samples directly to audio PWM. It doesn't play the audio even though there's a "reaction" when you tap the mic.

Is there a processing needed before outputing the samples, that would explain why it isn't working ?

Here is my main loop :

while (1) {

while (ADCChannelIsBusy(MyADCHandle)){}
ADCChannelRead  (MyADCHandle, Frame, ADC_BUFFER_SIZE);
while(OCPWMIsBusy(MyDACHandle));
OCPWMWrite  (MyDACHandle,Frame, OCPWM_FRAME_SIZE);

}

I wrote none of these functions, all were found in the dsPic starter Kit. Frame is the array where samples are stored.

This is unlikely to be a signal processing question. This is about the internals of the dsPic's firmware libraries and would be more suitable on another stackexchange site.

Things that could go wrong:

ADC input and PWM output might have a different data format (for example 12-bit right justified for the ADC, 8-bit for the PWM). You might need to scale / offset the data.
ADC_BUFFER_SIZE and OCPWM_FRAME_SIZE might be of different sizes, requiring you to use an intermediate buffer and counters to keep track of data pouring in/out of this buffer. Or, if the implementation allows it, you might instead read/write partial buffers of a smaller size - identical for both calls.


Why is my signal attenuated after applying a HighPass-filter? I am new to DSP, but I've been attempting apply a butterworth high-pass filter on a 1D array (Values peeking at 300~) to eliminate any noise, yet when applying a high pass filter, I'm ending up with an attenuated result?

So far I have attempted to approach the problem with different specs, but it seems that my lack of experience in Signal Processing is holding me back.

def butter_highpass(highcut, fs, order):
    nyq = 0.5 * fs
    high = highcut / nyq
    b, a = butter(order, high, btype="highpass")
    return b, a

fs      = 30.0            #Sampling rate Hz
b, a = butter_highpass(5, fs, 5)
filtered      = signal.lfilter(b, a, data)

Result:

My Goal : Reduce all noise to 0 while retaining the peak values.

Any suggestions or recommendations are more than welcome

The correct answer to this question has been provided by: Stanley Pawlukiewicz but has since been removed.

  Your filter attenuates the high frequencies of the total timeseries,
  signal and noise.

  Your results are expected.

  You really can’t perfectly remove just the noise.

  When the noise and signal overlap in time and frequency, you have to
  make a tradeoff between attenuating signal and attenuating noise, by
  adjusting the cutoff frequency.


How to perform this spectral decomposition in MATLAB? Given a filter $X(z)$ I want to find $G(z)$ such that it is stable, causal and minimum-phase, and it accomplishes that $$X(z)=K_0G(z)G^*(1/z^*)$$

where $K_0\in\mathbb{R}$. Of course, $G^*(1/z^*)$ would then be anticausal and maximum-phase. The filter $x[n]$ has complex coefficients.

I'm trying to achieve this in MATLAB but I'm finding it more difficult than I thought.

% This is my filter x[n]. You can try with any coefficients, it doesn't matter
x = dfilt.dffir(q_k + 1/(10^(SNR_MFB/10)));
% Here I find its zeros
zeros_x = zpk(x);
% And now I identify those who are inside and outside the unit circle
zeros_min = zeros_x(abs(zeros_x) < 1);
zeros_max = zeros_x(abs(zeros_x) > 1);
% Now I wanted to build the filter g[n] using this information but don't know how

I tried using spectralfact() but it doesn't support inputs with complex coefficients.

Does anyone have any idea on how to get the coefficients of $g[n]$?

The factorization

$$X(z)=G(z)G^*\left(\frac{1}{z^*}\right)\tag{2}$$

with minimum-phase $G(z)$ is only possible if $X(z)$ is a non-negative real function, i.e., $X(z)$ is non-negative and real on the unit circle $|z|=1$.

The following Octave/Matlab script is a simple example of how to find $G(z)$ from a given $X(z)$ (note that $x[n]$ is complex-valued):

N = 10;         % signal length is 2*N-1
x = randn(N,1) + 1i*randn(N,1);
x = conv(x,conj(x(end:-1:1)));
% x has conjugate symmetry and its spectrum is non-negative

r = roots(x);
I = find( abs(r) < 1 );
ri = r(I);
g = poly(ri);
K = real(x(N))/real(g(:)'*g(:));

x2 = K * conv(g,conj(g(end:-1:1)));    % same as x up to numerical errors


How to write time waveform of transmitted sequence with certain chip rate? I have been studying a book which describes how to write the waveform of  sequence that has certain chip rate $T_c$.

The sequence with chip rate $T_c$ is for example

$$X=[-1 -1 +1 +1 +1 +1 +1 +1 +1 -1 +1 -1 -1 +1 +1 -1 -1 -1 +1 +1 -1 -1 -1 -1 +1 -1 +1 -1 +1 -1 -1 +1
+1 +1 -1 -1 -1 -1 -1 -1 -1 +1 -1 +1 +1 -1 -1 +1 -1 -1 +1 +1 -1 -1 -1 -1 +1 -1 +1 -1 +1 -1 -1 +1
+1 +1 -1 -1 -1 -1 -1 -1 -1 +1 -1 +1 +1 -1 -1 +1 +1 +1 -1 -1 +1 +1 +1 +1 -1 +1 -1 +1 -1 +1 +1 -1
+1 +1 -1 -1 -1 -1 -1 -1 -1 +1 -1 +1 +1 -1 -1 +1 -1 -1 +1 +1 -1 -1 -1 -1 +1 -1 +1 -1 +1 -1 -1 +1]$$

Lets assume that $X$ is of length 128.

So the book explains that the waveform can be written as

$$r(nT_c)= X(n ) \underbrace{\exp( (j\pi \frac{n}{2})}_{???}\hspace{1cm} n=0,1,\cdots,127$$

The chip rate $Tc=0.57 ns$...

I don't understand is the part I have ??? under...

Any ideas.

Best Regards.

While a reference to where you got this from would be helpful, I have a theory on where the term in question comes from. Note its form for increasing $n$:

$$
\exp\left({j\frac{\pi n}{2}}\right) = [1, j, -1, -j, 1, j, -1, -j, \ldots ], \ n = [0, 1, 2, 3, 4, 5, 6, 7, \ldots]
$$

This is a discrete-time complex sinusoid whose frequency is one-fourth of the sample rate. In digital communications signal processing, this is a very common intermediate frequency selection when using real-valued sampling. It is often chosen because the center of the desired passband lies in the center of the single independent sideband that you get when using real samples. Shifting the resulting signal to complex baseband is often a subsequent processing step.

In your case, the equation that you gave may express that you have a chip sequence that is modulated onto a complex-valued carrier at one-fourth of the sample rate. Since you don't indicate what the values of the signal are at time instants other than $nT_c$, it's hard to say whether this is actually the case.


Confusion regarding application of time invariance to flip system? I am reading book signal processing first
Chapter 5

I am reading article of time invariance

Where author mentions example of a flip system as shown in attached photo ,i have drawn a red curve above confusion equation

I am confused, how he is doing this when he say"if we delay input and then reverse the order"

I think he is not properly reversing the order,as he reversed the sign of only "n" but he didn't reversed the sign of "n0"

Why this discrepancy in reversal/flipping?

The book is correct, there is no discrepancy. When we reverse a system in time, only the time-variable will get negated and not the shift. Time-reversal does not mean that the whole argument of $x[n]$ gets negated.

Take example of a sequence : $x[n] = {\hat{0},1,2,3,4,5}$. Shift this by 2 samples, so $x_{shift}[n] = {\hat{0},0,0,1,2,3,4,5}$. Now, reverse the time, what do you think should be expected sequence : $x_{final}[n] = 5,4,3,2,1,0,0, \hat{0}$.

You will get it only when $x_{final}[n] = x[-n-2]$ and not $x[-n+2]$. Check it out.
Whatever operation happens, time scaling, shifting and reversal, it happens on n and not on the complete argument.

That is why when the author, reverses the time first and then delays, the shift gets negated too.

Time-reversal : $x_{reverted}[n] = x[-n]$ and then shift of $n_o$, so again shift will happen to $n$ and not to $-n$, and hence, $x_{final}[n] = x[-(n-n_o)]$.

Lets take the same sequence as an example : $x[n] = \hat{0},1,2,3,4,5$.

Revert time first : $x_{reverted}[n] = 5,4,3,2,1,\hat{0}$

Now, when delay the sequence by 2 samples, the expected sequence should look like : $x_{final}[n] = 5,4,3,\hat{2},1,0$. Check that this would happen iff $x_{final}[n] = x[-n+2]$ and not if $x_{final}[n] = x[-n-2]$. Just the put the values of n and it will be clear.

It always helps to take an example sequence and do the shifting, delaying and reversal on that sequence and then see what should you do to the argument on $x$ to get he final expected sequence.


Designing a root-locus compensator to attain design requirements I am working with Matlab to design a compensator capable to make a system attain the design requirements ($\omega_n\geq 0.3$ rad/s & $\xi\geq 0.5$).

In order to do so I have a discrete-time transfer function with a sampling time of 1 second, which is the following:

$$G(z)=\dfrac{0.5z+0.5}{z^2-2z+1} $$

The functions has already a zero = -1 and a double pole = 1.

I am also using Matlab SISOtool to see where is the feasible region in which I should have the root-locus to fulfill the design requirements. Such region is the following (the white one):

The image shows the root-locus of the plant, G(z). Now I want to add a controller, C(z), that allows me to move the root locus into the white region (at least a part of it).

The controller should be in the form:

$$C(z)=\dfrac{(z-z_1)(z-z_2)...(z-z_n)}{(z-p_1)(z-p_2)...(z-p_m)}$$

And therefore, $C(z)G(z)$ should be evaluated to see that part of the root locus is located into the white area.

I want to find the value of the poles and zeros of the controller $C(z)$ that allow me to have part of the root locus into the white region provided in the image, and then being able to fulfill the design requirements.

This is driving me crazy because I have tried all the possible ways of adding poles and zeros and I am unable to obtain a region of the root-locus inside of the white area. I only need some idea on where to place the zeros/poles or if it is even feasible to attain this design requirements (it should be).

Matlab code:

n=[0.5 0.5];
d=[1 -2 1];
G=tf(n,d,1);

After that, to run SISOtool just type sisotool('rlocus',G) in Matlab console.

I don't have access to matlab, but I have played with this in R (see code below) and can confirm that if you're trying to transform your $G(s)$ into:
$$
H(s) = \frac{K G(s)} {1 + KG(s)}
$$
using only a scalar $K$, then you cannot get to the target pole locations (based on the values of $\xi$ and $\omega_n$).

In the diagram:

The black circles are the starting root (pole) positions.
The green circles are the target pole positions.
The red circles are the root locus for varying values of $K$.
The large black circle is the unit circle.

R Code Below

#30341

library("signal")

xi <- 0.5
omega_n <- 0.3

target_den <- c(1, 2*xi*omega_n, omega_n*omega_n)

target_roots <- roots(target_den)

num <- c(0, 0.5, 0.5)
den <- c(1, -2, 1)

start_roots <- roots(den)

par(pty="s")
plot(Re(start_roots),Im(start_roots),xlim=c(-5,5),ylim=c(-5,5), col="black", lwd=5)
points(Re(target_roots), Im(target_roots), col="green", lwd = 3)
for (k in seq(-100.0, 100.0, 0.1))
{
  rts <- roots(k*num + den)
  points(Re(rts), Im(rts), col="red")
}

theta <- seq(-1,1,.01)*pi
lines(sin(theta), cos(theta), col="black")


What are some suitable methods for removing low-frequency line noise from measured data prior to system identification/deconvolution? I am performing a system identification by striking an experimental model with an instrumented impact hammer and measuring the strain response.

The timebase for the recorded signals is 5 microseconds and the hammer pulse width is typically around 200 microseconds. The data are recorded for 20 milliseconds prior to the hammer hit and for 60 milliseconds after.

An example of the recorded data is provided here. A close-up of the pre-trigger strain signal is provided here.

Prior to the hammer hit, the strain signal is non-zero with both low and high frequency noise present. The low frequency noise is of most concern as this affects the quality of the system impulse response function that I subsequently compute.

Unfortunately when completing the tests I did not record a zero-reference prior to each hammer hit. Thus my questions are:

1) What methods are available for removing the low-frequency line noise from the entire signal given that the strain response is also of low-frequency.

2) Is is possible to use the first 20 milliseconds to construct a longer zero-reference signal that I may then subtract from the recorded data? If so, how?

Taking a quick look at your acquisition: this doesn't look too terrible in terms of noise, however, the acquisition window looks way too short. With a 60ms window your frequency resolution is only 160 Hz. From the looks of it, you will also get a massive truncation error. As a rough rule of thumb, you want the acquisition window long enough so that the impulse response visually disappears in the noise.

Specific answers (although they may not help given the truncation error)

If it's really line (like in AC supply voltage hum) noise you can measure the exact location and notch it out for the fundamental and the harmonics.
Sure. Take the mean and subtract. However, again, the window is way too short to get decent low frequency resolution.

A bunch of tips

Make the acquisition window long enough to capture the complete impulse response
Take an acquisition with no excitation. That's your noise floor
Take an acquisition with excitation. That's your signal
Compare signal spectrum to noise floor spectrum. Calculate SNR as a function of frequency
You want at least 10 dB SNR for a decent measurement. Frequencies must either meet that or be so low in level that they don't matter (< at least 20 dB of "the bulk")

If that can't be done there are more esoteric methods like Wiener filtering and/or reconstructive modelling, but it's certainly preferable to clean up the measurement process.


Why is the carrier suppressed by AM-DSB-SC modulation? Lets say I have a message that is a sine wave at 1KHz. As I understand it, I can perform amplitude modulation by multiplying this signal with a carrier, lets say 10KHz, and this process is essentially DSB-SC modulation. Looking at the resulting signal in a spectrum analyzer I see a spectral line at 9KHz and 11Khz. My question is, why is there no spectral line at 10KHz? If you look at the waveform the 10KHz carrier wave is still clearly visible, it simply periodically inverts its phase 180 degrees. Is the fact that no spectral line appears in a frequency analyzer simply because it is being hidden, due to a limitation of how the spectrum analyzer works, effectively canceling out the 10KHz signal over time? If so, then is the carrier in DSB-SC truly suppressed? What does it actually mean for the carrier to be suppressed?

Edit. For clarification, I already know that the math says there should only be frequency content above and bellow the carrier frequency. Perhaps my confusion will become more clear if instead of modulating the carrier with a sine wave we modulate it with a square wave. Bellow is a picture of such a signal in Audacity. According to the spectrum analyzer this signal contains nothing at the original carrier frequency. This is despite the fact that the signal is literally made up of sections of the otherwise unmodified carrier, simply shifted 180 degrees periodically. It's not intuitive that a waveform that can be constructed from sections of a single sine wave of a particular frequency would not contain that frequency at all.

It's a simple trigonometric identity:

$$\cos(\omega_mt)\cos(\omega_ct)=\frac12\left[\cos((\omega_m+\omega_c)t)+\cos((\omega_m-\omega_c))\right]\tag{1}$$

Multiplying two sinusoids of different frequencies results in the sum of two sinusoids with the sum and the difference of the two frequencies.

So this type of AM (DSB-SC) really results in a suppressed carrier, i.e., the modulated signal has no component at the carrier frequency. This is true as long as the message signal has no DC component.

What happens is that the spectrum of the message signal is simply shifted to $\pm \omega_c$, so DC gets mapped to the carrier frequency $\omega_c$, and as long as there is no DC component in the message there is also no carrier component in the modulated signal.

EDIT: Your square wave example is perfect. A square wave is a piecewise constant function, so you could say it is just switched DC, or, paraphrasing your question, it's "made up of sections of otherwise unmodified DC". However, if the duration and amplitude of the positive and the negative part of the square wave are equal, then the corresponding signal's DC component is exactly zero. If you accept that, you must also accept that the signal in your example has no component at the carrier frequency.


How to extract single side of signal spectrum in SIMULINK Signal spectrum have two side, positive and negative. I want to make these separate in two signal, by MATALB SIMULINK.
But how? I can't find it's block on DSP or communication toolbox.
I found a block that interactively generate filter, but this filters response is conjugate symmetric and that's not my desire.
As replica of Fat32 :
Hilbert approach not worked?!!

If you do not have a single side band separation block-diagram in Simulink, then you would try the following algorithm to separate the upper sideband from the lower sideband given that you have a block called Hilbert transformer.

Given that your original real valued baseband input signal is $x[n]$, then the following signal is called as the analytic signal:

$$  x^+[n] = x[n] + j \hat{x}[n] $$ where $$\hat{x}[n] = \mathcal{H} \{x[n] \} = \frac{1}{\pi n} \star x[n]$$ is called as the Hilbert transform of the signal $x[n]$.

It can be shown that the Fourier spectrum of the analytic signal is:
$$ X^+(\omega) = \begin{cases}
2 X(\omega) &, \text{for } ~~~  0 < \omega < \pi \\
0 &, \text{for } ~~~ -\pi < \omega < 0\\ \end{cases}
$$

Which is 2x the upper side band of the signal $x[n]$


Image processing functions for quality control I have to use image processing for quality control of a cookie manufacturer (as a part of my project). The program must be able to detect oversized, undersized, distorted cookies, cookies having cracks, etc. I am new to opencv. Can you please suggest which all functions would help me?

Your first task is to collect examples of various defects and narrow the requirements down from "find defective cookies" to "identify these particular precisely-defined categories of defects".  You will also want to define allowable failure rates (and how much you prefer false positives vs. false negatives) if you want to have a requirement that can actually be met.  Once you have that, you can start designing the optics and lighting you will need to acquire images that highlight the various defects.  This might require more than one imaging operation, as different lighting/optical conditions will highlight different features.

Once the above is done, you can start capturing images; once you have images, you can start to explore what vision tools will work for identifying the various defect categories you defined previously.  A few thoughts:

Identify cookie shape/size.  A back-light is often ideal for this kind of application, but you can sometimes get away with more direct lighting if you have a good even background that is easy to identify.  Your first goal is to separate the cookie from the background; next, you can use a tool like Blob Analysis to look at the distribution of the cookie's shape.
Look for cracks.  Here you will want to investigate edges.  Using the segmented image from above (so you know what part of the image is 'cookie'), you use various edge filters, like the Sobel operator, to identify areas with sharp transitions.  You would then want to experiment with various 'normal' and 'cracked' cookies to understand how the edges vary between them.  You might be able to look at a threshold for the total intensity of edges in the cookie to classify a defect, but this is something that will require more investigation before you can make a good design decision.

One side note: the process of designing lighting/optics, taking images, designing/implementing algorithms, and testing their results is, in practice, an iterative one.  You will learn from your experiments what areas you need to improve and will make adjustments to your earlier imaging system designs to better achieve the results you're after.  You will also likely discover new categories of defects or new ways that the system behavior can change that you did not anticipate.  A long-term systems approach to this kind of task involves a lot of ongoing monitoring and data gathering to characterize the performance of the system over time.  You might find interesting trends like an increase in failures between 1:00 and 3:00 pm, when external light comes through a window and happens to interfere with your processing, for instance.


Extracting filter coefficients from raised-cosine frequency response The frequency response of filter transfer function is given as

$$H(j\omega)=
\begin{cases}
0 &, \textrm{if }\quad (1+r)\frac{\pi}{2} < \lvert \omega \rvert < \pi  \\
 1 + \cos \left ( \frac{\pi}{2r}\left ( \frac{\lvert 2\omega \rvert}{\pi} + r-1 \right ) \right ) &, \textrm{if }\quad (1-r)\frac{\pi}{2} < \lvert \omega \rvert\leq  (1+r)\frac{\pi}{2}\\
 2&, \textrm{if }\quad \lvert \omega \rvert\leq (1-r)\frac{\pi}{2}
\end{cases}
$$

where $r = 0.1$.
I need to extract filter coefficients for the given FIR filter response or order $24$. How can I perform the task using MATLAB?

This is a well-known filter: it's called raised cosine filter. The best strategy is to explicitly compute the inverse Fourier transform of the given frequency response. The corresponding result can be found in the wikipedia article linked to above, and it's also given in this article, which has some additional information on pulse-shaping filters.

Since the corresponding impulse response is of infinite length, you need to truncate it symmetrically. Since the given filter order is $24$, you have 25 filter taps (assuming a common FIR solution). This means one tap at $n=0$, and $12$ taps for $n>0$ and $n<0$, respectively. For truncating the ideal impulse response, you can use any type of window. The most straightforward way is simple truncation, i.e., using a rectangular window.

You can use Matlab to evaluate the formula for the filter's impulse response, and for applying a window.


Derivation of the lowpass to bandpass transformation I have a basic question.
The "well known" lowpass to bandpass transformation is
$$
s \longmapsto \frac{\bar{s}^2 + \omega_1\omega_2}{\bar{s}(\omega_1 - \omega_2)},
$$
which gives a bandpass transfer function of
$$
\frac{1}{s + 1} \longmapsto \frac{\bar{s}(\omega_1 - \omega_2)}{\bar{s}^2 + \bar{s}(\omega_1 - \omega_2) + \omega_1 \omega_2}.
$$
My intuition is that a bandpass should be the product of a lowpass and a highpass.  However, this product gives a different transfer function:
$$
\frac{\omega_1}{s + \omega_1} \frac{s}{s + \omega_2} = \frac{\omega_1 s}{s^2 (\omega_1 + \omega_2) s + \omega_1 \omega_2},
$$
which indicates that the bandpass transformation does not give this cascade of lowpass and highpass.

My question is, how is the bandpass transformation designed, in terms of either combining lowpass filters or by pole placement?

Related question, but using a different derivation technique, and reference is made to the lowpass/highpass derivation, but it is not shown: How is the lowpass to bandpass transformation derived?

You're right that the multiplication of a low pass and a high pass filter results in a band pass filter, as long as the cut-off frequency of the low pass is higher than the cut-off frequency of the high pass. The problem with that approach is that low pass and high pass filters with magnitude responses that are optimal according to some chosen criterion (Butterworth, Chebyshev, Cauer) will not result in an optimum band pass filter.
On the other hand, mapping a single optimal filter will result in another optimal filter. Using $\omega_l\omega_u=\omega_0^2$, where $\omega_l$ and $\omega_u$ are lower and upper band edges, respectively, and $\omega_0$ is the center frequency of the band pass filter, and leaving out constants for the sake of simplicity, the transformation can be written as
$$s\longmapsto \frac{s^2+\omega_0^2}{s}\tag{1}$$
[Note that $\omega_l$ and $\omega_u$ are denoted as $\omega_1$ and $\omega_2$ in the OP, but they're used in a different way in the figure below.]
The mapping $(1)$ maps DC ($\omega=0$) to the desired center frequency $\omega_0$. Furthermore, $s=\pm\infty$ is mapped to $s=0$ and $s=\infty$. So the whole frequency axis of the low pass filter is mapped to the positive frequency axis of the band pass filter. (The same is true for the negative half axis of the band pass filter):

(from: Digital Filter Design by Parks and Burrus)


FFT: do I need to include gaps in data in my FFT domain? I am currently doing an FFT on a very large amount of vibration data, which I am going to sample with no pattern which I want to PSD eventually

My question is in essence can I go from my data

[1, 3, 5, 10, 2, 4, 3 , 23, 24]

then sample it to

A = [0, 0, 5, 10, 0, 0, 0, 23, 24]

If I fft(A),  will it give me the same answer as if I fft(B), where

B = [5, 10, 23, 24]

I think this is fine because if A[n] = 0 it provides no contribution to the sum in the DFT but this feels very wrong and I don't know anywhere near enough about this stuff to feel confident.

I can't access data set A because the way my data is sampled can only give me B.

More detail:

I have very large time series data for vibrations caused by air flow over an object, I am trying to parse that data to identify what happens to the vibrations at specific wind speeds and wind directions (this is the sampling I am talking about, and why its sporadic) .

However the data comes in a strange file format (.gcf) which I can only deal with using obspy, and the only way I can find to sample the data (because there is so much and my computer is rubbish) is using obspy's trim() function, which essentially deletes the rest of the not useful data. (I know this is stupid but I cant find a better way) And so this is what is reducing the dimension of the data .

And the end goal is to create a PSD plot but the problem is the fft really.

A = [0, 0, 5, 10, 0, 0, 0, 23, 24]

  If I fft A, will it give me the same answer as if I fft B, where

  B = [5, 10, 23, 24]

No.

The reason for this is not in what effect does the sample has to the amplitude but in that you would effectively be changing the phase. That is, the relative timing of the sample with respect to the rest.

The Discrete Fourier Transform for a sequence of samples is:

$$ X_k = \sum_{n=0}^{N-1} x_n \cdot \left(\cos\left(\frac{2 \pi kn}{N}\right)-i \sin\left(\frac{2 \pi kn}{N}\right)\right)$$

Where:

$x_n$ is the $n^{th}$ sample of your data sequence $x$ (your vibration data)
$N$ is the length of the sequence

Given the sequence:

$$A = [0, 0, 5, 10, 0, 0, 0, 23, 24]$$

Each sample occurs at a rate of $\frac{1}{N}$. Consequently, the $\cos, \sin$ functions "catch" a given sample at a precise timing which contributes to the weighting to the sum.

For example:

The value of $10$ occurs at $n=3$ which means that the argument of the trigonometric functions at that time instant would be $\frac{2 \pi k 3}{9} = \frac{2 \pi k}{3}$.
This means that during your summation, when you are assessing how much does that particular sample contribute to a sinusoid of frequency $k$, you will be catching that sinusoid at phase $k \cdot \frac{2 \pi}{3}$. The value $10$ occurs at $n=3$ and that sample at that amplitude and position in time contributes $\cos(1 \cdot \frac{2 \pi}{3})$ to a sinusoid of frequency $k=1$, $\cos(2 \cdot \frac{2 \pi}{3})$ to a sinusoid of frequency $k=2$ and so on (and of course, respectively to the sum corresponding to $\sin$).
By isolating certain samples from $A$ to construct a $B$ vector that does not respect the relevant timings, you are altering the time instant that a $\cos,\sin$ at some frequency $k$, "catches" the sample. Very simplistically, if the sample catches the $\cos, \sin$ at a peak, it contributes $1.0$, if it catches it at a "nodal point" it contributes $0$. If you start moving the samples around without respecting their original positions in the time series, your sums are not going to be "right".

So, what do we do?

Two things:

You need to become familiar with the concepts of normalised frequency and the DFTs "frequency resolution".

The fact that you have DFT's of different lengths on its own does not permit you from comparing the two frequency distributions. All that you have to do is evaluate the shorter DFT with more points (but with the same sampling frequency).

There is nothing stopping you from applying an "irregularly" sampled DFT as long as you also preserve the timing of the sample. That is, instead of creating a $B$ that is a vector, you need to create a $B$ that is a $N \times 2$ "vector" where the first column is the time instance that the sample occurs and the second column is the amplitude of the signal at that time instance.

Number 2, would be the most straightforward way to get to what you are after. Another thing that you can do is to take the irregular samples and "resample" them on to a regular grid. You can do this with a number of ways but it basically boils down to interpolation.

Hope this helps.


How to do people detection from bird eye view? I'm trying to do people detection from bird eye view.

Actually, it is not a normal bird eye view because I don't have RGB images but I have disparity map, depth map as well.

In these kind of images, people are shown like lines that represent a single side of the person. I try to re-explain it: these lines represent the front or the back of the person.

I searched for papers on internet but I got no results.

Anyone have any information about that?

PS: I'm working on Matlab.

The following image shows how people looking like from a bird-eye view depth map. You don't look at the big stain. The 3 red circles contain the bird-eye view depth map of (from the left to right): a street fountain, a person and a roadsign.

It appears that motion is your only cue in this case, assuming that you have a time-series of these images. I would start by trying to track these objects using a Kalman filter with a constant velocity motion model. The fountain and the street sign would be stationary. If your images contain other moving objects, e. g. cars, then you can try using speed to distinguish them from people.

Here is an example of how to track multiple objects using the Computer Vision System Toolbox for MATLAB. The example detects moving foreground objects using the Gaussian mixture models. You may want to replace that part or augment it using the depth information.


Why a DFT of two sinusoids is very noisy even with frequency sampling 5 times higher? I've set up this case to try to understand DFT implementing a real case in Excel

Frame Size $\;\color{blue}{(T = 5 \; s})$
Time Sampling $\;\color{blue}{(T_S = 0.1 \; s})$
Block Size $\;\color{blue}{(N = T/T_S = 50)}\;$  (It was 51 in my question)
Sampling Rate  $\;\color{blue}{(F_S = T/T_S+1 =10 Hz}$
Frequency Resolution $\;\color{blue}{(F_R = F_S/N = 0.2 Hz)}$

In time spectrum I have put just the sum of 2 sin functions

The first sine  has amplitude $\color{blue}{10}$, phase $ \color{blue}{\pi/3}$ and frequency  $\color{blue}{2}$ Hz.
The second sine  has amplitude $\color{blue}{5}$, phase $\color{blue}{0}$ and frequency $\color{blue}{1}$ Hz and for $ \color{blue}{\; n = 0\; to\; 49 \;}$ we have :

$$ \color{blue}{x_n = 10\, \sin{(2\pi t 2 + \pi/3)}+ 5\, \sin{(2\pi t 1)}}$$

So I have generated DFT values where $\color{blue}{\; k=-25\; to\; 24\, }$   :

$$\color{blue}{\quad X_k=a_k + ib_k}\;$$

$$\color{blue}{a_k = (1/N) \; \Sigma_{n=0,N-1} x_n \cos{(2 \pi n k /N)}} $$
$$\color{blue}{b_k = (1/N) \; \Sigma_{n=0,N-1} x_n \sin{(2 \pi n k /N)}} $$

After, I've calculated module $\color{blue}{A}$ and phase $\color{blue}{\Phi}$ for $ \color{blue}{\;  k = -25\; to \;24\,}$:

$$\color{blue}{A_k = |X_k| = \sqrt{(a_k)^2  + (b_k)^2}} $$
$$\color{blue}{\Phi_k = atan2(b_k, a_k)} $$

I consider in my analysis for $ \color{blue}{\;  k = -25\; to \;24\,}$:
$$\color{blue}{ Freq_k = k.F_S/N} $$

Old text: Now my doubt arise. In my understanding, the signal is concentrated in the values of frequencies $\color{blue}{-1}$ Hz, $\color{blue}{-2}$ Hz, $\color{blue}{1}$ Hz  and $\color{blue}{2}$ Hz. The negative frequencies I disregard, of course. However, in short, the frequency spectrum is noisy. For instance, $\color{blue}{2.2}$ Hz amplitude is around $\color{blue}{25%}$ of $\color{blue}{2}$ Hz amplitude. Is it normal? What's happening?

New text As Fat32 has showed me, I was using an additional and spurious point ($\color{blue}{51}$ in total) in my time sample that was distorting the whole calculation. If I use a specific time frame ($\color{blue}{5\;s}$ in this case), I need to suppress the final time point (time = $\color{blue}{5\;s}$)

$$\star\star\star$$

I'm working with frequency sampling $\color{blue}{10}$ Hz,  five times than maximum detected frequency ($\color{blue}{2}$ Hz).

To check if I had made a mistake, I had calculated the inverse DTF.

Old text:  The original signal was restored well, where $\color{blue}{\;n=0\; to\; 50}\;$ and imaginary part $\;\color{blue}{\sim 0}\;$.

New text:  The original signal was restored with no flaws, where $\color{blue}{\;n=0\; to\; 49}\;$ and imaginary part $\;\color{blue}{= 0}\;$.
The only positive frequences other than $\color{blue}{0}\;$ is $\color{blue}{1\;hz}\;$ and $\color{blue}{2\;hz}\;$ , respectively with value $\color{blue}{125}\;$ and $\color{blue}{250}\;$, $\color{blue}{25}\;$ times bigger than the real amplitude. I guess that the factor is $\color{blue}{25}\;$ because is the half the samples ($\color{blue}{N=50}\;$).

$$\star\star\star$$

$$\color{blue}{ Re(x_n) =(1/N)\Sigma_{k=-N/2,N/2-1}A_k \cos{\Phi_k} \cos{(2 \pi k n/N)} - A_k \sin{\Phi_k} \sin{(2 \pi k n/N)}}$$

$$\color{blue}{ Im(x_n) = -(1/N)\Sigma_{k=-N/2,N/2-1}A_k \sin{\Phi_k} \cos{(2 \pi k n /N)} + A_k \cos{\Phi_k}\sin{(2 \pi k n /N)}}$$

The 50 time data (it was 51 time date before the answer) is available from this link

The $\color{blue}{50}$ transformed data(Re(z), Im(z), |z|, $\Phi$(z))
(it was $\color{blue}{51}$ time date before the answer) is available from this link

As far as I understood your question, you are complaining about the windowing effect and its consequence of smearing of the spectrum.

So you have a mathematically defined pure sine wave that extends from minus infinity to plus infinity in time, whose theoretical spectrum, you expect to be, a single line in frequency. But you eventually represent that infinitely long ideal sine wave with a finite length record (the windowed version, aka truncated) in practice and take the DFT of the finite length version instead. Then the result will be a broadened single line into a sinc waveform mostly having nonzero amplitudes at all frequencies in general which is what you called noise?

That's the most fundamental manifestation of the effect of time domain windowing (aka aperture effect) on the frequency spectrum; the smearing of the single frequency line and also leakage when multiple tones (including DC) co-exist.

So those are not noise, but an irrecoverable side effect of the practical utilization of finite length DFT. Note however that you might be adding your computational error on top of that. Hence for best verification, you better test your algorithm with an established library from Matlab / Octave kind of platforms.

The following plots are generated by Matlab from your uploaded data, to test your problem and it confirms with the explanation; i.e., those nonzero values you see in the spectrum are due to the windowing effect...

You cannot get rid of it but the effect will minimize by using longer data, however is not the most practical approach. Instead by using special window types such as Hamming, Hanning, Blackman, Kaiser etc, you can tradeoff certain features of it. Incidentally you are using a rectangular window, whose peaks are sharpest but those tails (side lobes) are also largest. All other windows will supress the tails (side lobes) but widen the peaks...

And this is the result of Hamming window applied to your data :

Note that your signal $x[n]$ seems 1 sample longer than it would be...


How is Point Spread Function (PSF) related to Image Priors in Blind Deconvolution? We are researching for our thesis about enhancements in Blind Deconvolution Image Deblurring Algorithm Applied in Image Restoration.

What really is PSF? Is it one kind of image prior? What do these both terms have to do with Blind Deconvolution?

Please enlighten us. Thanks! xo

A "point spread function" denotes how a single point from the source, or a punctual object (mathematically a Dirac) would spread on the observed image through the imaging system.

In blind deconvolution, both image priors and PSF priors are useful to unmix those intricate systems (even more with noise). They are different, as the PSF is somehow invariant to the imaging system (there exist space variant PSFs or blurs, let us skip that for the moment), while different images can be acquired.

PSF priors may be related to some parametric shape (e.g. a Gaussian of unknown width), a bounded support in space or frequency, smoothness, etc, while data priors are often more generic to some image class, related to sparsity, bounded variation, etc.


High Speed Digital Signal Processors Can anybody suggest me the highest speed DSPs available in the market now.
Also, what would be the best for a filer design application a DSP or an FPGA?

This question is too broad and without some more information about the application and the constraints impossible to answer. For actual filter "design" a dedicated DSP is probably not a great choice and you are better off with a high end general purpose processor such as an Intel Xeon or then new-core I7s, Xeons also work well in multiprocessor setups.

Filter implementation really depends on whether you have high-complexity/low-clock or low-complexity/high-clock rate. Best solution for an 100k+ tap room response filter for audio, or a line EQ for a cable modem are totally different. In general for low complexity applications ASICs are faster than FPGAs are faster than general purpose DSP processors.