How do I evaluate the performance of a segmentation algorithm? I'm currently on a project of liver tumor segmentation. I segmented the liver using region growing, and I have to assess how accurate the result is.  I recently learnt that there are certain metrics for assessing the accuracy of the region growing algorithm segmentation, like Tanimoto Coefficient, Correlation, etc. But I don't know how to implement them in Matlab. Check out https://stackoverflow.com/questions/9553204/tanimoto-coefficient-using-matlab

Given that you are working only on Tanimoto Coefficient, i am trying to be more specific rather than giving generic answer with various different approaches.

The basic notation of Tanimoto Coefficient is as follows:

$$
T(A,B) = { N_{AB} \over N_A + N_B - N_{AB} }
$$

where $T$ is the desired result, over $A$ and $B$ images.

In this measure, we identify pixels either as belonging to a given segment i.e. it is segment pixel or it is background.
$N$ refers to number of pixels that are classified as segment pixel in the respective image. And $N_{AB}$ refers to the number of pixels which are classified as segment pixel in both images.

In this measure all pixels which qualifies as neither in A nor B is not calculated; only the pixels.

Also, both image must have same resolution and must have identical locations to for the segmented objects else even if the segmentation shape is right, the resultant overlap may not be right.

I am not getting into your MATLAB code, but here is the pseudo code that looks like.

Initialize N_A, N_B, N_AB;
for( all pixels @ x,y)
{
   if(image_A[x][y] == SEGMENT_CLASS_PIXEL)
    N_A += 1;

   if(image_B[x][y] == SEGMENT_CLASS_PIXEL)
    N_A += 1;

   if(image_A[x][y] == SEGMENT_CLASS_PIXEL
      && image_B[x][y] == SEGMENT_CLASS_PIXEL)
    N_AB += 1;

}

T = N_AB / (N_A + N_B - N_AB);


Scaling of the PWELCH function in MATLAB I am trying to compute the power spectral density of a random signal using the PWELCH function in MATLAB.

Since I think have not understood properly how Pwelch scales the PSD, I wrote a sample program, in which I generate a sum of two sinusoids with given amplitudes - A1 and A2 - and given frequencies - f1 and f2 - superimposed to random noise.

I computed the PSD using the FFT and then using Pwelch. It seems to me that the PSD computed via FFT is correct, because it returns the power corresponding to the sinusoids (peaks at $A1^2$ and $A2^2$).

However Pwelch seems to have a different scaling because the peak values depend on the number of samples I put in the PWELCH window. This should make sense up to a certain extent, but the variation seems too large to me. I think I haven't properly understoon either how to use PWELCH or how to scale the PSD returned by PWLECH.

Can you help me understand this? Thanks. The sample code is reported below.

Thank you
E.

clear all
close all
clc

rng default

Fs = 10000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 50000;             % Length of signal

t = (0:L-1)*T;

A1=2;
A2=3;

f1=5;
f2=20;

X=A1*sin(2*pi*f1*t) + A2*sin(2*pi*f2*t)+ 2*randn(size(t));

figure(1)
hold on
grid on
plot(t,X)
xlabel('time (s)')
ylabel('signal (Unit)')

%% Power Spectrum via FFT
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L

Y = fft(X);

P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = Fs*(0:(L/2))/L;

%% via FFT

figure(2)
hold on
grid on
plot(f,P1.^2)
title('Single-Sided Power Spectrum of X(t)')
xlim([0 50])
xlabel('f (Hz)')
ylabel('PSD (Unit^2/Hz)')
set(gca,'xscale','log')
% set(gca,'yscale','log')

%% via PWELCH

nfft= L;
window= rectwin(nfft);

[Pxx_f,f]=pwelch(X,window,[],nfft,Fs);

figure(2)
plot(f,Pxx_f,'r','linewidth',1.5);

legend('via FFT','via Pwelch','location','best')

I believe, there are two aspects of your question.

For the FFT calculation, the power density spectrum is the square of the FFT modulus multiplied by the sampling period (in your notation: T) and divided by the number of samples (in your notation:L):
T/L * abs(FFT(X))^2  or   T/L *abs(Y(f))^2

Here is now where I believe you have an error, you have to add the power of Y(f) and (Y(-f)) and not add first Y(f) and (Y(-f) and then perform the power calculation. By adding the amplitude first, in your notation,   P1(2:end-1) = 2*P1(2:end-1)  then after calculating the power, you have the double of the correct value.

A far as the Welch method for the calculation of the Power spectrum, it takes into account all the necessary normalization. You can chose between the power spectrum (option "power") or the power density spectrum (default option "psd")

Trying to use the option "power" for the welch calculation, it will gives as answer the amplitude of the sin squared divided by 2 (a sin as a power of 1/2): (A1^2)/2.

To calculate the power density spectrum, you have to divide the previous results by the applied frequency resolution, that is the sampling frequency by the number of samples, in your notation Fs/L. With the welch function, you have to use the option 'psd'.
Sorry for the long answer and hoping, this helps.


Tinnitus Masking with broadband Noise I have the problem of masking tinnitus. I have the following questions:

Can broadband noise mask tone tinnitus, or any kind of tinnitus (such as narrowband, or chirp)?
Assume an audiogram test has identified tinnitus as tonal with loudness level of $x$ dB and  frequency/pitch of $y$ Hz. How can I create broadband noise  using this information (note that we don't know the variance of the noise)? How could this be implemented in Matlab?

Thanks.

sedy

Read about tinnitus masking on Wiki. Yes, you can mask it with various kinds of noise and/or natural sounds that are "noisy" (ocean, rainfall, etc)

For the synthesis of masking noise, there is plenty of information via google & StackOverflow. Keep in mind that broadband noise (versus bandlimited noise) is present across the spectrum and thus, you only need to be approximately in the same bark frequency band of the tone you want to mask.

Here is a plot showing one measurement of Tinnitus spectrum (from http://www.thelancet.com/cms/attachment/2021743353/2041580240/gr1.jpg). I wasn't able to access the full thelancet article:


Kalman filter continues to update position even though speed is zero Background
I set up a conventional Kalman filter that makes use of smartphone GPS only (no inertial sensors). That is, it uses the position, doppler speed, and course, in order to create positions that represent a more accurate track.
When the device is stationary (i.e. the smartphone is placed on a table), the GPS drifts, and the doppler speed reports 0 m/s, as expected.
Problem and expected outcome
Given that the doppler speeds are all 0, I expect that the Kalman filter would "dismiss" subsequent position measurements and retain the original position measurement.
Details

I convert GPS coordinates from geodetic coordinates to NED (Cartesian) coordinates, $p_x, p_y$, relative to an initial coordinate
$\psi$ is course
$\sigma_h^2=h_{\tt accuracy}^2$, where $h_{\tt accuracy}$ is the horizontal accuracy, in meters, and represents 68% probability that the position is within its radius
$\sigma_s^2=s_{\tt accuracy}^2$, where $s_{\tt accuracy}$ is the speed accuracy, in meters/second, and represents 68% probability that the speed is within its range
$H=I$

Here are the Kalman filter equations, step by step:
$F=\begin{bmatrix}1&0&\Delta t&0\\0&1&0&\Delta t\\0&0&1&0\\0&0&0&1\end{bmatrix} \tag{1}$
$\hat x^-=F\hat x\tag{2}$
$Q=\begin{bmatrix}\sigma_h^2&0&\sigma_s^2&0\\0&\sigma_h^2&0&\sigma_s^2\\0&0&\sigma_s^2&0\\0&0&0&\sigma_s^2\end{bmatrix} \tag{3}$
$P=FPF^T+Q \tag{4}$
$R=\begin{bmatrix}\sigma_h^2&0&\sigma_s^2&0\\0&\sigma_h^2&0&\sigma_s^2\\0&0&\sigma_s^2&0\\0&0&0&\sigma_s^2\end{bmatrix} \tag{5}$
$K=PH^T(HPH^T+R)^{-1} \tag{6}$
$z=\begin{bmatrix}p_x\\p_y\\s\cos(\psi)\\s\sin(\psi)\end{bmatrix} \tag{7}$
$\hat x^+=\hat x^-+K(z-H\hat x^-) \tag{8}$
$P^+=P^--KHP^- \tag{9}$
Question
If $\sigma_h^2$ is large (poor accuracy) and $\sigma_s^2$ is low (good accuracy), I expect that the predicted position will not change, but it continues to be influenced by the measurement and drift. Why is that?

First, you err in using the square root of accuracy in your covariance matrices.  These should be squared.  If your positional accuracy is $\sigma_h$ and your velocity accuracy is $\sigma_s$, then your measurement matrix should be
$$\mathbf R = \begin{bmatrix}
\sigma_h^2 & 0 & 0 & 0 \\
0 & \sigma_h^2 & 0 & 0 \\
0 & 0 & \sigma_v^2 & 0 \\
0 & 0 & 0 & \sigma_v^2 \\
\end{bmatrix}$$
Second, your $\mathbf Q$ matrix is telling the Kalman filter that the position moves a distance with a deviation of $\sigma_h$ each step, and its velocity changes with a deviation of $\sigma_v$ each step.
So you're telling the filter "hey, this phone is moving around in this huge and unpredictable way every single time step".  As a consequence, at each time step the filter is going to put as much or more credence in the measurement as it does in its filtered position and velocity.
If you want the filter to track when you're moving and settle when you put it on a table, you either need GPS/IMU fusion, or you need to "tell" it (by setting the $\mathbf Q$ matrix very small) that it's sitting still when its on the table, and (by setting the $\mathbf Q$ matrix high) that it's moving when it's not on the table.
Figuring out how to do that, absent of full GPS/IMU fusion, I leave to you.


finding phase and group delay from frequency response (From Schaum's Outlines, DSP, second edition, problem 5.25, second part of problem)

What's the procedure to find the Phase and Group Delay of:

$$
H(e^{j\omega}) = e^{-j\theta}\left( \frac{e^{-j(\omega-\theta)}- |\alpha|
}{1-|\alpha|\ e^{-j(\omega - \theta)}} \right)
$$

Book says that group delay should be:

$$
\tau(\omega) = \tau_0(\omega - \theta) = \frac{1-\alpha^2}{|1-\alpha\ e^{-j(\omega-\theta)}|^2}
$$

I want to say the phase is $\phi(\omega) = (\omega - \theta)$ and group delay is: $\tau(\omega) = - \phi ' (\omega)= 1$... but i know that's not correct because the book gives a more complicated group delay.  it looks like book version has something to do with "$\tau_0$", which i'm guessing is the initial group delay at time zero??  actually, i have no idea... its not explained in the book. (Any ideas how to do this same problem using the method in second part of problem?)

I feel there's a 50% chance i need to do this the complicated algebra and calculus way using this expresss as start:

$$
\phi(\omega) = \arg\Big\{ H_R(e^{j\omega})+jH_I(e^{j\omega}) \Big\}
$$

but then, i'm not really sure if there's a shortcut to just convert it to polar form instead of slugging it out in complex Cartesian form alegbra. like it seems the phase is $\phi(\omega) = (\omega - \theta)$...

with a lot of algebra I obtained this result:

$$
\tau(\omega) = \frac{(-\alpha^2_I - \alpha^2_R + 1)}{1 + \alpha_{R}^{2} + \alpha_{I}^{2} -2 \alpha_{R} \cos(\omega) - 2 \alpha_{I} \sin(\omega)}
$$

still would be nice if there were a simpler way..like the way the second part of the problem suggests...but doesn't really explain...

Fancy factoring to match book result:

$$1 + \alpha_{R}^{2} + \alpha_{I}^{2} -2 \alpha_{R} \cos(\omega) - 2 \alpha_{I} \sin(\omega)$$

$$= (1-\alpha e^{-j\omega}) (1- \alpha^{*} e^{j\omega})$$

$$=|1-\alpha\ e^{-j\omega}|^2$$

other factor:

$$|z|^2 = z_I^2 + z_R^2 $$

$$(-\alpha^2_I - \alpha^2_R + 1)$$

$$=(1 -(\alpha^2_I + \alpha^2_R) )$$

$$=(1 - |\alpha|^2)$$

can remove absolute value as long as $\alpha$ remains squared.

$$=(1 - \alpha^2)$$

Define a frequency response

$$\tilde{H}(e^{j\omega})=\frac{e^{-j\omega}-|\alpha|}{1-|\alpha|e^{-j\omega}}\tag{1}$$

and note that the group delay $\tilde{\tau}(\omega)$ corresponding to $(1)$ is related to the group delay $\tau(\omega)$ of the original frequency response $H(e^{j\omega})$ by

$$\tau(\omega)=\tilde{\tau}(\omega-\theta)\tag{2}$$

because $(1)$ is just a scaled and frequency-shifted version of the original $H(e^{j\omega})$.

$\tilde{H}(e^{j\omega})$ is the frequency response of a first-order all-pass filter with a real-valued pole at $z_{\infty}=|\alpha|<1$, and a real-valued zero at $z_0=1/|\alpha|$. The original frequency response $H(e^{j\omega})$ describes a first-order all-pass filter with a complex-valued pole at $z_{\infty}=\alpha=|\alpha|e^{j\theta}$ ($|\alpha|<1$), and a zero at $z_0=1/\alpha^*$.

Now write

$$\tilde{H}(e^{j\omega})=\frac{B(e^{j\omega})}{A(e^{j\omega})}\tag{3}$$

and note that

$$\tilde{\phi}(\omega)=\phi_B(\omega)-\phi_A(\omega)\tag{4}$$

where $\tilde{\phi}(\omega)$ is the phase of $\tilde{H}(e^{j\omega})$, and $\phi_A(\omega)$ and $\phi_B(\omega)$ are the phases of $A(e^{j\omega})$ and $B(e^{j\omega})$, respectively. Consequently,

$$\tilde{\tau}(\omega)=\tau_B(\omega)-\tau_A(\omega)\tag{5}$$

where $\tau_A(\omega)$ and $\tau_B(\omega)$ are the group delays of $A(e^{j\omega})$ and $B(e^{j\omega})$, respectively.

We have

$$\phi_B(\omega)=-\arctan\left(\frac{\sin(\omega)}{\cos(\omega)-|\alpha|}\right)\tag{6}$$

and

$$\phi_A(\omega)=\arctan\left(\frac{|\alpha|\sin(\omega)}{1-|\alpha|\cos(\omega)}\right)\tag{7}$$

where I've ignored any ambiguities by $\pm\pi$ due to the principal value of $\arctan()$ because they are irrelevant for computing the group delay.

Taking the negative derivative of $(6)$ and $(7)$ with respect to $\omega$ gives

$$\tau_B(\omega)=\frac{1-|\alpha|\cos(\omega)}{1-2|\alpha|\cos(\omega)+|\alpha|^2}\tag{8}$$

and

$$\tau_A(\omega)=\frac{|\alpha|^2-|\alpha|\cos(\omega)}{1-2|\alpha|\cos(\omega)+|\alpha|^2}\tag{9}$$

From $(5)$ we get

$$\tilde{\tau}(\omega)=\frac{1-|\alpha|^2}{1-2|\alpha|\cos(\omega)+|\alpha|^2}=\frac{1-|\alpha|^2}{\big|1-|\alpha|e^{-j\omega}\big|^2}\tag{10}$$

From $(2)$ the group delay of the complex first-order all-pass is given by

$$\tau(\omega)=\frac{1-|\alpha|^2}{1-2|\alpha|\cos(\omega-\theta)+|\alpha|^2}=\frac{1-|\alpha|^2}{\big|1-|\alpha|e^{-j(\omega-\theta)}\big|^2}\tag{11}$$

The difference with the given solution in Schaum's outline is the fact that in $(11)$ you have $|\alpha|$ instead of the complex-valued $\alpha$. Note that $(11)$ is correct and Schaum's formula is wrong because a complex-valued group delay doesn't make much sense.


Why does the Hamming window attenuate the FFT? I am working on a very simple voice recognition system using the msp430. In my research I've read that for real life audio signals its necessary to apply a window to the FFT data because the audio signals won't be periodic. In my code, I've applied the hamming window and plotted it against the non-windowed FFT spectrum. It appears that the hamming attenuates the signal in the lower valued N data points.

I think this makes sense based upon the equation of the hamming window, but why would I want this attenuate for my signal? It seems that identifying the spectral characteristics would be more difficult after the window. Any thoughts?

You can figure the change in energy applied to your signal by analyzing the energy in your window.  The Hamming window's energy (or scaling) is well known and Google will tell you what it is.  You are doing something different. Can you post your code for further help?

From your code, ans your image, it really appears that you applied the window to your transformed data.  The window (in this case) is applied to the time data before the transform.  The window will reduce sidelobes from strong frequency components.  Take note that taper windows cause the so-called "spectral leakage" where energy leaks into other frequency bins


Which image transformations are considered continuous transformations? In the article J. Matas, O. Chum, M. Urban, T. Pajdla: Robust Wide Baseline Stereo from Maximally Stable Extremal Regions it is stated that:

  The set (of extremal regions) is closed under continuous (and thus perspective) transformation of image coordinates, ...

I know what perspective transformations (and consequently, affine transformations) are, but could anyone give examples of other transformations that fall under the category of continuous transformations of image coordinates? How would such transformations be characterized?

I think "continuous" is meant in the strict mathematical sense here, i.e. no abrupt changes in the coordinate mapping.

Most non-rigid image alignment algorithms (try to) produce continuous mappings. Distortion filters in popular photo-editing applications (like Photoshop's Warp effect) are usually continuous mappings, too.

ADD: Here are two examples of transformations that are continuous but not affine or perspective:

How to flatten image label on a food jar

Flatten a sudoku grid


Mean Square Error and Gibbs oscillations While studying the convergence of Fourier transform, I got to know two conditions.

$$\sum_{n=-\infty}^{\infty}|x(n)|<\infty$$
$$\sum|x(n)|^{2} \leq [\sum|x(n)|]^{2}$$

While I was reading the text, I found this paragraph quite confusing. I didn't understood this.

  If a sequence is not absolutely summable but has finite energy, one may employ a type of convergence in which the series converges so the mean square error is 0. The "attendant Gibbs oscillations" at a discontinuity are of practical significance in filter design.

I don't know what Gibbs oscillations mean here. So if someone could please explain the meaning of this paragraph (taken for DSP by Alan V. Oppenheim).

The first condition mentioned in your question is absolute summability, which is sufficient for the discrete-time Fourier transform (DTFT) to exist. In this case, the sum given by the DTFT of a sequence converges uniformly. The other condition you probably mean is square summability:

$$\sum_{n=-\infty}^{\infty}|x[n]|^2\lt\infty\tag{1}$$

in which case the DTFT exists if we relax the condition of uniform convergence. That other type of convergence is called mean-square convergence. In this case, the result of the sum oscillates around points of discontinuities, and those oscillations do not decrease with increasing number of elements in the sum. You can find a more in-depth explanation of this so-called Gibbs phenomenon in this article.

The importance of the Gibbs phenomenon in digital filter design is that you might want to approximate an ideal frequency response (with discontinuities) by a filter of finite length (finite impulse response, FIR), and this results in large errors close to the discontinuities. This problem can be alleviated by introducing transition bands and/or by (non-rectangular) windowing instead of simple truncation of the sum.


How to model the noise in Kalman Filter? Background:

I am a newbie in DSP. I am implementing a simple Kalman Filter that estimates the heading direction of a robot. The robot is equipped with a compass and a gyroscope.

My Understanding:

I am thinking about representing my state as a 2D vector $(x, \dot{x})$, where $x$ is the current heading direction and  $\dot{x}$ is the rotation rate reported by the gyroscope.

Questions:

If my understanding is correct, there will be no control term, $u$ in my filter. Is it true? What if I take the state as a 1D vector $(x)$? Then does my $\dot{x}$becomes the control term $u$? Will these two methods yield different results?
As we know, the main noise source comes from the compass when the compass is in a distorted magnetic field. Here, I suppose the Gaussian noise is less significant. But the magnetic distortion is totally unpredictable. How do we model it in the Kalman Filter?
In Kalman Filter, is the assumption that "all the noises are white" necessary? Say, if my noise distribution is actually a Laplacian distribution, can I still use a Kalman Filter? Or I have to switch to another filter, like Extended Kalman Filter?

You are most likely not having a control vector $u$. Maybe you can somehow model it into the filter, but it is not required in your case and uncommon. Furthermore $\overset{.}{x}$ is not your gyro input but it is the Kalman filtered rotation rate calculated after using all your measurements.
Usually, the main effects on a magnetometer are hard iron and soft iron effects. There is a lot of articles out there that describe the cause and how to mitigate those. I would probably do this in a preprocessing step and leave it out of the Kalman filter entirely. Furthermore, if you have some means to detect if you are in a strongly distorted field (e.g. because the absolute strength is going up), you can increase the measurement noise you report to the Kalman filter for those measurements. This will indicate to the filter to not trust the values coming from the magnetometer anymore. So if you know the total strength of the earth's magnetic field and you suddenly measure a much stronger field, you could use this technique. This is also commonly done with accelerometers to put less trust on values measured when the vehicle is actually accelerating and hence, the measurement is not pointing "down". See here for an example
Your noise has to be white and it has to be Gaussian or near Gaussian. You can only have correlation between the different measurements of one time step but not between multiple time steps. If you search around in the academic world, there are many papers presenting adapted Kalman filters for other noise types that work more or less well (see here for example, haven't tested this, though). EKF and UKF approaches don't help here, either. You use those, when you have a nonlinear system model. As already said, if your noise is very non-white and non-gaussian that you get big estimation errors, then a particle filter is an option.


How to compute 2D displacement vector for binary image registration? My inputs are several binary images like:

They have globally the same content but may be unregistered as they are produced from a hand camera. What I would like to compute is the 2D displacement vector from the first image to the second. I use OpenCV and my first attempt was to compute feature points (SURF algorithm) and the affine transform between both images. But of course, feature description are kinda poor on a binary image, therefore the matching is difficult and the matching map very inaccurate.

Does anyone has an idea how I could do this?

I am somewhat surprised that feature points don't work that well. I have had success registering shapes like yours using either

Harris points, this is a corner detector, in combination with the RANSAC algorithm. See the wiki or Peter Kovesi his site
Using a feature detector like SURF or SIFT in combination with an edge map of the image prior to feature detection followed by some form of robust matching.

EDIT

I gave it a go with matlab and tried some variations on this theme. Currently I use

SIFT features from the edge map to ascertain rotation and scaling, although they are very small between the images you provided
RANSAC for robust matching
Cross correlation to ascertain the translation between the two images

Source on github.

An alternative for RANSAC could be the Hough transform voting/binning approach that was proposed by the inventor of SIFT, David Lowe.


Obtaining an IR from a sine sweep without the dry sweep I have a test vinyl record that includes an exponential sine sweep as one of the tracks. I have a recording of the sweep through my turntable/mixer, and now I want to calculate the system IR. The problem is I don't have a copy of the clean sweep signal. What's the best way to approach this?
My first thought is to take a best guess at the sweep start and end frequencies and the $\Delta t$ between them, then use that to figure out the sweep rate and manually reconstruct the original sweep. But this seems kind of clumsy, is there a better way?

If this is used for channel sounding then you aren't concerned with the actual accuracy of your turntable (causing a frequency shift) in which case I suggest measuring the response in a static condition (non-changing channel) and determining the least square linear fit of the frequency sweep (I would assume but don't know that the test pattern would be a linear sweep) and then from that reconstruct a reference sweep to then use for channel estimation.
Measuring the frequency involves using some for of frequency discriminator process, most of which are also sensitive to amplitude variation, so I recommend hard-limiting the signal prior to the discriminator. Also pay attention to the relationship between instantaneous frequency and the underlying frequency of your ramp as detailed in this "DSP Quiz" question:  Simulation of a Frequency ramp


How can I design Nyquist interpolation filters with the Parks-McClellan algorithm? We can easily design interpolation filters that obey certain frequency-domain constraints using the Parks-McClellan algorithm.  However, it's not immediately clear how to enforce time-domain constraints; in particular, I'm interested in generating Nyquist filters.  So if I'm oversampling by a factor of N, I want the filter to have zero-crossings at kN, for non-zero integer k (this ensures that the input samples to my interpolator will appear in the output sequence).

I've seen Harris1 talk about a technique for designing half-band filters, i.e. the special case where N=2.  Is there a general solution for this?  (I know that we can easily design filters with the window method, but that doesn't give us the same control.)

[1] Multirate Signal Processing for Communication Systems, pp. 208-209

One design method, albeit one that is limited to powers of two, would be to start with one halfband filter, insert zeros at every other (creates a spectral replica), then convolve it with a second halfband filter having a wider transition band.  Repeat the process until you get to the required power of 2.

Here's an example that creates a lowpass filter with Fc=fs/8 and zero crossings every 4 samples:

b0=remez(34,[0 .45 .55 1],[1 1 0 0])';
b1=remez(6,[0 .25 .75 1],[1 1 0 0])';
b0up = zeros(1,2*length(b0)-1);
b0up(1:2:end) = b0;
B0up=freqz(b0up);
b2 = conv(b0up,b1);  % length = 34*2+1 + 6 = 75 coefficients


How to flatten the image of a label on a food jar? I'd like to take pictures of labels on a jar of food, and be able to transform them so the label is flat, with the right and left side resized to be even with the center of the image.

Ideally, I'd like to use the contrast between the label and the background in order to find the edges and apply the correction.  Otherwise, I can ask the user to somehow identify the corners and sides of the image.

I'm looking for general techniques and algorithms to take an image that is skewed spherically (cylindrically in my case) and can flatten the image.  Currently the image of a label that is wrapped around a jar or bottle, will have features and text that shrinks as it recedes to the right or left of the image.  Also the lines that denote the edge of the label, will only be parallel in the center of the image, and will skew towards each-other on the right and left extreme of the label.

After manipulating the image, I would like to be left with an almost perfect rectangle where the text and features are uniformly sized, as if I took a picture of the label when it was not on the jar or bottle.

Also, I would like it if the technique could automatically detect the edges of the label, in order to apply the suitable correction.  Otherwise I would have to ask my user to indicate the label boundaries.

I've already Googled and found articles like this one:  flattening curved documents, but I am looking for something a bit simpler, as my needs are for labels with a simple curve.

A similar question was asked on Mathematica.Stackexchange. My answer over there evolved and got quite long in the end, so I'll summarize the algorithm here.

Abstract

The basic idea is:

Find the label.
Find the borders of the label
Find a mapping that maps image coordinates to cylinder coordinates so that it maps the pixels along the top border of the label to ([anything] / 0), the pixels along the right border to (1 / [anything]) and so on.
Transform the image using this mapping

The algorithm only works for images where:

the label is brighter than the background (this is needed for the label detection)
the label is rectangular (this is used to measure the quality of a mapping)
the jar is (almost) vertical (this is used to keep the mapping function simple)
the jar is cylindrical (this is used to keep the mapping function simple)

However, the algorithm is modular. At least in principle, you could write your own label detection that does not require a dark background, or you could write your own quality measurement function that can cope with elliptical or octagonal labels.

Results

These images were processed fully automatically, i.e. the algorithm takes the source image, works for a few seconds, then shows the mapping (left) and the un-distorted image (right):

The next images were processed with a modified version of the algorithm, were the user selects the left and right borders of the jar (not the label), because the curvature of the label cannot be estimated from the image in a frontal shot (i.e. the fully automatic algorithm would return images that are slightly distorted):

Implementation:

1. Find the label

The label is bright in front of a dark background, so I can find it easily using binarization:

src = Import["http://i.stack.imgur.com/rfNu7.png"];
binary = FillingTransform[DeleteBorderComponents[Binarize[src]]]

I simply pick the largest connected component and assume that's the label:

labelMask = Image[SortBy[ComponentMeasurements[binary, {"Area", "Mask"}][[All, 2]], First][[-1, 2]]]

2. Find the borders of the label

Next step: find the top/bottom/left/right borders using simple derivative convolution masks:

topBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{1}, {-1}}]];
bottomBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{-1}, {1}}]];
leftBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{1, -1}}]];
rightBorder = DeleteSmallComponents[ImageConvolve[labelMask, {{-1, 1}}]];

This is a little helper function that finds all white pixels in one of these four images and converts the indices to coordinates (Position returns indices, and indices are 1-based {y,x}-tuples, where y=1 is at the top of the image. But all the image processing functions expect coordinates, which are 0-based {x,y}-tuples, where y=0 is the bottom of the image):

{w, h} = ImageDimensions[topBorder];
maskToPoints = Function[mask, {#[[2]]-1, h - #[[1]]+1} & /@ Position[ImageData[mask], 1.]];

3. Find a mapping from image to cylinder coordinates

Now I have four separate lists of coordinates of the top, bottom, left, right borders of the label. I define a mapping from image coordinates to cylinder coordinates:

arcSinSeries = Normal[Series[ArcSin[\[Alpha]], {\[Alpha], 0, 10}]]
Clear[mapping];
mapping[{x_, y_}] :=
   {
    c1 + c2*(arcSinSeries /. \[Alpha] -> (x - cx)/r) + c3*y + c4*x*y,
    top + y*height + tilt1*Sqrt[Clip[r^2 - (x - cx)^2, {0.01, \[Infinity]}]] + tilt2*y*Sqrt[Clip[r^2 - (x - cx)^2, {0.01, \[Infinity]}]]
   }

This is a cylindrical mapping, that maps X/Y-coordinates in the source image to cylindrical coordinates. The mapping has 10 degrees of freedom for height/radius/center/perspective/tilt. I used the Taylor series to approximate the arc sine, because I couldn't get the optimization working with ArcSin directly. The Clip calls are my ad-hoc attempt to prevent complex numbers during the optimization. There's a trade-off here: On the one hand, the function should be as close to an exact cylindrical mapping as possible, to give the lowest possible distortion. On the other hand, if it's to complicated, it gets much harder to find optimal values for the degrees of freedom automatically. (The nice thing about doing image processing with Mathematica is that you can play around with mathematical models like this very easily, introduce additional terms for different distortions and use the same optimization functions to get final results. I've never been able to do anything like that using OpenCV or Matlab. But I never tried the symbolic toolbox for Matlab, maybe that makes it more useful.)

Next I define an "error function" that measures the quality of a image -> cylinder coordinate mapping. It's just the sum of squared errors for the border pixels:

errorFunction =
  Flatten[{
    (mapping[#][[1]])^2 & /@ maskToPoints[leftBorder],
    (mapping[#][[1]] - 1)^2 & /@ maskToPoints[rightBorder],
    (mapping[#][[2]] - 1)^2 & /@ maskToPoints[topBorder],
    (mapping[#][[2]])^2 & /@ maskToPoints[bottomBorder]
    }];

This error function measures the "quality" of a mapping: It's lowest if the points on the left border are mapped to (0 / [anything]), pixels on the top border are mapped to ([anything] / 0) and so on.

Now I can tell Mathematica to find coefficients that minimize this error function. I can make "educated guesses" about some of the coefficients (e.g. the radius and center of the jar in the image). I use these as starting points of the optimization:

leftMean = Mean[maskToPoints[leftBorder]][[1]];
rightMean = Mean[maskToPoints[rightBorder]][[1]];
topMean = Mean[maskToPoints[topBorder]][[2]];
bottomMean = Mean[maskToPoints[bottomBorder]][[2]];
solution =
 FindMinimum[
   Total[errorFunction],
    {{c1, 0}, {c2, rightMean - leftMean}, {c3, 0}, {c4, 0},
     {cx, (leftMean + rightMean)/2},
     {top, topMean},
     {r, rightMean - leftMean},
     {height, bottomMean - topMean},
     {tilt1, 0}, {tilt2, 0}}][[2]]

FindMinimum finds values for the 10 degrees of freedom of my mapping function that minimize the error function. Combine the generic mapping and this solution and I get a mapping from X/Y image coordinates, that fits the label area. I can visualize this mapping using Mathematica's ContourPlot function:

Show[src,
 ContourPlot[mapping[{x, y}][[1]] /. solution, {x, 0, w}, {y, 0, h},
  ContourShading -> None, ContourStyle -> Red,
  Contours -> Range[0, 1, 0.1],
  RegionFunction -> Function[{x, y}, 0 <= (mapping[{x, y}][[2]] /. solution) <= 1]],
 ContourPlot[mapping[{x, y}][[2]] /. solution, {x, 0, w}, {y, 0, h},
  ContourShading -> None, ContourStyle -> Red,
  Contours -> Range[0, 1, 0.2],
  RegionFunction -> Function[{x, y}, 0 <= (mapping[{x, y}][[1]] /. solution) <= 1]]]

4. Transform the image

Finally, I use Mathematica's ImageForwardTransform function to distort the image according to this mapping:

ImageForwardTransformation[src, mapping[#] /. solution &, {400, 300}, DataRange -> Full, PlotRange -> {{0, 1}, {0, 1}}]

That gives the results as shown above.

Manually assisted version

The algorithm above is full-automatic. No adjustments required. It works reasonably well as long as the picture is taken from above or below. But if it's a frontal shot, the radius of the jar can not be estimated from the shape of the label. In these cases, I get much better results if I let the user enter the left/right borders of the jar manually, and set the corresponding degrees of freedom in the mapping explicitly.

This code lets the user select the left/right borders:

LocatorPane[Dynamic[{{xLeft, y1}, {xRight, y2}}],
 Dynamic[Show[src,
   Graphics[{Red, Line[{{xLeft, 0}, {xLeft, h}}],
     Line[{{xRight, 0}, {xRight, h}}]}]]]]

This is the alternative optimization code, where the center&radius are given explicitly.

manualAdjustments = {cx -> (xLeft + xRight)/2, r -> (xRight - xLeft)/2};
solution =
  FindMinimum[
   Total[minimize /. manualAdjustments],
    {{c1, 0}, {c2, rightMean - leftMean}, {c3, 0}, {c4, 0},
     {top, topMean},
     {height, bottomMean - topMean},
     {tilt1, 0}, {tilt2, 0}}][[2]]
solution = Join[solution, manualAdjustments]


Any basic paper on biomedical signal processing to practise with? I am quite interested in biomedical signal processing and I am looking for some biomedical signal processing papers in order to reproduce the results that they obtained. The idea is to redo what they did so I can learn about the topic and improve my research skills.

The papers should fit the following requisites:

The data that they use must be public.
Most of the work could be done in Matlab or similar software.

Papers about biomedical image processing are also welcomed.

Thanks!

I would suggest you go to PLOS Biology and do a search for "Signaling MatLab".  There are many scientific papers there on various biology research topics.  Some even give the Matlab commands used.


Intuitive explanation of the Fourier Transform for some of the functions Does anyone have a mechanism to understand intuitively  (and automatically) why the Fourier Transform of certain functions have certain shapes (at least for some functions, not necessarily for all)? I know what kind of operator the Fourier Transform is and what it does to a function but somehow I can't see intuitively and automatically why why the Fourier Transform of certain functions have certain shapes. For example, is there a intuitive reason why the Fourier Transform of a pulse (box function) is in the $\mathbb{sinc}$?

Thanks in advance.

What might seem intuitive differs greatly between individuals.  But let's start with a few basic things about a Fourier transform that people who have studied it might know about it.  Intuitively or not.

Basic concept 1:
Something symmetric around t=0 in the time domain is strictly real in the frequency domain (as only cosine functions and DC are purely symmetric in the time domain).

Basic concept 2:
The duality principle of the FT implies that if any transform looks intuitive in one direction, the transform in the opposite direction should also.

Transform 1:
A DC level in the time domain corresponds to an impulse at 0 in the frequency domain.

Transform !B:
Getting away from impulses and infinities, but staying with the same picture, a strictly real and really skinny Gaussian pulse in the frequency domain at DC must convert to a really really wide Gaussian envelope in the time domain, also symmetric around t=0.  And vice versa.  Going in between those two (and depending on your graph scaling), a medium width Gaussian in one domain should correspond to a medium width Gaussian in its transform domain.

Possible intuitive transform concept 2 (if you have studied signal modulation concepts):
Adding strictly real and symmetric sidebands to a carrier in the frequency domain corresponds to AM modulation with a cosine sinusoid in the time domain.  The farther away the sidebands from the carrier center, the higher the frequency of the modulation.  The closer the sidebands, the lower the frequency of the modulating sinusoid.

Thus if we add symmetric wide sidebands to a wide Gaussian, that will correspond to a narrow Gaussian that is high frequency AM modulated in the other domain.  If we add only narrow sidebands close in to a narrow Gaussian, that will correspond to a really wide Gaussian with a low frequency AM modulation in the other domain.

Possible intuitive transform concept 3 (if you have done much filter design):
If you have had to design  FIR filters, you might know that the sharper the cut-off and the flatter the pass-band desired in the frequency domain, the longer the FIR filter becomes to meet those requirements, corresponding to a longer and wider impulse response in the time domain.  Or, in the other direction, the sharper-edged and flatter-topped a pulse becomes in the time domain, the wider the bulk of its frequency response spans in the frequency domain.

So, let's take that Gaussian pulse, and add some similar shaped side-bands within the bulk of its main lobe so that its top become flatter. The Gaussian in the other domain will get amplitude (AM) modulated with sinusoidal modulators, mostly of frequencies where the tops of the added sidebands are much higher that the original Gaussian roll-off. (To make it even flatter, perhaps remove the some or all of carrier (original Gaussian) to produce a modulation over 100%, e.g. that crosses zero.)

Let's also chop off the tails of the Gaussian with a brick wall filter and increase the flatness of the pass-band top even more to “square it up” in one domain.  In the other domain, that will cause lots of added ripples extending wider than before, or a “fattening up” of the tails of the original Gaussian.

So as we “square up” a Gaussian in one domain, it gets amplitude modulated plus fatter tails in the other domain.

Thus, a skinny “squared up” Gaussian pulse in one domain might corresponds to a wide slowly AC modulated pulse with fat tails in the other domain,  And a wide but “squared up” Gaussian might correspond to something like a narrow pulse with fast AC modulation, but also with fat tails, in the other domain.

An AC modulated hump with fat tails starts to approach looking more similar to a Sinc waveform than the original Gaussian waveform, and with a width inversely proportional to the width of the "squared up" looking waveform in the other domain.

Intuitive?  Who knows…


A clarinet has no even harmonics. What would produce no odd harmonics? According to this link, the waveforms of clarinets do not have even-numbered components in their harmonic series:

  A closed cylindrical air column will produce resonant standing waves at a fundamental frequency and at odd harmonics... As can be seen from a sample waveform, the even harmonics missing from the tone

Is there an example of a system that would produce no odd harmonics?

Producing "only" even harmonics is equivalent to producing all harmonics.

Say your original signal contains all harmonics:

1f, 2f, 3f, 4f, 5f, ...

Now remove all the odd harmonics:

2f, 4f, 6f, 8f, 10f, ...

The original 1f has been removed, so now you will hear 2f as the fundamental, and 4f will sound like the 2nd harmonic, 6f will sound like the 3rd harmonic, etc.  You're back to 1f, 2f, 3f, 4f, 5f, ... except an octave higher.  I think this might be something that happens in instruments, like "overblowing"?  But I'm not sure.  Anyway it will just produce a harmonic sound an octave higher, but with a different timbre.

Somewhat related: https://dsp.stackexchange.com/a/6178/29


Calculating penalty to SNR for non-matched filter If a matched filter is optimal for detecting a known pulse in an AWGN channel, any other filter must have a degraded SNR. Given a discrete pulse shape and an arbitrary non-matched filter, is it possible to calculate how much worse the SNR will be, relative to the matched filter case?

Let's consider a received signal

$$Y(t)=Ap(t)+N(t)\tag{1}$$

where $A$ is the information symbol (modeled as a random variable), $p(t)$ is the transmit pulse, and $N(t)$ is additive white Gaussian noise (AWGN, modeled as a random process). For the sake of simplicity let's assume that all signals are real-valued (baseband case).

Filtering with an LTI system with impulse response $h(t)$ and sampling at $t=0$ gives

$$Y_0=\int_{-\infty}^{\infty}Y(t-\tau)h(\tau)d\tau\Big |_{t=0}=\int_{-\infty}^{\infty}Y(-\tau)h(\tau)d\tau\tag{2}$$

With $(1)$ you get

$$Y_0=A\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau+\int_{-\infty}^{\infty}N(-\tau)h(\tau)d\tau\tag{3}$$

The power of the desired signal (the first term in $(3)$) is

$$P_S=E[A^2]\left(\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau\right)^2\tag{4}$$

where $E[\cdot]$ is the expected value, and $E[A^2]$ is the power of the information symbol $A$. From the second term in $(3)$, the noise power is

$$P_N=\frac{N_0}{2}\int_{-\infty}^{\infty}h^2(\tau)d\tau\tag{5}$$

where $N_0$ is the one-sided power spectral density of the noise. From $(4)$ and $(5)$ the SNR is

$$SNR=\frac{P_S}{P_N}=\frac{E[A^2]\displaystyle\left(\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau\right)^2}{\displaystyle\frac{N_0}{2}\int_{-\infty}^{\infty}h^2(\tau)d\tau}\tag{6}$$

We know that $(6)$ can be maximized by a matched filter:

$$h_{opt}(t)=Cp(-t)\tag{7}$$

with some constant $C$. From $(6)$ and $(7)$, the optimal (maximum) SNR is given by

$$SNR_{opt}=\frac{E[A^2]}{\frac{N_0}{2}}\tag{8}$$

The actual SNR for an arbitrary receive filter impulse response $h(t)$ can be written in terms of $SNR_{opt}$:

$$SNR=SNR_{opt}\underbrace{\frac{\displaystyle\left(\int_{-\infty}^{\infty}p(-\tau)h(\tau)d\tau\right)^2}{\displaystyle\int_{-\infty}^{\infty}h^2(\tau)d\tau}}_{\le 1}\tag{9}$$

where the last term in $(9)$ (or its inverse) could be called the "SNR penalty".


MATLAB: Decimation image I'm trying to decimate given picture

Vertically.
Horizontally.
Horizontally + Vertically.

I know that I should use low-pass filter* and then downsample** it not to see any aliasing in image.
To this my code is below:

Horizontally filter:

filter_hor    = [0.25,0.5,0.25];
filtered_hor  = imfilter(big_shirt,filter_hor);%% big_shirt is shirt image given above
decimated_hor = filtered_hor(1:2:end,:);%%decimated_hor = downsample(filtered_hor,2);
imshow(decimated_hor)
title('Horizontally filtered')
axis on

Vertically filter:

filter_ver    = [0.25;0.5;0.25];
filtered_ver  = imfilter(big_shirt,filter_ver);
decimated_ver = filtered_ver(:,1:2:end);
imshow(decimated_ver)
title('Vertically filtered')
axis on

Horizontally + Vertically filter:

filtered_2D  = imfilter(decimated_hor,filter_ver);
decimated_2D = filtered_2D(:,1:2:end);
imshow(decimated_2D)
title('2D filtered')

What I get is :

*My low pass filter is: [0.25,0.5,0.25]

**By a factor of 2

My expectation is:

When I use horizontal filter, I just get rid of aliasing on lines of shirt collar.

When I use vertical filter, I just get rid of aliasing in vertically lines.

Vertical+horizontal will be nice picture like what I get.

However in first two cases are not satisfied. What can be the problem. My resolution on computer is 1366X768 and image size is 800X600. When I use imshow, I get:

Warning: Image is too big to fit on screen; displaying at 67%

Can be problem related to this or am I doing something wrong?

After my trials, I realize that problem is directly about Warning: Image is too big to fit on screen; displaying at 67% warning. I use imwrite to my matrix and convert it to jpg file. After looking the image using jpg file , there is no problem.


FIR filter designing in Matlab I am studying digital signal processing and our professor gave us a task of modulation and demodulation and designing an appropriate FIR low-pass filter in Matlab.

I have done the modulation and demodulation but I am stuck on the filter designing. I want to know how the coefficients of filer are determined.

Have a look at the examples for fir1, fir2, or firpm. The first two are straightforward, while firpm is a little more complicated as it derives the minimum filter length needed to match the design criteria, but it's handy to know.

For a simple lowpass filter this should be enough: Feel free to play around with N or look into using firpm to get the minimum order of N.

fc = 100; fs = 1024; N = 256;
b = fir1(N, 2*fc/fs);
freqz(b, 1, 65536, fs)


Baseline correlation for template matching In normal cross-correlation or convolution between a signal and a template(say a small portion of the signal) the template is moved along length of the signal without changing the amplitude offset of the template. But is there any function (in MATLAB or general) that alters the amplitude offset of the template according to the changes in the amplitude of the signal during cross-correlation or convolution?

For example, consider the following ECG signal which is has a large amount of baseline wander due to motion artifacts. The sharp peaks and trough are the heartbeat signals. The beat signal contained in the red box is taken as a template to perform cross-correlation with the full ECG signal.

Due to the large baseline wander, after doing normal cross-correlation, the heartbeat peaks are getting de-emphasized.But, there is strong possibility that the heart-beats will be emphasized if the amplitude offset(offset along y axis) of the template follows the signal amplitude during correlation. In other words, the template must be more or less alighted with the heart-beat in both time and amplitude during convolution or xcross operation. I hope the above illustration adds clarity to the question and the requirement. If there is any alternate method for template matching that will give better results for the above scenario, please do suggest the same.

This would be a cumbersome way to detect heart beats (or the QRS complex), if that is what you are trying to do ultimately.

A little bit about what you are trying to do currently:

Your observations are correct and to these I would like to add that no two heart beats are the same and therefore, strictly speaking, your template will be aligning just with similar heart beats. What if the subject started running in the mean time? That would be quite a different heart beat template to compare against.

If you absolutely must use a template matching technique and apply a hard threshold on the peaks of the cross correlation, here is a set of points that might help you:

Apply a high pass filter in the region of 0-5Hz. This will remove the slow varying components and what you observe as "drift".
Select a number of heart beats, at random, from the filtered signal and create an average "beat". Use this as your template.

A little bit about how is this task usually done:

Picking up the heart beats from an ECG signal usually involves rectification of the signal and sliding window integration with the possible addition of some filtering to reduce the effect of artefacts. The "ideal" result of this is a time series that looks like a square pulse with the positive transition of the pulse almost aligned with the start of the QRS complex and the negative transition of the pulse aligned with the end of the QRS complex. But in practice, you might find that this signal is more "floating" due to noise and artefacts which means that you might also have to use an adaptive threshold technique to improve the detection performance.

For more information please see this link. For much more information, including signals, please see this link.

Hope this helps.


What is obtained from the cross correlation plot? Let’s assume that we have two audio signals, x(t) and y(t) affected by the noise as shown below.
And we would like to cross-correlate these two signals and the cross-correlation plot is shown as below.

In this correlation plot there is a peak value around -11 msec.
I am trying to understand how we interpret this peak in this plot? What does it mean?
Please also explain me what do we obtain from the Fourier transform of the correlation function.

Thanks in advance!

I may answer things that you already know,but i will do my best anyway (i had similar problems once..).
When you do the Correlation of the two sequences you simply shift the one relatively to the other (as you may already know) so each time you shift it you get lets say one "dot" for your plot over there,one single point.This points in a row construct the plot you are getting there.Now,in order to get each one of them you perform one operation,lets say better one calculation which is:the number of equal/same bits minus the number of non equal/same bits,after you perform that you shift again and calculate again....Now when it come to the point that two identical sequences are perfectly aligned,you can understand that the number of "non equal bits" will be zero (cause apparently will be matching perfectly and the result will be a massive number (compared to the rest) that will be YOUR CORRELATION PEAK,meaning that you have reached correlation.As i see there that is something similar to BOC for satellite coms that i have been using in the past,is that case that time   you need to shift the waveform in order to reach to the peak is the time that you need to uses to calculate the distance in your GPS system!
I hope that helps...


Signal in frequency domain with OpenCV dft I am experimenting with cv::dft: a 1HZ sinus signal is generated, and displayed in the frequency domain. But for some reason it hasn't got the maximum component at 1Hz. My code is the following:
       const int FRAME_RATE = 20;             //!< sampling rate in [Hz]
       const int WINDOW_SIZE = 256;

       double len = double(WINDOW_SIZE)/double(FRAME_RATE);   // signal length in seconds
       double Fb = 1./len;                                     // frequency bin in Hz

       // Constructing frequency vector
       std::vector<double> f;
       double freq_step = 0;
       for (int i = 0; i < WINDOW_SIZE; ++i)
       {
          f.push_back(freq_step);
          freq_step += Fb;
       }

   // Create time vector
   std::vector<double> t;
   double time_step = 0;
   for(int i = 0; i<WINDOW_SIZE; ++i)
   {
       t.push_back(time_step);
       time_step += 1./double(FRAME_RATE);
   }

   // Creating sin signal with 1Hz period
   std::vector<double> y;
   for(auto val : t)
   {
       y.push_back(sin(1*FRAME_RATE*val));
   }

   // Compute DFT
   cv::Mat fd;
   cv::dft(y, fd, cv::DFT_REAL_OUTPUT);
   fd = cv::abs(fd);

If I plot the signal in time and frequency domain: plot(t, y); plot(f, fd) the result is the following:

The time signal is good, but the frequency signal has maximum around 6HZ instead of 1HZ.
Where did I take the mistake?

Your sine wave should use 2pi instead of frame rate.  Also, your frequency data should have window size/2 points and go from 0 to frame rate/2, not frame rate.


Demonstrating the effect of aliasing How does the signal look when we don't use the Nyquist rate to remove aliasing from a signal during sampling?

Let's suppose the signal is sinusoidal, with a frequency of 500 Hz and an amplitude of 2.

signal = 2*cos(2*pi*500*t)

If I sample it, (replacing t=nTs , Ts = sampling period and n represent number of samples) and plotting the sampled signals with a different sampling period using the subplot command in MATLAB, how could I identify the aliasing in a sampled signal?

Here is the example code that plotted two signals, one at the Nyquist rate while the other less than the Nyquist rate:

A = 2;
Fmax = 10;
Fs = 2*Fmax;
n = 0:1/Fs:1;
Cont = A*sin(2*pi*(Fmax/Fs)*n);
Cont1 = A*sin(2*pi*(Fmax/18)*n);
subplot(2,1,1)
stem(n,Cont)
hold on
stem(n,Cont1)

and here is the waveform:

I wasn't able to identify the aliasing. How did it affect the signal when Nyquist rate didn't use?

You can't identify aliasing with a simple sinusoid at a specific frequency and that in a way is the whole point about trying to avoid it. You can't know if the sinusoid you are "looking at" is $Q$ Hz or $2Fs-Q$Hz.

A single aliased sinusoidal component looks just like a non-aliased sinusoid. If you want to experience aliasing, you have to attempt it with either a more complex waveform or a sinusoid that is changing in time.

One way to "experience aliasing" is to undersample a chirp in the following way:

Fs = 8000;t=0:(1./Fs):(5-1./Fs);p=2.*pi.*t; %Sampling Frequency, Time Vector, Phase Vector
y1 = chirp(t,0,5,Fs/2.0); %Create a chirp that goes from DC to Fs/2.0
spectrogram(y1); %Have a look at it through spectrogram, please pay attention at the axis labels. This is basically going to be a "line" increasing with time.
soundsc(y1,Fs); %Listen to it...It clearly "goes up" in frequency
y2 = chirp(t,0,5,Fs); %Now create a chirp that goes from DC to Fs
spectrogram(y2); %Have a look at it through spectrogram
soundsc(y2,Fs); %Listen to it...Do you "get" the folding of the spectrum?

In general, you can think of Sampling as Modulation because this is what is effectively happening at the sample and hold element of an ADC converter.

This will enable you to understand more easily concepts like undersampling for example (and applications where it is perfectly OK to sample at lower than the Nyquist rates). But also, you can load a WAV file in MATLAB (with 'wavread') containing some more complex signal and prior to listening to it with 'soundsc', simply multiply it with a "square" wave* (you might want to lookup function 'square') at some frequency lower than the WAV file's $F_s$. This will effectively introduce the key (unwanted) characteristic of aliasing which is this folding of the spectrum. The result is not very pleasant so you might want to keep your speakers volume down.

I hope this helps.

*EDIT: Obviously, "square" returns a square with an amplitude in the interval [-1,1] so prior to multiplying it with your signal it would be better to rescale it as:

aSquareWave = (square(100.*p)+1.0)/2.0 % Where p is the phase vector and here we are producing a square wave at 100Hz (given an Fs of 8kHz as above). aSquareWave's amplitude is now in the interval [0,1]


Is $y[n] = n x[n]$ an LTI system? How can I test if $y[n] = n x[n]$ is an LTI system? And any other system for that matter?

For example, how come $y[n] = \left( \frac{1}{2} \right)^n u[n]$, where $u[n]$ is the unit step, is an LTI system?

Here instead of a general proof method,  let us use the method known as the counter example which relies on the fact that a single counter-example that violates the claim is sufficient (enough) to prove that the claim is wrong.

Assume that the system was LTI,(the claim is that system is LTI). Then consider the following inputs $x_1[n] = \delta[n]$ and $x_2[n] = \delta[n-1]$ with corresponding outputs $y_1[n] = 0 \cdot \delta[n] = 0$ for all $n$ and $y_2[n] = 1.\delta[n-1]=\delta[n-1]$.

Now if the system is actually LTI, then its output should satisfy $y_2[n] = y_1[n-1]$. However one can see that this is not the case, as $$ y_2[n]= \delta[n-1] \neq y_1[n-1] = 0.$$ Hence this counter-example has proved that the claim was wrong; i.e., the system is not LTI. In fact, what we have actually shown is that the system is not time-invariant, but then the fact that it's also not LTI follows as a consequence.


RoC and Stability of a Rectangular Signal If we have a system with an impulse defined as: $$h(t)=u(t)-u(t-2)$$
Then the Laplace Transform of h(t) would be the transfer function: $$H(s)=\frac{1}{s}-\frac{e^{-2s}}{s}, \quad Re(s)>0$$
We also know that a stable system must have a RoC region that passes through the imaginary axis. In our case for $h(t)$, $\textrm{Re}(s)>0$ so it doesn't pass the imaginary axis, indicating that the system is not stable.
However, the system is stable because $$\int_{-\infty}^{\infty}|h(t)| \leq M$$
So my question is how is this possible? Am I missing something in understanding the stability condition in Laplace Transform?

Note that for $s\to 0$, the transfer function $H(s)$ doesn't have a pole because also the numerator approaches zero. In fact, you have
$$\lim_{s\to 0}H(s)=2\tag{1}$$
The ROC is the whole $s$-plane, except for $\textrm{Re}(s)\to -\infty$.


how to stop matlab from rounding my values and keep them at same precision during calculations? here are my IIR coefficients:

CoefIIR=[1.630078934734 1.267893338385 1.112390120011]

and here are the result given to me by matlab commands:

CoefIIR =

1.6301    1.2679    1.1124

How can I keep the precision for data storage and during calculations?

The precision is kept. It is just the display values that are truncated.

Type format long to change the displayed precision.


FIR Filter Design: Window vs Parks McClellan and Least Squares Are there any advantages to use a window approach over Parks-McClellan (further abbreviated here as PMcC) or Least Squares algorithms for FIR filter design of a low pass filter?  Assume with today's computational power that the complexity of the algorithms themselves is not a factor.

This question is not comparing PMcC to Least Squares but specifically if there is any reason to use any window FIR design technique instead of those algorithms, or were windowing techniques to filter design obsoleted by those algorithms and relegated to didactic purposes?

Below is one comparison where I had compared a Hamming window to my favored design approach with Least-Squared, using the same number of taps. I widened the passband in the Least Squared approach to closely match that of the Hamming Window, and in this case it was quite clear that the Least-Squared would outperform (offering significantly more stop band rejection). I have not done this with all windows, which leads me to the question if you could ever out-perform PMcC and least-squares, or if there are other applications for a FIR low pass filter where a windowing approach would be preferred?

I agree that the windowing filter design method is not one of the most important design methods anymore, and it might indeed be the case that it is overrepresented in traditional textbooks, probably due to historical reasons.

However, I think that its use can be justified in certain situations. I do not agree that computational complexity is no issue anymore. This depends on the platform. Sitting at our desktop computer and designing a filter, we indeed don't need to worry about complexity. However, on specific platforms and in situations where the design needs to be done in quasi-realtime, computational complexity is an issue, and a simple suboptimal design technique will be preferred over an optimal technique that is much more complex. As an example, I once worked on a system for beamforming where the filter (beamformer) would need to be re-designed on the fly, and so computational complexity was indeed an issue.

I'm also convinced that in many practical situations we don't need to worry about the difference between the optimal and the suboptimal design. This becomes even more true if we need to use fixed-point arithmetic with quantized coefficients and quantized results of arithmetic operations.

Another issue is the numerical stability of the optimal filter design methods and their implementations. I've come across several cases where the Parks-McClellan algorithm (I should say, the implementation I used) did simply not converge. This will happen if the specification doesn't make much sense, but it can also happen with totally reasonable specs. The same is true for the least squares design method where a system of linear equations needs to be solved, which can become an ill-conditioned problem. Under these circumstances, the windowing method will never let you down.

A remark about your comparison between the window method and the least squares design: I do not think that this comparison shows any general superiority of the least squares method over the windowing method. First, you seem to look at stop band attenuation, which is no design goal for either of the two methods. The windowing method is not optimal in any sense, and the least squares design minimizes the stop band energy, and doesn't care at all about stop band ripple size. What can be seen is that the pass band edge of the window design is larger than the one of the least squares design, whereas the stop band edge is smaller. Consequently, the transition band width of the filter designed by windowing is smaller which will result in higher stop band ripples. The difference in transition band width may be small, but filter properties are very sensitive to this parameter. There is no doubt that the least squares filter outperforms the other filter when it comes to stop band energy, but that's not as easy to see as ripple size. And the question remains if that difference would actually make a difference in a practical application.

Let me show you that such comparisons can often be made to look the way one would like them to look. In the figure below I compare a least squares optimal low pass filter designed with the Matlab/Octave function firls.m (blue) to a low pass filter designed with the window method using a Kaiser window (red).

From the figure, one could even conclude that the filter designed by windowing is slightly better than the least squares optimal filter. This is of course non-sense because we didn't even define "better", and the least squares filter must have a smaller mean squared approximation error. However, you don't see that directly in the figure. Anyway, this is just to support my claim that one must be very careful and clear when doing such comparisons.

In sum, apart from being useful to learn for DSP students for purely didactical reasons, I think that despite the technological advances since the 1970's the use of the windowing method can be justified in certain practical scenarios, and I don't think that that will change very soon.


Linearity and time-shifting of $\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}$ To preface, this is not a homework related question but purely for self-study purposes.

Hi there, I try to calculate $\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}$ by using the properties of Discrete time Fourier transform.

Based on the linearity property, we have:

\begin{align}
\mathcal{F}\{0.8^nu[n]\cos(0.1πn)\}&=\mathcal{F}\{0.8^nu[n]\}\cdot\mathcal{F}\{\cos(0.1πn)\}
\end{align}

I got

\begin{align}
\frac{1} {1-0.8e^{-jw}}\pi( \delta(w+0.1\pi)+\delta(w-0.1\pi))\tag1
\end{align}

But if I applied the time-shifting property:

\begin{align}
\mathcal{F}\{0.8^n\cos(0.1πn)u[n]\}
&= \mathcal{F}\left\{0.8^n\frac{e^{j0.1\pi n}+e^{j0.1\pi n}}{2} u[n]\right\}\\
&= 0.5\left(\frac{1}{1-0.8e^{-j(w+0.1\pi) }}+\frac{1}{1-0.8e^{-j(w-0.1\pi) }}\right) \tag2
\end{align}
I am not quite sure if equation (1) and (2) are equivalent to each other. It looks like we have $\frac{\pi} {1-0.8e^{-jw}}$    at (1) when $w$= $\pm0.1\pi$, but I cannot get the same value at (2).

Why I get the different result here? Am I doing something wrong?

$\displaystyle \hbox{ DTFT}(\alpha x_1 + \beta x_2) = \alpha\cdot \hbox{ DTFT}(x_1) + \beta \cdot\hbox{ DTFT}(x_2)$

Linearity property cannot applied here.


What is the exact meaning of unstable system in DSP? In physical systems I understand what is the meaning of stability or unstability. An operational amplifier for example, if working in positive feedback will either saturate or start osscilating (i.e will not have any stable state). thats clear to me.

But I am unable to understand what exactly we mean when we say an IIR filter (or any other digital system), for example, can become very unstable.

What exactly inside the Digital Signal Proccessor  happens, what happens to the output physically?
What exactly we mean by unstable system in this context?

Unstable typically means and unbounded output for a bounded input. In other words the output of, say , a filter can get infinitely large although the input is perfectly okay and of "normal" size.
A simple example would be the difference equation $y[n] =  x[n] + y[n-1]$. If we calculate the step response, i.e. $x[n] = u[n]$, we get y[0] = 1, y[1] = 2, y[2] = 3 ...
The output grows infinitely even though the input is a perfectly well behaved signal, bounded by 1.


detect impact from digital audio I've started a hobby project which is to detect and count the number of impacts (for example someone clapping, bouncing a ball or hitting a boxing bag).

I've been reading and playing around (in java) with a few variations of total amplitude, FFT, IIR etc but with so much information and concepts to grasp (I'm normally a web/database developer) I can't figure out which is the right path to follow.

practically speaking I know i'll need a high frequency sample rate in order to distinguish between peaks and i think i'll need some form of filtered amplitude to distinguish between the crisp leading edge of the impact and the rumble noise that may follow.

I'm hoping someone here can point me in the right direction.

Thanks for your time!

If your environment's noise level allows it, you can use a quiet no-brainer solution. Just make sure that the impacts have greater peak amplitude than the noise floor in every case.

Specify a threshold level based on previous experiments. This level has to be great enough for all the impacts you generate and you want to detect. If your signal crosses the level, you can assume, that an impact has happened.
If a detection has occurred, you can start a timer, that prevents unwanted re-detection for the same impact because you can't expect that your impact will be below the threshold for the next sound sample. For this reason, the timer has to run longer than the maximum detectable length of any impacts. Until the timer runs, you disable the threshold detection.
After the timer has fired, you can enable the threshold detection again.

This method requires previous experiment and knowledge about the detectable impacts. You have to determine the minimum sampling rate that allows you to catch at least a few sample of the impact. Standard audio quality sampling rate would be good enough.

You can of course apply signal conditioning by filtering the signal with a digital filter, but this necessity is highly depends on your environment and noise floor.

Note that this is possibly the most no-brainer solution. There are many more sophisticated methods for this particular problem.


Why are polynomial filters for gradients in the opposite sign? I was looking at the filters, and the ones for 1st derivatives, gradient http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter#Tables_of_selected_convolution_coefficients
by testing the 9-points one:

D=[-4;-3;-2;-1;0;1;2;3;4]/60;
y=sin((0:0.2:6)');
yd = conv(D,y);
plot([repmat(0,128,1) y yd(5:end-4)])

we see that actually [4;3;2;1;0;-1;-2;-3;-4]/60 would give the right signs, so why is it like that?

If you look at the formula above the table you linked to, you see that the output values are computed like

$$Y^{\prime}_j=\sum_{n=-4}^{4}D_ny_{j+n}$$

In Matlab this could simply be implemented by a vector product (for one output value $Y^{\prime}_j$, where the vector $y$ contains the input values around the current index $j$):

yd = D'*y;

On the other hand, convolution of these sequences is defined by

$$Y^{\prime}_j=\sum_{n=-4}^{4}D_ny_{j-n}$$

(note the '-' sign in the index of $y_j$!). So if you want to use the 'conv()' function in Matlab, you need to flip the coefficients in the vector $D$.


How to convert a low-pass Gaussian filter to a high pass filter I have the following code for a applying a Gaussian filter to an image. From what I understand this is a low pass filter. Now I have to convert this into a high-pass filter, and from what we were told from the instructions, one difference between  Gaussian low and high pass filters is that for a high-pass, the sum of the elements in the filter kernel should sum up to zero as opposed to one like for the-low pass filter. How should I go about doing this?

%%horizontal and vertical coordinates
coordinates = -floor(N/2) : floor(N/2);
[X Y] = meshgrid(coordinates, coordinates);
N = 12;
sigma=2;

gfilter = exp(-(X.^2 + Y.^2)/(2*sigma.^2));
gfilter = gfilter/sum(gfilter(:));

I = mean(double(imread('image.png')),3);
GaussBlur = conv2(I,gfilter);
%%plot images side by side
subplot(1,2,1);
imshow(I, []);
subplot(1,2,2);
imshow(GaussBlur, []);

The simplest option is to remove the low-pass signal from the data. Just subtract the outcome of your code from the original image. This amounts to computing:

$$I_\textrm{high} = \delta \ast I - g\ast I  = (\delta -g) \ast I\,.$$

The term $\delta \ast I$ corresponds to the convolution of a Dirac delta to the image, which is an invariant operation.
As the sum of the coefficients of the $\delta$ filter is one, like $g$, the resulting filter $\delta -g$ has a zero sum.
Warning though: a filter whose coefficients sum to one is not necessarily an high-pass filter. It is a low-cut, and can be for instance a band-pass.

A more convincing method consists in using the modulation property:  you can multiply the coefficients of the low-pass filter pointwise by alternating $1$ and $-1$ (alternating signs are discussed here: why the sum of masks should be one?). The following code does it in 1D:

b = [1 4 6 4 1];b=b/sum(b) ;a=1;
c = b.* (-1).^(1:length(b));
[hb,w] = freqz(b,a,1024);
[hc,w] = freqz(c,a,1024);
figure(1);plot(w,abs([hb,(hc)]));axis tight

The adaptation in 2D follows:

N =13;
%%horizontal and vertical coordinates
coordinates = -floor(N/2) : floor(N/2);
[X Y] = meshgrid(coordinates, coordinates);
sigma=2;

gfilter = exp(-(X.^2 + Y.^2)/(2*sigma.^2));
gfilter = gfilter/sum(gfilter(:));

gridInv = (-1).^(X+Y);
hfilter = gfilter.*gridInv;

subplot(2,2,1);
mesh(gfilter);
subplot(2,2,2);
mesh(hfilter);

subplot(2,2,3);
imagesc(abs(fftshift(fft2(gfilter))));
subplot(2,2,4);
imagesc(abs(fftshift(fft2(hfilter))));


The difference between DFT and DFS In the literature, I've found that DFS and DFT are one and the same. If they are one and the same why to use two different names for them? If there is really a difference what is it and what is the significance of discrete Fourier series?

I think part of the problem is an awkward and inconsistent naming convention. There are 4 flavors of Fourier Transforms depending on which domain is continuous or discrete (which maps to being aperiodic or perodic in the other domain). So we have

         Name                  Time                  Frequency
Fourier Transform    continous/aperiodic     continous/aperiodic
Fourier Series       continous/periodic      discrete/aperiodic
Discrete Time FT     discrete/aperiodic      continous/periodic
DFT or DFS           discrete/periodic       discrete/periodic

A better naming would have been

         Name                  Time                  Frequency
Fourier Transform    continous/aperiodic     continous/aperiodic
Fourier Series       continous/periodic      discrete/aperiodic
Discrete FT          discrete/aperiodic      continous/periodic
Discrete FS          discrete/periodic       discrete/periodic

so that discrete refers to "discrete in time and periodic in frequency" and "series" refers "discrete in frequency and periodic in time". In other words "series" means sums and "transform" means integrals. Discrete mean sums and continuous means integrals.


Detecting scene transitions in a video How would you reliably detect scene transitions in a video? The simplest case wouldn't be too hard, but sometimes it might be a fast moving scene, or there might be an effect like lightning that would throw it off.

Scenes generally exhibit a fade-to-black transition. You could capture those frames with simple image processing tools, which are found in many libraries, such as OpenCV.

If you want to rely on the change of the mis-an-scene to be robust against sudden changes, then of cours algorithms for detecting temporal differences is more appropriate. For that, you might want to consider:

http://vis.uky.edu/~cheung/courses/ee639_fall04/readings/spie99.pdf
http://www.mitsubishi-electric-itce.fr/uk-rce/pubdocs/VIL07-D071.pdf

If you want a simple algorithm do as described here:
https://stackoverflow.com/questions/4801053/video-scene-detection-implementation

For harder cases, you might even try GMM-like background models to test if the transition was a temporarily quick shot or a complete change.


Sample and hold - LTI filter I've found a material explaining the sample and hold operation pretty well, but there's one thing I can't understand.

http://www.springer.com/cda/content/document/cda_downloaddocument/9781848001183-c1.pdf?SGWID=0-0-45-484012-p173780368

Next, we would like to take this ideally sampled signal and hold the values
between samples. As we have a train of impulses with the correct areas, we
need a “block” that takes an impulse with area A, transforms it into a rectangular
pulse of height A that starts at the time at which the delta function is
input to the block, and persists for exactly Ts seconds. A little bit of thought
shows that what we need is a linear, time-invariant (LTI) filter whose impulse
response, h(t), is 1 between t = 0 and t = Ts and is zero elsewhere.

How do we know this is the impulse response that will give us a pulse of HEIGHT A, not area A? Could anyone explain this "a little bit of thought" part? Thanks in advance.

An impulse of "area" $A$ means $A\delta(t)$. The nomenclature arises because
$\displaystyle \int_{-\infty}^\infty A\delta(t) \,\mathrm dt = A,$
that is, the "area under the curve" is $A$.

Since the system is linear and
time-invariant, and the response to input $\delta(t)$ is a pulse of height $1$
and duration from $0$ to $T_s$, the response to $A\delta(t)$ is a pulse of
height $A$ and duration from $0$ to $T_s$.  Its "area" is $AT_s$, that is, $A$
multiplied by the area $T_s$ of the response to the unit impulse.


Inverse DTFT of phase shifted complex exponential I have been working on this problem for a few days now and I think this is the closest I have gotten.
I am getting an Answer of zero and I would like to know if that is correct and if someone could check my work please? Thanks for your time, here is what I have so far (also Im quite sure evaluating at n = -10 isnt correct, could someone please give me the correct process?)
Not looking for an answer, just how to start/continue to process. Thanks
I put all my steps in the attatched picture. (sorry, Im new to Stack Exchange and couldnt figure out how to format correctly)
Thanks for your time :)
-Dom

HINTS:

$e^{j\pi/2}=j$, so $X(e^{j\omega})=j\omega e^{-j10\omega}$
The term $e^{-j10\omega}$ corresponds to a delay of $10$ samples in the time domain, so you can just ignore it for now, and then add a delay of $10$ samples. So you compute the IDTFT of $j\omega$, and then replace the variable $n$ of the result by $n-10$.
The IDTFT of $j\omega$ is $$\frac{1}{2\pi}\int_{-\pi}^{\pi}j\omega e^{jn\omega}d\omega$$ That integral is most easily computed using integration by parts.

Note that $X(e^{j\omega})$ is the frequency response of an ideal discrete-time differentiator with a delay of $10$ samples.


Synchronizing data sent over sound beeps with Goertzel Detection, C# I'm trying to send some text over sound. For example, I try to send the words "Hello John". I do it this way: Every char is converted into its binary ascii code. Then I construct a wav file with no gaps by NAudio library in this manner:

    StartDelimiter = 12KHz,
    1's = 13kHz,
    0's = 14kHz(for the message),
    End Delimiter = 10kHz.
    The durations of whole beeps are similar: 100ms.

I'm only at the beginning of my work and I Record the series of beeps and after recording I analyze the record(in the future I want it to be near Real Time) by:

Find the start delimiter's index with Goertzel filter(max power along the recorded buffer).
Skipping a known duration(100ms) in a buffer to decode the 1 OR 0 by using another 2 Goertzel filters. If the End delimiter is found I break the search and show the decoded message to a box.

My problem that I've encountered so are:

After sending more than 5 chars I get for example: "Hello Pohny". The first chars are fine but then it start to rubbish.
Some noise defects the series.

Are there some strategies which I can use to improve the skills of my the above algorithm?
Thanks
I work with C# on WPF

You can't simply skip ahead some amount of time, and won't be able to skip ahead when decoding in real time.

You need to simply decode as the values come in, and switch your output depending on what comes in.  Change states only on input changes.
So, like this:

Feed input signal into the bank of goertzel filters,  the detected tone is the filter output that goes over a defined limit.
When detected tone changes to 12kHz, set message start flag
When detected tone changes to 13khz, send output bit 1
When detected tone changes to 14kHz, send output bit 0
When detected tone changes to 10kHz, end of message.
Repeat from step 1 until end of message.

Working that way, you should stay synchronized and have fewer errors.  You should also be able to shorten the beep lengths, though you may need to work with the length of the Goertzel filters to make it work with shorter beeps.

You will find that you will have to make a compromise between the length of the beeps, the length of the goertzel filter, and the separation between the individual tones in order to get fast, reliable communication.

To get faster communication, you will need shorter filters which in turn require larger separation between the tones to get reliable detection.

What you are doing is called frequency shift keying (FSK,) and has been used to transmit data over audio paths for a long time.

Noise will always be a problem.  There are two ways to fight it, and you might have to use them both.
First, you make your goertzel filters as long as possible, but still short enough to react within your desired bit rate.
Second, most systems encode the data with checksums.  Checksum codes can be used that allow you to recover from bit errors - they don't just say "bad data," they can also say "bad data, bit X should be value Z."  There are of course limits.  The better the checksum and error correction, the more extra information you have to send, so the slower your effective data rate.

A third thing you could do is to lower the data rate and extend the length of the goertzel filters, but you are already running a very low rate.


Low pass FIR filter for EEG signal I'm trying to do Analysis of EEG Signal for the Detection  of epileptic seizure so in the pre-processing of the EEG signal I use a low pass FIR filter

d1=designfilt('lowpassfir','PassbandFrequency',0.15,'StopbandFrequency',0.2,'PassbandRipple',3,'StopbandAttenuation',60,'DesignMethod','equiripple');
a = filtfilt(d1,deriv1);

I want to know if the values above for PassbandFrequency and StopbandFrequency are right.

Assuming MATLAB functions to be used, passband frequency 0.15 means $\omega_p = 0.15 \pi$ whereas stopband frequency 0.2 means $\omega_s = 0.2 \pi$ and since $\omega_s > \omega_p$ these settings qualify for a proper lowpass filter.

Also passband ripple being $3dB$ and stopband attenuation being $60 dB$ are valid choices (although more exact specifications requires the knowledge application requirements)

The equiripple option designs a filter in which the designed filter magnitude error is distributed accross the pass and/or stop bands in equal amounts so as to reduce the filter length to a minimum while still meeting the maximum of the allowed peak error criteria.


OpenCV camshift algorithm I got a lab report to make in an image processing course in my university, and I was asked a question about the OpenCV camshift algorithm (for video tracking):

the camshift algorithm works on the hue component of the HSV color space.  I need to explain why working on this component should give better results than working in the RGB color space.

I must admit that I'm not familiar with OpenCV at all, and basically I'm leaning all the color model theory specifically for this course, so please consider that...thanks.

According to SR who posted here http://answers.opencv.org/question/18311/why-use-hue-channel-of-hsv-for-mean-shift-tracking/ the H channel in HSV

  "... is not influenced by the image brightness and saturation. You get the "pure" color disregarding the actual lighting (which might change over time)."


Most natural interpolation to simulate increasing distance? Background: I'm generating multiple datasets with varying scales of resolution to simulate taking pictures from varying distances and then comparing image classification performances. For example original images at 32x32, then downsampling to 16x16 and measuring the decrease in classification performance.

I'm using Matlab's imresize function to accomplish this. The docs say by default that it uses bicubic interpolation. Now I'm not sure if my following question even makes sense, but what I'm wondering is, is there a type of interpolation method that best mimics the decrease in resolution of an object that would happen naturally if you moved a camera farther away from said object? There's many different choices for the interpolation method parameter, and I'm wondering if I should choose a different one. Thanks!

One typically interpolates going up, 16 $\times$16 to 32 $\times$ 32.  Going down is low pass filtering (anti aliasing) and decimation.  As you have acknowledged, this isn’t the same as the effect of distance.  As long you understand that your original image was captured with a camera with a specific MTF, and sensor resolution, and are going to neglect scattering (which would require something like MODTRAN),  your original image has been acquired with some angular resolution, $\theta_0$, your image has a linear spatial resolution $\Delta_0$ such that $$
\Delta_0 = r_0 \theta_0
$$ where $r_0$ is the distance that the image was acquired at.  At some greater distance $r_g$, $\theta_0$ can be considered the same so,
$$ \Delta_g= \theta_0 r_g. $$
Given some coordinate system associated with the original system, you could use something like a Gaussian blur with a standard deviation related to $\Delta_g$, and resample. This doesn’t really help with a large depth of field.


Is zero forcing Optimization technique or Equalization technique? I got some problem in the equalization techniques. I am new to communication simulations but I have done some data analysis labs in my college. I read some where on Internet that in a research paper that
zero forcing is Optimization technique but on Wikipedia it says it is Equalization technique?

I am confused in equalization and optimization techniques in communication. I am doing a OFDM simulation project in MATLAB and I am stuck in terms of these. I have read on Google but they both terms are kind confusing

Can some body please explain the different equalization techniques and the use of zero forcing?
Why to use zero forcing and what are basic possible equalization techniques that can be used with zero forcing.

Zero-forcing is not an optimiziation technique and neither an equalization technique. It is simply one of several possible criteria to design (or optimize) an equalizer. A zero-forcing equalizer is designed to completely remove intersymbol interference, but does not take noise into account. Another, more practical criterion is the MSE (mean-square error) criterion, which results in an equalizer that minimizes intersymbol interference and reduces noise power at the same time.

I would suggest you consult one of the many good books on digital communication to learn more about this quite complex topic. Here are just two which I like:

Digital Communication by Lee and Messerschmitt
Digital Communications by Proakis


In-cylinder combustion pressure I have pressure data with respect to crank angle and I want to convert it into Sound pressure level with respect to frequency. Please suggest any known procedure .Thank you

If you assume that the engine speed $N$ is constant, given in RPM (Revolution per minute), with an angular sampling of $\alpha$, you should convert it into a sampling interval:

$$ \Delta t = \frac{\alpha}{6N}\,.$$

The sampling frequency in Hertz is the inverse.

If the pressure data $P$ is given in bars, you should convert it to the SI unit in Pascals (about $10^{−5}$ bars) to get an absolute sound pressure. The Sound Pressure Level (SPL) is a logarithmic ratio of the absolute sound pressure and a reference level, depending on the application.

Be aware that crank angle sampling is often far for uniform, and at certain speeds the angular sampling is too high to capture phenomena without aliasing.


MATLAB: Implementing Least Squares Estimator for a Given Model The formula to estimate $\mathbf{h}$ is then
 $$\hat{\mathbf{h}} = (X^T X)^{-1} X^T \vec{y}\tag{2}$$

I think this can be implemented in Matlab using hat_h = pinv(X)*X*y

Question 1 : What is the lag of the model?

Question 2: I don't know how to create the design matrix X in order to estimate h. Can somebody please provide the complete implementation? Thank you

The equation you're trying to solve is
$$
\mathbf{y}=\mathbf{X}\mathbf{h},
$$
where $\mathbf{h}$ is your unknown. The matrix $\mathbf{X}$ is going to have a time-shifted structure that reflects the convolution operator. If we assume that the $\mathbf{y}$ vector starts with y(3) i.e. ignores the first two zeroed out elements of y,  then the corresponding $\mathbf{X}$ matrix is given by:

$$
\mathbf{X}=\left[\begin{matrix}
x(3) & x(2) & x(1)\\
x(4) & x(3) & x(2) \\
x(5) & x(4) & x(3) \\
... & ... & ...
\end{matrix}  \right]
$$

You can add as many rows as you are have observations to support them. You can then solve for $\mathbf{h}$ by
$$
\mathbf{h}_{est}=\text{pinv}(\mathbf{X})\mathbf{y}
$$


sampling a signal below nyquist rate This question is for people working in multi rate signal processing. Using compressive sensing techniques, we can still recover a signal by sampling below nyquist rate(by decimating the signal). How is that possible, can you explain?

There are some answers available on the web , but I want to hear it from SE group.

phani tej

No, recovering an "arbitrary" signal from its samples taken below its Nyquist rate (assuming it is bandlimited) is not possible in general.

Compressive sensing is based on the assumption that if you know a-priori that a signal will be sparse in some of its support domains, then those bands will not be carrying information about the signal and can be discarted in principle, therefore, reducing the pre-computed Nyquist sampling rate, which was computed based solely on the bandlimitedness (Shannon-Nyquist sampling theorem) of the signal.

The trick here is on the "a-priori" knowledge which also contains information, and this must be taken into account when computing the total information present in a signal. This a-priori knowledge can be the result of imposed structures on the signal generation mechanism, or well known models of the signal source.

If you can impose no constraints on the signal other than its bandlimit, you have to sample it above its Nyquist rate, so that exact reconstruction (using unimplementable ideal lowpass reconstruction filters) is possible at least in principle.

The topic of sampling is bound by information theory. There is no way to walk around it. Why would you ever sample a DC signal, or a single sine wave, or a periodic signal? It is based what is known and what is not.


Lower frequencies is not visible in FFT plot when adding into the original signal I plotted a simple Cosine wave and After taking its FFT the frequency can be seen easily but after modifying the original signal to add more lower and higher frequencies, the higher ones are visible but the lower ones is not showing in the FFT plot,

Original Signal,

Amplitude = 4
x=A*cos(2*pi*f*n);

Frequency of Wave: 100Hz

Amplitude = 4

Sampling Frequency = 10000Hz

The 100Hz frequency (original signal frequency) is supposed to be plotted at point 1 of FFT, which is fine,

When adding another higher frequency into the original signal,

x= 1*cos(2*pi*500*n) + x ;

500HZ frequency is added into the original signal and supposed to be plotted at point 5 of the FFt.

Results in,

Which is fine, but when i am trying to add lower frequencies than higher they are not visible.

Adding 50Hz frequency to the original 100Hz Cosine wave,

what i think is 50Hz should be plotted at the 0.5 of the x-axis as it is half of the original frequency.

x= 1*cos(2*pi*50*n) + x ;

Results in,

Where is the 50Hz frequency ?

What you're assuming about how FFT works is actually correct for Fourier transform, but – despite its name – FFT is not such a continuous all-$\mathbb{R}$ transformation. It is actually just a DFT, discrete Fourier transform, i.e. both the frequency and time spaces are quantised. In your example, you selected the quantisation in frequency space to be $100\:\mathrm{Hz}$, so there exists no $50\: \mathrm{Hz}$ bin.


How to construct a preemptive digital filter to neutralize an unwanted filter I'm trying to think of how to construct a filter based on the following scenario. I have some time trace with a certain known power spectral density (I can verify this with a periodogram of the time trace). The time trace then goes into a machine that realizes it, which does so with some level of imperfection. Treating the device as a black box, it is essentially a filter that is applied to the time trace, and produces a slightly adjusted power spectral density. This power spectral density I can then measure with a spectrum analyzer.

So I'm wondering, can I apply a preemptive filter to my time trace, before it goes into the black box, so that I get my desired power spectral density? I suppose I should have all the information I need as I know the power spectral density before and after the black box, which should somehow allow me to construct a type of inverse filter. How exactly I am not sure of however, and that is my main question.

I should also note that I tried a much less clean method. I took the power spectral density after the black box, and simply mirrored it, and used this as a preemptive filter. This actually seems to work quite well in the sense that the resulting PSD has the correct shape (it was sloping downwards before while it should be straight, which this method achieves), but the amplitudes are no longer correct and rescaling the data is problematic. So I figured a more analytic approach would be required.

If $S_x(\omega)$ is the desired input power spectral density (PSD), and $S_y(\omega)$ is the PSD at the output of the black box, which has the (unknown) frequency response $H(\omega)$, then the following relation holds:

$$S_y(\omega)=S_x(\omega)|H(\omega)|^2\tag{1}$$

Now you want to cascade the black box with a filter to recover the original PSD $S_x(\omega)$ from the distorted PSD $S_y(\omega)$:

$$S_y(\omega)|G(\omega)|^2=S_x(\omega)\tag{2}$$

From $(2)$ it immediately follows that

$$|G(\omega)|=\sqrt{\frac{S_x(\omega)}{S_y(\omega)}},\quad S_y(\omega)\neq 0\tag{3}$$

You need to design a filter with magnitude response $|G(\omega)|$ given by $(3)$. Note that for this application, the filter's phase response is irrelevant.


QAM constellation slicer with unknown carrier phase As shown in the image, I have a 16QAM constellation that is misaligned due to a phase rotation. In this case you can see that rotation amount is approximately $\theta = \pi/4$, but this won't be the case in general. For real-world data the phase could be a slowly varying function of time, $\theta(t)$, so that it is not enough to apply some fixed correction factor.

I am aware of differential mapping schemes that solve the phase ambiguity problem due to the constellation having $\pi/2$ symmetry, but it seems $\theta$ must still be known to perform the slicing.

One suggested solution was to try to map the received constellation point to the nearest QAM constellation point and feed a phase-locked-loop with the result, but it is not clear how this would perform when $\theta$ varies over time.

What techniques exist to recover the symbols? I have already tried various carrier recovery schemes based on feedback loops, with no success, and am interested in decision directed approaches that may avoid having to find the phase.

What you need is carrier phase synchronization. This is a complicated topic with many different approaches. The approach that you'll choose could depend on things like:

Data-aided versus blind: Does the underlying sequence contain any known data (e.g. a training or sync sequence of some kind) that you can use to divine the phase offset? Or, do you have to synchronize with no knowledge of the modulating symbols?

Blind approaches are more general, but you can get better phase estimation performance if you're aided by the data. Also, blind approaches often introduce ambiguities into the recovered phase (i.e. there are multiple phase offset solutions that equally satisfy the blind optimization criterion).

I haven't come across any blind solutions that work equally well across all modulation types; for instance, most blind phase estimators are best applied to PSK signals. They can be made to work, suboptimally, on QAM. Various blind carrier phase estimators for QAM are documented in the literature, but I don't have any great recommendations for QAM-specific estimators. The most commonly-used blind estimator for PSK-type signals is the power-law detector. Here is an example of a paper that talks about its application to QAM.
Coherent versus differential: As you noted, one way to get around having to synchronize to the carrier phase directly is to use differential modulation. In that case, information is carried by the difference in phase between consecutive symbols. Since the carrier phase is likely to be approximately constant across a two-symbol time period, the carrier component cancels out. This makes synchronization easier, but there is a symbol error rate performance loss on the order of 1-2 dB for differentially-coherent modulations versus fully-coherent operation.
Feedforward versus feedback: Are you processing a continuous, indefinitely-long stream of symbols, or do you have a finitely-sized batch? A feedback approach like a phase-locked loop might be appropriate for the former, while feedforward techniques that estimate the bulk phase offset of a block of symbols at a time are best for the latter (feedback techniques have some acquisition period during which you won't get good output; if you only have short blocks of data at a time, this can be a problem).

If you're looking for a book reference, my go-to is Synchronization Techniques for Digital Receivers by Mengali. It's expensive and hard to find a copy, but I find it very thorough.


Can we do cyclic prefix with zeros in OFDM rather than with samples from end of the ofdm symbol? considering an OFDM symbol of x(0)x(1)x(2)x(3)x(4)....x(n-2)x(n-1)x(n).
To perform cyclic prefix we prefix some samples from end of this OFDM symbol to the beginning of the OFDM symbol,  like

x(n-3)x(n-2)x(n-1)x(n)x(0)x(1)x(2)x(3)........ x(n-4)x(n-3)x(n-2)x(n-1)x(n).

Instead of taking x(n-3)x(n-2)x(n-1)x(n), if we prefix with zeros what will happen?

There are two major approaches for cyclic extension in OFDM systems - CP (cyclic prefix) and ZP (zero padding, also called Trailing Zeros, TZ). Generally they show the same performance, I mean in AWGN or Fading channel. CP method is the simplest one so it is preferred in the most cases. ZP approach leads to slightly less transmission power level for obvious reason. Doppler performance is the same. The main idea of using ZP as I understood is the existence of equalization techniques for combating spectral nulls that is the major problem of OFDM equalizers. I can't explain this technique because I haven't modeled or implemented it, but you can search some articles in the Web. In one of the article it stands that for CP the simple equalizer exists but symbol recovery isn't guaranteed. But for ZP symbol recovery is guaranteed with complexity of equalizer increasing. I think it is the ckear idea. But I suppose it leads to some matrix subspace algorithms e.g. SVD, so it can be too difficult to implement. I've heard ZP is popular in hydroacoustic OFDM techniques of equalization so you can also refer to this topic.

Hope it helps.