Capturing the sampling rate from a basketball being dribbled I'm taking an online image and signals class and the structure is absolutely horrible.

I'm wondering if you can provide me hints, or tips to solving the below question.  I'm not looking for an answer, because I'd actually like to learn something!  The structure of the class is just not conducive to learning as they just point us to online resources without any real "teaching."

  Consider a basketball being dribbled. If the height of the basketball can be described by a sine wave of maximum height 2h, average height h and minimum height 0, and the ball hits the ground once per second, how fast would a video camera have to sample the dribbling to extract its frequency? What happens to the frequency estimate if the sampling rate is too low?

So just thinking out loud... the ball is bouncing at a rate of $1\frac{bounce}{sec}$, which from my reading tells me the frequency is just the inverse of this, so $1\frac{sec}{bounce}$.

I don't understand, though, how to extract how fast a video camera would need to sample this to extract its frequency as it seems like I already have the frequency.

The information you are given is $f_\textrm{bounce} = 1\,\textrm{Hz}$ - all the technical details about the sine wave and its amplitude is irrelevant. Please note that frequency is always measured in something happening per second, i.e. the unit is always $\frac{1}{\textrm{s}}$.

Since you do not want to have an answer, I try to give you a hint:

Imagine you had a camera that would capture an image once a second. Assume that the moment the image is taken is always when the basketball exactly in the middle between the floor and the player's hand. What can you say about the motion of the basketball after having watched this video?
Imagine your camera would capture the scene every $0.5\,\textrm{s}$, i.e. the frame rate would be $2f_\textrm{bounce} = 2\,\textrm{Hz}$, again exactly when the basketball is between the player's hand and the floor. What would be the difference between this video and the first video? What would you learn about the motion of the basketball by watching this video?

If you want to get more information, you could search for the Nyquist–Shannon sampling theorem.


Derivative of signal with missing samples I have software that tracks an object moving (in the x-dimension only) across a video shot from a stationary camera. I need to find the velocity and acceleration of the object as functions of time. This calculation can be done offline, so there are no real-time constraints. The motion of the object is also pretty regular, decelerating approximately constantly throughout the duration of motion.

Normally I would design an FIR differentiator to do the job, however the object in the video becomes occluded for up to a few frames at a time on a regular basis throughout it's motion, so I don't have a point for it for several frames. In other words, my data looks like this:

frame_index = [1 2 3 6 7 8 9 13 14 15 16 17 18 21 22 ... ]
x_position  = [0 0.2 0.39 1.14 1.22 ... ]

I've thought about interpolating and then using an FIR differentiator but am unsure what a good scheme for interpolating irregular data like this would be. I've also thought about using a formulation based on Lagrange interpolating polynomials to calculate derivatives directly, but do not know enough about them to understand trade-offs. My goal is accuracy.

So, what would be the best way to go about calculating approximations of the first and second derivatives of the data?

See Smooth noise-robust differentiators. While Pavel doesn't address your specific problem there, he and I corresponded about a similar "missing samples" problem earlier this year, and he has some good ideas that you can try. In particular, he recommends polynomial approximation in the region of the missing samples. In addition I'll mention that if you use Neville's Algorithm you can guarantee that the polynomial passes through your "good" sample values. You can find Pavel's contact info in various places on his Website.

Greg Berchin


Power Spectral Density units I am trying to compare power spectral density values (obtained from dsp.SpectrumAnalyzer) with PSD values obtained from hardware which has a VSA.
MATLAB's documentation says that " The power spectral density is the magnitude squared of the spectrum normalized to a bandwidth of 1 hertz." and expressed in dBm/Hz

whereas the hardware expresses PSD in units  dBm/RBW per 100 kHz    So,I'm assuming that it is normalized w.r.t 100Khz bandwidth
This is an example snippet of hardware's vsa( not my actual output)

Do I need to do any kind of conversion or something in terms of PSD ?
For ex:  dBm-->mW--> multiply by 100e3  something like that to properly compare? Kind of confused with these units.
MATLAB's output is off by approx 20 dBm

Here is the output comparison

EDIT>>>

0)  RBW is set to 100KHz on hardware and in MATLAB too
 1) I'm guessing there is a factor 100 which I need to multiply to PSD in mW to make it similar, but not sure why I need to do that

2)The reference load for matlab is set to 1 ohm by default, whereas  hardware (litepoint ) reference impedance is 50 ohms

3) Updating basic specAnalyzer code in matlab if that helps in understanding

nSpectralAverages = 3;
specAnalyzer = dsp.SpectrumAnalyzer('SpectrumType','Power density',...
                              'SampleRate',fs);
specAnalyzer.FrequencySpan = 'Span and center frequency';
specAnalyzer.Span = psdSpanLpVSA;
specAnalyzer.SpectralAverages = nSpectralAverages;
% Resolution Bandwidth
specAnalyzer.RBWSource = 'Property';
specAnalyzer.RBW = resolutionBandwidthHz;

specAnalyzer.AveragingMethod = 'Exponential';
specAnalyzer.ForgettingFactor = 1;
% Run the spec An and get collective result
specAnalyzer(WiFi_iqData_bw20MHz);
specData = specAnalyzer.getFramework.Visual;

I'm using exponential averaging method as its pretty clean and similar to litepoint's output

4) My other question is: Can I even compare outputs directly like this? Or should I just add the difference in spectrum peaks and then compare (after adding 50 Ohms to the MATLAB's output)

EDIT 5)
 I have converted dB/Hz (which matlab provides) to dBm/100Khz for comparison and also calculated power values for 50 ohm

Pwr in dBm/Hz --> convert to mW --> multiply by 10000/50 (100e3 for 100KHz and 50 for 50 ohms conversion)  --> then reconvert to dBm

Here is the output

Any comments /suggestions

Thanks

The division should be by 100000, but I doubt that hardware is giving the result in actual milliwatts, maybe a normalization is required in the power being reported by the hardware as well followed by a subsequent normalization in bandwidth. This would require knowing how the hardware is generating this number,

considerations like, what is amount of duration over which hardware measures this spectrum, does it average out over those durations or does it provide cumulative value which needs to be normalized.


Why use Wavelet Denosing instead of Band pass Filtering What is the important advantage of adopting a wavelet based denosing scheme such as SURE, SUREshrink, etc. than using approaches involving low pass, high pass, band pass, band reject filter for signal denosing?

More precisely, what is the advantage of using wavelet denosing on EEG signals than other methods.

Thank you!

If you use frequency selective filters for denoising you implicitly assume that the desired signal and the noise occupy different frequency bands. You only remove noise in frequency regions where there is no desired signal, whereas in the passbands of the filter the noise is left unchanged. For some applications the assumption that noise and signal occupy different frequency regions is justified, for some it isn't. Wavelet based denoising methods compose a noisy signal into different scales and attempt to remove the noise while preserving the signal, regardless of the frequency content.

In practice you have to look very carefully at the characteristics of the noise. It is always advisable to apply a frequency selective filter to suppress all out-of-band noise. For narrow band noise sources a notch filter is often sufficient. In cases where you can get a reference signal of the noise (i.e. an independent recording of the actual noise source) the best strategy is using an adaptive noise cancellation algorithm. Only if these more basic methods do not lead to a satisfactory result, and if signal and noise largely overlap in frequency and time, you should look into more sophisticated methods such as wavelet denoising, or blind source separation, which is discussed in the context of your application in this paper.


How to determine variance/std deviation of Gaussian noise from measured data If I have a signal measured in discrete time intervals of a process containing some signal and some Gaussian noise how do I go about measuring the variance of this noise [i.e. the $\sigma^2$ term in $\mathcal{N}(0, \sigma^2)$].

For example if I generate some pure Gaussian noise and plot the time trace, the histogram of counts and the power spectral density like so:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import scipy.signal

sigma = 0.1
sample_freq = 1e6 # 1 MS/s
dt = 1/sample_freq
r = np.random.normal(0, sigma, 1000000)
t = np.arange(0, 1000000*dt, dt)

# Create 2x2 sub plots
gs = gridspec.GridSpec(2, 2)

fig = plt.figure(figsize=[10, 7])
ax1 = fig.add_subplot(gs[0, 0]) # row 0, col 0
ax1.plot(t, r)

ax2 = fig.add_subplot(gs[0, 1]) # row 0, col 1
ax2.hist(r, bins=100, orientation="horizontal")
for tick in ax2.get_xticklabels():
    tick.set_rotation(270)

ax3 = fig.add_subplot(gs[1, :]) # row 1, span all columns
freqs, PSD = scipy.signal.welch(r, sample_freq, nperseg=1000)
PSD = PSD[freqs.argsort()]
freqs.sort()
ax3.semilogy(freqs, PSD)
ax3.set_ylabel("$S_{xx}$ ($V^2/Hz$)")
ax3.set_xlabel("$f$ (Hz)")

I could fit the equation for the probability density of a Gaussian [i.e. $f(x|\mu, \sigma^2) = \dfrac{1}{\sqrt{2\pi\sigma^2}}e^{\frac{(\mu-x)^2}{2\sigma^2}}$] to the histogram to get the std deviation $\sigma$ of the noise.

However, in the case where we have some signal with some Gaussian noise we get something like that which I have below for a pure 100kHz sinusoidal signal:

rs = r + 1*np.sin(2*np.pi*100e3*t)

# Create 2x2 sub plots
gs = gridspec.GridSpec(2, 2)

fig = plt.figure(figsize=[10, 7])
ax1 = fig.add_subplot(gs[0, 0]) # row 0, col 0
ax1.plot(t, rs)

ax2 = fig.add_subplot(gs[0, 1]) # row 0, col 1
ax2.hist(rs, bins=100, orientation="horizontal")
for tick in ax2.get_xticklabels():
    tick.set_rotation(270)

ax3 = fig.add_subplot(gs[1, :]) # row 1, span all columns
freqs, PSD = scipy.signal.welch(rs, sample_freq, nperseg=1000)
PSD = PSD[freqs.argsort()]
freqs.sort()
ax3.semilogy(freqs, PSD)
ax3.set_ylabel("$S_{xx}$ ($V^2/Hz$)")
ax3.set_xlabel("$f$ (Hz)")

We can now no longer fit to the histogram without deconvolving the sinusoidal signal some how (which may not be so trivial for a more complex signal). However the noise level in the power spectral density remains the same, so there should be some way of getting the std deviation $\sigma$ of the Gaussian noise from this.

I understand that the equation for the power spectral density in discrete time is:

$$S_{xx}(\omega) = \dfrac{(\Delta t)^2}{T}|\sum_{n=1}^{N}x_n e^{-i\omega n}|^2$$

and the equation for the probability density of $x$ (where $x$ is drawn from $\mathcal{N}(0, \sigma^2)$) is

$$f(x|\mu, \sigma^2) = \dfrac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(\mu-x)^2}{2\sigma^2}}$$

But I am unsure how to derive the relationship between the level of the noise in the PSD and the std deviation $\sigma$ of the noise in order to extract the std deviation $\sigma$ in the noise on a measured signal.

What is typically done is to apply an orthogonal transform (or a close to tight frame) that

preserves the Gaussianity,
better separates the signal and the noise.

Fourier for simple sums of sines, wavelets for piecewise polynomial signals. Then, identify or parse some coefficient bins (bands, wavelet subbands), expecting that signal coefficients are sparse enough. Then, compute a robust scale estimator (like the Median Absolute Deviation, MAD) on the absolute values of the coefficients. And apply some factor to the result. This could work also with band-pass filters, derivative operators, etc.

You can find implementations in Matlab:

wnoisest: Estimate noise of 1-D wavelet coefficients
adapted for Ricean noise: Robust rician noise estimation for 3D MR

References:

Nonlinear wavelet thresholding: A recursive method to determine the optimal denoising threshold


Gravity center of a discrete signal One of the problems of my home assignment was this.

If i use the definition to solve this(I'm not seeing other way) i have an expression like $\sum\sum nx[n]e^{-j\omega n}$ on the numerator and $\sum\sum x[n]e^{-j\omega n}$ on denominator. How do i simplify those expressions? The double sum is making me confuse.

Thanks.

Your double sums do not make sense. The numerator and the denominator of the expression for the center of gravity are numbers, and it's pointless to take the discrete-time Fourier transform (DTFT) of these numbers. What you are supposed to do is express the center of gravity in terms of the DTFT of $x[n]$.

Since this is a homework problem, instead of solving it for you I'll give you a few hints that should help you figure out the solution by yourself:

The discrete-time Fourier transform (DTFT) of the sequence $nx[n]$ is given by $$\text{DTFT}\{nx[n]\}=j\frac{dX(\omega)}{d\omega}$$
The sum over a sequence (if it exists) is equal to the DTFT of that sequence evaluated at $\omega=0$: $$\sum_{n=-\infty}^{\infty}x[n]=X(0)$$


Chromatogram peak detection - bunching vs others? Reading presentations from existing tool providers I noticed that in order to detect peaks they first use bunching (average N points) and then use slope and curvature to detect peaks. I'm guessing they use bunching to get rid of high-frequency noise.

Aren't there better alternatives like removing high frequencies with FT? I'm concerned that bunching requires peak width to be known to define Bunching Factor. Which means we first need to detect peak width in some other way before using curvature. But I've seen more than 1 presentation talking about bunching, so there have to be reasons to use this approach instead of others.

peak(s) determination tends to be application specific.  There are many heuristics and one usually needs to try and test a lot.

It helps a lot if there is a reliable statistical model for your data.

One trick is to note that peak is larger than its adjacent neighbors.  Peaks on boundaries can be determined by reflection.

Another issue is the set that you are finding peaks on.  Some data consists of multiple scans that can be averaged to reduce noise. One can also do what you refer to as bunching and bunching has numerous variations such as summing, weighted summing, taking a median and combinations.

Another concern may be strong peaks bleeding into weak peaks or masking them entirely.  Many applications have this bin bias problem and there are ways to approach that which are beyond what can be answered in a simple answer given here.

One can also use thresholds to isolate peaks and again that has many variations.

Usually, a community has conventions on how to approach problems like finding peaks and there is probably some documents that constitute a standard. You might want to look for those.


Rectify signal basic question When rectifying a signal by squaring it, is it necessary to upsample by a factor of 2?

Why is this?

The only reason I can think of may be to preserve the frequencies over which original signal had spectrum (before squaring). When you square a signal $x[n]$, whose spectrum $X(e^{j\omega})$ extends from  $-\omega_0 \le \omega \le \omega_0$, you are multiplying by itself. So in frequency domain, the effect is to (periodically) convolve $X(e^{j\omega})$ by itself.
$$
S(e^{j\omega}) = \int_0^{2\pi} X(e^{j\theta})X(e^{j(\omega - \theta)})d\theta
$$
The convolution will result in $S(e^{j\omega})$ non-zero over $-2\omega_0 \le \omega \le 2\omega_0$. If you had upsampled before squaring, the spectrum before squaring would have been non-zero only from $-\omega_0/2 \le \omega \le \omega_0/2$. So after squaring (convolution in frequency domain), its spectrum is still $-\omega_0 \le \omega \le \omega_0$.


How to find difference equation, given poles and zeros? I'm given poles at ${1+i}$ and ${1-i}$, and zero at $0$.

I have to find the difference equation and find out whether the system is stable.

Now I found in a similar question, which structure I tried to follow ( How to find system function, H(z) in the z-domain, given zero-pole plot of the system? ), that we can express $H(z)$ as $$H(z)=G_0z^{-1}\frac{(z-z_0)(z-z_1)\ldots}{(z-p_0)(z-p_1)\ldots} $$

It's clear to me that the numerator of H(z) is simply = z, but I can't wrap my head around the denominator.

I tried to go for: (z- (1-i))*(z- (1+i)) - my prof looked over my shoulder, nodded and told me to simply multiply it, adding "there will be no i in the polynomial". Now I really don't find any combination of tricks to achieve this, and I've been stuck for several hours right now.

I'd very much appreciate a hint!

We have a zero at $z=0$ and poles at $z=1+i$ and $z=1-i$

Therefore the Transfer function is:

$$H(z) = k \dfrac{z}{[z-(1+i)][z-(1-i)]}$$

$$\Rightarrow H(z) = k\dfrac{z}{(z-1)^2+1}=\dfrac{kz}{z^2-2z+2}$$

$$\Rightarrow \dfrac{Y(z)}{X(z)} = \dfrac{kz}{z^2-2z+2}$$

$$\therefore Y(z)[z^2-2z+2] = kzX(z)$$

This gives:
$$z^2Y(z)-2zY(z)+2Y(z) =  kzX(z)$$
Taking inverse Z transform:
$$y(n+2)-2y(n+1)+2y(n) = kx(n+1)$$


What phase rotation occurs when you take the derivative of an audio signal? If you take the derivative of an audio signal, it provides a 6 dB/oct upward sloping filter (increasing high frequencies / cutting low frequencies) all the way across the spectrum.
What is the result in terms of the phase? Is there phase rotation? Is it uniform across the spectrum or does it vary with frequency? What would the plot of this phase rotation look like?
Thanks

Simple Fourier Transform propery for differentiation is
$$\mathcal{F}\{\frac{\partial}{\partial t}x(t)\}=j\omega \cdot\mathcal{F}\{x(t)\}$$
Differentiation in time corresponds to multiplication with $j\omega$ in frequency. Hence the +6dB/octave slope. The phase shift is a constant 90 degrees for all frequencies.


Causal systems - Paradox? Confusion of definition? For a system such that $y[n]=\cos(x[n])$, I'd like to determine whether the system is causal or not.

Definition 1: A system is causal if and only if the output $y[n]$ is a function of a linear combination of inputs $x[n-k]$ such that $k \ge0 $.

We can see that the system is causal by this definition because $y[n]$ only ever depends on $x[n]$.
Definition 2: A system is causal if and only if the impulse response $h[n]=0$ for all $n<0$.

The response of the system to $\delta[n]$ is $h[n]=\cos(\delta[n])$.

$$h[-1]=\cos(0)=1 \ne 0 \quad \therefore\text{The system is not causal.}$$

Have I hit a paradox? Am I missing something here?

You have your definition of causality wrong. It's actually much simpler and much more intuitive. A causal system is a system in which the output does not depend on future values of the input. This property is not exclusive linear systems and can apply to systems in general. Here are a few examples to illustrate the point:
Causal linear time-invariant system:
$$y[n] = x[n] - 2x[n-1] + 0.5y[n-1]$$
Here $y[n]$ only depends on current and previous values of $x$ and $y$.
Non-causal linear time-invariant system:
$$y[n] = \frac{1}{2}x[n-1] + \frac{1}{2}x[n+1]$$
Also known as a central moving average, this function is non-causal because for output $y[n]$, the $x[n+1]$ term peeks into the future of our input.
Causal non-linear time-invariant system:
$$y[n] = \cos(x[n])$$
This is the example you provided. As you may have guessed by now, since we never look at a future input, this system is causal. In fact, it's even more special than that. Because the current output is a function only of the current input (not past or future inputs) this system is called a memoryless system.
Causal linear time-varying system:
$$y[n] = (n+1)x[n]$$
I made this one up to be somewhat of a brainteaser. At the first glance, it may seem like this system is non-causal because there's this $(n+1)$ term. But this doesn't mater because it's not a time index of $x$. We still only look at the current value of the input and are not peeking ahead. This is also an example of a memoryless system.
A tricky causal non-linear time-varying system:
$$y[n] = e^nx[n] + \ln\left(\left|x[n-1]\right|+1\right) - \pi y[n+1]$$
Clearly, clearly this is a non-causal system because of the $y[n+1]$ term, right? Wrong! This is a classic definition of a recursive causal system with a few terms rearranged. We proceed to bring it into a more standard form in three steps:
$$y[n] = e^nx[n] + \ln\left(\left|x[n-1]\right|+1\right) - \pi y[n+1]$$
$$\pi y[n+1] = e^nx[n] + \ln\left(\left|x[n-1]\right|+1\right) - y[n]$$
$$\pi y[k] = e^nx[k-1] + \ln\left(\left|x[k-2]\right|+1\right) - y[k-1]$$
$$y[k] = \frac{e^nx[k-1] + \ln\left(\left|x[k-2]\right|+1\right) - y[k-1]}{\pi}$$
Second to last line was achieved by substituting $k=n+1$. The trick here is that there is a relationships between outputs at different times. But once you work it out, no output ever depends on future value of an input relative to itself. Non-linearity and time-variance were thrown in to make it more fun and challenging. Make sure you understand what this one says.
As you can see causality can be a property of all sorts of systems and there are many more fun and bizarre examples one can come up with.

Now let's get to resolving your paradox. The key lies in realizing that your definitions of causality and linearity (and perhaps time-invariance, too) a bit entangled and confused. Solving the paradox amusingly is as easy as adding the word linear into both of your definitions (time-invariance will also play a subtle role). Here's how.

Definition 1: A linear system is causal if and only if the output $y[n]$ is a
function of a linear combination of inputs $x[n-k]$ such that $k \ge0$.

This is because not all systems are linear combinations of inputs. Linear systems are. A causal system depends only on past and current inputs, therefore a causal linear system is a linear combination of current and previous inputs. Actually, to be accurate with the definition, linear systems are linear combinations of both inputs and outputs, and in causal linear systems outputs at time $n$ cannot depend on an input at time $m \gt n$.
Also,

Definition 2: A linear time-invariant system is causal if and only if the impulse response
$h[n]=0$ for all $n<0$.

This one is also quite intuitive. All entries of $h[n-k]$ for $k \ge 0$ are the coefficients with which you multiply current and past values of $x$ to get the current output (whether the system is recursive is irrelevant for this case). Notice that this definition only makes sense for linear time-invariant systems! This is because convolution only exists for those. Let's see why linearity is necessary. If you have $h[n+1] \ne 0$, then $x[n+1]$ is contributing to $y[n]$ and that's what makes it non-causal. This is true because in linear systems, if nothing goes in then nothing goes out. This is not generally true for non-linear systems (like in the example you gave), so this definition wouldn't apply.
The system also needs to be time-invariant because it cannot be completely defined by its impulse response unless it's full on LTI. If you run a impulse function $\delta [n]$ through the system
$$y[n] = x[n] + (n+1)x[n+1],$$
you're going to get an output which is an impulse itself (causal... right?). However, the system is clearly non-causal. This is why time-invariance is important. When you run a shifted impulse function $\delta[n-m]$ for $m \ne 0$, its non-causal nature starts manifesting itself.
So, all talking aside, your system is perfectly causal, but your definitions only apply to linear systems, while your system is non-linear. The correct definition of a causal system is that any output $y[n]$ cannot depend on input $x[n+k]$ where $k \gt 0$.


Why can Sum of Squared Differences (SSD) seen as cross-correlation? In this slides, on page 4,  the SSD can be seen as cross-correlation:

http://www.lira.dist.unige.it/teaching/SINA/slides-current/pyramids.pdf
http://www.cse.yorku.ca/~kosta/CompVis_Notes/ssd_cross-correlation.pdf

The SSD is defined as:

$$SSD(m,n) = \sum_i \sum_j [g(i,j) - t(i-m,j-n)]^2,$$

and expand the SSD. It looks like this:

$$SSD(m,n) = \sum_i \sum_j [g(i,j)^2 + t(i-m,j-n)^2 - 2g(i,j)t(i-m,j-n)], $$

where the squared terms can be assumed as constant. As a result, the cross-correlation remains:

$$R(m,n)= \sum_i \sum_j [g(i,j) t(i-m, j-n)].$$

Question: Why can we assume that the squared terms are constant?

Short Answer:
The statement "the the squared terms are constant" is true only if some further theoretical assumptions are made. In practice, the terms are often not constant (i.e. they are dependent on $m$ and $n$). For many discussions however, it's ok to treat the terms as constant.
Details:
The term $\sum_{i,j}\left[g(i, j)\right]^2$ is not interesting since it is clearly independent on $m$ and $n$. The term $\sum_{i,j} \left[t(i-m, j-n)\right]^2$, however, formally depends on $m$ and $n$, as $m$ and $n$ appear in the summed expression. Nevertheless, it may be effectively constant depending on how the bounds for the running variables $i$ and $j$ are chosen and whether $t$'s support is bounded or not.
In theoretical discussions, the image function $t(x,y)$ is usually assumed to have bounded support. This means that there are minimum and maximum values of both $x$ and $y$, outside of which all image pixels can be assumed to be  zero. This assumption is plausible because a real digital image only has a finite number of pixels, e.g., $1024$x$768$ pixels.
The sum is often assumed to cover the entire space, that is $i$ and $j$ go from negative to positive infinity: $\sum_{i,j}$ really means $\sum_{i=-\infty}^{+\infty}\sum_{j=-\infty}^{+\infty}$.
Only if these two assumptions (bounded image and infinite sum) are silently made, then the squared terms are independent of $n$ and $m$. No matter what $n$ and $m$ are you will always "hit" all the pixels of image $t$:
$\sum_{i,j} \left[t(i-m, j-n)\right]^2 = \sum_{i,j} \left[t(i, j)\right]^2, \textrm{for all } n, m$.
In practical template matching, the sum will often cover only the pixels of the template, for performance reasons. This means that, in the (larger) image, only the pixels are regarded which are currently covered by the template given the current $m$ and $n$. In these cases, the squared terms are not independent if $n$ and $m$.
As a last note, in your squared difference implementation you can choose to omit calculating the squared terms and only regard the mixed term. This way, you get the same result as if you would have taken the infinite sum (up to a constant), and thus the theoretical considerations will hold again.


Multichannel sampling with aliasing if i have 2 sensors, sensor1 and sensor2, that sample a signal on complementary points so that sensor2 samples always between sensor1s sampling points. can I achieve the doubled sampling rate with the combination of both ?
so with merging both samples, although they are aliased individually.

Yes; what you describe doesn't mathematically look any different than a single sampler at twice the rate. But they sensors need to sense the same, band-limited thing!

Actually, that's how some very high-speed ADCs work.

If you're not actually doing GS/s, getting a faster ADC would probably be easier than staggering multiple ADCs.


Matched Filter BER vs SNR curve does not match the theoretical curve I am creating a matched filter using filter command. The coefficients are also generated in another file using sigma delta modulator these are ternary coefficients (+1,-1,0). I have generated 1 and 0 with equal probability, mapped the bits with a rectangular pulse then added noise and used filter command to filter the received signal with ternary coefficients. After retrieving of the bits when I try to calculate the bit error, the plot comes out to be flat, it does not resemble the theoretical curve and the error rate gives strange values. I'll really appreciate it if someone helps me, please.

The Pb values does not seem to be right. There is not much difference that is why the plot is flat

I dont know, if this is exactly what you wanted to do, but here's a script that creates your theoretic curve. I used your mirrored coefficients as the transmit pulse, and the coefficients as the filter coefficients. This way, the pair becomes a matched filter pair.

N=10^4;        %No: of bits
 A=10;          %amplitude of S(t)
 T=length(v);          %duration of S(t)
 S=fliplr(v)*A;       %the transmit pulse
 E=norm(S)^2;         %S(t) Energy
 input = rand(1,N)>0.5;   % generating 0,1 with equal probability
 bits=(2*input-1);      %mapping the bits with S(t)
 x=kron(bits,S);

 SNR_db=0:0.25:10;
 n_var=(E/2)*10.^(-SNR_db/10);

 %Noise that gives SNR=0dB to SNR=7dB

 for i = 1:length(n_var)

  n=randn(1,length(x))*sqrt(n_var(i));   %noise
   r=x+n; %received signal

  %v=ternary coefficients, filtering the received signal and coefficients
   filtered=filter(v,1,r);

   %sampling at T & using thresholding operation
   z=sign(filtered);
   z=(z(T:T:end));

  %recovering the sent bits
  received_bits=(z+1)/2;

  %calculating the bit error
  Pb(i)=mean(abs((input)-(received_bits)));

 end

 %Plotting Pb Vs SNR
  figure(1);
 SNR=10.^(SNR_db/10);
 P_theory=(0.5)*erfc(sqrt(2.*SNR)./sqrt(2));
 semilogy(SNR_db,Pb,'o',SNR_db,P_theory,'r-')
 title ('SNR Vs Probabolity of Error');
 xlabel ('SNR (dB)');
 ylabel ('Probability of Error');
 legend('practical curve','theoratical curve');
 grid on;


Classifing audio signals to detect fault The basic idea i'm working on is automatic fault detection using audio signals captured from Motor. I have like set of sample audio signal which are recorded when there is no fault and with fault at different rpm. While hearing the audio I can clearly identify and differentiate the audio signal. But how can I differentiate that technically or programtically. Basically im from computer science background , so im having a hard time to find the method of implementation whether it can be found by audio analysis itself or if I want to implement Machine learning on audio signals

Sample signals are here

1.3500rpm-clean
2.3500rpm-fault

You might want to look at:

Detection of Abrupt Changes:
Theory and Application1
Michele Basseville `
IRISA/CNRS
Rennes, France
Igor V. Nikiforov
Institute of Control Sciences
Moscow, Russia

with the free pdf at:

ftp://ftp.irisa.fr/local/as/mb/k11.pdf

It's a little old, but was advanced when it came out.

There are a lot of Google results for this topic.


Convolution in frequency space with fft and ifft -> num output samples? Let's say you had two 1 second sound samples at 44.1khz, and you wanted to convolve them in frequency space.

You FFT each one to get two arrays of complex numbers that each have 44100 items in them.

You then multiply each bin...
for i = 0 to 44099
c[i] = A[i] * B[i]

Then it's time to turn it back into time domain samples.  If you do an ifft you'll get 44100 time samples out, but when you convolve these two audio samples in time space, it'll give you 88199 samples out (44100+44100-1).

Obviously getting 44100 samples out is incorrect.  At what part in the process did I go wrong? (:

Thank you!!

What you do is called circular convolution but you want "linear" convolution: Say you have two signals of length $N$ and $M$. Append zeros to both time-domain signals, so that their length is $L=N+M-1$, then calculate the $L$-DFT of both signals, multiply them in frequency-domain and convert the result back to time-domain.

In your case, you should append 44099 zeros to each signal before transforming it to the frequency domain. The result will then be 88199 samples long.

If the available FFT is shorter than 88199 you have to divide the signal into shorter blocks and apply the overlap add method.


Why is it necessary to have two state variables I am learning about control theory.
Let's consider this system.

$$
m a(t) + b v(t) + k x(t) = f(t)
$$

$a$ is acceleration
$v$ is velocity
$x$ is displacement
$f$ is external force

In my textbook, in chapter about "state space model" , two state variables ($x$ and $v$) are necessary
to completely describe this system.
However, I think $x$ can completely describe this mass's motion by itself.
If we know $x$, $v$ can derived from $x$.
So I have two questions.

Why two variables are needed.
Why "$x$ and $v$" ?  Can "$v$ and $a$" or "$x$ and $a$" also describe this system ?

Let's rewrite your system as

$$
m \ddot{x}(t) + b \dot{x}(t) + k x(t) = f(t)
$$

then you can see what you're saying: Everything I need to know is in $x(t)$!

But is it really?

State space systems are predicated on having a single, first order differential equation to solve.  As you can see in the rewritten equation above, there is a double-derivative term.

In order to write that equation as a first order equation, we need to have two state variables:
$$
X(t) = \left[ \begin{array}{c} x\\ \dot{x} \end{array} \right]
$$
and so
$$
A\dot{X}(t) + B X(t) = f(t)
$$
where $A=[b\ m]$ and $B=[k\ 0]$.

Let's try to do the same thing using $x$ and $\ddot{x}$.

$$
X'(t) = \left[ \begin{array}{c} x\\ \ddot{x} \end{array} \right]
$$
and so
$$
A'\dot{X'}(t) + B' X'(t) = f(t)
$$
where $A'=[b\ 0]$ and $B'=[k\ m]$.

So I suppose you could chose the state variables to be $x$ and $\ddot{x}$... but it would mean an unused $\dddot{x}$ variable.

Some examples of first order systems can be found in this document.


duration for a signal after filtering I have a signal $s(t)$ of duration $T$. This signal doesn't have a limited bandwidth. To solve this problem, I filter it with a pulse shaping filter (raised cosine filter).

My question:
Does filtering change the duration of the signal?

A classical result states that a (non-trivial)  time-limited signal does not have a  limited bandwidth (in terms of support). The converse happens as well: a bandwidth-limited signal cannot be of finite length.

And there is a similar rule-of-thumb principle:  when you try to reduce the support of a signal in one domain (time or frequency), you generally spread it in the other domain.


Finding the transfer function of a discrete signal described by two equations A discrete time system is described by the following system of equations.

$$q[n] = \big(x[n]-\frac k4q[n-1]\big)$$
$$y[n] = \big(q[n]-\frac k3q[n-1]\big)$$

Find the systen function and then find the values of $k$ for which the system is stable. Also find the values of $k$ for which the system is of mimimum-phase.

Can anybody give me some guidelines as to how to handle this system of equations in order to find the system function?

This is an unsolved exercise given by my professor.

HINT:

(because it's a homework problem)

Apply the $\mathcal{Z}$-transform to both equations. Use the first equation to express $Q(z)$ in terms of $X(z)$, and plug that into the second equation to obtain $Y(z)$ in terms of $X(z)$. From this you can easily get the transfer function $H(z)$.
Compute the poles and zeros of $H(z)$ and determine the constant $k$ such that the system is stable and/or minimum-phase. You should know what stability and minimum-phase mean in terms of poles and zeros (if not, it should be easy to find that online or in a textbook).


How to determine if the system is invertible Is there any systematic way to determine if the system is invertible?

My general approach is first trying to find the inverse system by using mathematical method; that is, solving for the output in terms of input. If this is difficult, my guess would be that the system is non-invertible, and the goal is then to sit and think of the two inputs $x_{1}(t)$ and $x_{2}(t)$ that generate the same output $y(t)$.

However, I was just wondering if this method is 100% correct:

If we can find the input $x(t)$ in terms of output $y(t)$ (by solving mathematically), does that always mean the system is invertible?
Is there any possibility that can lead to the failure in reasoning?

For example, determine if the following system is invertible:
$$ y(t)=\int_{-\infty}^{t}e^{-(t-\tau)}x(\tau)d\tau $$

Firstly
\begin{align}
y(t)&=e^{-t}\int_{-\infty}^{t}e^{\tau}x(\tau)d\tau\\
e^{t}y(t)&=\int_{-\infty}^{t}e^{\tau}x(\tau)d\tau\\
\frac{d}{dt}\left(e^{t}y(t)\right) &= e^{t}x(t)\\
x(t)&=\frac{1}{e^{t}}\frac{d}{dt}\left(e^{t}y(t)\right)
\end{align}
So the inverse system is:
$$ y^{-1}(t)=\frac{1}{e^{t}}\frac{d}{dt}\left(e^{t}x(t)\right) $$

It depends on what exactly you mean by "invertible". In system theory, what is often meant is if there is a causal and stable system that can invert a given system, because otherwise there might be an inverse system but you can't implement it.

For linear time-invariant systems there is a straightforward method, as mentioned in the comments by Robert Bristow-Johnson. First, note that the input-output relationship of an LTI system is given by the convolution integral

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

Comparing $(1)$ to the system given in your question shows that the system is LTI with impulse response

$$h(t)=e^{-t}u(t)\tag{2}$$

where $u(t)$ is the unit step function. The corresponding transfer function is

$$H(s)=\frac{1}{1+s}\tag{3}$$

This system is causal and stable. However, its inverse system

$$G(s)=\frac{1}{H(s)}=1+s\tag{4}$$

is not stable. It has a pole at infinity. The input-output relation of the inverse system is

$$y(t)=x(t)+x'(t)\tag{5}$$

where $x'(t)$ is the derivative of $x(t)$. So in order to implement the inverse system you need a differentiator, which is not stable in the bounded-input-bounded-output (BIBO) sense. I.e., bounded input signals can lead to unbounded output signals (just use a rectangular input function as an example).


Viterbi equalizer in a GSM system I'm about to tackle implementing a Viterbi equalizer fo a GSM system in gnuradio.
I went through the internet and ieeexplore and I think that this document explains it best:

  Application of Viterbi Equalizer in GSM
  System

My question concerns the Matched Filter in the block diagram.
Normal Burst in GSM has 8 possible midambles. First of all, the Matched Filter block should determine which midable has been transmited. Midambles have highly-peaked auto-correlation function, so the first thing to do is to compute the modulated version of the 8 midambles and then convolve each of them with the received midamble and the one with the most peak-like is the most likely to have been sent. Once I have that done, it's easy to compute the channel response estimate.

I don't have much background in DSP and I was hoping someone could let me know if I understand it correctly.

Thanks for the answer. You are correct - the training sequence is known by the time you establish any communication.

In case if somebody stumbles upon this in the future, here's how to get this done:

The synchronization burst contains an extended training sequence. It's longer than the regular one and there is only one sequence. You are supposed to use that to equalize, demodulate and decode the SCH burst. Among other useful information, it contains the color code of the BTS, which identifies the training sequence used by that BTS.

The procedure is nicely explained in this book:

Sajal Kumar Das: "Mobile Handset Design"


Conversion from laplace transform to z-transform I would like to know if

$$
  \text {Z-Transform ( }  G(s)H(s) \text{ )} = \text {Z-Transform (}G(s) \text{)} \text { Z-Transform (} H(s) \text{) } = G(z)H(z)
$$

where G(s), H(s) are the Laplace transform representations of g and h, and G(z) and H(z) are the Z-transform representation of g and h.

Is this relation true?

Your question makes no sense. Z transform is performed on a discrete signal/series.

Since $H(s)$ is a continuous function, you can't just calculate a Z-transform of $H(s)$ without first sampling it, to make it discrete. Also, it doesn't make much sense to do a time->spectrum transform (such as a Z-transform) on a spectral representation ($H(s)$)

I'm assuming that when you write "$Z-Transform(H(s))$", what you really want to do is to convert $H(s)\to H(z)$, meaning to calculate the Z-transform of $h[nT]$, where $h[nT]$ is $h(t)$, sampled at intervals of $T$, and $h(t)$ is the inverse Laplace transform of $H(s)$.

If I'm correct in my assumption, the transformation you are seeking is known as "star transform", which would provide a transform function $H^{*}(s)$, in terms of $e^{sT}$, which may be easily converted to $H(z)$ by way of the substitution $z=e^{sT}$.

Edit - some elaboration on the conversion process: what you need to do is calculate $H^{*}(s)$ from $H(s)$ using one of the two relations described in "Relation to Laplace transform" in the Wikipedia article, then do the substitution $z=e^{sT}$, to get $H(z)$.


Difference between frequency components and harmonic components - Fourier What is the difference between frequency components and harmonic components? The first concern the continuous domain of frequency, while the second concern the discrete domain of frequency ($f_{k}=kf_{0}$ with $f_{0}$ foundamental frequency)?

Thank you for your time.

Frequency components are $f_k$ for $k=1,\ldots,K$.  There is no relationship between them.

Harmonic components are $f_p$ for $p=1,\ldots,P$ and
$$
f_p = p f_1
$$


Regarding the Nyquist Criterion I recently read that with compressive Sensing it possible to sample the signal at a rate lesser than that suggested by Nyquist Criterion. However, I am still not getting how this is possible. Can someone suggest an intuitive way of understanding this concept of sampling at lower rates? And if this is possible, do modern communication techniques use this concept?

The Nyquist criteria refers not to the frequency, but to the bandwidth, which is related to information density in a signal.  A very high frequency signal, of approximately known frequency, with a sufficiently small bandwidth, will still be aliased or folded down with baseband frequencies by undersampling.  But if the bandwidth (or other known characteristics) of the signal is known to rule out other aliased frequencies (such as the non-existence of baseband spectrum content), that knowledge in conjunction with the undersampled samples may still provide enough information to reconstruct the signal.


How does resizing an image affect the intrinsic camera matrix? I have a camera matrix (I know both intrinsic and extrinsic parameters) known for image of size HxW. (I use this matrix for some calculations I need).

I want to use a smaller image, say: $\frac{H}{2}\times \frac{W}{2}$ (half the original).
What changes do I need to make to the matrix, in order to keep the same relation ?

I have, $K$ as the intrinsic parameters, ($R$,$T$ rotation and translation)

$$\text{cam} = K \cdot [R T]$$

$$K = \left( \begin{array}&a_x &0 &u_0\\0 &a_y &v_0 \\ 0 &0 &1\end{array} \right)$$

$K$ is 3*3, I thought on multiplying $a_x$, $a_y$, $u_0$, and $v_0$ by 0.5 (the factor the image was resized) , but I'm not sure.

Note: That depends on what coordinates you use in the resized image. I am assuming that you are using zero-based system (like C, unlike Matlab) and 0 is transformed to 0. Also, I am assuming that you have no skew between coordinates. If you do have a skew, it should be multiplied as well

Short answer: Assuming that you are using a coordinate system in which $u' = \frac{u}{2} , v' = \frac{v}{2}$, yes, you should multiply $a_x,a_y,u_0,v_0$ by 0.5.

Detailed answer The function that converts a point $P$ in world coordinates to camera coordinates  $(x,y,z,1)->(u,v,S)$ is:

$$ \left( \begin{array}{ccc}
a_x & 0 & u_0 \\
0 & a_y & v_0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
R_{11}  & R_{12} & R_{13} & T_x \\
R_{21}  & R_{22} & R_{23} & T_y  \\
R_{31}  & R_{32} & R_{33} & T_z \\
0  & 0& 0 & 1
\end{array} \right)
\left( \begin{array}{ccc}
x \\
y  \\
z \\
 1
\end{array} \right)
$$

Where $(u,v,S)->(u/S,v/S,1)$, since the coordinates are homogenous.

In short this can be written as
$ u= \frac{m_1 P}{m_3 P} , v = \frac{m_2 P}{m_3 P}$
where $M$ is the product of the two matrixes mentioned above, and $m_i$ is the i'th row of the matrix $M$. (The product is scalar product).

Re-sizing the image can be thought of:

$$ u'=u/2, v'=v/2 $$

Thus

$$ u' = (1/2) \frac {M_1 P} {M_3 P}  \\
v' = (1/2) \frac {M_2 P} {M_3 P}
$$

Converting back to matrix form gives us:

$$
\left( \begin{array}{ccc}
0.5 & 0 & 0 \\
0 & 0.5 & 0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
a_x & 0 & u_0 \\
0 & a_y & v_0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
R_{11}  & R_{12} & R_{13} & T_x \\
R_{21}  & R_{22} & R_{23} & T_y  \\
R_{31}  & R_{32} & R_{33} & T_z \\
0  & 0& 0 & 1
\end{array} \right)
\left( \begin{array}{ccc}
x \\
y  \\
z \\
 1
\end{array} \right)
$$

Which is equal to

$$ \left( \begin{array}{ccc}
0.5 a_x & 0 & 0.5 u_0 \\
0 & 0.5 a_y & 0.5 v_0 \\
0 & 0 & 1 \end{array} \right)
\left( \begin{array}{ccc}
R_{11}  & R_{12} & R_{13} & T_x \\
R_{21}  & R_{22} & R_{23} & T_y  \\
R_{31}  & R_{32} & R_{33} & T_z \\
0  & 0& 0 & 1
\end{array} \right)
\left( \begin{array}{ccc}
x \\
y  \\
z \\
 1
\end{array} \right)
$$

For additional information, refer to Forsyth, chapter 3 - Geometric camera calibration.


How to remove historical values from exponential filter Suppose I have two time series $A:=[0,0,0,4,5,6]$ and $B:=[1,2,3,4,5,6]$. I implement the following filter (with initial value at the first element):

$$ f(t) = f(t-1) + \alpha*(f(t)-f(t-1))$$

with $\alpha \in [0,1]$.

Is there any way I can use filtered values of $B$ to recover the filtered values of $A$? For example suppose I can store all historical filtered values of $B$ for any choice of $\alpha$ in memory, is it possible to reconstruct what the filtered values for $A$ would have been for some choice of $\alpha$?

Intuitively this would involve subtracting the some filtered value(s) of $B$ at time step $3$ from the filtered value(s) of $B$ at time step $6$, to recover an approximation of $A$.

Yes. Your filter $f(t)$ is an LTI system, so you can use that to your advantage. Let:

$$
C(t) = B(t) - A(t)
$$

Then, if $y_A(t)$, $y_B(t)$, and $y_C(t)$ are the system's responses to $A(t)$, $B(t)$, and $C(t)$ as inputs, then you can straightforwardly do what you want:

$$
A(t) = B(t) - C(t) \rightarrow y_A(t) = y_B(t) - y_C(t)
$$

So, pass $B(t)$ through the filter, then pass $C(t) = [1, 2, 3, 0, 0, 0]$ through the filter and find the difference between the two filter outputs. That's exactly what you would get if you passed $A(t)$ through the filter on its own.


Is the path loss exponent $\alpha$ is always 2 in friis equation model? In the friis equation,
the reference distance is 1m at the log distance path loss.

Then the formula would be:

$$P_r = P_tA_rA_t(\lambda/(4\pi))^2 /d^\alpha$$

where $d$ is distance between terminals, $P_r$ is received power, $P_t$ is transmission power, $A_r, A_t$ are the rx, tx antenna gains, respectively, and
$\lambda$ is wavelength.

It seems that the path loss exponent $\alpha$ is not mentioned.

But I haven't seen a case where this $\alpha$ isn't 2.

Is this $\alpha$ always 2 in the friis equation?

And can the simulation environment be used even if this $\alpha$ is not 2?

Path loss exponent is not always equal to 2, only in the free space scenario is this true. In general the received power $P_r$ is proportional to $d^{-\alpha}$, where $\alpha$ is somewhere around 2-6 depending on the situation. The choice of the path loss exponent should match the scenario you're trying to simulate.

For example, for communication at a high altitude with little obstructions then $\alpha$ is probably close to 2 since the situation is approximately free space condition.

For communication in a building, within the same floor, then $\alpha$ may be around 2-3. This is because walls other obstructions attenuate the signal power further than the free space condition.

Things can get even worse if we again consider communication in a building but now if the transmitter is on the first floor and receiver is on the third floor, then $\alpha$ could be as high as 4-6 depending on construction of the building, etc.


how to detect center of a blurry circle with opencv I have got the following image:

There are curvs on the picture. i would like to find center of the circles containing curves.

i tried opencv and hough circle transform but had no results.

You'll need to improve the contrast of the image first, then filter it strongly to remove the noise. since the circles are 'thick' (blurry), you can filter quite a bit without destroying the circles structure.

I would then apply some edge detection algorithm to get a binary edge image that can be processed by the Circular Hough Transform.

I get the following edge image from your image:

Using the following MATLAB commands:

 % x is the input grayscale image. First we adaptively improve the contrast over the image
 y= adapthisteq(x);

 % next we use the Canny edge detector with a strong Gaussian lowpass filter
 ee=edge(y, 'canny', [], 5);


How to verify if two images have the same object/content? I am truly new to computer vision field, but it is fascinating me! I have now a challenge in my hands and I am seeking for mentors/advisers to give me some follow up.

My project is: From a picture of a video-game cover, search that pic through a video-game cover database, and if there is a very good match the app will return a string with the name of the videogame and the platform.

Problem example:

Take a photo of a cover similar to this:

Cover match to this in the database:

App gives string: "Fifa 12 Playstation 2"

In my preliminary research I found that I should save in my cover database, the name of the game, the platform, the URL for the cover, and the image features (keypoints and descriptors).

I am using SURF features detector/extractor.

There are some concepts that I am still confused about... I am not looking for similarity right? I just need to observe if there are some sort of good keypoints-paired/matched, right? Because in the examples images mentioned before I get "img1 - 1087 features (query image), img2 - 1755 features - 30 % - SIMILARITY - 321/333 inliers/matched" What are the inliers? My calculation of similarity seems wrong to me... I would say that these two images are like 70 % look alike...

ps: I am using Python so...

Results so far:

Thanks for your time and help. Sorry If I could not explain any better my problems/concerns.

I am no expert, but I would like to give some pointers which could help you to solve.

Regarding the method you have mentioned/used, is for image registering/stitching. You are interpreting the results wrong. Calculating the percentage of feature match is not the intention of inliers, rather computing corresponding feature points in both images which could be used for image registering/stitching. Inlier is the feature point for which the algorithm has found a corresponding feature point in the other image. Feature points for which the algorithm has not founda match are called outliers. More inliers, better the match. Keep in mind, the algorithm can sometimes match some feature points of completely unrelated images(try running your method on two unrelated images)
Can you use the method you have mentioned to solve your problem ?.
Yes,You can use number of inliers for each image in the database and decide the one that has the highest number of inliers as the match. You can have a minimum count of inliers to be present, in case the picture of a video-game doesnt exist in your database.
Other ways/methods which you could explore to solve this.

Explore more on the features you can use. SURF does well for some applications, for some it doesnt work well. You have to choose based on your data (image database). You can check with SIFT, and also color histogram as features.
You can check bag of words for image search. (check this)
Whenever you add a new game to database, some algorithm(eg - bag of words) need the code book to be recomputed everytime you add a new image to database, so choose the algorithm wisely , coz as the image database gets bigger, the cost of computing  increases.
Also think in lines where scale of the images differ(i.e image is database and search query image) Some features like SURF/SIFT are scale invariant, if you choose to use any other, keep this in mind.

These points might help you in deciding and I hope I did not complicate your life :)


Nyquist Frequency Phase Shift The figure below shows in dashed lines sinusoidal signals of the same frequency at three different phase shifts.  The signals are then sampled such that the sinusoidal frequency is exactly a half of the sampling frequency, i.e. the frequency of all the sinusoids is the Nyquist frequency.  The samples taken from this signal are represented by the circles.

From the figure, it seems that the amplitude of the digital sinusoidal signals is dependent upon the sampling rate and instant of sampling.  In fact if the sampling times coincide with the zero-crossings of the sinusoid, then no signal will be detected at all.

I had initially thought that sampling a bandlimited analog signal at the appropriate sampling frequency would enable perfect reconstruction, but this counter-example has left me stumped.  It seems that this sinusoid will generally not be reconstructed properly if digitized and then reconstructed at this rate.  Have I gone wrong in my understanding, and if so, can someone please point me in the right direction?

Cheers!

The sampling theorem states that $f_\mathrm{S} \geq 2f_\mathrm{max}$, where $f_\mathrm{S}$ and $f_\mathrm{max}$ are the sampling and maximum signal freuqency, respectively. But there's an additional condition: The equal sign only holds if the signal spectrum does not contain a dirac impulse at $f_\mathrm{S}/2$ which is clearly the case in your example. Therefore $f_\mathrm{S} > 2f_\mathrm{max}$ has to be fulfilled.
A reference for this statement can be found here, for example:

If a signal contains a component at exactly B hertz, then samples spaced at exactly 1/(2B) seconds do not completely determine the signal, Shannon's statement notwithstanding. This sufficient condition can be weakened, as discussed at Sampling of non-baseband signals below.
More recent statements of the theorem are sometimes careful to exclude the equality condition; that is, the condition is if x(t) contains no frequencies higher than or equal to B; this condition is equivalent to Shannon's except when the function includes a steady sinusoidal component at exactly frequency B.

The above is an excerpt of an earlier version of the according Wikipedia article.


When and how can multiple EQ filters be merged into one? As I understand it, it should be possible to merge similar filters, e.g. if I have two shelving filters with the same Q and frequency but different gains, I can just merge them into one by adding the gains, and it will have the same effect as two. Is this correct?

In what other cases is this possible?

This concept generalizes to any series cascade of linear filters. If I have a collection of linear filters with frequency responses $H_1(f), H_2(f), H_3(f), \ldots$ that are applied in series (i.e. an input signal is applied to $H_1$, its output is applied as the input to $H_2$, and so on), that series cascade could be replaced by a single filter whose overall frequency response is equal to $H(f) = H_1(f)H_2(f)H_3(f) \ldots$. This follows directly from the properties of LTI systems.

With that said, I'm not too familiar with some of the audio EQ terminology that you're using, so I'm not sure. If you have two filters with the same frequency response apart from a constant gain factor, like:

$$
H_1(f) = k_1 S(f)
$$

$$
H_2(f) = k_2 S(f)
$$

where $S(f)$ defines the shape of the frequency response, then the frequency response of the cascade of the two filters is:

$$
H(f) = H_1(f) H_2(f) = k_1 k_2 S^2(f)
$$

Note the squaring of $S(f)$; this is an important distinction. Applying two filters in series with the same frequency response shape (apart from a constant) will have a different effect than if you just created one similar filter with a different gain. Specifically, due to the squaring of $S(f)$, the passband will experience twice the ripple and the stopband will experience twice the attenuation than you would observe had you just used one of the two filters.


pole/zero plot & impulse response in a mock exam of dsp was asked this question, and even in knowing how to compute DFT's of the plot, I wondered if there was a fast way of pairing these plots to their respective impulse response?

Thanks for reading!

It's pretty simple if you know a few basic things:

a single pole on the real axis corresponds to an exponential sequence. The elements of that sequence either have the same sign, or their signs alternate, depending on the angle of the pole ($0$ or $\pi$). This should allow you to pair figures $(c)$ and $(e)$ with figures $(1)$ and $(6)$.
a complex conjugate pole pair corresponds to a sinusoid with exponentially decaying or increasing amplitude, depending on whether the poles are inside or outside the unit circle. This should help you pair figure $(f)$ with the corresponding sequence.
if there are only zeros in the PZ-diagram (and poles at the origin or at infinity), then the corresponding sequence has finite length. Note that since all shown sequences are causal, there must be the same number of poles at the origin as there are zeros. Those poles are not shown in the figures, which could be considered bad practice.
linear phase sequences have certain symmetry properties in the time domain as well as in their PZ-diagrams.
the formula for the geometric series helps with figure $(5)$.

Also take a look at this related question and its answers.

  a5,b2,c6,d4,e1,f3


Difference Between `fft(signal, nfft)` and `fft(signal)` I want to know the difference between these following commands (matlab) and why I get two different outputs?

fft(signal);

fft(signal, nfft);

I found that the output arrays are not the same, but I don't know why!
how it works with the signal of any length without padding the original signal with zeros!

signal length is't a power of 2
nfft is the next power of 2

Assuming a signal x[n] of length L, The Matlab function fft(x) computes an L-point DFT of the signal, whereas fft(x,N) computes an N point DFT of x[n]. For certain reasons N must be greater or equal to L, otherwise a wrong (aliased) DFT result would be computed, so that one cannot reconstruct x[n] back from such a DFT via inverse DFT.

The two calls would result in two DFTs which are of different lengths and which ,in general, take different values. This is because they are in fact L and N samples of DTFT $X(e^{j\omega})$ of the signal x[n].


AM Modulator Demodulator in Python I am trying to make a AM modem using Python Numpy and Matplotlib. I am successful in generating the AM Signal but I cannot demodulate it using a lowpass filter function that I made. The demodulation equation that I am using is;

$$A_{c}m(t)\cos(\omega_{c}t)\cdot 2\cos(\omega_{c}t)$$

My modulated output is this;

The demodulated signal is however this;

When it should have been this;

The lowpass filter that I am using is;
$$h(t) = \frac{\sin(2\pi Bt)}{\pi t}$$

Finally, this is my code;

import matplotlib.pyplot as plt
import numpy as np

t = np.arange(0, 100, 0.1)

message_freq = 50
message = np.sin(6.28*message_freq*t)
carrier_freq = 1000

class Filters():
    def LowPassFilter(self, cutoff):
        f = np.sin(6.28*cutoff*t)/(3.14*t)
        return f

class AmplitudeModem:
    def modulate(self, message, fc):
        c = np.cos(6.28*fc*t)

        mod = message*c
        return mod

    def demodulate(self, received_array, fc, fm):
        c = np.cos(6.28*fc*t)

        demod = 2*c*received_array
        baseband = demod #- self.modulate(Ac, received_array, 2*fc)
        return np.convolve(Filters().LowPassFilter(message_freq), baseband)

modulated = AmplitudeModem().modulate(message, carrier_freq)
demodulated = AmplitudeModem().demodulate(modulated, carrier_freq, message_freq)

Another problem that I am facing is that, I cannot plot the convolution (used in demodulate() function) with respect to time, due to their different dimensions.

This is wrong
Because if your carrier frequency is less than your message frequency, then you will have trouble in filtering .


I do not understand this simple bandstop filter example I have a periodic signal like this one with period $T=120$:

I would like to apply a bandstop filter to the signal, the stop band should go from 90% of $\frac{1}{T}$ to 110% of $\frac{1}{T}$.

So I run this code in Octave

A=-ones(1,35);
B=ones(1,45);
C=-ones(1,20);
D=ones(1,20);
E=[A B C D];
x=repmat(E,[1 6]);

figure
plot(x,'b');
axis([1 numel(x) -1.5 1.5]);
title('Input signal');
print -dpng input.png

figure
plot(x,'b');
axis([1 numel(x) -1.5 1.5]);

signal_period=numel(E);
signal_frequency=1/signal_period;

sampling_period=1;
sampling_frequency=1/sampling_period;

% "The Nyquist frequency is half the sample rate",
% from http://it.mathworks.com/help/signal/ref/fir1.html#inputarg_Wn
Nyquist_frequency=sampling_frequency/2;

signal_normalized_frequency = signal_frequency/Nyquist_frequency;

b = fir1(40,[0.9*signal_normalized_frequency 1.1*signal_normalized_frequency],'stop');
y = filter(b,1,x);

hold on
plot(y,'r');
title('Input (blue), output (red)');

print -dpng filter_response.png

and I get this filtered signal:

I would expect an attenuation of the signal but I only get a slightly "delayed" signal with some "ringing".

What am I missing?

My ideas:

I made a mistake in fir1 with the definition of the lower and
upper cutoff frequencies, perhaps I do not understand the normalized
frequency concept?
Maybe the cutoff frequencies are too low? The frequency response of the filter is the following: and I cannot see the classical "deep" stop band of a notch filter in the magnitude graph.

fir1 designs your filter using a Hamming window one bin larger than the degree you specify. 41 bins is too small to give you the frequency resolution you want. Try a much larger degree, say 300-600.


How can a signal be both periodic and random? Do any examples of such signals exist where the signal is both periodic and random? Because as I see it, if a signal is periodic then the randomness kinda goes away because you know what the signal will look like after some time.

An example you run typically across in a text book (Papoulis as an example) is the sine with random phase
$$
x(t)=\sin(2\pi f + \phi)
$$
where $\phi$ is a random variable, distributed uniformly, over $0$ to $2\pi$.

Any realization will have $\phi$ take on a particular value, but it’s random, just like a 6 on a dice after a throw. You could not predict it would be 6 before the throw or in this case you could not predict $\phi$ before you received the signal.

At the other extreme, you can construct a very “random” looking waveform using something like randn() in Matlab. It will pass several statistical tests for randomness but if you use the same seed, the waveform is perfectly repeatable (deterministic).

Every waveform has a story ( a context ) and the story is how you judge if something is random or not.


Cross correlate a 2D array I have a 2D array of eeg data with shape (64,512) - 64 electrodes, 512 timepoints

I want to calculate the maximum cross correlation (irrespective of lag/time shift) between every single electrode, so I end up with a 64x64 array containing max cross correlation values between all pairs

Is there an efficient way of doing this in python/numpy/scipy without iterating through all pairs of electrodes?

I tried using scipy.signal.correlate2d but I'm not sure its doing what I think its doing as I end up with a 2D array of size 127x1023 rather than 64x64:

from scipy import signal
import numpy as np

data = np.random.randint(1,100,(64,512))

xcorr = signal.correlate2d(data,data)

What you have (conceptually) is not a 2D array but a collection of 1D arrays. correlate2D is designed to perform a 2D correlation calculation, so that's not what you need.
Iterating through all pairs is not a big ask really - you can still use numpy to perform the cross correlation, you'll just need to have two loops (nested) to determine which signals to perform the calculation on. With only 64 signals that shouldn't take long.

for first in range(n_signals):
    for second in range(first + 1, n_signals):
        corr = numpy.correlate(data[first], data[second], mode='full')
        max_corr[first, second] = numpy.max(corr)


How to detect the fuzzy edge in this image? I am trying to process a lot of images of bacteria. As you can see from this image, when several bacterias are very close, the boundary tends to be very fuzzy.

I try to binary the image or detect the edge use canny operator.But the result is unsatisfying.

So is there a right way to distinguish every cell in this situation?I want both get the right count and position of every cell. Thanks to all of you.

bacteria segmentation is a broad field with many approaches. here's an example. For your spesific image, I think you can do the following tricks:

d=double(imread('bac.png'));
d=d./max(d(:));

ws=7; % window size
mean_d=imfilter(d,fspecial('average',ws),'replicate'); % this is a local mean operation
d2=d-mean_d;

d3 = conv2(d2.*(d2>0), fspecial('laplacian',0),'same');
subplot(3,1,1);imagesc(d);colorbar
subplot(3,1,2);imagesc(d2);colorbar
subplot(3,1,3);imagesc(d3);colorbar

you can take it from there...


Is sampling at double the desired reproduction frequency accurate? This is probably a general principles question, though I'm thinking specifically in relation to sound, which is commonly sampled at a rate of 44.1 kHz in part because the maximum frequency the average person can hear tops out at around 20 kHz.

Given that the sampling theorem says the sampling frequency must be greater than twice the maximum frequency one wishes to reproduce, that means you need a sample rate of at least 40 kHz if you're trying to reproduce audio.  However, isn't it the case that the theorem is establishing only the lower-bound at which accurate reproduction becomes possible, rather than the sampling rate at which accurate reproduction is guaranteed?

By way of illustration, here are some examples showing a sine wave with samples taken at twice its frequency (in each example the source data and the interval between samples is the same, and the only thing that varies is the timing of the first sample in the interval):

In the 'In Phase' example, your samples occur exactly in time with when the wave is at a peak/trough, and you get an accurate capture that could reproduce the source wave (or something reasonably like it).  I assume this represents the ideal case, and rarely happens in practice.

In the 'Out of Phase' example, however, a terrible stroke of luck has the samples occurring at exactly the same time that the source data is crossing '0' amplitude.  This results in a capture of [0,0,0,0], which I assume can only reproduce as silence.  If so, that's not accurate at all.  I assume this represents the worst-possible case, and also rarely happens in practice.

In the 'Mid-Phase' example, the samples occur midway between the source data's peaks and troughs.  This captures a sine-wave pattern ([1,-1,1,-1]), but at half the amplitude of the source signal.  That doesn't seem like adequate accuracy for proper reproduction of the source.  I assume some variant of this scenario (where the sample rate is neither perfectly in-phase nor perfectly out-of-phase with the source) is the most common practical case, by far.

Thus it seems like if you sample at double the target reproduction frequency, most of the time your recorded data will have the wrong amplitude.  It may occasionally be accurate, but it's also just as likely to occasionally miss the signal completely.

Is this correct?  And if so, what is the sample rate one needs to use to make accurate reproduction guaranteed rather than merely possible?

If that can only be answered by imposing some constraints upon the definition of "accurate reproduction", let's say there must be zero chance of capturing silence wherever a signal is actually present, and that the recorded amplitude must always be accurate to within 5% of the actual source.

There is not only one sampling rate that theoretically allows exact reproduction. The sampling theorem states that it is sufficient for perfect reconstruction if the sampling frequency $f_s$ is greater than twice the highest signal frequency:

$$f_s>f_{max}\tag{1}$$

Without further restrictions, Eq. $(1)$ is also necessary for baseband signals occupying the whole frequency range $[0,f_{max}]$. Some people like to discuss whether in $(1)$ there should be a "greater or equal" sign or just a "greater" sign. In practice, this question is totally irrelevant because you always have to sample at a frequency twice the maximum frequency plus some extra tolerance to allow for non-ideal filters.

As you've noticed, sampling a sinusoid at exactly twice its frequency will generally not allow you to reconstruct it. As soon as you sample at a rate slightly larger than twice its frequency you will be able to reconstruct it perfectly.

Note that there is band pass sampling for band pass signals where the sampling rate can be considerably smaller than twice the highest signal frequency for perfect reconstruction to be possible. However, it's generally not true that sampling at twice the signal's bandwidth is sufficient, even thought that may be the case for specific signals.

These are only the basics of (uniform) sampling of 1-dimensional signals. There is much more to be learned. Browse the "sampling" tag on this site, and read the overview article Sampling - 50 Years After Shannon by M. Unser.


What is the source of the vibrations apparent the FFT I have a system which is composed of 3 Rotating components
•   Motor i.e. Shaft (running at around ~7rev/s)
•   3 Gears of x16 teeth mounted in series
•   shaft hold Bearings (number of rolling elements unknown)

I am doing an FFT using Matlab on the accelerometer readings, mounted over the shaft holding bearing.

i stack these FFTs and plot them in time(called sample in the image) and get the image below

Now my questions are 2:

This is the main question: what is the source of vibration indicated in frequency range 200-250Hz, since i know the gear meshing frequency ~120Hz (clear in the FFTs image), i can eliminate it being the gears. But maybe it is the gear in the middle meshing with the 2 other gears on each side? or is it just noise ? any help on this one be much appreciated
second question: Why do the frequencies drift to the right, i.e. the meshing frequency seems to be increasing. my take on this is the motor heating up after running and so drawing less current, therefore less torque is applied and then the motor runs at higher speed ?

First, I find the 3D waterfall style plots pretty hard to interpret. My preferred method is a 2D color map (e.g. pcolor() in matlab).  It's a lot easier to see your frequencies.

To your first question, it's going to be hard to answer that without knowing a lot more about your system.  Vibration at twice gear mesh frequency is certainly possible.  It's also possible that you have a structural mode.  The first thing to do to investigate this is to do a frequency sweep of your drive speed. Look to see whether this higher frequency tracks with the drive speed.  If it is synchronous with the drive speed, then you know it is definitely related to something rotating.

Your second question is essentially impossible to answer without knowing a ton more information about your system.  e.g. what kind of motor is it?  permanent magnet DC? 3 phase AC induction?  How are you controlling it?  Is there any load on the gearbox or is it spinning freely?  You'd probably need to post all of the electric schematics, pictures of the entire setup, and more data.  I doubt it's the kind of question that would be easily answered on this kind of forum.


Checking linearity/time-invariance etc. for recursive system I'm familiar with the general methodology for checking if a system is linear:

Define output $y_1[n]$ for input $x_1[n]$, define output $y_2[n]$ for input $x_2[n]$
Define $x_3[n] = ax_1[n] + bx_2[n]$
Compute output $y_3[n]$ for input $x_3[n]$ and verify whether $y_3[n]=ay_1[n]+by_2[n]$. If it is, the system is linear.

However, I am having trouble applying this method to systems that are recursive. Some recursive systems can be converted to an iterative form by inspection. For example:
$y[n] - y[n-1] = x[n] \implies y[n] = x[n] + y[n-1] \implies y[n] = \sum\limits_{k=-\infty}^n x[k]$

For others, I am having more trouble, such as $y[n]+y[n+1]=x[n]$

First, I do a change of variables:

$k = n+1 \implies y[k]+y[k-1]=x[k-1]$

Following the general approach:

$y_1[k]+y_1[k-1] = x_1[k-1]$

$y_2[k]+y_2[k-1] = x_2[k-1]$

$x_3[k] = ax_1[k]+bx_2[k]$

$y_3[k]+y_3[k-1] = x_3[k-1] = ax_1[k-1]+bx_2[k-1]$

Here is where I got stuck. I was anticipating that the system should be linear, but due to the recursive term, I do not see how the approach I took could yield $y_3[n] = ay_1[n]+by_2[n]$.

Does anyone know where I made a mistake, or if there is a different approach I should be taking? (For instance, is there a general process for going from a recursive expression to an iterative one that I could use?) I'm a student using Oppenheim and Schafer for reference.

You made no mistake. Following the hint in Dilip Sarwate's comment, you will be able to show the following input-output relation of the system:

$$y[n]=x[n-1]-x[n-2]+x[n-3]-+\ldots=\sum_{k=-\infty}^{n-1}(-1)^{n-k-1}x[k]\tag{1}$$

Using this relation, showing linearity will be much easier.

Furthermore, note that any system described by a linear difference equation with constant coefficients is linear and time-invariant (if the initial conditions are zero!). And your system obviously belongs to this class.


Which of these signals have same Nyquist rate as $x(t)$? The options are:

$x^2(t)$
$x(2t)$
Time derivative of $x(t)$
Convolution of $x(t)$ with itself

I guess, the Nyquist rate remains the same as long as bandwidth of $X(w)$ remains the same, as Nyquist rate = 2 * bandwidth.

The Fourier transform of the fourth option is $X^2(w)$. I think this has the same bandwidth as $X{w}$, because the zero terms outside of the bandwidth give zero squares. So this is correct

The third option: Derivative leads to multiplication by jw. So, again, jw(X(w)) is zero outside of the bandwidth of $X(w)$. So this is correct

Second option: i think it's incorrect, because scaling time domain also scales w domain, which changes badnwidth.

First option: This leads to convolution in the w domain. I think convolution does not preserve bandwidth. Because, say, for some k within the bandwidth, the term $X(k)X(w-k)$ goes outisde the bandwidth.

Is my work right?

Your understanding is right.

With the usual Fourier transform pair notation :

$$x(t) \longleftrightarrow X(\omega)$$

and assuming bandwidth of $x(t)$ is $W$, then the four cases will be

1-) $x^2(t) \longleftrightarrow X(\omega) \star X(\omega) \implies Bandwidth = 2W$
2-) $x(2t) \longleftrightarrow \frac{1}{2}X(\omega/2) ~~~~~ ~\implies Bandwidth = 2W$
3-) $ x'(t) \longleftrightarrow j\omega X(\omega) ~~~~~ ~~~~~\implies Bandwidth = W$
4-) $ x(t)\star x(t) \longleftrightarrow X^2(\omega) ~~~\implies Bandwidth = W$

You can see that, cases 3 and 4 has the same bandwidth as $x(t)$ and therefore their Nyquist sampling frequencies will also be the same.


Morphological transformations for images with discrete values I have a series of images in which the pixels take integer values corresponding to a classification. I'd like to try some morphological transformations (e.g. erosion, dilation) to remove noise and smooth across pixels that are classified differently than its neighbors.

However the only morphological transformation functions I know in Python operate on continuous values. This is obviously a problem for discrete values.

I'm sure other people have wanted to apply the techniques to images with classification values. Does anyone know of functions that exist for this sort of thing?

Morphological operators are primarily defined over binary images where values 0,1 could correspond to areas identified as foreground and background respectively.

In your case, you can create a binary image of the region you want to apply the operator on by assigning the "islands" of noisy classified regions to the "foreground" (e.g 4  now becomes the value 1) and the broad area around them that is consistently classified in a stable way to the "foreground" (e.g value 3). Apply erosion to this image and then substitute the pixels back from the new binary image to the segmentation mask.

Alternatively, if you want to eliminate spurious classifications so that you end up with consistent areas, you could apply nonlinear filters such as median, max and min over some sliding m by n window.

For more information, please see this link and this link

As far as Python is concerned please see this or this and this.  Another excellent python module for this job is this one.

(Also, please note that you could apply any kind of filter that could even produce a Real type output (e.g 2.43, 1.76, etc) but either round or truncate its final output to make it to an integer. Not sure these operators would be very helpful though in this specific case.)

Hope this helps.


Frequency response of a system given its block diagram Given the following block diagram, I am asked to find the frequency response $H(f)$ of the system.

This is what I have done: The output of the first block is $j2\pi fX(f)$, and after going through $t_0$, we have $j2\pi fX(f)e^{-j2\pi t_0 f}$. Then, the input of the last block is $j2\pi fX(f)(1-e^{-j2\pi t_0 f})$ and the resulting output is $\frac{j2\pi fX(f)(1-e^{-j2\pi t_0 f})}{j2\pi f}=X(f)(1-e^{-j2\pi t_0 f})=Y(f)$. Finally, we have that $H(f)=\frac{Y(f)}{X(f)}=1-e^{-j2\pi t_0 f}$.

I know that my procedure can be written in a more formal way but, in general terms, is it correct? And if I was asked to plot $|H(f)|$, one way could be computing the square root of $|H(f)|^2=(1-e^{-j2\pi t_0 f})(1-e^{j2\pi t_0 f})=2-2\cos(2\pi t_0f)$, am I wrong?

Your math looks correct, thanks for including what you have done. Also for such a block diagram of a linear system, you can rearrange each of the three blocks in any order at the points where the nodes come together (can't break loops) - For example, you can move the derivative to the end without changing the overall result. This may make it even more obvious that the derivative cancels with the integrator and we are left with the simple delay and subtract.

It helps to have an intuitive view of this filter when confirming the math. The delay and subtract is referred to as a comb filter.  A delay line (the $t_o$ block) does not change the amplitude of the signal but will change the phase proportional to the frequency of the signal coming in: just picture the difference between a very low frequency sine wave and a high frequency sine wave, shift the two the same amount in time and it should be clear that the phase shift is directly proportional to frequency.

The delay and subtract comb filter will have a null at DC or f=0 (this one is obvious as DC will not have any phase shift when it goes through the delay line, and a constant subtracted from a constant is 0). As the frequency increases to the point where the phase is 180°, the signal will have gain of two since the direct path and the delayed path will be in phase.  This occurs when $f=\frac{1}{2t_o}$. As the frequency increases further the phase will reach 360° creating a null again at $f=\frac{1}{t_o}$, and this pattern repeats for every additional multiply of 180°.

One you take the square root of your result, you will see that this "comb" response is the absolute value of a sine where the pattern extends to $f = +/-\infty$:


Matlab `filtfilt` provides excessive transient I observed that filtfilt suffers from an undesired behavior when I provide IIR bandpass filters having steep transition bands. Specifically, the output signal exhibits excessive transient response; nevertheless, this behaviour does not emerge if I use a 'home-made' version of zero phase filtering based on filter. I guess that it is not a numerical issue involving the filter coefficients or structure, nor the filter function.

% design bandpass filter having transition bandwidth of 200 Hz (Fs = 8000)
bp = designfilt('bandpassiir', 'StopbandFrequency1', 50, 'PassbandFrequency1', 250,...
    'PassbandFrequency2', 3600, 'StopbandFrequency2', 3800, 'StopbandAttenuation1', 30,...
    'PassbandRipple', 0.1, 'StopbandAttenuation2', 30, 'SampleRate', 8000, 'DesignMethod', 'cheby2');
% check stability
assert(isstable(bp),'Unstable filter');
% apply filtfilt to a random (white) long input signal; output signal shows an undesirable transient
x = randn(2^20,1);
y = filtfilt(bp,x);
% apply 'home-made' filtfilt to the same input; output signal shows a more acceptable transient
y2 = flipud(filter(bp,flipud(filter(bp,x))));
% compare effects
figure; semilogy(abs(y-y2));

As a rule of thumb, the effect grows as the transition bands get narrower, while it tends to vanish as they get broader.

Where is the problem? Have I missed some recommendations or hints in the function's help?

It must have to do with the initial conditions used by the function filtfilt.m. The idea is to match initial conditions in a way such that startup and end transients are minimized. This, however, doesn't always seem to work, and it appears that for your filter specifications it actually does more harm than good. As far as I know there is no way to tell filtfilt.m not to mess with initial conditions. What remains in this case is to simply implement a straightforward version, just as you did.

Also take a look at this related answer.

EDIT:

I tried to reproduce your results with Octave, but I don't get any of the effects that you describe. I use a $10^{th}$ order Chebyshev 2 filter (which is the lowest order satisfying your specs), and I used a modified version of filtfilt allowing me to both switch on and off the use of initial conditions and the extrapolation of the input signal. There are two explanations for this: either your version of Matlab's filtfilt has a bug, or your filter is substantially different from mine and it introduces numerical problems (due to poles very close to the unit circle), which happen to manifest themselves when the input signal is extrapolated, and/or when special initial conditions are used (as is the case with filtfilt). If you provide your filter coefficients we can probably solve this mystery.

EDIT 2:

Tendero kindly reproduced the filter coefficients using the Matlab commands provided in the question:

b = [0.682585947577945 -0.012906667355225 -3.354515884719406 0.025370108447479 6.651816094305592 -0.000000000000000 -6.651816094305592 -0.025370108447479 3.354515884719406 0.012906667355225 -0.682585947577945];
a = [1.000000000000000 -0.194376467633997 -4.144930396202897 0.592472615834755 7.013126532356848 -0.691483043645839 -6.033921428993976 0.364443043103532 2.633924737223088 -0.072944768631229 -0.465923573748158];

I used these coefficients in Octave with the function filtfilt.m with a random input signal x = randn(2^20,1) just like in the question. As mentioned above, I used a modified version of filtfilt.m to be able to switch on and off the different measures for reducing transients:

no extrapolation of the input signal, zero initial conditions
no extrapolation of the input signal, optimal initial conditions (as implemented in the standard version of filtfilt.m)
extrapolation of the input signal, zero initial conditions
extrapolation of the input signal, optimal initial conditions

The last combination is the processing implemented in the standard version of filtfilt.m.

The plot below shows the first and last $50$ samples of the $4$ different output signals. For all other samples the signals are virtually identical.

As expected, one can see slight differences during startup and at the end, but these differences are very small, and the problem of the OP cannot be reproduced using the original filter coefficients and the Octave function filtfilt.m. There must be either a bug in the Matlab version of filtfilt.m used by the OP, or another problem that is neither related to the filter coefficients nor to the function filtfilt.m.


Signal processing techniques for an accelerometer signal? I am running some tests where I am recording accelerometer measurements.  I am looking to use elements of signal processing on this signal, but I am unsure about where to begin, or what my approach should be.

My ultimate goal is to be able to monitor the acceleration readings in real-time, and then display a notification when the event occurs.  As you can see around the 150,000 sample time, an event occurs.

If I am monitoring this data in real time, what sort of signal processing techniques could be implemented to react to this event?
Would a Short-Time Fourier Transform (STFT) be an option?

I am monitoring my data in Python, and they have a decent STFT function.

The arguments of this function are as follows:

scipy.signal.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None,
detrend=False, return_onesided=True, boundary='zeros', padded=True, axis=-1)

How do I determine optimal parameters to use to process this signal?
Are there any other methods that you folks think may help me in identifying when the event occurs in real-time (as opposed to just using the magnitude of the acceleration)?

EDIT 1:

My STFT has been added above.

I'm wondering why the STFT pops out. To me, wouldn't a simple threshold on the signal itself or on its envelope do better / just as well, after removal of the g offset?
Once you decide what "measure" is best to detect your event, you can apply the work of Basseville and Nikiforov, that I answered here.

The classic reference for that problem is Detection of Abrupt Changes - Theory and Application by Basseville and Nikiforov.  The whole book is available as a PDF download.
My recommendation is that you read Chapter 2.2 on the CUSUM (cumulative sum) algorithm.


Concept of Orthogonality and Orthonormality Please explain me the concept of Orthogonality and Orthonormality. I understand that the mathematical definition of orthogonality between two vectors is that their dot product is zero. It means that there is no correlation between the two vectors. And in case of Orthonormality the vector length would be 1.

I am unable to understand the concepts clearly. Please explain me the concept of Orthogonality and Orthonormality in terms of Image Processing.

Orthogonal means that the inner product is zero. For example, in the case of using dot product as your inner product, two perpendicular vectors are orthogonal. Orthonormal means these vectors have been normalized so that their length is 1.

Orthogonal vectors are useful for creating a basis for a space. This is because every point in the space can be represented as a (linear) combination of the vectors. So for example in 3D space, x=[0,0,1]  y=[0,1,0] and z=[1,0,0] form an orthornormal basis. No component of x can be represented with a component of the other vectors. This is because they are linearly independent

This is a different type of independence to statistical dependence however.

Correlation is a different notion. There are different ways of representing the correlation of two vectors (random variables X and Y).

$corr(X,Y) = \frac{cov(X,Y)}{\sigma_x \sigma_y} = \frac{E[X - \mu_x]E[Y - \mu_y]}{\sigma_x \sigma_y}$

Covariance is a measure of how two variables change together.

Another measure of correlation is cross-correlation. This is measure of the similarity of two functions (or vectors).

These concepts ares used in many different ways in image processing.

For example cross-correlation is used in template matching. When looking for a small image inside a much larger image, the small template image is 'slided' over the large image and the cross-correlation is computed for each position. Locations with a high cross-correlation are likely to contain the image we are looking for.

The concept of Linearly independent component vectors is used in principal component analysis (PCA). PCA takes a cloud of points in space and calculates a set of orthogonal basis vectors to represent them.

Independent component analysis (ICA) is used for separating mixed signals. ie: separating a speaker at a noisy cocktail party. ICA uses properties of the signal correlation to separate the signals.


Sample Frequency - Recommended Practical Nyquist, Lower Used by majorities appications Doubt1: What is the Lowest Sampling Frequency used by majorities sound applications. Is there limited frequency sampling lower than 8KHz In C/C++?

Comment In Java I haven't detected any problem, but using VCL program into one wav sound file (with lower Freq to 8K) was not played, producing Error!

Doubt2: Can I to capture sound using arbitrary Frequency Sampling in C/C++? This process/function/method depends on specific Library? Can I to do this using pure C/C++ Code (without specialized Libraries), where?

Maybe this is a foolish question. I have sound with the 100Hz-500Hz Spectrum.

Doubt3: What is the recommended lowest Nyquist Sampling Frequency (I Know that teorical is two), but is there another practical Criterion? Example: Is recommended to use 2.5, 3 or 4.
In my sound the TEORICAL Nyquist sampling Frequency must be 1KHz.

Thank you

44.1 kHz and 48 kHz are sampling frequencies that are almost always supported by software and drivers but if you need to go lower than that then stop at 8 kHz, as it is commonly supported.


How to know if a audio file is real lossless using its spectrogram? I need to know if a certain audio file is real lossless and not a transcode (i.e. lossy converted in lossless, but lossy at the end). I was told there was a way to do this by creating a spectrogram of the file and then checking something else. The problem is that, since I am not well versed in audio processing, I don't know what else to check.

I would like to learn how to identify this fake and true lossless files using a spectrogram (since, as a Linux user, I found a program which creates them easily). If you want to give me a hand with the file in question, here it is:

P.D.: I'm not completely sure if this question belongs to this site. I put it here since is a site about audio, but if not just tell me and I will migrate it were it belongs.

Eventhough there could be opposing theoretical answers to your question, practically speaking it may have a solution which depends on the lossy codec you are using. In that case what you eventually need is a practical knowledge of codec under concern that is used for the audio clip. However, considering the number of different audio codecs and the fact that there could be any number of nonstandard implementations to exist, your problem may not have a guarantied solution. It is because for every different lossy auido codec a different procedure (provided in theory it does exist) might be necessary to test whether a given piece of audio belongs to its output or not. And since for a nonstandard/proprietary/unknown audio codec, you will not have such a procedure at hand you may not even be able to test to see it.

For example, mp3 (MPEG-I layer 3) standard defines an optional low-pass filter about 16khz, before encoding the raw pcm audio, hence if you find some signal acitivity above 16 khz inside an audio clip, then it is either not an mp3 file or the mp3 profile used to encode that clip, did not activate the lowpass prefilter (as it can be one of the many optional features of the codec) Note that in either case you cannot say that it is a raw file just because it's not an mp3 as it might come from another lossy codec... That is the problem...

Finally, if a lossy codec will include a distinctive feature as a mandatory (nonoptional) requirement of the standard, then you can check for that property to at least decide whether the given piece may belong to its output or not.


Detecting lines with houghlines in matlab I have This image which is basically made of sloped lines (almost horizontal and parallel)

I tried exactly the code in matlab documentation to detect those lines (except the rotation part) but I got a very strange result (vertical lines detected and no single horizontal line!):

I really can't find explanation of this output. I tried the same code on a different image and it was ok. I tried to change all the possible parameters in houghlines and houghpeaks but I always get vertical lines even though my image has no vertical lines or anything similar to be vertical!

I need to detect the parallel lines in the image or somehow get their slope. If you think there's something wrong in my image or additional modification I should do to the example code. or if you have any better idea that would be great!

The biggest issue is that the pixels for each line are too few (I will explain with more details below). I would consider to stretch, and then dilate your raw image a little bit:
file='http://i.stack.imgur.com/LmIJJ.png';
I=imread(file);
I=imcrop(I,[1 206 size(I,2) size(I,1)]);
I=imresize(I,[size(I,1) 256]);
I=imdilate(I,strel('line',1,0));

I assume you are using the sample code from Mathworks on hough, houghlines, and houghpeaks. The problem is that there are not very many pixels along one slope in your image. You feel like it is a line, but it is actually a curve that consists of many short lines with different orientation. I think one parameter that you need to tune is the 'FillGap' which will connect the very short segments. I tried:
lines = houghlines(BW,T,R,P,'FillGap',30,'MinLength',15);

And the image I obtained:

There is still plenty of room for improvement, yet at least most lines are detected.
Note that the theta-axis in hough peaks indicates the slopes of the lines, you may also check them from lines, T, and R. But since the image has already been stretched along the 2nd dimension of raw data, you may need to scale them back (with some tangent transformations) if you plan to get the slopes of those detected lines.
Hope it helps, thank you.


How to calculate power of a signal at one specific frequency? Given a signal, $s(t)$, sampled at frequency $f_s$, how does one calculate the power of the signal at one specific frequency $f_q$?  Since I am only interested in the power at $f_q$, computing the FFT seems like over kill.
I can identify the following pre-existing content on stackexchange:

This post effectively suggests Goetzel's algorithm.  There is a pre-existing MATLAB function which appears to provide what's needed.

This post suggest band-pass filtering the signal then taking the sum of the squared signal.  This seems to make sense to me but I wonder how selective the band-pass filter can actually be given the sync-like frequency content of the windowing function.

This post deals with power in a specific signal band with suggestions of using MUSIC, goertzel, etc.

My interest is not in a frequency band but the power at one specific frequency.  My thought is to convolve a sinusoid whose frequency is fq with the signal and then summing the power in the resultant signal (a matched filter).
I would have thought there would be a more straightforward solution but I struggle to find the words to start an effective literature review.

Ultimately if your noise is white then the matched filter would be the best approach (multiply by one frequency and sum-- if it is complex you would multiply by the complex conjugate- both a form of correlation). This is equivalent to computing a single bin in the DFT (assuming your frequency is on bin center, otherwise equivalently the operations in the DFT can be modified to be centered on your particular frequency). This approach has bandwidth which for white noise will have the equivalent power of a similar brickwall filter with a bandwidth of $1/T$ where $T$ is the duration of your correlation; but importantly the frequency response is a Sinc function (or Dirichlet Kernel if discrete as a sampled Sinc) in frequency which means this approach is sensitive to non-white noise at other frequency locations (a strong tone very relatively far away would show up in your power measurement); the main lobes of the Sinc only roll-off at 6 dB per octave in frequency.  So for white noise, this approach is optimum, for other applications with colored noise, applying a properly designed filter to the signal to attenuate "out of band" noise before measuring the power would be recommended.
For measuring a single tone with the equivalent of a single DFT bin, the Goertzel is a well known choice that is efficient at doing this very thing.


Sampling rate issue in microcontroller I am using a 3-axis LIS2D12H accelerometer interfaced with NRF51832 through the SPI interface.

The sensor has a sampling-rate from 1 Hz to 5.376 KHz which is happening at sensor only.
I am taking sensor data at 400 Hz by setting register values in the sensor.
I am implementing an FFT algorithm in microcontroller over a raw sensor data received by the 3 axes, for vibration analysis.

I can have resolution of FFT from 64 to 2048.

Q: My confusion is how should I sample it at microcontroller for FFT? There are FFT functions in ARM cortex which I am using.

Shall I take raw samples in buffer and sample them or shall I take each reading from sensor and sample them individually?
If that is the case how shall I select sampling frequency at microcontroller level and/or at what rate microcontroller shall read data?

I understand the sampling process is completely handled by the sensor itself, i.e. it delivers digital data. Then you just collect the desired number of samples from the sensor and calculate its FFT. The frequency range depends only on the sampling rate ($f_{\text{max}}=f_s/2)$, the frequency resolution depends only on the number of samples $N$, in your case $N\in[64,2048]$.

There is no sampling done on the uC, as it already receives sampled data. However, what you will want to do is employ a lowpass filter before calculating the FFT, with a cutoff frequency slightly lower than $f_{\text{max}}$, to get rid of high frequency content, that would otherwise cause aliasing in your spectrum.