Time invariance of a System I have this small question about the time invariance of a system. Which is:

  If the current output is multiplied by the current input (see both are variables) will the system be time variant or time invariant?

To be precise, the question I am dealing with is:
$$ y(n-3) + y(n-1) + y(n)x(n) = x(n-3) $$

I would be glad if anyone could help.

Yes, this system is time-invariant assuming $x[n]\neq0 $ for all $n $.

The test for this is
$$x_{1}[n]=x[n-n_{0}] $$
$$y_{1}[n]=y[n-n_{0}] $$

So let's first express in terms of just $y[n]$.
$$ y[n]=\frac{x[n-3]-y[n-3]-y[n-1]}{x[n]} $$
Next, we go through the test.
$$y[n-n_{0}]=\frac{x[n-3-n_{0}]-y[n-3-n_{0}]-y[n-1-n_{0}]}{x[n-n_{0}]} $$
Substituting $y_{1}[n] $
$$y_{1}[n]=\frac{x[n-3-n_{0}]-y[n-3-n_{0}]-y[n-1-n_{0}]}{x[n-n_{0}]} $$
$$y_{1}[n]=\frac{x[n-3-n_{0}]-y_{1}[n-3]-y_{1}[n-1]}{x[n-n_{0}]} $$
And finally substituting $x_{1}[n] $
$$y_{1}[n]=\frac{x_{1}[n-3]-y_{1}[n-3]-y_{1}[n-1]}{x_{1}[n]} $$
Because the shifted sequence has the exact same relationship, it is said to be time-invariant. Most systems that don't alter $n$ or $t$ meet this. Be wary of anything multiplying or otherwise messing with the vector arguments beyond simple delays.


Are there frequency detection methods that are subject to less latency than the Goertzel algorithm? I am trying to detect oscillations in a small range around 1.5Hz as early as possible. I am currently using the Goertzel algorithm with bins 1.0Hz, 1.2Hz,..., 1.8Hz, and 2.0Hz; with block size N = 500 and sampling rate Fs = 100 (the sampling rate, in my circumstance, cannot be modified).

I am applying the Goertzel algorithm on two signals, and checking for three conditions that, if met, result in the oscillation I'm concerned with being detected:

amplitude of signal 1 > threshold_1
amplitude of signal 2 > threshold_2
phase difference between signals > threshold_3

If all conditions are met, a flag is raised, and that is the purpose of this whole oscillation detection endeavor.

This all works fine, except that the Goertzel algorithm is inherently subject to latency in accordance with N. I am wondering if there are any other methods that do better than the Goertzel algorithm in terms of latency? Or if there's any way to reduce the latency of this implementation while still confidently detecting oscillations around 1.5Hz. I suppose I could reduce N but I'm not sure I want to be increasing the bin width.

I might be mentioning this too often, but if you need to detect frequencies finer than $f_\text{sample}/N$, then superresolution frequency estimators might be the way to go.

In your case, an algorithm called MUSIC would sound promising. $N$ would in that case be the amount of samples you use to estimate the autocorrelation matrix that algorithm works on, and can be as small as $n_\text{tones}+1$, but realistically, don't go too low, and average the estimate a bit.


Phase response distortion due to IQ demodulation I have a bandpass signal, I do IQ demodulation meaning:

first step:

$IQ(t) = x(t)*exp(-2\pi f_ct)$

second step:

low pass using $filtfilt$ function in MATLAB. (in order to avoid phase distortion)

Here, $F_c=5 MHz$ and $F_S=40 MHz$

For low-pass filter I use this:

[b,a] = butter(5,2*Fc/Fs) ;
IQ(t) = filtfilt(b,a,IQ(t))*2;

The plot shows frequency response of $x(t)$, $real(IQ(t))$ a.k.a $I$, $imag(IQ(t))$ a.k.a $Q$ and $IQ(t)$, respectively.

When I reconstruct this signal again to the original bandpass signal:

$x'(t) = IQ(t)*exp(+2\pi f_ct)$

The magnitude, phase response and power spectrum of the original and reconstructed signal is a follows:

note: To plot phase response I use matlab's $unwrap$ function.

The frequency spectrum of original and reconstructed signal:

The phase response of the $x'(t)$ seems distorted at $2f_c = 10MHz$ .

Question1: Is this normal?

Question2: Does this happen because low-pass filter is not reversible?

Question3: Can I recover that phase response?

I’ll hazard a guess on this that David’s comment is spot on.  You could help by providing some specifications for the low pass filter attenuation at the frequencies at which you are observing the phase distortion, or perhaps provide plots of your magnitude/phase responses for the input and output. If this is what you’re observing then:

1). Yes. For a DFT or FFT of a signal, there is likely to be phase distortion at low relative amplitudes based on different implementation details. This has to do with the precision/resolution of the data on your machine.

2). Sort of.  The filter could technically be reversible in the analytical sense (ie on paper) but quantization error from the limited precision of the data will make this impossible with many many applications.

3). Recover?  No, not really.  Assuming that you don’t have any errors in your implementation, and it sounds like you don’t, then that is simply the nature of the beast.  However, you are using filtfilt, which I understand works by applying the filter forwards and backwards, resulting in a net zero phase response.  As such, it would follow that your phase response should be unchanged, which I assume was your original hypothesis. You could just use the phase response from the original signal, if you have it.  This may or may not be practical for your application, and could be done more efficiently if it were an option.


How to produce a high-pass filter from a low-pass one? I already have a low-pass filter and want to remove the DC offset of a signal with a high-pass one.  I guess you can subtract a low-pass filtered signal from the original signal and the result is a high-pass filtered signal?

My low-pass filter:

x = input, raw signal
y = output, low-pass filtered signal
n = current sample
n - 1 = previous sample

y(n) = y(n - 1) + (x(n) - y(n - 1)) / alpha

High-pass filter:

x = input, raw signal
y = input, low-pass filtered signal
z = output, high-pass filtered signal
n = current sample

z(n) = x(n) - y(n)

Is this last equation correct?

The result will indeed be a high pass filter. From your difference equation, the transfer function of the low pass filter is

$$H_l(z)=\frac{\beta}{1-(1-\beta)z^{-1}}\tag{1}$$

with $\beta=1/\alpha$.
Note that this is actually a leaky integrator, not a classic low pass filter, because its frequency response does not have a zero at Nyquist.

The high pass filter has the transfer function

$$H_h(z)=1-H_l(z)=(1-\beta)\frac{1-z^{-1}}{1-(1-\beta)z^{-1}}\tag{2}$$

which has a zero at DC (i.e. at $z=1$), as it should. However, it is scaled, i.e. its value at Nyquist ($z=-1$) is not $1$, but $2(1-\beta)/(2-\beta)$. So in order to have a 0 dB gain at Nyquist, you need to scale your high pass filter with the inverse of that factor:

$$H'_h(z)=\frac{2-\beta}{2(1-\beta)}H_h(z)=\left(1-\frac{\beta}{2}\right)\frac{1-z^{-1}}{1-(1-\beta)z^{-1}}\tag{3}$$

EDIT:

You can compute the 3dB cut-off frequencies of the low pass filter (1) and of the high pass filter (2) by solving

$$|H_l(e^{j\omega_c})|^2=\frac12\quad\text{and}\quad |H'_h(e^{j\omega_c})|^2=\frac12$$

After some algebra this gives for the low pass filter

$$\cos\omega_c=1-\frac{\beta^2}{2(1-\beta)}$$

and for the high pass filter

$$\cos\omega_c=\frac{1}{1+\frac{\beta^2}{2(1-\beta)}}$$

where $\omega$ is the normalized frequency in radians:

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

with the sampling frequency $f_s$.


Does our brain constantly ifft to hear? From my understanding our ears have hairs/cilia in the cochlea that resonate at frequencies within our hearing range.  To me this means we are hearing in the frequency domain as opposed to the time domain. But does that mean we are constantly ifft-ing the sound in order to process it?

No. Why do you think it would? First of all, the human brain works very different then any human constructed computer (to date); so the assumption that it runs mathematical "algorithms" is somewhat dicey..

Here is roughly how it works:

Sound wiggles the air drum
That vibration is transferred by the ossicles to the cochlea. The ossicles  act mainly as a mechanical transformer.
Inside the cochlea the vibration induces a bending wave on the basilar membrane. Due to its shape certain frequencies wiggle preferably at different spots. That's how the frequency selectivity is created. In essence this is a mechanical way of doing a short term Fourier transform.
There are special type of neurons (nerve cells) called hair cells or cilia in the cochlea. When the basilar membrane wiggles it shears the hairs and induces the neurons to fire electrical impulses.
Interestingly enough about 20% of the neurons go from the brain back to the cochlea. This implements an active feedback mechanisms that greatly enhances the frequency selectivity. If that feedback loop goes out of whack you can "hear your ears ringing". That's one of the main causes of tinnitus.
Neurons are inherently somewhat digital. They are amplitude discrete (i.e. they fire or not) but they are time continuous. The more stimulus is there, the more frequent they fire, but the size of the impulse stays always the same.
From there on upward it all stays in the neural domain. Everything that can be observed are firing patterns of active cells. Every layer of the brain is connected to the next layer by a huge number of connections which transform one firing pattern into a different one. Certain firing patterns mean that you actually "hear" something and are related to conscious perception.
The human brain has a huge number of active neurons and glia. There are more than 100 billion cells in your head that can process information. Needless to say this is a very complicated process.
However, there is no need for the brain to reconstruct the original time domain waveform. It actually has no good way of representing a fast analog waveform. Neurons are pretty slow with a maximum firing rate of about 100 times per second. That's really not a good match for sound signals.


RF mask (power spectrum measurement bandwidth) I'm reading an IEEE standard and it says that the first channel needs to be 50dB below the highest power in a 6 MHz operating channel in 100 KHz bandwidth. It was mentioned several times that the power spectrum measurement needs to be done over a 100 KHz bandwidth. Isn't it sufficient to just say that the power of the first channel should be 50dB below the highest power. I guess my question is what is the significance of the measurement bandwidth, such as the 100 KHz. Why do we need to specify a measurement bandwidth?.

The simple answer: when you take a power measurement using some piece of test equipment (such as a spectrum analyzer or an RF power meter), the reported value indicates the amount of power measured within some detection bandwidth of the instrument. For some devices, that might be a fixed value, determined by the filtering present in its RF front end. Other types of equipment, like a spectrum analyzer, will typically have adjustable bandwidth (which in this case would usually be called the instrument's resolution bandwidth).

To see how the bandwidth setting affects the measurement, try injecting a modulated signal into a spectrum analyzer and see how the power measurements vary as you change the resolution bandwidth: you can get very different answers. Take a common $\frac{\sin{x}}{x}$ spectral envelope as an example:

If you use a very small measurement bandwidth, then you're taking power measurements over very small swaths of frequency. This allows you to observe the detail in the spectrum's shape.
As you increase the resolution bandwidth, your instrument will take more and more of the signal's spectrum into account for the power measurement, and the reported values will become more of an overall measure of the total received power.

As was indicated in the comments already, the specification calls out a specific bandwith to ensure repeatability. Specifications are used to validate implementations as conforming to an existing standard. This particular scheme for measuring the emission mask of a transmitter implementation is likely in use because that was the one that was chosen when the standard was designed. Due to the imprecise definitions of concepts like SNR and bandwidth, there are many ways that you could define such a requirement, but this is the one that the authors chose. That doesn't make it right or wrong; it just is.


Determine causality and stability from given filter structure I have the diagram above. I found the transfer function below from it;

The question asks me to find out if the system is causal and stable, but didn't it have to indicate whether it was left-sided or right-sided?

It looks like a Homework problem so I will only provide small hints:

Check from the block diagram implementation of the system whether the current output sample $y[n]$ depends on only current input sample $x[n]$ and past input & output samples $x[n-n_0]$ and $y[n-m_0]$. If so, the system is causal, else it will be non-causal.

You can use the information from the previous point (causal or non-causal) along with the pole location of $H(z)$ to check if the system is BIBO stable or not.

Any Causal System with poles inside unit-circle are BIBO stable.
You can refer to this Pole-Zero Analysis for better understanding.


Software-Emulated Amplifier: Cracking while change Volume I'm developing a software emulated analog synthesizer in the style of Euro Rack Modules. I have a problem with the volume amplifier.
Whenever I control the volume button via mouse/hand/manually rapidly (fast low and up, or fast from low to middle), I hear lots of little cracks while changing the volume. The confusing fact is, when controlling the button via automatism, like via LFO (Low-Frequency Oscillator), all is fine.
I already checked if all values are within the valid ranges (floating between -1 to 1).
Does anyone have an idea why this occurs and what I could do to avoid this?
Here's a short demonstration:
https://youtu.be/5PXcoTd32RM

I’ve heard it called popcorn noise before.  Think of the output signal as the product of the input signal and the volume signal. A discontinuity in the volume signal causes a discontinuity in the output signal. The solution is to limit the bandwidth of the volume signal.  The most common implementation I’ve seen is to limit how fast the volume can change.  Another option would be to feed the volume signal through a low pass filter.  Note that the output of an LFO is already (pretty much) band limited, which explains why you wouldn’t hear any noise in that case.


Are there any template matching / optical flow algorithms for large image displacement? I had some success with feature-based methods for image registration, but these methods are sometimes unsuitable (e.g. images containing repetitive patterns, very low contrast, low texture, no edges/corners).

I am therefore looking for a direct (pixel-based) global alignment alternative.

So far I looked on algorithms like Horn-Schunck, Lucas-Kanade, ECC but all the examples show only images with very small displacement (few pixels at most).

Is there any direct (pixel-based) method for registering images with a displacement of 200 pixels or more (in any direction)?

So far I am using simple template matching with exhaustive search over possible image shifts and Normalised Cross-Correlation error metric, which works, but it is extremely slow (basically $O(n^4)$) and does not allow for rotation.

Are there any faster methods allowing large image shifts/rotations?

I personally like this one. It is nicely designed for occlusions and large displacements.

But, recent trends in deep learning lead to better results in optical flow such as this and this.


How to plot Full spectrum with negative frequencies I'm trying to do frequency analysis on signal of a rotating shaft. The signal is generated in two measurement planes i,e X-position and Y-Position. I was able to use fft and transform into frequency domain for individual planes. The plot I get for each place is only Half spectra. What I'm looking for is a full spectra which could combine both the planes. I did some research and found the following description online.

Full spectrum: The full spectrum is an additional diagnostic tool and is also called the spectrum of an orbit. It shows the same information as an orbit but in a different format. It helps to determine the degree of ellipticity (or flattening) associated with the various machinery conditions along with the precessional direction for all the frequency components present. To obtain the full spectrum, the orthogonal X and Y transducer signals are fed into the direct and quadrature parts of the FFT input. The positive and negative vibration components for each frequency are obtained. Positive is defined to be the forward precession and the negative component as the reverse precession. These components yield the following ellipticity and precessional information for a given orbit of any particular frequency (1× or 2× or …):

• The sum of two components, forward and reverse, is the length of the orbit major axis. • The difference between the two components is the length of the orbit minor axis. • The larger of the two components, positive or negative, determines the direction of precession that is forward or reverse.

One of the possible applications of full spectrum is analysis of the rotor runout caused by mechanical, electrical or magnetic irregularities. Depending on the periodicity of such irregularities observed by the X–Y proximity probes, different combinations of forward and reverse components are observed. The method forms the basis for many useful machinery diagnostics.

The full spectrum (just like the normal FFT) can be obtained in a steady-state analysis (a single FFT or waterfall) and even in transient analysis, which would then be called the full spectrum cascade (Figure attached).

Could anyone help how the 2 half spectra is converted into full spectra plot. Does anyone have any idea about this? Any hints or suggestions please!!

Take the fft of x+i*y.

The phase angles in the picture you linked https://drive.google.com/file/d/0Byib6yrCQ9vnRUM2Yi1sMklXeFk/view are from the result of this fft.

The image linked have an error, from the article it came from, where the phases of the negative frequencies of the fft result are Beta. The image uses Alpha in the negative and positive results. In the end, the image just shows mathematical prove that the amplitudes of the fft are derived from the amplitudes of the individual signals.

The fullspectrum process is simple as the fft of x+i*y. The majoraxis of an ellipsis of one frequency will be the sum of the R+ and R-, and the direction will be the (Alpha+Beta)/2.

This paper could clarify it more "Use of directional spectra of vibration signals for diagnosis of misalignment in rotating machinery"


Sampling vs repetition code I'm working on a multiuser wireless system simulation and have run into something regarding sampling rates. From everything I've known for a matched filter receiver, the correlation step is performed followed by symbol rate sampling. Say the symbol rate is $T_s$. My question is: do we ever sample at higher than the symbol rate? What if I was able to sample at $\frac{T_s}{5}$? It seems to me then I'd have five samples of the same symbol, which is in the same spirit of repetition coding? Is there a difference between repetition coding and just sampling faster?

You are mixing up notions of matched filtering (e.g. using an LTI system) with a correlation receiver implementation of the same.  The correlation receiver is an "integrate and dump" device -- a time-varying circuit -- whose output resets to zero (that's the dumping part) immediately after the output is sampled at the correct time (once each symbol interval) so that the correlator starts with initial condition $0$ for each symbol correlation. Note that the correlation receiver output does not equal the matched filter output except at the sampling instants.  A simple illustration of this point is available towards the end of my answer to
Understanding the Matched Filter

So, if you sample the correlation receiver output multiple times (five times in your case) during one symbol interval, the question arises as to what happens to the correlator after each sampling?  Is it dumped after each sample, or does it continue to accumulate the correlation until the final sample is taken and the correlator is dumped right afterwards?  In the former case, the samples are independent random variables while in the latter case, they are dependent.  Also, how do you propose to process these samples?  In the former case, the optimum thing to do would be to add up the five samples (which has the same result as running the correlator over the entire symbol interval and sampling just once at the end of the interval) while in the latter case, the optimum thing to do would be be discard the first four samples and just keep the last one.  Here, optimum is with respect to the signals being received in white Gaussian noise; if your noise is something else, then you need to specify what it is, and to what extent is the concept of matched filtering or correlation even justified as the right thing to do.


kernel to calculate second order partial derivative of digital image im working on image stacks, and I need to calculate second order partial derivatiatives of it.
I already know how to calculate derivative on x and y axis, using coefficients calculator: https://web.media.mit.edu/~crtaylor/calculator.html
So for example to calculate five point second order derivative on x axis we have formula like:

So I already know how to calculate derivatives on x and y axes, but I also need to get formulas how to calculate first and second order partial derivatives for f'xy and f''xy. I know that I could simply calculate derivatives on x axis first, then calculate from it on y axis, but I would like to obtain simplified formula that uses x and y pixels in image to calculate f'xy/f''xy.
Could anyone explain me how to obtain them?
Update:
I found somewhere on the web the formula for second order f''xy derivative for images that looks like that:
(inX[2]-1inX[0]+inY[2]-1inY[0])/4;
Would like also to get an explanation whether its correct, and how its acquired

Your five-point derivative kernel is a 1D kernel. Applied along the x axis it gives the partial derivative for x, applied along the y axis it gives the partial derivative for y.
The $\frac{\partial^2}{\partial x \partial y}$ derivative would need a 2D square kernel. It is more efficient to compute this by applying two first order partial derivatives.
$$\frac{\partial^2}{\partial x \partial y} f = \left( \frac{\partial}{\partial x} \frac{\partial}{\partial x} \right) f = \frac{\partial}{\partial x} \left(  \frac{\partial}{\partial x} f \right)$$
Since $\frac{\partial}{\partial x} f$ is a convolution,
$$\frac{\partial}{\partial x} f = k_x \ast f$$
we can write the earlier equation as
$$\frac{\partial^2}{\partial x \partial y} f = \left( k_x \ast k_y \right) \ast f = k_x \ast \left(  k_y \ast f \right)$$
So we can say that the 2D square convolution kernel that you are looking for is $k_x \ast k_y$ (with $k_x$ the kernel that computes the x-derivative and $k_y$ the one that computes the y-derivative). Note that, if you think of your 1D kernels as vectors, this convolution is equivalent to the outer product of the two vectors.
However, it is more expensive to compute a convolution with one 5x5 kernel (25 multiplication and additions per pixel) than to compute two convolutions with 5x1 and 1x5 kernels (10 multiplications and additions per pixel). I would compute two first order derivatives in sequence.


Is there a linear filter whose step response is an S-curve? Is there a linear filter whose step-response (or ramp-response) is an S-curve?

I think the properties I'm looking for here are (1) time-reversal symmetry, and (2) at least a continuous first derivative.

Sums of harmonic cosines

The Hann window is a sum of a rectangular function and a truncated cosine. As an impulse response it has the qualities you asked for:

Figure 1. Hann window as impulse response, and the corresponding step response. The step response is valued $x - \frac{\sin(2\pi x)}{2\pi}$ over $0 \le x \le 1$, and is constant outside that range, because the impulse response is zero outside that range. The step response and its first two derivatives are continuous, because the impulse response and its first derivative are continuous.

A filter with Hann impulse response can be constructed as the (normalized) sum of a truncated integrator filter and the real part of the lowest-non-zero-frequency bin of a sliding discrete Fourier transform, both of which can be calculated recursively as truncated IIR filters.

The filters have poles on the z-plane unit circle (they are marginally stable) and require periodical refreshing of their state to cut accumulation of round-off errors. For fixed point input the truncated integrator can be implemented with extended bit depth so that no rounding errors accumulate and refreshing its state is not necessary.

Octave source for the composite filter:

r = 20;
w = 2*pi/r;
b0 = zeros(1, r+1);
b0(1) = 1;
b0(r+1) = -1;
a0 = [1, -1];
b1 = zeros(1, r+2);
b1(1) = -1; b1(2) = cos(w); b1(r+1) = 1; b1(r+2) = -cos(w);
a1 = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
ir = filter(b0, a0, impulse) + filter(b1, a1, impulse);
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/r #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");

Figure 2. Hann impulse response from the first Octave script.

It is possible to get rid of the unnecessary one-sample delay without additional cost by multiplying the complex transfer function of the sliding DFT bin by $e^{-i\omega}$ before taking the real part:

r = 20;
w = 2*pi/r;
b0 = zeros(1, r+1);
b0(1) = 1;
b0(r+1) = -1;
a0 = [1, -1];
b1 = zeros(1, r+2);
b1(1) = -cos(w); b1(2) = 1; b1(r+1) = cos(w); b1(r+2) = -1;
a1 = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
ir = filter(b0, a0, impulse) + filter(b1, a1, impulse);
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/r #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");

Figure 3. Undelayed Hann impulse response from the second Octave script.

There are other window functions with harmonic cosine terms, and filters with them as impulse responses could be implemented similarly.

B-splines

Piece-wise polynomials such as B-splines are another alternative for the impulse response shape. Cascaded integrator–comb (CIC) filters can be used to implement them. A triangle function would be the minimal piece-wise polynomial function (and B-spline) to meet your demands. The implementation is very simple: A truncated integrator in cascade with an identical truncated integrator. This comes from the favorable property of B-spline functions that they are each a convolution of an impulse by rectangular function multiple times. The number of continuous derivatives is incremented by each convolution. The earlier notes about stability apply.

Figure 4. A triangle impulse response and the corresponding step response.

Octave source:

r = 20;
b = zeros(1, r+1);
b(1) = 1;
b(r+1) = -1;
a = [1, -1];
length = 2*r+10;
impulse = eye(1, length);
ir = filter(b, a, filter(b, a, impulse));
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/(r*r) #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");

Figure 4. Triangular impulse response from the Octave script.


How to removed ripples cause by zero-padding in fft? Desired outcome: Interpolation in frequency domain
Attempts: Zero padding when doing fft
Problem: Ripples in frequency domain appeared after fft with zero padding
Code to reproduce the problem:
L0 = 1024;
t = linspace(-10*pi,10*pi,L0);
yt = sin(t)+cos(2*t)+0.1*randi(10,[1 L0]);

for i = 0:4
    mul = 2^i;                                 %zero padding multiplier
    L = L0*mul;
    T0 = (max(t)-min(t))/(L-1);                %period
    f = (0:L/2)/T0/L;                          %frequency axis

    yc = fft(yt,L)/L;
    yc = yc(1:L/2+1);

    figure
    plot(f,abs(yc));
end

Thanks in advance !

Windowing the time domain function prior to zero padding will reduce the ripples caused by zero padding (due to the otherwise sharp truncation of the time domain function). Without additional windowing you are equivalently multiplying your longer time domain function with a rectangular window, which Fourier Transform is a Sinc function. To reproduce this exactly in frequency the frequency must extend to infinity, which is doesn't in the discrete system so aliasing of these tails in frequency occurs, resulting in ripple. Choose a window whose Fourier Transform has tails that are much lower at the aliasing frequencies and the ripple will be significantly reduced! (I personally like the Kaiser window as you can adjust the trade between width of the main lobe (frequency resolution) and depth of the tails (dynamic range) by adjusting the Beta factor in the window construction.


Least Squares with Non Zero Mean Noise What happens if the noise has no zero mean? I mean, if the exercise is something like:
$$y(k) = A x + \eta(k)$$

When I have zero mean, I start from:
$$y = A x$$
$$\Rightarrow \hat{y} = A \hat{x}$$
Using algebra, I to get to this equation:
$$A^H y = A^H A \hat{x}$$

But what happens if it has no zero mean? I have to use the following inner product, I suppose:
$$\langle x,y \rangle = \mathbb{E}[(x-\mathbb{E}[x])^H(x-\mathbb{E}[x])]$$
But I can't see how to get to an equation from that.

Thanks!

Since this is a linear model if you add noise which isn't centered (Non zero mean noise) your estimation will be good up to a bias term.

The easy way to do so is to remove the bias from $ y $ and solve.

In case you know the mean of the added noise, just remove it from your measurement and your model becomes the classic Linear Least Squares.

If you don't know it, you have to estimate it. You can solve the problem either in 2 steps, namely estimate, remove, solve classic Linear Least Squares, or build a model which estimates the noise and solve the least squares at once.


Time Domain Example of Nyquist/Shannon? I looked at an old thread, Nyquist Frequency Phase Shift, but was left wondering still about oversampling.

I learned Shannon 40 years ago, then worked 30 years in (s/w) engineering, but now retired I have only begun using/configuring complete digital audio setups. I find myself baffled by exactly the sort of hands-on question that was posed in the referenced thread. Yet most of the relatively simple (IMO) questions I've thrown at the search engines have landed on generic "Nyquist says" posts, and almost all graphs have been frequency domain.

Like the author of the referenced post, if I draw a simple sine wave and "sample" it with something just over twice its frequency, then vary the phasing of the sample points (or the wave) a few times, I end up with very different squarewave patterns. Especially for an illustratively short burst, it's impossible for me to deduce how even the most perfect linear-phase brickwall filter could regenerate the exact analog input (amplitude plus phase), which would be critical in preserving harmonic structure for any instrument. Certainly Nyquist had little to draw upon except filters as a physical 'decoder'.

It is very easy to picture that, for a tiny increment over 2*Fm, something very similar to the zero-output example given in the referenced post would be generated for a fairly long burst. It's equally hard to see how any procedure would regenerate a perfect image of the input burst, intermediate oversampling etc. notwithstanding, whereas for a much lower frequency (same Fs) it looks feasible.

Since available hw/sw has improved so radically since CDs and DAT first showed up, has anyone seen a genuine graphed experiment? Something like a 20khz pure tone riding on a (reasonably short) 2khz tone burst (or any such two-tone combination to give a gauge for zero phase shift), sampled at a rate just higher (say 5% e.g.) than twice the higher frequency. I would love to see the results, in time domain, input vs. output. I'd also love for them to be exact over multiple sample/decode tests. Then I could get over the feeling that not just oversampling but huge oversampling might be a good idea after all.

Or does insistence on a readable burst of short length throw such a monkey wrench into the spectrum being sampled that Fmax would in reality be way, way higher?

A short burst of frequency f contains a ton of higher frequency spectra energy which is needed to make the burst short, or rectangular windowed, rather than infinitely long.  Thus, a finite burst no longer meets the Nyquust criteria for sampling at just a bit over 2 times f.

The transform of a rectangular window is a Sinc, which has infinite support in the frequency domain.  But a really really long window has a Sinc that decays faster, and thus may go below your noise floor soon enough.

If you want to sample at just a bit over 2 times f and not end up with aliasing, the signal may need to be stationary for a very very long time to allow reconstruction.

Or you can sample at a higher sample rate to acquire a sufficient amount of the spectra of the window that limits the length of a much shorter time domain burst.

So there's a trade-off between how close you can get to Nyquist frequencies versus the length of the signal.


Cross correlation with short or non periodic signals? I am a total newbe to signal processing and i have a problem which may be  quite trivial to most of you guys.

I want to align two short sensor signals (speaking of max. 1000 samples recorded with 20Hz). The signals are sampled with the exact same sample rate and are not necessarily periodic. I use octave to make first experiments and i recognized some confusing behaviour of the xcorr() function in the 'signal' package.
If i correlate two short signals the resulting lag value is not really correct. The longer one of the signals is, the better the lag value calculated by the xcorr() function.

Here is what i tested:

pkg load signal

clear all

%generate two short signals
x = -pi:0.05:pi;

sigA = sin(x);
sigB = cos(x);

%plot the original signals
figure 1
subplot(2,1,1);
plot((0:length(sigA)-1), sigA, 'r');
title("sin(x)");
xlabel("Sample index");
grid on;
subplot(2,1,2);
plot((0:length(sigB)-1), sigB, 'b');
title("cos(x)");
xlabel("Sample index");
grid on;

%calculate the cross correlation between signal B and signal A
[c, lag] = xcorr(sigB, sigA);
[~,i] = max(abs(c));
shift = lag(i) %the shifting offset between the two signals

%align the signals
sigA = sigA(abs(shift)+1:end);

maxLength = max(length(sigA), length(sigB)); %for sample index axis range

%plot the aligned signals
figure 2
subplot(2,1,1);
plot((0:length(sigA)-1), sigA, 'r');
title("sin(x) aligned");
xlabel("Sample index");
axis([0, maxLength, -1, 1]);
grid on;
subplot(2,1,2);
plot((0:length(sigB)-1), sigB, 'b');
title("cos(x) aligned");
xlabel("Sample index");
axis([0, maxLength, -1, 1]);
grid on;

The code produces the following figures and calculates a shift of -27 samples:

As you can clearly see, the signal A (sin(x)) is not perfectly aligned!
I would have to change the alignment call to

%align the signals
sigA = sigA(abs(shift)+6:end); %where does this offset (1+5) come from?

to align them perfectly:

Now if i generate a longer signal for sigA, this offset disappears completely:

x = -100:0.1:100;
y = -pi:0.05:pi;

sigA = sin(x);
sigB = cos(y);

Then the situation looks like this:

As you can see in the image above, the two signals are perfectly aligned without shifting signal A manually as i had to with the shorter version of signal A.

Now here is an actual signal from my sensor (where signal A is a stored signal and signal B is a subset of signal A, the x axis shows the time in ms):

Same Problem appears here. I would have to shift signal A by 29 samples to the right to align them perfectly.

Why is that? Is cross correlation not the right thing to do for short or non periodic finite signals? Or is there a way to calculate the additional offset to the signal  shift?

Given two signals with the same finite duration and number of samples $N$, we can calculate the periodic  or cyclic crosscorrelation function which has $N$ values, or the aperiodic crosscorrelation function which has $2N-1$ values. For signals with $N$ and $M$ samples respectively ($N > M$), the aperiodic crosscorrelation function has $N+M-1$ values.

When the two signals are the same except that one is a cyclic shift of the other, the periodic crosscorrelation function has a maximum at perfect alignment. However, the maximum of the aperiodic crosscorrelation function is not necessarily at perfect alignment because the part that is "sticking out" and not overlapping with the other signal is missing, which lowers the peak.

When one signal is a truncated version of the other, and the longer signal is periodic, then the aperiodic crosscorrelation function will have multiple maxima all of the same height. Look at the array that
xcorr gives you as th eoutput, and check how the max function works in MATLAB when there are multiple maxima. Does max report the leftmost
maximum if there are multiple maxima of the same height?  (I am not a
MATLABi and cannot help you on this last point).


Question on the definition of one-sided PSD in AWGN channel I am currently simulating the following in a telecommunications exercise in Python:

  Consider a binary sequence $m[n]$ that is being modulated with QPSK
  (Gray-mapped) that is being transmitted via an AWGN (additive white
  gaussian noise) channel that adds complex noise of the form $$n = x +
 jy$$ such that $x,y$ are independent with $x  \sim \mathcal N(0,
> \sigma_x^2), y \sim \mathcal N (0, \sigma_y^2)$. The variances should
  be selected in a way that the one-sided PSD of the white noise is $N_0
 / 2$.

I am a little bit confused by the definition of one-sided PSD:

Since the one-sided PSD is $N_0 / 2$ the two sided PSD must be

$$S_{n} (f) = \frac {N_0} 4$$

The inverse-DTFT transform gives us the discrete autocorellation function $$R_n [m]
 = \frac {N_0} 4 \delta [m]$$ which means that $\mathbb E[n^2] = N_0 / 4$

$$\mathbb E[n^2] = \mathbb E[x^2 + 2jxy + y^2] = \sigma_x^2 + \sigma_y^2 =  \frac {N_0} 4$$

So we can select $$\sigma_x^2 = \frac {N_0} 4, \sigma_y^2 = \frac {N_0} 4$$
for simulating $n = x + jy$

My question is: Have I understood the meaning of the one-sided PSD correctly (and therefore are my assumptions correct)?

The relation between one-sided and two-sided noise PSD is correct. But note that the power of a complex random variable $n$ is given by $E\{|n|^2\}$, and not by $E\{n^2\}$. So in your case you obtain $E\{|n|^2\}=\sigma^2_x+\sigma^2_y$.


Is there a name for the procedure of taking the FT over separate consecutive small time-blocks? Suppose we have a continuous time-interval $I=[a,b]$, and a signal $x \colon I \to \mathbb{R}$.

A procedure that is sometimes carried out (e.g. when doing bispectral analysis) is to partition $I$ into $k$ small sub-intervals $[a\!=\!a_0,a_1],[a_1,a_2],\ldots,[a_{k-2},a_{k-1}],[a_{k-1},a_k\!=\!b]$ and take the Fourier transform of $x$ over each of these intervals,
$$ \hat{x}(n,\xi) \ := \ \int_{a_{n-1}}^{a_n} x(t)e^{-2\pi i \xi t} \, dt \, , \hspace{8mm} n \in \{1,\ldots,k\}. $$

Is there a name for $\hat{x}(n,\xi)$ (e.g. something like "time-segmented Fourier transform")?

The Short Term Fourier Transform (STFT)

https://en.m.wikipedia.org/wiki/Short-time_Fourier_transform

.

For your question, the windows are finite and non overlapping.

The continuous time case is probably more of a mathematical device that could be approximated by diffraction over some periodic, channelized structure.

The discrete time case is very commonly used, although as RJB pointed out, usually with overlap.

If your signal consisted of flat noise and a set of sinusoids, for large DFT segments, the bins would be approximately statistically independent in time and frequency for a non overlapping window, so there are applications where overlapping is not used.


Why do waveforms that are symmetrical above and below their horizontal centerlines contain no even-numbered harmonics? All About Circuits site states that waveforms that are symmetrical above and below their horizontal centerlines contain no even-numbered harmonics. Can somebody explain this mathematically, or point to a resource? I do not know much about dsp, but understand what a Fourier transform is.

The complex Fourier coefficients of a $T$-periodic function $f(t)$ are given by

$$c_n=\frac{1}{T}\int_{0}^Tf(t)e^{-j2\pi nt/T}dt\tag{1}$$

with

$$f(t)=\sum_{n=-\infty}^{\infty}c_ne^{j2\pi nt/T}\tag{2}$$

The coefficients with even indices are

$$\begin{align}c_{2n}&=\frac{1}{T}\int_{0}^Tf(t)e^{-j4\pi nt/T}dt\\&=\frac{1}{T}\left[\int_{0}^{T/2}f(t)e^{-j4\pi nt/T}dt+\int_{T/2}^{T}f(t)e^{-j4\pi nt/T}dt\right]\\&=\frac{1}{T}\left[\int_{0}^{T/2}f(t)e^{-j4\pi nt/T}dt+\int_{0}^{T/2}f(t+T/2)e^{-j4\pi nt/T}e^{-j2\pi n}dt\right]\\&=\frac{1}{T}\int_{0}^{T/2}\left[f(t)+f(t+T/2)\right]e^{-j4\pi nt/T}dt\\&=\frac12\frac{2}{T}\int_{0}^{T/2}\left[f(t)+f(t+T/2)\right]e^{-j2\pi nt/(T/2)}dt\\&=\frac12 d_n\tag{3}\end{align}$$

where $d_n$ are the complex Fourier coefficients of the $T/2$-periodic function $g(t)=f(t)+f(t+T/2)$. These coefficients can only be zero if $g(t)=0$, i.e., if $f(t)=-f(t+T/2)$. This latter condition is exactly the symmetry condition you mention in your question. Consequently, the even Fourier coefficients are zero if (and only if) $f(t)$ satisfies

$$f(t)=-f(t+T/2)\tag{4}$$

Or, in other words, a $T$-periodic function $f(t)$ has only odd harmonics if it is $T/2$-antiperiodic.
I.e., if you shift the function by half its period and flip it across the horizontal axis, it must look the same as before.


Demodulation of BFSK signal with frequency hopping using Goertzel Algorithm I am implementing a BFSK algorithm with frequency hopping on an embedded system and I am using the Goertezel algorithm in my demodulation approach. Basically, I use the goertzel alorithm to determine the powers of the signal at the next expected frequencies for both both bits 0 and 1  and whichever gives a higher power (generally in order of 10000) is the bit received. My implementation of the algorithm in C is as follows :

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, int* modData)
{
int k,i;
float   floatnumSamples;
float   omega,sine,cosine,coeff,q0,q1,q2,result;

floatnumSamples = (float) numSamples;
k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
omega = (2.0 * 3.1416 * k) / floatnumSamples;
sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;

q0=0;
q1=0;
q2=0;
result=0;

for(i=0; i<numSamples; i++)
{
    q0 = coeff * q1 - q2 + modData[i];
    q2 = q1;
    q1 = q0;
}

result = sqrtf(q1 * q1 + q2 * q2 - q1 * q2 * coeff);
return result;
}

The data at the receiver end is always in the form of a 32-bit integer.

The problem is that my received signal is a baseband signal and the band width is spread from -2KHz to 2KHz. Hence half of the hopping frequencies are in the negative side of the spectrum and goertzel algorithm is not working correctly for the negative frequencies. So I tried to use the frequency bins in the next set i.e. at (sampling frequency + hopping frequency), just as in fft. However, my algorithm gives the same output power for both the frequencies, current hopping frequency and (sampling frequency + hopping frequency).
for example, my sampling frequency is 18000 and the current hopping is -2000. So the power should be higher at -2000Hz and 16000Hz than at 2000Hz, but in my case the power is the same at 2000 and 16000Hz, which defeats the demodulation purpose.
This is not the case when I use the goertzel algorithm defined in Matlab. The Matlab function gives the powers distinctly for both current hopping frequency and (sampling frequency + hopping frequency).

Is there something that I am missing from my C implementation which can solve this issue? Is there any other method (other than using DFT, I have already tried it) that I can use for the demodulation procedure.

--Anshu

Guys I am really sorry about this stupid question. Just shifting the center frequency of the spectrum from 0Hz to some higher value ( I used 3000Hz as the center frequency) and then calculating the power using goertzel algorithm at the shifted frequencies, solved the problem. I forgot about the 2-sided nature of the frequency spectrum.


Algorithm(s) to mix audio signals without clipping I'd like to mix two or more PCM audio channels (eg recorded samples) digitally in an acoustically-faithful manner, preferably in near-real-time (meaning little or no peek-ahead).

The physically "correct" way to do this is summing the samples. However when you add two arbitrary samples, the resulting value could be up to twice the maximum value.

For example, if your samples are 16-bit values, the result will be up to 65536*2. This results in clipping.

The naive solution here is to divide by N, where N is the number of channels being mixed. However, this results in each sample being 1/Nth as loud, which is completely unrealistic. In the real world, when two instruments play simultaneously, each instrument does not become half as loud.

From reading around, a common method of mixing is: result = A + B - AB, where A and B are the two normalized samples being mixed, and AB is a term to ensure louder sounds are increasingly "soft-clipped".

However, this introduces a distortion of the signal. Is this level of distortion acceptable in high-quality audio synthesis?

What other methods are there to solve this problem? I'm interested in efficient lesser-quality algorithms as well as less-efficient high-quality algorithms.

I'm asking my question in the context of digital music synthesis, for the purpose of mixing multiple instrument tracks together. The tracks could be synthesised audio, pre-recorded samples, or real-time microphone input.

The physically "correct" way to do this is summing the samples. However when you add two arbitrary samples, the resulting value could be up to twice the maximum value. ... The naive solution here is to divide by N, where N is the number of channels being mixed.

That's not the "naive" solution, its the only solution.  That's what every analog and digital mixer does, because it's what the air does, and it's what your brain does.

Unfortunately, this appears to be a common misconception, as demonstrated by these other incorrect non-linear "mixing" (distortion) algorithms:

Mixing digital audio (the wrong way)
A quick-and-dirty audio sample mixing technique to avoid clipping (don't do this)

The "dividing by N" is called headroom; the extra room for peaks that's allocated above the RMS level of the waveform.  The amount of headroom required for a signal is determined by the signal's crest factor.  (Misunderstanding of digital signal levels and headroom is probably partially to blame for the Loudness war and Elephunk.)

In analog hardware, the headroom is maybe 20 dB.  In a hardware DSP, fixed-point is often used, with a fixed headroom; AD's SigmaDSP, for instance, has 24 dB of headroom.  In computer software, the audio processing is usually performed in 32 bit floating point, so the headroom is enormous.

Ideally, you wouldn't need to divide by N at all, you'd just sum the signals together, because your signals wouldn't be generated at 0 dBFS in the first place.

Note that most signals are not correlated to each other, anyway, so it's uncommon for all the channels of a mixer to constructively interfere at the same moment.  Yes, mixing 10 identical, in-phase sine waves would increase the peak level by 10 times (20 dB), but mixing 10 non-coherent noise sources will only increase the peak level by 3.2 times (10 dB).  For real signals, the value will be between these extremes.

In order to get the mixed signal out of a DAC without clipping, you simply reduce the gain of the mix.  If you want to keep the RMS level of the mix high without hard clipping, you will need to apply some type of compression to limit the peaks of the waveform, but this is not part of mixing, it's a separate step.  You mix first, with plenty of headroom, and then put it through dynamic range compression later, if desired.


How to Apply the SSIM Measure on RGB Images? I want to compare the similarity of two images by looking at its structure and colors. I need a measure that takes both structural and color fidelity into account. When I checked the formulation of original paper, it is said only x and y mean and variance is measured. Does SSIM account for structure only or does it take color similarity into account too?

What usually can be done is calculate it per channel and have the average (Or weighted average).
So if you have $ {R}_{1}, {G}_{1}, {B}_{1} $ and $ {R}_{2}, {G}_{2}, {B}_{2} $ you would do:
$$ {SSIM}_{r} \left( {R}_{1}, {R}_{2} \right), {SSIM}_{g} \left( {G}_{1}, {G}_{2} \right), {SSIM}_{b} \left( {B}_{1}, {B}_{2} \right) $$
Then $ {SSIM}_{rgb} = {w}_{r} {SSIM}_{r} + {w}_{g} {SSIM}_{g} + {w}_{b} {SSIM}_{b} $.
For instance, MATLAB's ssim() does the above with $ {w}_{r} = {w}_{g} = {w}_{b} = \frac{1}{3} $.
SciKit Image's compare_ssim() also apply it channel wise using the multichannel flag:


What's the meaning of a complex zero/pole? I have been studying signal processing and control for a while now, and I use Laplace and Fourier transforms almost everyday. Also another tools such as Nyquist or Bode plots.

However, I had never thought of this until today: what is the physical meaning of a complex number when dealing with frequencies?

This may sound silly, but I was asked this question and I didn't know what to answer. Why do we talk about $j\omega$ and not just $\omega$ in, for example, Fourier transforms and Bode or Nyquist plots? What is the physical sense of the real and imaginary part of a zero or a pole in the Laplace domain?

We usually talk of $j\omega$ when we're also interested in the Laplace transform of a signal / system, but want to just talk about the frequency response.

The physical meaning of the imaginary part is that it refers to purely sinusoidal signals and are constant "amplitude". The real part refers to signals for which the "amplitude" decays or grows exponentially.


Explanation on Haarlets Would someone be able to give me some info or a link etc ... regarding haarlets aka Haar wavelet-like features. I'm reading several papers for my master dissertation and several of these papers briefly mention haarlets. I can't seem to really get what they are and how they are used. I apologize if this is a basic question. But i would really appreciate some clear, structured info on the subject.

Have you seen Robust Real-time Object Detection by Viola and Jones?  This is probably the most widely used face detection algorithm, and also the most famous example of the use of the Haar wavelet-like features.


Frequency of Signal of a point tracked in a video? I have used MATLAB to track a feature point (corner point) using Lucas-Kanade tracker and plotted its movement through each frame now I need to find its frequency content, how should I find it? (I have only taken $y$-axis movement of the corner point)

You could use Fourier transform by using FFT then find the magnitude of each frquency or simply find the power spectrum by functions like pwelch.


Is wavelet analysis useful for 1D signals? Wavelets seem to be very useful for image processing. Assuming that I'm only ever going to study signals of time, i.e. 1D signals, should I still take a course in wavelet analysis? Are they applicable to 1D signals?

They absolutely are applicable. They are used for signal compression, denoising, feature extraction, filter-bank signal processing and many more things. In fact, a course of wavelet analysis will (should) stress one-dimensional analysis, which 2D is an extension of.


Frequency Domain Equivalent to Zero-Phase Filtering (`filtfilt`) Related to this question by @Molem7b5, what would be the frequency domain equivalent (filter implementation through overlap-add of DFT's) for a zero-phase filter as implemented with MATLAB/Octave/Python scipy.signal filtfilt command?

In essence filtfilt does the the following
$$y[n] = x[n] \ast h[n] \ast h[-n]$$
It's convolution with the impulse response and then the time reversed impulse response (in any order since the convolution is commutative).
So in the frequency domain it's simply
$$Y(z) = X(z) \cdot H(z) \cdot H'(z) = X(z) \cdot |H(z)|^2 $$
which obviously is zero-phase since $|H(z)|$ is real.
For an FIR filter that's easy enough: If your filter length is $N$, the cascade length is $2*N-1$ and your FFT length needs to be at least $4N-3$. Take your impulse response, zero pad it to the FFT size, take the magnitude squared and then mutliply with $e^{-j2\pi n/N} $ to get it circulated to the proper time location. Then proceed as with any other overlap add.
Alternatively you can just calculate the aggregate impules respones as
$$ h_{agg}[n] = h[n] \ast h[-n] $$
IIR filters can be done using the same approach but it requires turning the IIR into an FIR by suitable truncation of the impulse response. filtfilt() does that implicitly as well but uses the length of the signal as the truncation length.


Relation between histogram equalization and 'Auto Levels' I need to enhance an image using only global linear adjustment of intensity, i.e.:

$$I'=aI+b,\qquad a,b\in \mathbb{R}$$

where $I$ and $I'$ are input and output (enhanced) image.

I already know how to determine best values of $a,b$ through histogram stretching given the minimum and maximum intensity within image. However, this will not work on the following image. I have added black and white pixel (in red circle):

The image already contains values 0 and MAX, further histogram stretching is not possible.

However, I tried "Auto Levels" feature in Paint.NET software  with the following result:

Looking on the dialog window, there is some extra factor that affects overall histogram shape (number in the middle right part):

Is this an additional scaling factor? Or it is some parameter of non-linear transform?

How is is determined?

From the paint.NET manual, the parameter in the middle is a gamma correction, which can be used to enhance the contrast in the dark tones or high tones. The gamma-curve is the simplest non-linear level transfer curve. A more sophisticated non-linear technique for automatic level adjustment is histogram equalization, which consists in applying a monotonic, non-linear map to the image, so that the CDF of the resulting image is linear.


Why DFT-SPREAD OFDM has lower PAPR than OFDM? I just read 4G LTE/LTE-Advanced
for Mobile Broadband, and found the content as below,

  If the DFT size M equals the IDFT size N, the cascaded DFT/IDFT
  processing would obviously completely cancel each other out. However,
  if M is smaller than N and the remaining inputs to the IDFT are set to
  zero, the output of the IDFT will be a signal with “single-carrier”
  properties – that is, a signal with low power variations, and with a
  bandwidth that depends on M.

My question is that why DFT-SPREAD OFDM has lower variations in comparison with OFDM? The picture is DFT-SPREAD OFDM modulation scheme.

Thank you.

It is because the time domain symbols have lower PAPR than the frequency domain symbols.

The PAPR of OFDM is due to the 2 reasons. One is due to the IDFT(IFFT) and the other is due to the symbol type. The IDFT converts the time domain symbols to frequency domain symbols that have very high PAPR. The PAPR of 16QAM is also higher than that of QPSK. However, even when these symbol types are the same, the PAPR of DFT-SPREAD OFDM is still lower than that of regular OFDM. The reason is the (N point) IDFT effect is mostly cancelled out by the (M point) DFT. If the M=N (no zero), they are completely canceled out and the symbols would be just a time domain symbols. Even though N is much larger than M, the resulting symbol has more like time domain symbols, then the PAPR is lower. LTE up-link uses this type of OFDM that is called SC-FDMA. SC means Single Carrier because the symbols look like single carrier time domain symbols even though it has an OFDM structure. They use SC-FDMA for the up-link LTE for the lower PAPR seeking the longer battery life.


What is doppler tolerance Can anybody explain the meaning of Doppler tolerance for the waveforms(Linear frequency modulated, Barker code, Polyphase codes) used for Pulse compression.

The Doppler frequency shift $f_{d}$ due to a moving target of velocity $v_{r}$ is given by

$f_{d} = 2*v_{r}/\lambda$ where $\lambda$ is the wavelength of the transmitted signal.

Polyphase codes are known to be having high Doppler tolerance compared to LFM waveforms.
Can anybody explain this fact in simple terms.

In pulse radars you perform matched filtering in receiver where you compute cross-correlation of the received waveform and reference one. Due to Doppler shift matched filter performance  will decrease. It's like decreasing of correlation function with increasing of time delay between waveforms. Doppler tolerance defines the maximum available Doppler shift for given waveform such that you can still achieve correlation peak bigger than some threshold (SNR is assumed high enough to produce such a correlation).

There is fundamental property of all waveforms used in Radar (and other apps): ambiguity function. Look https://en.wikipedia.org/wiki/Ambiguity_function. In a few words this function determines how correlation between 2 identical waveforms changes with time delay and also with doppler shift. Some waveforms have very narrow ambiguity function, while others - more spreading.

For example maximum length sequences have narrow AF both in time and frequency shifts. LFM has slightly wide AF especially in frequency shift. They are actually very tolerable to Doppler shift. COZAC sequences (used in LTE) have narrow AF in time axe but wide in frequency one.

So if you are choosing waveform for your application you should refer to its AF at first.


Designing an allpass IIR filter with arbitrary phase response could you tell me step by step the algorithm to design an arbitrary phase response IIR allpass filter design ?

I know Matlab uses a "brute force" algorithm to find an arbitrary group delay response but how do they do that?

Depends a bit on what exactly you are trying to do. Generally, strict all-pass filters have reciprocal zeros and poles. This will always result in 180 phase shift per filter order so an arbitrary phase response that's strictly allpass cannot be designed. Example: let's say you want something that's in phase at DC,  sits at -90 degrees in the band and then goes slowly back to zero at Nyquist: this can't be implemented using strict all-passes.

Of course you can always approximate but the best method tends to depend on the specifics of your design problems and what types of errors you can tolerate.


How can I detect potholes using an accelerometer mounted in a vehicle? I am currently studying DSP and FFT, I am very new to this and have been doing electronics for a long time with Arduino and hobby projects. Recently, I am doing a project with the goal to map potholes along your daily commute in the car, measuring road roughness. This is not the road profile, but the subjected roughness of the ride the driver feels during the trip. I have an acceleration mounted to measure Z-axis "vertical" acceleration of the car, when a pothole is struck the shocks and springs dampen the force as per the quarter car model for simplicity.

Basically I want to be able to create a detector that will detect potholes by using FFT to find the dominating frequency of the vertical acceleration, then I can maybe match it to pattern to detect door slamming of the car door, or a man-hole cover from a pothole maybe in the future.

I'm not sure if FFT is the way to go, any input or advice and ideas are greatly appreciated, I have made a low pass filter to filter out the higher frequency vibration from the drive train and engine.

Any ideas on the subjected "roughness" felt by the passengers/driver in the car? I was thinking of doing road tests on a smooth road, gravel road, brick road, bad road, and coming up with a threshold of the magnitude of the vertical acceleration and making a scale somehow on what is considered "rough"

The potholes detected would be also logged with GPS their locations, I'm trying to use the Rasberry Pi for this, with data storage and its speed as a good project for it.

Thanks any help or ideas is appreciated, books, tutorials, wisdom, etc

UPDATE:

Carry over additional data here:

https://electronics.stackexchange.com/questions/56238/accelerometer-data-smoothing-filtering-pothole-detection

You are not really interested in the frequency information in your signal so much as detectable features - as such I would suggest that there isn't that much utility in an FFT here.

What you really want to do is correlate the continuous 'car signal' recorded in the car with a different signal which is the 'signature' of the pothole bump. You will need to acquire this signature by recording it in a car going over a pothole (perhaps recording many and using them to produce a generic signal which is a good representation of a pothole). Once you have it you can cross correlate it with your car signal. You may want to do this in the frequency domain in which case you should cross correlate the FFT of your signature with a sliding FFT of the car signal.

You may want to process the car signal before you use it to remove, for example, high frequency noise (with a lowpass filter). You can do this pretty effectively by eye - play around with various filters on car signals to see which removes the most extraneous information while leaving the important pothole information intact.


ouput of low pass filter $$y[n]=(-1)^n w[n]+w[n]$$

$n$ ranges from $-\infty$ to $\infty$

where $w[n]$ is output of  frequency response block with cutoff frequency = $\pi/2$, and input is unit impulse function.

I know output of frequency response block is $w[n]$ is $h[n]$
because $\delta[n]$ conv with  $h[n]$ =  $h[n]$  since delta exists only at zero

May be we can solve this in frequency domain?

The half-band lowpass filter, has the following impulse response :
$$ h[n] = \frac{\sin(\frac{\pi}{2} n) }{ \pi n} $$

And for $x[n] = \delta[n]$ ,  $w[n] = h[n]$.

Then the output $y[n]$ will be :

$$y[n] = w[n] \big( 1 + (-1)^n \big) $$

which can be simplified by expanding the parenthesis as:

$$
\begin{align}
y[n] &= \frac{\sin(\frac{\pi}{2} n) }{ \pi n} \cdot \big( e^{j 2 \pi n} + e^{j \pi n} \big) \\ \\
y[n] &= \frac{\sin(\frac{\pi}{2} n) }{ \pi n} \cdot  e^{j \frac{3 \pi}{2} n} \cdot  \big( e^{j \frac{\pi}{2} n} + e^{-j \frac{\pi}{2} n} \big)\\ \\
y[n] &= \frac{ \sin(\frac{\pi}{2} n) }{ \pi n} \cdot 2 \cdot  \cos( \frac{\pi}{2} n) \cdot  e^{j \frac{3 \pi}{2}  n} \\ \\
y[n] &= \frac{\sin(\pi n) }{ \pi n}  \cdot  e^{j \frac{3 \pi}{2} n}
\end{align}
$$
Which can be shown to be :

$$y[n] =  \begin{cases}{ 1 ~~~~, n = 0 \\ 0 ~~~~ , n \neq 0  } \end{cases} $$

Therefore $y[n] = \delta[n]$. The system is a zero-phase all pass filter. And finally, the required quantity is
$$ M = \sum_{n=-\infty}^{\infty} y[n] = \sum_{n=-\infty}^{\infty} \delta[n] = 1 $$


Why is this function BIBO stable? I have a continuous system described as:
$$y(t)=2x(t)+0.5x(t-2)+0.25x(t+2)$$
and I'm trying to understand why the system is stable. I know that a system is unstable when you provide an input and as a result the output is not bounded. So, when I looked at the system I set $t=\infty$ and assumed that the output goes to infinity because I'd have
$$ y(\infty)=2x(\infty)+0.5x(\infty-2)+0.25x(\infty+2)$$
making it unbounded. Looking at this equation, isn't the output going to infinity resulting in an unstable system? The solutions state it's stable and I can't see why. Thanks

Triangle inequality is your friend: if $y = a_1x_1+  a_2x_2+ a_3x_3 $, then $$|y| \le |a_1||x_1|+  |a_2||x_2|+ |a_3||x_3|$$
Hence if $x_{.}$ is bounded, the sum is bounded as well.


Linear convolution of discrete signals with defined lengths What is the maximum count of non-zero elements, that can a linear convolution of discrete signals of "lengths" 5 and 7 have?

When I label the length of signal $x[n]$ as $M$, and the length of
signal $h[n]$ as $N$, then the length of their convolution is the signal $y[n]=M+N-1$. However this is not that thing I was looking for. So how can I find the maximum count of non-zero elements?

It seems like you have already the correct answer, but try to visualize what's going on

First understand that signals of length $n_0$ are really infinite length, but have nonzero values at $n = 0$ and $n = n_0 - 1$. The values in between can be anything, but for the purposes of this problem take them to be nonzero as well.

Now perform the discrete convolution by literally shifting the length-5 signal and dot multiplying it with the length-7 signal. Your result will also be an infinite length signal with nonzero values only where the two signals overlap (when they dont overlap, you should find the convolution to be zero). In this case, there are 11 points where the signals actually overlap.

If some parts within the signal are zero, it is possible that you get fewer nonzero values in the result. However, in the max case where the full signal is nonzero you get this max, $11 = 7 + 5 - 1$ samples


Does zero-forcing equalizer need known channel impulse response? I'm studying some of the basic equalizing structres and I understand how Zero-forcing works, but it seems to me that a known channel impulse response is needed. Am I right? If so, what's the point? I mean, you're not likely going to know how the channel is, so how Zero-forcing is in any way useful?

If the received signal can be written as

$$\mathbf{y} = \mathbf{H}\,\mathbf{x} + \mathbf{n}$$

where $\mathbf{H}$ is the channel matrix, $\mathbf{x}$ is the transmitted vector, and $\mathbf{n}$ is the AWGN of the channel, then a zero forcing equalizer is simply (assuming that the channel matrix is square, and it's estimated perfectly at the receiver)

$$\mathbf{H}^{-1}\mathbf{y} = \mathbf{x}+\mathbf{H}^{-1}\mathbf{n}$$

Obviously, you need the channel impulse response, which is captured in the channel matrix. This channel matrix is estimated in practice using any channel estimation technique, but the estimation is usually not perfect, and thus the aforementioned ZF equalizer serves as a theoretical limit.


On the signal recovery in OFDM after estimation Given a modulated signal X which is transmitted based on OFDM with N = 256 sub-carriers through channel h = [0.7 0.3 0 0 0 0.4] of length L = 6; following the steps : [X --> ifft(X,N) --> adding CP --> ... ext] and then we have Y representing the received signal representing the output of FFT at receiver as below :

X = The modulated signal of length 256;

h = The channel of length 6;

Y = the received signal after taking the FFT of length 256 too.

Assuming the channel is estimated at receiver perfectly and h_est' = [0.7 0.3 0 0 0 0.4];

I's asking how can we get the sign X_est back ? Is it by X_est = Y. / fft(h_est, 256);    however when I try that way in Matlab, it doesn't give the estimated transmitted signal !

her is the code as reference:

clear; close all; clc

n_subc = 2^8;                          % # of subcarriers
n_ofdm_sym = 1;                        % # of OFDM symbols
n_data_sym = n_subc*n_ofdm_sym;         % # of data symbols to transmit
CP = n_subc/4;                          % CP length

M = 4;                                  %Modulation order setting

%Set channel
h = [0.7 0.5 0 0 0 0.2];

iter = 10^2;    %% Iteration
FFT = dftmtx(n_subc);           %% FFT matrix

SNR = [0:2:40];

for i = 1:length(SNR)
    i
    for j = 1:iter                                % Monte Carlo iterations

    d_symb = randi([0 3],1,n_data_sym);            %generate the data to transmit
    x_mod=qammod(d_symb,M) / sqrt(2);           %modulate

    x_mod_s2p = reshape(x_mod,n_subc,[]);        % serial to parallel

    x_ifft = FFT'*x_mod_s2p;             %Ifft conversion

     x_cp = [x_ifft(end-CP+1:end,:); x_ifft];     % add CP

    x_ifft_p2s = reshape(x_cp,1,[]);             % parallel to serial

    % --------------------> Channel mixing and adding noise
    y = conv(x_ifft_p2s,h,'same');                          % Linear convolution
   % y = awgn(y,SNR(i),'measured');                          % add noise
    % --------------------> Receiver

    y_s2p = reshape(y,n_subc+CP,[]);             % serial to parallel conversion

    y_cp = y_s2p(CP+1:end,:);                    % Removing CP

    y_fft = FFT*y_cp;                            % FFT

    y_p2s = reshape(y_fft,1,[]);                 % Parallel to serial

    X_est = y_p2s ./ fft(h,n_subc);
    y_demod= qamdemod(X_est,M);

    [~,ber(j)] = biterr(d_symb,y_demod);
    end
    BER(i)=mean(ber);
end

figure
semilogy(SNR,BER)

There is a mistake int the conv function which you are using in your code. In ofdm, the channel must be convoluted with every symbol.
the link provided in the above comments are ok, but you need to modify them according to your parameters.  so replace the command of y = conv(x_ifft_p2s,h,'same'); by below command:

for jj = 1 : n_ofdm_sym
         y(jj,:) = conv(h(jj,:),x_cp(jj,:));
end

then add the noise into y.

NP: be careful if the n_ofdm_sym  is more than one, you need to build your channel $h$ as a matrix of size $nofdmsym$ x $nsubc$ .

Good luck


Time-invariance, causality and stability of $h(t)$ of four given systems Question:

  The impulse response functions of four linear systems $S_1,\ S_2,\ S_3,\ S_4$ are given respectively by

  \begin{align}
h_1(t)&=1\\
h_2(t)&=u(t)\\
h_3(t)&=\frac{u(t)}{(t+1)}\\
h_4(t)&=\exp(-3t)*u(t)
\end{align}

  where $u(t)$ is the unit step function. Which of these systems is time-invariant, causal and stable?

  A) S1
  B) S2
  C) S3
  D) S4

My Approach:

we know if $h(t)$ is not equal to $0$ for $t<0$ then the system is non-causal; so option $A$ is non-causal. We also know if the impulse response function of the system is finite then the system is stable; so all the systems are stable. Now what is the condition to check Time-invariancy??

You should remember the necessary and sufficient condition for (BIBO)-stability of a system described by an impulse response $h(t)$:

$$\int_{-\infty}^{\infty}|h(t)|dt<\infty\tag{1}$$

If you use condition $(1)$ you will find out that actually only one of the four systems is stable.

The question about time-invariance does not make much sense, because if we assume that the systems are characterized by their impulse response then we already imply that the systems are linear and time-invariant. For linear time-varying systems you would need a two-dimensional impulse response, and non-linear systems cannot be characterized by an impulse response.


Scale and Rotation invariant Template Matching I'm looking for a method for scale and rotation invariant Template matching. I already tried some, but they didn't work so good for my examples or took for ever to execute
.
SIFT and SURF Feature detection failed totally. I also tried to implement a Log-Polar Template Matching function, but I never finished (didn't know exactly how to).

In these Articles (the first is in german)

http://cvpr.uni-muenster.de/teaching/ss08/seminarSS08/downloads/Wentker-Vortrag.pdf

http://www.jprr.org/index.php/jprr/article/viewFile/355/148

I read about that method. Mapping the polar coordinates worked, but I don't know if it's right. The Images look like this.

source_log_polar.png http://www.shareimages.com/images/pics/0/0/3/62394-pZSfl5WenZysnpyVnKg-source_log_polar.png

And after matching these 2 images with OpenCV's Template Matching function i got that result

Now I don't how to go on.

My templates are always simple symbols in building blueprints and the blueprints itself. The symbols can differ in size and orientation.

For example my simple blueprint:

And my template

In this example there is only one template, but in the blueprints it should find all occurrences, even the ones with sizes and/or orientations.

Does anyone have an approach how I could solve this?

Edit:

An addition to Andrey's approach. The distance capturing algorithm for a radial profile. (Using EmguCV)

private float[] getRadialProfile( Image<Gray, byte> image, Point center, int resolution )
 {

 var roi = image.ROI;

 if ( !roi.Contains( center ) )
  {
   return null;
  }

 var steps = resolution;
 var degreeSteps = 360 / (double)resolution;
 var data = image.Data;
 var peak = 0.0f;
 var bottom = double.MaxValue;
 var bottomIndex = 0;
 var width = roi.Width;
 var height = roi.Height;
 var minX = roi.X;
 var minY = roi.Y;

 float[] distances = new float[resolution];
 for ( var i = 0; i < steps; i++ )
  {
   var degree = i * degreeSteps;
   var radial = degree * Math.PI / 180.0;
   var dy = Math.Sin( radial );
   var dx = Math.Cos( radial );

   var x = (double)center.X;
   var y = (double)center.Y;

   while ( true )
    {
    x += dx;
    y += dy;
    if ( x >= minX + width || y >= minY + height || x <= minX || y <= minY )
     {
      x = -1;
      y = -1;
      break;
     }
    var pixel = data[(int)y, (int)x, 0];
    if ( pixel == 0 )
     {
      break;
     }
    }

    float distance = 0.0f;
    if ( x != -1 && y != -1 )
    {
      distance = (float)Math.Sqrt( Math.Pow( (center.X - x), 2 ) + Math.Pow( (center.Y - y), 2 ) );
    }

    distances[i] = distance;
    if ( distance > peak )
    {
      peak = distance;
    }
    if ( distance < bottom )
    {
      bottom = distance;
      bottomIndex = i;
    }
   }

    // Scale invariance. Divide by peak
   for ( var i = 0; i < distances.Length; i++ )
   {
     distances[i] /= peak;
   }

    // rotation invariance, shift to lowest value
   for ( var i = 0; i < bottomIndex; i++ )
   {
     distances.ShiftLeft(); // Just rotates the array nothing special
   }

   return distances;
}

I think that you can solve you problem in a much easier way. Considering that you are dealing with blueprints, you should not worry about edge connectivity, noise, and many other things that SIFT and SURF were built to accommodate for. Your template is a hollow shape with specific edge shapes.

Thus, My recommendation is:

Walk around the perimeter and find a profile of distances of the edges around the center of the template. This is the radial profile of the template. Divide by largest distance, to be scale invariant. Rotate the vector so that the smallest distance is the first, to be rotation invariant. (If your template has no dominant distance, you can change step 2 later)

Find blobs in the image. Compute the radial profile described at part (1), and compare the two vectors by normalized correlation. If your template has no dominant distance, correlation becomes normalized cross-correlation, and selecting maximum). Those who pass some threshold are considered matches.

Here is some Matlab code for you to start with - I wrote the part that finds distance profile for a specific blob and calculated it for the template:

function Doors
    im = imread('http://i.stack.imgur.com/Tf8EV.png');
    im = im(:,:,1);
    template = imread('http://i.stack.imgur.com/PlP4i.png');
    template = template(:,:,1);

    blobs = regionprops(template>0,'Area','Image');
    largestBlob = GetLargestBlob(blobs);
    [prof,edgeImage] = GetBlobRadialProfile(largestBlob);

    figure;
    subplot(1,2,1);plot(prof); title('Radial profile')
    subplot(1,2,2);imshow(edgeImage); title('Template');

end

function [prof,edgeImage] = GetBlobRadialProfile(blob)
    paddedImage = padarray( blob.Image,[8 8]);
    erodedImage = imerode(paddedImage,strel('disk',1));
    edgeImage = xor(erodedImage,paddedImage);

    c = regionprops(paddedImage,'Centroid');
    cx  = c.Centroid(1);
    cy  = c.Centroid(2);

    [y,x] = find(edgeImage);
    rad = (x(:)-cx).^2 + (y(:)-cy).^2;
    [~,minIndex] = min(rad);
    contour = bwtraceboundary(edgeImage, [y(minIndex), x(minIndex)],'N');
    prof = (contour(:,2)-cx).^2 + (contour(:,1)-cy).^2;
    prof = prof./max(prof);
end

function largestBlob = GetLargestBlob(blobs)
    area = [blobs.Area];
    [~,index] = max(area);
    largestBlob = blobs(index);
end


Variance of a filtered signal I'm using a very simple 1st order Butterworth Filter shown in Matlab code:

order = 1;
cutOff = 0.1;
[b, a] = butter(order, (2*cutOff)/SampleRate, 'high');

So from the Matlab documentation of butter, it is not hard to know the transfer function is:

$$H(z) = \frac{b_1 + b_2z^{-1}}{a_1 + a_2z^{-1}}$$

The question is, if I have a signal $x$ with variance $v$, after the signal $x$ is filtered by the above high-pass filter, what is the variance of the output signal? Will it change? How can I calculate it? And is the filter linear?

Your question (as expanded in the comments) is asking if we start with
$$
x(t)  = x_{\rm true}(t) + n_1(t)
$$
and filtering it using a filter $H(\omega)$ to get
$$
x_{\rm hp}(t) = h(t) * x(t) = x_{\rm true}(t) + n_2(t)
$$
where the variance of $n_1$ is $v_1$, then what is $v_2$, the variance of $n_2$?

That seems ill-posed because the filter $H$ will change $x_{\rm true}$ as well as the noise.

This page looks at the problem when $x_{\rm true}(t) = 0,\ \forall t$.

Here it is assumed that $n_1$ is white noise.

Where they get it is that
$$
r_{x_{\rm hp}}(l) = v_1 \cdot (h * h)(l)
$$
and $v_2$ is just this evaluated at $l=0$ i.e.
$$
v_2 = r_{x_{\rm hp}}(0)
$$

So you need to find the impulse response of your filter, convolve it with itself and find the zero-lag value of this.

As I said above, this doesn't quite answer your question but it might be a start.


Detecting changes across two images I am trying to calculate the change in pixels across two images.
Visually looking at them, there seem to be too many changes but calculating the distance between pixels gives me around about the same percent change values.

I am trying to calculate the distance here by adding the squared differences for each pixel and its square root which gives me the distance of value at rgb level for each pixel. i.e how much has the pixel value changed (sqrt((img1.r - img2.r)**2+(img1.g - img2.g)**2+(img1.b - img2.b)**2).

Then I am trying to check how many such pixels have changed by a threshold of 10.

I am not able to figure out, why does it show about the same difference in any two pair of images.

Is there a problem in the program or have I been thinking incorrectly

import pylab,skimage
import imageio,math,numpy as np

img1 = imageio.imread("029.jpg")
img2 = imageio.imread("028.jpg")
img3 = imageio.imread("030.jpg")
wh = img1.shape[0]*img1.shape[1]
squared_diff = ((img1 - img2)**2)
squared_diff1 = ((img3 - img2)**2)
squared_diff2 = ((img1 - img3)**2)
print float((np.sqrt(squared_diff[:,:,0] + squared_diff[:,:,1] + squared_diff[:,:,2])>10).sum())/wh*100
print float((np.sqrt(squared_diff1[:,:,0] + squared_diff1[:,:,1] + squared_diff1[:,:,2])>10).sum())/wh*100
print float((np.sqrt(squared_diff2[:,:,0] + squared_diff2[:,:,1] + squared_diff2[:,:,2])>10).sum())/wh*100

Output :

45.2615017361
46.1771918403
46.2858072917

For higher thresholds like 20 it gives 0 . np.max gives max distance as 15.x  everytime.

Is there any other way to find a feature like number of pixel changed with a threshold?

Any help is appreciated.

The value of :

img1 = imageio.imread(".jpg")

are rgb value ? i mean are they ranged between 0-255 or between 0 and 1.
Cause if they where all between 0 and 1 then:

((img1 - img2)**2)

will have always value below 10 indeed the maximum value of this operations could be only 1 and therefore the maximum value obtained here :

np.sqrt(squared_diff[:,:,0] + squared_diff[:,:,1] + squared_diff[:,:,2])

could be np.sqrt(3) and that's always minor than 10 and ">10" operator will return an array with all False values


1/n octave smoothing Given a frequency response obtained with FFT, I would like to apply a 1/n octave smoothing. What filter should I be using and how? Maybe someone could point to a good reference (a paper or book on the subject).

Typically "smoothing" means "replace the current value with average over the neighboring ones". Most common is energy smoothing, where the smoothing results in the energy average over the smoothing interval and the phase information is lost. Complex smoothing can be done as well but it's tricky business because of phase wrapping.

Energy smoothing can be expressed as
$$Y(k)=\sqrt{\frac{1}{N}\cdot \sum_{i=0}^{N-1}X(i)\cdot X^{*}(i)\cdot W_{k}(i)}$$

where $W_{k}(i)$ is some suitable window function. In the case of, say, third-octave smoothing this could be derived as the the magnitude squared of the transfer function of a third octave band pass filter around frequency k. This also means that for a, say, 1024 point FFT you need to design 1024 different bandpass filters, so that's a fair bit of work.

Things can be simplified if the exact shape of the smoothing filter is flexible. Rectangular smoothing can be done as
$$Y(k)=\sqrt{\frac{1}{b-a+1}\cdot \sum_{i=a}^{b}X(i)\cdot X^{*}(i)}$$

where
$$a = round(k*2^{-\frac{1}{2\cdot n}}), b = round(k*2^{\frac{1}{2\cdot n}})$$

are simply the indices of the band edges for $n^{th}$ octave smoothing.

There are a few more methods that are in between in the arbitrary window and the rectangular one in terms of complexity.


Best way to morph between waveforms I am making an low-frequency oscillator (LFO) that should be able to output different kinds of waveforms.
If I, say, want to be able to create a smooth transition from a sine wave to a square wave, what would the best approach be?
I have thought about creating several hardcoded lookup-tables and "scroll" between these, but I am worried about the transition not being smooth enough, unless there are many of these tables.

If anyone knows about the best way/usual way to solve this problem, I would be happy to hear about it.

this is what we sometimes call "Vector Synthesis" and other times call "Wavetable Synthesis".  linear crossfading between equal-sized circular lookup tables usually suffices.

but make sure the waveforms are aligned.  you can circularly delay (or "rotate") one lookup table relative to the other.  choose a relative rotation so that the cross-correlation between the two waveforms is maximum.  if you somehow crossfade between a sine wave and a square wave where the fundamental of the square wave is outa phase with the sine, you'll get a glitch in the crossfade.


How to group similar patterns together that are shifted along the x axis? I have a set of time series plots that are shifted along the x-axis. I am looking for an algorithm using which I can group them together. My end goal is to average them.

While reviewing the literature on this topic, I found out that Fourier Transforms might help. Since all the patterns have a similar shape, I might be able to detect initial changes in frequency for each pattern and then group them all together. I was hoping that someone could elaborate on whether my approach is correct. I am also open to other approaches that could tackle this problem.

I'm a little confused about your plots.  You say it is a time series plot, yet there are several places where your rises and falls "bend backward" so they aren't functions at all.  I assume time is not going backward at those moments, so are these complex plane plots?

Disregarding that for a moment.  You have a characteristic shape, so I think your best approach is to extract some characteristic values and use those for your comparisons.  Prior to doing that, I would also recommend that your "smooth" your shapes to reduce the impact on any noise/jaggedness in your data.

Assuming they are a time series, x[n], with evenly in time sampled points, I would run a differentiating filter on them and locate the three points the define the top of the "M".  It also looks like you have a local minimum ahead of the initial rise and one following the terminal fall.  Locating those will give you five calibration points for each graph to align them along the x axis.  The best displacement value can be found by doing a least squares fit.  If they need to be rescaled, it get a little bit more complicated, but still doable.

Same is true for vertical rescaling.

For the smoothing and differentiating, my technique of choice is to do a forward exponential smoothing, a backward exponential smoothing, then average the two for a smoothed function, and subtract the forward from the backward to get a smoothed "derivative".  It's not exactly the derivative, but close.  Details in this article of mine:

Exponential Smoothing with a Wrinkle

Basically, any smoothing or differentiating technique you choose should work as you are looking for characteristic values and as long as you apply the same approach to each graph you will get consistent compare-able results.

If you are looking to find parametric equations that define your characteristic shapes, then want to compare parameters, that's a different topic.  If you do, I would align them as much as possible on the location of characteristic points first for better results.  Fourier series would be good for that.

That should get you started.

Ced


How does this represent a low-pass filter? I am faced with the following code snippet to get an understanding of a F0 estimation algorithm (DIO). Have a quick look below:
static void DesignLowCutFilter(int N, int fft_size, double *low_cut_filter) {
  for (int i = 1; i <= N; ++i)
    low_cut_filter[i - 1] = 0.5 - 0.5 * cos(i * 2.0 * world::kPi / (N + 1));
  for (int i = N; i < fft_size; ++i)
    low_cut_filter[i] = 0.0;
  double sum_of_amplitude = 0.0;
  for (int i = 0; i < N; ++i)
    sum_of_amplitude += low_cut_filter[i];
  for (int i = 0; i < N; ++i)
    low_cut_filter[i] = -low_cut_filter[i] / sum_of_amplitude;
  for (int i = 0; i < (N - 1) / 2; ++i)
    low_cut_filter[fft_size - (N - 1) / 2 + i] = low_cut_filter[i];
  for (int i = 0; i < N; ++i)
    low_cut_filter[i] = low_cut_filter[i + (N - 1) / 2];
  low_cut_filter[0] += 1.0;
}

I struggle to understand how this realises a low-pass filter.

The arguments: N seems to be the cut-off frequency in samples, fft_size is the number of bins in our FFT, and low cut filter seems to be a pointer to an array that we should populate with our low-pass filter coefficients
The first line seems to be a Hann function with a window length of N.
Then a cumulative amplitude is calculated to negate and normalise the filter. I guess this is to ensure all numbers below are below -1.
Then we seem to do something which looks like copying the signal to the other end of the spectrum.
Then there is another copying part, which I don't get.
Then we add one to the first bin.

I realise the main confusion for me is that I would start making a low pass filter in a completely different way and I don't see the connections with the above algorithm.
For an n-th order Butterworth low pass filter, I would make a transfer function:
$$ G^{2}(\omega) = \frac{G_0^2}{1 + (\frac{j \omega}{j \omega_c})^{2n}} $$
I would then multiply this transfer function together with the FFT of my signal, that would be the easiest. I can see that this is not a Butterworth filter, but I don't know what kind of filter the above is.
Could you help me in providing some pointers to find the missing concepts in my understanding?

This is one tortured piece of code. I strongly recommend to toss it and write whatever you need from scratch.
Arguments

N filter order or some design parameter
fft_size: size of the final vector
low_cut: the actual filter impulse Response padded to FFT size

Code

First loop creates a cosine window of length N and bias it with 0.5
2nd loop zero pads to FFT length
3rd loop calcuates the sum of the amplitudes
4th loop divides by the sum and multiplies with -1.
5th loop seems to time reverse the filter (pointless since it's symmetrical) AND copies towards the end of FFT buffer AND also writes past the end of the FFT buffer in the process.
6th loop shift N samples forward by N/2-1 overwriting what was previously in there

What happens exactly here depends on values of N and fft_size and there is some minor low cut behavior in there, but overall it seems  like and extremely convoluted, imprecise and dangerous way of doing a low cut


Absolute maximum number of possible PN sequences for given length I previously asked if it was possible to use encryption-grade randomness as a CDMA spreading sequence and learned that the codes must be orthogonal, not necessarily random.

I've tried different ways of generating codes and very few can coexist without mutual harm. Everything I have seen on the subject says that you can have up to n orthogonal possibilities for an n-bit spreading code. Is there any type of PN sequence that can produce more than n orthogonal CDMA codes?

No, that isn't possible. Think about it in the geometric sense, where each bit in the sequence is a dimension. In that context, you can recast your question as:

  In an $n$-dimensional space, are there any sets of $M$ vectors (where $M > n$) that are mutually orthogonal?

For example, in three-dimensional space, can you construct a set of four vectors that are mutually orthogonal to one another? You can't. Stated in linear algebraic terms, you can only have as many basis vectors as the dimensionality of the space.


Visual Tracking / Object Tracking: Estimate the Direction the Whale Is Moving I am trying to find out which way the whale is moving. So I have an image of the whale swimming in the water and I want to find out where is it located and which direction does it swim (with respect to X/Y coordinates of the image).

For example here the angle is approximately 30 degrees in respect to X coordinate.

So my idea was to use canny edge detection to extract edges, construct a matrix from the points of the edges, calculate SVD and use the U matrix from SVD as a rotation matrix to get the angle.

So I have done something like this in python:

Where the first image is the starting image, the second one is the edges, and the third one is the scatterplot with points from the edges and the red arrow is a direction of the angle I got from the SVD. The angles I got

Angle 146.73
Angle 135.08

are not even close to real angles (which are approximately 30 and 85 degrees).

I have tested that I get a correct angle when I SVD the gaussian matrix, which I rotated by some degree.

This is my first post, so if you need some clarifications (code that I used, explanation why I think this approach will work or anything else), please let me know.

So am I on the right track to find where the whale is going (and just made some silly mistake) or should I use another idea?

I wouldn't rely on edges as they are not connected and foam/waves also generate them.

I would rely on the contrast between the colors of skin, water and foam, to get a solid blob.

Possibly a two step process: discard the more saturated pixels (the sea), then among the achromatic keep the dark ones. After experimenting I saw that erosion can be useful to get rid of unwanted residues.

In the end you should get a big blob, and determine its dominant orientation (equivalent ellipse of inertia, Feret diameter, most distant points...)

Unsaturated pixels:

Dark pixels:

Combination:

Eroded:

Blobs and ellipse of inertia:

As there is a risk of blob fragmentation, you can consider merging the largest blobs.


Analysis frequency of the DFT? So I'm reading Richard G. Lyons' book "Understanding Digital Signal Processing" and I've just started to make my way through the DFT chapter. While most of his examples and explanations make sense, I am a bit confused about one particular term he defined. This term is the "analysis frequency" defined as

$$f_\text{a} = \frac{mf_\text{s}} N$$

where $f_\text{s}$ is the sampling frequency, $m$ is the index of the current DFT term we're calculating, and $N$ is the length of the discrete, input signal $x[n]$.

So in the book it's stated that $f_\text{a}$ tells us, for each value of $m$, which specific frequency we're currently searching $x[n]$ for. And for each value of $m$ we've got a specific complex sinusoid that we're matching up with $x[n]$ point by point to calculate a single value for $X[m]$. The frequency of that complex sinusoid is said to be $f_\text{a}$ but when the sinusoid is written out explicitly, it does not seem to have the frequency one would expect.

Here are the specifics of the example in the book.

$$x_\text{c}(t) = \sin(2\pi \cdot 1000 \cdot t) + 0.5\sin(2\pi \cdot t + \frac{3\pi}4)$$
$$x[n] = x_\text{c}(n \, T_\text{s}) = \sin(2\pi \cdot 1000 \cdot n \, T_\text{s}) + 0.5\sin(2\pi n \, T_\text{s} + \frac{3\pi}4)$$
$$f_\text{s} \triangleq \frac{1}{T_\text{s}} = 8000 \text{ Hz}$$
$$N=8$$
$$f_\text{a} = \frac{1 \cdot 8000}8 = 1 \text{ kHz}$$
$$X[1]= \sum_{n=0}^7{x[n] \cos(2\pi \cdot n/8) - jx[n] \sin(2\pi \cdot n/8)}$$

So based on my understanding of the text, the $X[1]$ term is searching $x[n]$ to find frequency content at a frequency of $f_\text{a} = 1\text{ kHz}$. This is done by summing a point for point product of $x[n]$ with a complex sinusiod that has a frequency of $f_\text{a} = 1\text{ kHz}$. However, the sinusoid I wrote above doesn't appear to have that frequency.

So I guess my question is two fold. One, how do we come up with the idea of $f_\text{a}$ and how does it relate to the DFT? Second, why doesn't the complex sinusoid have the expected frequency? Is it because we're dealing with a discrete complex sinusoid?

If I have made any mistakes in understanding the problem or the concepts please let me know. I appreciate any help you guys can provide and if this post is too incoherent or "busy" then let me know and I'll try and fix it. Thanks in advance!

I will provide two interpretations of the DFT samples $X[k]$;

The first interpretation relies on the concept of an inner product, $\langle x(t), \phi_k(t) \rangle $, between the signal $x(t)$ and a base function $\phi_k(t)$ which measures and quantifies the amount of the base function $\phi_k(t)$ in the signal $x(t)$. This operation is essentially performing an analysis of $x(t)$.

When this interpretation is applied to the DFT, those discrete functions $\phi_k[n] = e^{j\frac{2 \pi}{N} k n}$ become the base functions (as complex exponentials of a specific frequency) of an orthonormal basis and those values $X[k]$ become the inner products between $x[n]$ and the complex exponentials; computing the amount of that particular base function in the signal $x[n]$. Note that the analysis is performed in the discrete-time sampled frequency $\omega_k$ and its equivalent continuous-time frequency $\Omega_k$ can be found from the  sampling operation and frequency normalization, which provides the second interpretation below.

Consider a continuous-time bandlimited signal $x_\text{c}(t)$ which is sampled by $f_\text{s}$ samples per second, or equivalently by a sampling period of $T_\text{s} \triangleq \tfrac{1}{f_\text{s}}$ seconds per sample, without aliasing, to produce the discrete-time samples $x[n]$ of length $N$;

$$x[n] = x_\text{c}(n T_\text{s})$$

An $N$-point DFT of this sequence be defined as:

$$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k } ~~~ , ~~~\text{for} ~~~ k=0,1,...,N-1$$

These DFT values, $X[k]$, can also be considered as the samples of the DTFT $X(e^{j \omega})$ associated with $x[n]$ as;

$$X[k] = X(e^{j \omega}) \bigg|_{\omega_k = \frac{2 \pi}{N}k}$$

The DTFT itself is a frequency scaled and repeated version of the original continuous time signal's continuous Fourier transform $X_\text{c}(j\Omega)$ as

$$X(e^{j \omega}) = \frac{1}{T_\text{s}} X_\text{c}\left(j \tfrac{\omega}{T_\text{s}} \right)$$

in the base period $-\pi \le \omega < \pi$

Therefore the corresponding continuous-time frequency (in Hz.) that a given DFT sample $X[k]$ resides can be found by the relation

$$f_k = \frac{f_\text{s}}{N} k$$

for $ -\tfrac{N}{2} \le k < \tfrac{N}{2}$. Note the use of negative $k$ values.

So your equation provides this frequency $f_k$ (saying $f_\text{a}$ to indicate its analysis nature) as an equivalent continuous-time frequency that's being analysed by the complex exponential base function $\phi_k[n]=e^{j\frac{2\pi}{N}k n}$ with a discrete-time frequency $\omega_k = \frac{2 \pi}{N} k$.

Note that the value of $f_\text{a}$ is obtained by considering the sampling operation applied on the continuous time-signal $x_\text{c}(t)$ that's used to obtain the samples $x[n]$


Intro Question to Signal Processing (Low-Pass Filter) I have a noisy signal file in Matlab and I have to denoise the polluted signal using a discrete Fourier transform.

I'm asked to perform the fourier transform, then take its absolute value. Then study/examine the absolute values to then implement a low-pass filter for the actual sound (and corresponding high-pass filter for the background noise) in Matlab.

Any ideas on general approaches? This is an intro course so suggestions shouldn't be too formal. Also this is homework so only hints please.

Based on your description it sounds like the idea is to use the FFT to figure out where your signal energy is in terms of frequencies.  Once you know that you know what the passband of your filter should be, and from that you can decide on a reasonable cutoff frequency.  You want to cutoff as much of the noise energy as possible which means you want to have the cutoff frequency as close to the passband frequency as possible, but you need some space for the transition band to make the filter implementable with a reasonable number of taps.

By the way, I often like to look at the signal energy in the log scale because it helps to get a better understanding of the weaker portions of the signal.  In Matlab you do it like so:

plot(20*log10(abs(fft(signal))))