How to simplify the Fourier Series Using an Approximation? I have a signal, $f(t)$.  I know a function that can be used to generate this signal, such that I can determine its Fourier series.  I want to express this Fourier series in simpler terms so that the waveform can be rendered by a GPU in parallel, with the added ability to change the lowpass cutoff frequency of the signal dynamically.

Take for example a saw wave with a fundamental frequency of 1 Hz, represented by the function:

$$f(t) = {-2\tan^{-1}(\cot(\pi t)) \over \pi}$$

I'm using this form of the saw wave because it's cyclic and easy for a GPU to calculate.  Its Fourier series is rather difficult to solve, but I have the series on hand:

$$g(t, n) = {-2 \over \pi} \sum_{k=1}^n {\sin(k \pi t) \over k}$$

Now, I could technically use additive synthesis to generate the waveform using this series, but my project has real-time constraints and this simply won't do.  A single saw wave might have hundreds of partials up to the Nyquist frequency.  This takes a lot of time to compute, and on a GPU it's a branching nightmare.  Since I'm rendering all of the samples in parallel, there is no history buffer for existing filtering algorithms to operate on, and I can't do that either.

What I'm trying to do is approximate this series as a function of $n \in Z$, without the summation.  That way, it doesn't matter what cut-off I use, I get consistent performance across the board that's most likely a huge improvement over additive synthesis.  Right now I'm trying to fit the curve manually, but this is a laborious process.  Is there a way to find $h(t, n) \approx g(t, n)$, where $h(t, n)$ is the approximation of the series?

Sounds like bandlimited synthesis.

  with the added ability to change the lowpass cutoff frequency of the signal dynamically.

That makes it even more complicated.

The simplest way to do this right is additive synthesis, as you know.   That way you can adjust each harmonic's frequency and phase and amplitude on the fly. Yes, there are a lot of calculations if the fundamental frequency is low, but you can use approximations of sin() if that speeds things up.  The 2-MAC resonant filter/oscillator (1, 2), for instance:

or various polynomial approximations with varying amounts of distortion:

If the waveform is fixed and a submultiple of the sample rate, you could use IFFT synthesis, but for anything else I imagine this would be difficult because the harmonics cover multiple FFT bins.  Maybe with overlapping and fading between chunks?

Otherwise the best solution depends on what exact waveform you're trying to synthesize.

If you only need a "sawtooth-like" waveform, for instance, you can use the DSF directly, which produces waveforms with harmonics that fall off linearly:

but if you want an actual square wave you use the DSF to generate a "bipolar" band-limited impulse train (BLIT):

and then leaky-integrate the sum. 2

These methods can't be used for arbitrary waveforms, however, and can't be used for swept-frequency chirps or swept-frequency LPFs (as you want?), since the harmonics suddenly disappearing from one sample to the next produces a pop or click which can't be filtered out after the fact.

You could probably use this method for the low harmonics and then combine with additive synthesis of the top harmonics and then fade those in and out to prevent the clicks.  In other words, as a harmonic disappears from one sample, you seamlessly replace it with an additive-synthesis harmonic at the same frequency and phase, which you can then fade out more smoothly.

Wavetable synthesis, fading between wavetables as harmonics disappear?
Then there are other methods.  :)

https://ccrma.stanford.edu/~stilti/papers/blit.pdf


Can I create a transformation matrix from rotation/translation vectors? I'm trying to deskew an image that contains an element of known size. Here's a test image:

I can use aruco::estimatePoseBoard which returns rotation and translation vectors. Is there a way to use that information to deskew everything that's in the same plane as the marker board? (Presumably by using the rotation and translation vectors to create a transformation matrix to pass to warpPerspective.)

I know how to deskew the marker board itself. What I want to be able to do is deskew the other things (in this case, the cloud-shaped object) that are in the same plane as the marker board. I'm trying to determine whether or not that's possible and, if so, how to do it.

I can already put four individual markers around the object I want to deskew and use the detected corners along with the known distance between them as input to getPerspectiveTransform. That works very well. But for our real-world application it may be difficult for the user to place markers exactly. It would be much easier if they could place a single marker board in the frame and have the software deskew the other objects.

I was stuck on the assumption that the destination points in the call to getPerspectiveTransform had to be the corners of the output image. Once it dawned on me that the destination points could be somewhere within the output image I had my answer.

float boardX = 1240;
float boardY = 1570;
float boardWidth = 1730;
float boardHeight = 1400;

vector<Point2f> destinationCorners;
destinationCorners(Point2f(boardX+boardWidth, boardY));
destinationCorners(Point2f(boardX+boardWidth, boardY+boardHeight));
destinationCorners(Point2f(boardX, boardY+boardHeight));
destinationCorners(Point2f(boardX, boardY));

Mat h = getPerspectiveTransform(detectedCorners, destinationCorners);

Mat bigImage(image.size() * 3, image.type(), Scalar(0, 50, 50));

warpPerspective(image, bigImage, h, bigImage.size());

This fixed the perspective of the board and everything in its plane. (The waviness of the board is due to the fact that the paper wasn't lying flat in the original photo.)


PSNR of two images of different size in matlab I performed bicubic interpolation on a 256*256 image(img)

 dest = interp2(img,'bicubic')

and  i got a 511 * 511 image.I want to compute PSNR of a 512 * 512 image(original) and the 'dest' image as follows

original = double(original);
dest = double(dest);
[M N] = size(original);
error = original - dest;
MSE = sum(sum(error .* error)) / (M * N);

if(MSE > 0)
PSNR = 10*log(255*255/MSE) / log(10);
disp(['PSNR = ', num2str(PSNR)])
else
PSNR = 99;
disp(['PSNR = ', num2str(PSNR)])
end

But I'm getting error due to different matrix dimensions.How to avoid this error.Is it possible to calculate PSNR of images with different size?Please help

When interpolating, why not use the same size as the image to compare with:

% some random inputs
original = rand(435,782);
img = rand(100,200);
% position where to interpolate
[x,y] = meshgrid(linspace(1,size(img,2),size(original,2)),linspace(1,size(img,1),size(original,1)));
% interpolate
dest = interp2(img,x,y,'bicubic');
% display size of result
disp(size(dest));


What is Finite Rate of Innovation Signal? I have read about Finite Rate of Innovation signal by Martin Vetterli in here. But i do not understand several basic things

The paper said that finite rate of innovation is the number of degree of freedom of a signal. But i cannot imagine a degree of freedom in signal term. If i have a robot that move lateral in X, Y, and Z, hence i can say it has 3 degree of freedom. How can i imagine it in signal term?
A standard sampling is just sample a signal for every time interval, but in FRI sampling, there is a smoothing kernel, what is the purpose of this?
what is annihilation filter?

Thank you

If a signal can be exactly represented by $N$ real numbers per time interval, then its number of degrees of freedom for that time interval equals $N$. The most well-known example are band-limited signals, which can be represented by $2B$ samples per second, if their bandwidth is not greater than $B$. I.e., such signals have $2B$ degrees of freedom per second. This is equivalent to saying that these signals have a finite rate of innovation, and that rate of innovation equals $2B$.

A signal doesn't need to be band-limited in order to have a finite rate of innovation (i.e., a finite number of degrees of freedom). You can think of any signal that can be parameterized with a finite number of parameters per time unit. A simple example would be a signal of the form

$$x(t)=\sum_kw_kg(t-kT)\tag{1}$$

where $g(t)$ is a rectangular impulse with a value of $1$ in the interval $[0,T]$ and zero otherwise. You need exactly one sample per $T$ seconds to represent $x(t)$, even though $x(t)$ is clearly not band-limited.

A smoothing kernel is just an optional filter that is applied to the continuous-time signal before sampling. For sampling a (quasi-)band-limited signal, the anti-aliasing filter would be an instance of such a smoothing filter.

An annihilating filter for a specific signal is an LTI system that produces an output of zero when that specific signal is applied at its input. E.g., for a sinusoidal signal with frequency $\omega_0$, a filter with a notch at $\omega_0$ is an annihilating filter.


Does the instantaneous frequency contradict the uncertainty principle? The uncertainty principle states that there is a trade off between time and frequency. So, finding frequency components at specific time is impossible. However, the instantaneous frequency measure the frequency as a function of time. Which means using the instantaneous frequency, the frequency components could be found for a signal at a specific time. How can you interpret this? Why don't we use the instantaneous frequency for time-frequency analysis?

The uncertainty principle works in the presence of (an uncertain amount of) noise or other signals (including possible harmonics), corrupting the exact phase, and thus the rate of change of phase of the signal of interest.  If the phase is corrupt, or mixed with the phase of other signals, then deriving an instantaneous frequency from the 1st derivative of that phase might produce nonsense.  Time-frequency analysis might be one way to (statistically?) separate information about the signal of interest out of these potentially existing "corrupting" influences.

Whereas the phase of a perfectly analytic signal including zero additive noise is better defined.


Image Processing - What Is a Gradient Domain? What is the gradient domain? Does it relate to the Laplace Pyramids?

I found only articles assuming you know what it is and not define it, can someone please define what is it and what it is used for?

The reference used is Seamless Image Stitching in the Gradient Domain

You may want to read this lecture's slides on blending--They also describe Laplace pyramids as another blending tool:

Image Blending and Compositing (15-463: Computational Photography, Alexei Efros, CMU, Fall 2011)

The "gradient domain" is the name often given to a class of processing techniques which operate not on the value of the pixels in images but on differences between the values of neighboring pixels, ie. local spatial transitions.

It means that in the core computations of the algorithm, the image is represented by its first spatial derivatives (horizontal and vertical)--remember the first terms of the Taylor serie of a 1D function $f(x)$ around a point $x_0$, $f(x_0) + (x-x_0) * f'(x_0)$.

For example in the following image, for each pixel inside the white mask the image is not represented by an intensity value, but by a vector of a given direction and norm:

In the context of the stitching paper you mention, the stitched result in the overlapping region consists of a blend of the two original images which tries to respect (that is, be very close to) the observed spatial derivatives in the two images (and not the observed pixel values).


The two types of stability and "Why exponential" When I was learning about LTI systems,
I noticed that
LTI system is said to be BIBO stable
if and only if its square sum of impulse response is finite.
This expression is found on many textbook of signal processing.

I continued study and found this expression
LTI system is stable if and only if
its impulse response is exponentially bounded.
I could not understand this at first,
but by further research,
I elaborately acquire this interpretation.
However, I could not find explicit description about this.
So my question is "is these things true ?",

There are two types of stability (of the system), internal stability and external stability.
BIBO stability is equivalent to external stability.
If the system is controllable and observable, internal and external stability is equivalent.
If the system is linear, "exponentially bounded" is equivalent to internal stability.

and if 4. is true,
why the function must be "exponential decaying" ?
How is this restricted to be "exponential" ?

Let's consider discrete-time systems. The generalization to continuous-time systems is straightforward.

First note that BIBO stability does not require that the impulse response be square summable. It is necessary and sufficient that the impulse response be absolutely summable:

$$\sum_{n=-\infty}^{\infty}\left|h[n]\right|<\infty\tag{1}$$

In the context of stability of a causal LTI system, an exponentially bounded impulse response must satisfy

$$|h[n]|\le a\rho^{n},\quad n\ge 0\tag{2}$$

for some positive constants $a$ and $\rho$, with $0<\rho<1$.

If a causal LTI system is described by a rational transfer function, it can be shown that condition (1) for BIBO stability is equivalent to the condition that all poles of the transfer function lie inside the unit circle of the $z$-plane. This is equivalent with the system's impulse response being a weighted sum of decaying exponentials, and, consequently, condition (2) is always met by choosing $\rho$ equal to the largest pole radius, and by selecting an appropriate value for $a$.

In sum, condition (2) is a special case of general BIBO stability for LTI systems with rational transfer functions.

Also take a look at this answer, which also mentions non-linear systems.


Why does the periodic signal in time always give a discrete frequency spectrum I would like to know, why does the periodic signal in time always give a discrete frequency spectrum in FT?

I know the equations, but I simply dont understand why is it so.

Thanks!

Here's an intuitive explanation if the convolution theorem is taken for granted:

Since the time-domain signal is periodic, one can say that it can be built by "copying and pasting" the same block of signal every period: your periodic signal can be expressed as a little block of signal (spanning one period) convolved with a dirac comb.

Thus, its Fourier transform will be the Fourier transform of the little block multiplied by the Fourier transform of a Dirac comb (which is another Dirac comb). Multiplying a continuous signal by a Dirac comb yields a discrete signal.

The same reasoning is also true the other way round (discrete in time implies periodic in frequency).


Is there a way to obtain the transfer function from a bode plot on Python? (I know that it is possible on Matlab) Quite simply, I have a bode plot obtained from a source signal.
Now I wish to obtain the transfer function.

I know it is possible with Matlab: http://www.mathworks.com/help/ident/examples/frequency-domain-identification-estimating-models-using-frequency-domain-data.html?prodcode=ID&language=en

So, I'm asking if anyone knows if there is a way to do the same with Python and where can I search for it?

Thanks in advance

Yes, use the control toolbox that is available in python. You can install it using pip install control

From here you can implement transfer functions, bode plots, step responses etc similar to what you would be accustomed to doing with the Matlab control toolbox.

Here are further details:

https://github.com/python-control/python-control


Timing syncronisation using polyphase decimation filter This question is regarding application of polyphase decimation filter in symbol time correction. How does choosing a perticular "polyphase branch" in a polyphase decimation filter result in timing adjustment ?

Suppose I have a 256:1 decimation filter with a sampling frequency of $256\textrm{ kHz}$ with total of $256\cdot 8$ taps (or $256$ filter banks with $8$ tap each). These are linear phase FIR filters, each filter bank having a data rate of $1\textrm{ kHz}$. Consider following two scenario:

I start with polyphase filter bank index $1$ and carry on convolving input data and roll over polyphase bank index at $256$.
I start with polyphase filter bank index $2$ and roll over polyphase bank index at $256$.

How does this result in timing adjustment ?

The phases of a  polyphase filter can be interpreted as a fractional delays. Consider your filter with $8\cdot 256$ taps in the up-sampled domain. In order to avoid aliasing the filter needs to be a low pass that doesn't allow any content above $1/256$ of the sample rate.

Since there is no energy above $1/256\cdot f_s$ in the transfer function you can actually down-sample the impulse response by a factor of $256$ without aliasing. That means that all phases have the same magnitude response in the down-sampled domain. The only difference between the phases is a linear phase that corresponds to a fractional delay. If phase $0$ has a delay of "$0$" samples, then phase $1$ will have a delay of $1/256$ samples etc. in the down sampled domain.


Using white noise to test filter freq. response? If I want to test a digital filter, is using white noise a good way? I mean input white noise to the filter, and visualize the frequency out from the output.

Is that procedure the same as testing a sine wave freq. by freq., to plot the freq. response.

Hope the question makes sense.

Yes you can also identify a system (more spcifically an LTI system) by inputing white noise to it and measuring the corresponding output. Yet this will be a statistical approach compared to the sine-wave measurement which is deterministic.

First, you have to be able to generate a WSS (wide-sense stationary) white noise of power $\sigma_x^2$ as the input $x[n]$ for testing.

Then you have to measure the output signal $y[n]$ and estimate it's power spectral density (PSD), assuming it's also WSS (which is indeed). So you will use some statistical approaches to obtain those quantities.

Then the relation between the input PSD of the white noise of power $\sigma_x^2$ and output PSD of the (colored) noise through the LTI is:

$$ S_{yy}(e^{j\omega}) = |H(e^{j\omega})|^2 S_{xx}(e^{j\omega}) $$

where $S_{xx}(e^{j\omega})$, and $S_{yy}(e^{j\omega})$ are the discrete-time input and output PSDs respectively. Since $S_{xx}(e^{j\omega}) = \sigma_x^2$ the output PSD will be:
$$ S_{yy}(e^{j\omega}) = |H(e^{j\omega})|^2 \sigma_x^2 $$

However, what is worth noting here is that,  this method yields not the frequency response $H(e^{j\omega})$, but its magnitude square. So the phase information is lost.

A better approach, which yields the frquency response $H(e^{j\omega})$ directly is to use cross-correlation or cross-power spectral densities instead.
The relation between input noise $x[n]$, output noise $y[n]$ can be expressed as

$$ S_{xy}(\omega) = H(\omega) S_{xx}(\omega) $$

Which shows a relation between input PSD and cross-PSD, and reveals $H(\omega)$ directly when input PSD is known as in the case with white-noise input PSD of $S_{xx}(e^{j\omega}) = \sigma_x^2$:

$$ S_{xy}(\omega) = \sigma_x^2 H(\omega). $$

Here the cross power spectral density, $S_{xy}(\omega)$ between the input $x[n]$ and the output $y[n]$ can be obtained through the DTFT of the cross-correlation between $x[n]$ and $y[n]$ from

$$\phi_{xy}[m] = \mathcal{E}\{x[n]y[n+m]^*\}$$

Which is indeed equivalent to the impulse response of the LTI system to be identified :

$$\phi_{xy}[m] = \sigma_x^2 h[m]$$


Samples per symbol and number of symbols for QAM I am trying to understand the terminologies in QAM. This is the information that I have:

number of symbols -
bits per symbol -
samples per symbol -

I understand bits per symbol but I am confused with number of symbols and samples per symbol. I read that samples per symbol is sampling rate but what is the purpose of it.

Any link or explanation is greatly appreciated.

Thanks.

Number of symbols

A baseband signal model that is used is: $x(t)=\sum_{k=1}^N a_k p(t-kT)$, where $a_k$ are the symbols (there are many types, PSK, QAM, PAM for example), $p(t)$ is the pulse shaping filter, $T$ is the symbol period, and $N$ is the number of symbols. The pulse shaping filter limits the bandwidth of $x(t)$.

  Samples per symbol

Having multiple samples per symbol is called oversampling. The samples per symbol, $\text{sps}$, is usually not explicitly chosen but is determined from the symbol rate of the signal, $F$, and sample rate, $F_s$, on a particular platform: $\text{sps} = \frac{F_s}{F}$. Your statement that "samples per symbol is sampling rate" is not correct. Visually, having more samples per symbol "fills in" the symbol and will make it look smoother. There are many reasons, besides the pulse shaping filter, to oversample. For example, some algorithms in the receiver are designed to operate on the oversampled received signal, while others expect to work on the received symbols.

Code: https://github.com/B-William/DSPSE/blob/master/spsScript.m


Why does the representation of the filter matter? How explain it to people outside signal processing theory? Recently I was able to solve a problem on concerning filtering by changing the representation from zero-pole-gain to 2nd-order-section when dealing with higher sampling rate data.

But: What do representations do to the function of the filter? Why did this work?

The short answer is that each for each structure there is a different sensitivity to small perturbations of the coefficients. The locus of poles and zeros change under those perturbations.

The different forms correspond to different factorizations of the polynomials that correspond to the filters transfer function.

So with infinite precision arithmetic the different forms are exactly equivalent but the factorizations are not equivalent with respect to small perturbations.

Similar ideas are found in numerical linear algebra. A Cholesky factorization is more accurate than the unfactored matrix the factorization corresponds to.  This is why rank one updates are favored in various algorithms


Why correlation property of $\mathcal Z$-transform contains a time reversal operation I'm reading through Digital Signal Processing, Proakis and Manolakis, third edition.  I've reached section 3.2: Properties of $\mathcal Z$-transform.

One property is the convolution:

$$x(n) = x_1(n)\star x_2(n) \longleftrightarrow X(z) = X_1(z)X_2(z)$$

Other property is correlation:

$$r_{x_1,\ x_2}(l) = \sum_{n=-\infty}^{\infty}x_1(n)x_2(n-l) \longleftrightarrow R_{x_1,\ x_2}(z)=X_1(z)X_2\left(z^{-1}\right)$$

I'm confused by the time reversal (folding) which appears in the correlation property at $X_2\left(z^{-1}\right)$. I know that the folding step occurs on convolution but not on cross/auto/correlation. Why is this and what am I missing?

To clarify the "time reversal" expression I'm referring to another property of $\mathcal Z$-transform taken from the same source:

If $$x(n) \longleftrightarrow X(z)$$
Then
$$x(-n) \longleftrightarrow X(z^{-1})$$

For the two answers: Matt's one clarifies the matter using a mathematical perspective while Peter's one shows in an intuitive way how the folding occures in polynomial multiplication. So thanks, I wish SE had a way to mark two answers as accepted.

Note that (discrete-time) convolution is defined as

$$x_1[n]\star x_2[n]=\sum_kx_1[k]x_2[n-k]\tag{1}$$

and correlation is defined as

$$r_{x_1,x_2}[n]=\sum_kx_1[k]x_2[k-n]\tag{2}$$

Comparing $(1)$ and $(2)$ we see that correlation can be written as the following convolution:

$$r_{x_1,x_2}[n]=x_1[n]\star x_2[-n]\tag{3}$$

The $\mathcal{Z}$-transform of $x_2[-n]$ is given by

$$\sum_nx_2[-n]z^{-n}=\sum_nx_2[n]z^n=X_2\left(\frac{1}{z}\right)\tag{4}$$

Consequently, from $(3)$ the $\mathcal{Z}$-transform of $r_{x_1,x_2}[n]$ must be

$$R_{x_1,x_2}(z)= X_1(z)X_2\left(\frac{1}{z}\right)\tag{5}$$


Summation limits in DFT Assume a discrete time signal $(x_n)$ is given. Some texts define the DFT as
$$ X[k] = \sum_{n=-N}^N x_n\exp\left(\frac{-2\pi j k n}{N}\right) $$
while others define it as
$$X[k] = \sum_{n=0}^{N-1} x_n\exp\left(\frac{-2\pi j k n}{N}\right) $$
or
$$X[k] = \sum_{n=1}^N x_n\exp\left(\frac{-2\pi j k n}{N}\right)$$
or many other seemingly random variants with different indices. So what is the true index limits for DFT? Or is there no unique consensus and any random author is allowed to introduce a new one?

If $N$ is the DFT length, the summation is always over $N$ points. So under this assumption, the first formula doesn't make any sense. Note that the term $e^{-j2\pi nk/N}$ is $N$-periodic in $k$ and $n$, so summing over one period always gives the same value, regardless of the starting index. With the additional assumption that $x[n]$ is periodically continued, i.e., $x[n]=x[n+N]$, you can compute the given sum over any interval of length $N$. A notation emphasizing this fact is
$$X[k]=\sum_{<N>}x[n]e^{-j2\pi nk/N}\tag{1}$$


What are some of the differences between DFT and FFT that make FFT so fast? I'm trying to understand FFTs, here's what I have so far:

In order to find the magnitude of frequencies in a waveform, one must probe for them by multiplying the wave by the frequency they are searching for, in two different phases (sin and cos) and averaging each.  The phase is found by its relation to the two, and the code for that is something like this:

//simple pseudocode
var wave = [...];                //an array of floats representing amplitude of wave
var numSamples = wave.length;
var spectrum = [1,2,3,4,5,6...]  //all frequencies being tested for.

function getMagnitudesOfSpectrum() {
   var magnitudesOut = [];
   var phasesOut = [];

   for(freq in spectrum) {
       var magnitudeSin = 0;
       var magnitudeCos = 0;

       for(sample in numSamples) {
          magnitudeSin += amplitudeSinAt(sample, freq) * wave[sample];
          magnitudeCos += amplitudeCosAt(sample, freq) * wave[sample];
       }

       magnitudesOut[freq] = (magnitudeSin + magnitudeCos)/numSamples;
       phasesOut[freq] = //based off magnitudeSin and magnitudeCos
   }

   return magnitudesOut and phasesOut;
}

In order to do this for very many frequencies very quickly, FFTs use many tricks.

What are some of the tricks used to make FFTs so much faster than DFT?

P.S. I have tried looking at completed FFT algorithms on the web, but all the tricks tend to be condensed into one beautiful piece of code without much explanation.  What I need first, before I can understand the entire thing, is some introduction to each of these efficient changes as concepts.

Thank you.

The naive implementation of an $N$-point DFT is basically a multiplication by a $N \times N$ matrix. This results in a complexity of $\mathcal{O}(N^2)$.

One of the most common Fast Fourier Transform (FFT) algorithm is the radix-2 Cooley-Tukey Decimation-in-Time FFT algorithm. This is a basic divide and conquer approach.

First define the "twiddle factor" as:
$$W_N \triangleq e^{-j\frac{2\pi}{N}}$$
where $j \triangleq \sqrt{-1}$ is the imaginary unit, then the DFT $X[k]$ of $x[n]$ is given by
$$X[k] = \sum_{n=0}^{N-1} x[n] \, W_N^{kn} \, .$$
If $N$ is even (and $\tfrac{N}{2}$ is an integer), the sum can then be divided in two sums as follows
$$X[k] = \sum_{n=0}^{N/2-1} x[2n]W_N^{2kn} + \sum_{n=0}^{N/2-1} x[2n+1]W_N^{k(2n+1)}$$
where the first summation deals with the even samples of $x[n]$ and the second with the odd samples of $x[n]$. Defining $x_e[n] \triangleq x[2n]$ and $x_o[n] \triangleq x[2n+1]$ and using the fact that

$W_N^{k(2n+1)} = W_N^{2kn}W_N^k$, and
$W_N^{2kn} = W_{N/2}^{kn}$,

this can be re-written as
$$
\begin{align}
X[k] &= \sum_{n=0}^{N/2-1} x_e[n] W_{N/2}^{kn} + W_N^k\sum_{n=0}^{N/2-1} x_o[n]W_{N/2}^{kn} \\
& = X_e[k] + W_N^k X_o[k]
\end{align}
$$
where $X_e[k]$ and $X_o[k]$ are the $\tfrac{N}{2}$-point DFT transforms of the even and odd samples of $x[n]$ respectively. So we just transformed a single $N$-point DFT into two smaller $\tfrac{N}{2}$-point DFTs.  This reduces computational cost because
$$ 2 \left( \frac{N}{2} \right)^2 + N < N^2 $$
when $N > 2$.

We can then re-iterate the same process on those two smaller DFTs. This divide-and-conquer approach allows to reach complexity of $\mathcal{O}(N\log N)$, which is way better than the $\mathcal{O}(N^2)$ we had with the naive DFT implementation (as is greatly illustrated by leftaroundabout's answer).


Estimation of occupied frequency bins (location of non-zero fourier coefficients) I'm working in circuits fields and I am not very familiar with spectrum sensing techniques.

Is there a method to identify location of non-zero Fourier coefficients of a signal (just locations, not values) with minimum computational cost (less that computations of FFT)?

Suppose my signal is a vector of length $N$ associated with some white noise.Suppose the demanded frequency accuracy is
$$\frac{f_\text{sample}}{10N}$$

and imagine my signal is a 1000 samples of a signal  and we have quantized this signal through a low resolution ADC. And signal is sparse in frequency domain, meaning that and only a fraction of its DFT coefficients are non-zero (like less than 10% of coefficents). Currently I am using a method called Periodogram which is based on FFT.

Other than FFT, do you think using multiple band-pass filters work? I was thinking to Band-pass filter, but it seems to be more computationally costly.
I mean I can pass the signal through multiple band-pass (like 100) filters which altogether cover the whole spectrum and then measure the energy of the each frequency portion of the signal?

Is there a method to identify location of non-zero Fourier coefficients of a signal (just locations, not values) with minimum computational cost (less that computations of FFT)?

So, first of all, computation of FFT:

Remember, the $10N$-point FFT, which you need for the resolution you demand:

  Suppose my signal is a vector of length $N$ associated with some white noise.Suppose the demanded frequency accuracy is
  $$\frac{f_\text{sample}}{10N}$$

has complexity of $\mathcal O\left(10N\log(10N)\right)$.

You've got a sparsity of ca. $\frac1{20}$:

  comment clarify I give an example, imagin we have signal with lengh of 1000 or 2000. We take its DFT and there are less than 50 coefficient non-zero.

which means that it's not that sparse, and all algorithms exploiting sparsity by throwing linear algebra at the problem can't really benefit much from the situation; a $K$/$L$-rank matrix/matrix operation typically has best-case complexity $K\cdot L$, so considering $K=N$  this can only work out¹ in the favour of the sparse operation if $L<10\log(10N)$.  Now you gave $N\approx 1000$, and let's chose a very modest logarithm base here – $\log_{10}$ – so, you'd end up with the condition $L < 10\log_{10} 10000 = 40$, and as you can see, with $L\approx 50$, this is not generally the case – this can really only work out if your sparsity-exploiting "tone finder" is very, very, very efficient.

But, you'll say, there's superefficient parametric estimators such as ESPRIT, that can get a much higher resolution with less samples!

Yes, that's what I'd have suggested in your situation, in general, too, but you've got to realize that these subspace/Eigenvalue based approaches work on a autocorrelation signal model – i.e. they need an estimate of the auto-correlation-matrix, which you typically get by calculating the dyadic product of the $N$-sample vector – a operation with complexity $\frac 12 N^2$, i.e. much much more than the complexity of the FFT.

So, to answer your question: No, there's to my knowledge no estimator with lower computational cost, but there's more mathematically efficient estimators than the DFT; however, you might actually get a $\frac1{10}$ resolution with less than $N$ samples, but really, you need to have a rank of your matrix estimate sufficiently higher than the number of tones you're looking for, so I'm afraid this won't change that much.

¹ I'm well aware that's not really how Landau symbols work. However, FFT implementations are really highly optimized, and so is LAPack, and thus, my gut feeling says I can really use the complexity limits as "constants".


What is averaging and how can it be done? I'm studying about the analysis of (mainly) fMRI and EEG data. Multiple times it's mentioned that to reduce noise, you can use averaging, but nothing more detailed than that.

There never is literally stated that you can use this or that technique to average.
So, what exactly is averaging? What techniques can be used to achieve it?

In general, EEG, fMRI (and also MEG, SPECT and PET datasets -the so called "functional modalities") are obtained from a subject (e.g. a human being) while it is engaged in one or more "activities".

This "activity" could even be "Stay at rest with your eyes closed" or "try to solve this equation" or "simply look at the images on this screen".

Usually, more than one subjects will be engaged in a study and these subjects may be divided into groups according to some criterion to compare and contrast their responses (e.g. "Controls" (subjects that are "normal" or do not suffer from a disease) and "Diseased" (subjects that have been verified to suffer from some disease)).

Both of these modalities, in addition to sampling time they also sample space AND each fMRI / EEG dataset is itself a sample from a population of subjects mentioned above.

Having described this, you can almost see the dimensions across which you would apply averaging and why:

You can do (the common) temporal averaging (like a moving average or median filter) over each channel (for the EEG) or voxel time series (for fMRI) to reduce the effect of noise.
You can do Spatial averaging in a number of different ways and for different purposes but all that spatial averaging entails is summing together the responses from $N$ individual channels ($N$ individual vectors) and diving them by $N$ to get a mean response.

For example, perhaps you would like to obtain the average response corresponding to some brain region (occipital, parietal, frontal, other). Then you can select the $N$ vectors (channels) corresponding to those regions and obtain ONE vector by averaging them (for example, if your channels are in X(channelN, timeSampleM) you can obtain ONE average by doing mean(X) which will return one vector of size 1 x timeSampleM with the mean of each column (that corresponds to the sample taken by each channel at some time instant $M$).

Another thing you might want to do, is obtain the average response of ONE subject depending on the activity they were carrying out. For example, you would like to see the differences in the average response of a subject's brain between the conditions of REST and LISTENINGTOMUSIC. In this case, you need to know which datasets correspond to which activity and you simply obtain the average of those $X_{k,l}$ (where $k$ denotes the dataset and $l$ denotes the activity).
You can do spatial averaging over a population (and / or "activity" they were carrying out). So for example, you might want to see what is the average response of the brain at condition REST between two subject groups (Group A, Group B). In this case you select all datasets from all subjects in Group A and obtain an average dataset and similarly for Group B. You can do the same (but changing the membership) to obtain the average response of ONE group between different conditions.

(And of course, you could apply any combination of spatial / temporal / group / condition averaging depending on what you are trying to observe)

Finally, another area where averaging is inherent to the data processing is that of Evoked Potentials.

For example, imagine that you would like to examine the response of the brain to a burst of audible noise. In this case, you provide a stimulus to the subject and record its response for a number of times. After that, you average all recordings (alligned at the stimulus instance) and produce ONE signal that represents the response of the brain. A famous such signal is the P300.

As you can see, although the operation of averaging is the same (from a mathematical point of view) there are many different ways to apply it depending on the objective. Therefore, I would recommend going over a specific paper carefully to understand

a) why they apply averaging
and b) how is averaging applied.

Furthermore, i would like to recommend to you the following resources that you might find useful and which cotain far more detailed information about EEG and fMRI analysis (incudling averaging) in general:

Electroencephalograhy by Lopes Da Silva (obviously for EEG, key
chapters to focus on are 2,4,14,53,58,59)
Statistical Parametric Mapping by Friston et al (for a very
thorough treatment of the analysis of fMRI images. Earlier versions
of that book (and of the also excellent "Human Brain Function" are
available for free from here (you need to have a close look at
chapters 3-5 and especially 5 to understand the experimental side))
Scholarpedia articles on Neuroimaging (This is a very good
resource in general)


Shock filtering using structure tensor my name is niladri, I am new to image processing(actually this is my first code). I want to implement shock filter using structure tensor. I have rough idea of what structure tensor is and implemented in MATLAB. But to design a shock I need to calculate the sign using the dominating eigenvector. But as for my understanding dominating eigenvector will be number(real and complex) for any matrix and structure tensor is basically a 2X2 matrix for each point(x,y). please correct me if my understanding is wrong. My code so far is
Edit:11/08/16: I have completed the Code but still not getting the correct response. I am unable to understand my misconception

 clc;
 clear all;
 I = imread('C:\Users\Niladri\Desktop\miramarp2gs.bmp');

 [m,n] = size(I);
 struct_scale = 3;
 %inte_scale = 5;

 x = -2*struct_scale:2*struct_scale;
 g  = exp(-0.5*(x/struct_scale).^2);
 g  = g/sum(g);
 gd = -x.*g/(struct_scale^2);
 Ix = conv2(g', conv2(gd, double(I)));

 Iy = conv2(g, conv2(gd', double(I)));

 sx = Ix^2;
 sxy = Ix*Iy;
 sy = Iy^2;

 %fileID = fopen('D:\IP\test.txt','w');
 fileID1 = fopen('D:\IP\test2.txt','w');

 for i = 1:512
    for j = 1:512
    s = [sx(i,j) sxy(i,j);sxy(i,j) sy(i,j)];
        %[v ,d] = sort(eigs(s));
        [v1, d1] = sort(diag(real(eig(s))));
        %[u , s1, d] = svd(s);
        v_perp = [0 -1;1 0 ]*v1(:,1);
        t = abs(dot(v1(:,1),v_perp));

        if v1(dot(v_perp, u(:,1)))>1E-9
            disp(t);
        end

        c  = [v1(:,1),v_perp];
        temp = det(c);

        if temp > 0
            res = 1;
        elseif temp < 0
            res = 0;
        end

        fprintf(fileID1,'%1d',res);
   end
   fprintf(fileID1,'\n');
end
fclose(fileID1);

I think that the question here is "how is the dominant direction of the local slope actually estimated?" with the outlook of using it in a shock filter.

A Shock Filter is applied iteratively and each time it propagates grayscale values along the direction of an edge but not across the edge, therefore, progressively, the edge is preserved but the area towards the broad direction of the structure is "filled". This is similar to applying dilation or erosion but taking into account the local contrast gradient. The edge can be defined at various scale levels, which results to the average direction of an edge defined over a bigger area surrounding a pixel of interest $(x,y)$.

The code seems alright, in that it sets up the normalised Gaussian g and then derives its first order derivative in gd. It uses that to obtain the gradients Ix,Iy but sx, sxy should be pointwise, rather than matrix, operations. Otherwise, Ix, Iy, "filter" each other. By controlling for the variance parameter of g (and by extension, gd too), we control how big the area around the pixel of interest $(x,y)$ is (or, in terms of dilation and erosion, how big the disk -or structuring element- is).

At this point, Ix, Iy could be used to estimate the local contrast gradient but if we come across areas with roughly opposite slopes, then averaging their components would result to 0 and the (false) perception that the area is flat. Instead, the Structure Tensor is used and more specifically the eigenvector that corresponds to the largest eigenvalue (and there is one of those for each $(x,y)$ pair). For more information please see this link and specifically sections 2, 3 and 4.

So, with that vector indicating the direction of the gradient of an edge, we find the perpendicular vector (with that sign flipping) and propagate values in that direction.

Hope this helps.


Identification of sampling frequency I have the following vector:

%create vector of sampling
t=0:pi/100:2*pi;
x=sin(t);

What is the sampling frequency (of this vector) ?

Thanks

A discrete sine wave is defined as
$$\sin (2\pi  f n)$$
Now your index is n=2*pi*(0:1/200:1). This means your sampling frequency is 200Hz.


Forward Euler Discretization I don't understand why the substitution $s=\frac{z-1}{T}$ allows us to discretize a transfer function from laplace to z-transform through Forward Euler Discretization.
Can you explain to me ?

Multiplication with $s$ in the Laplace transform domain equals differentiation in the time domain. In the discrete-time domain we can approximate differentiation by the equation

$$y[n]=\frac{x[n+1]-x[n]}{T}\tag{1}$$

where $T$ is the sampling interval.

In the Z-transform domain, Eq. $(1)$ becomes

$$Y(z)=X(z)\frac{z-1}{T}\tag{2}$$

I.e., the transfer function

$$H(z)=\frac{z-1}{T}\tag{3}$$

approximates differentiation, and replacing $s$ in a continuous-time transfer function by $H(z)$ is thus a way (usually not the best one) to approximate a continuous-time system by a discrete-time system.

One of the problems with the forward Euler method is that transforming a stable continuous-time system could result in an unstable discrete-time system.


Confusion regarding Sobel operator I am taking computer vision course, where I learned about Sobel operator for detecting edges. I somewhat understand how it works as described here. But I have a basic confusion. The operator can be divided along $x$ and $y$ axis (as $G_x$ and $G_y$). But these divided operators are finding derivative along one direction, then why are they 2-dimensional?

Actually the answer is already provided in the wiki page as follows:

"Since the Sobel kernels can be decomposed as the products of an averaging and a differentiation kernel, they compute the gradient with smoothing..."

Indeed those 1-D differentiation and averaging kernels are pre-multiplied to form the 2D final kernels to compute directional derivaties.

Hence, even though as you correctly expected the row and column derivates can be computed via 1D kernels, the averaging function being utilized into, renders them to be 2D.


Remove background from image I want to remove the background from image. Let's say I have this background:

.

Then I put pen drive and tippex on it.

.

Is it possible to remove the background image from second image? If so, then how can I do it?

I'm using the AForge.NET Framework.

UPDATE

I want to remove background like in this research paper. Is there any library for this kind of work?.

Since you already know your background image beforehand, it should be simple. I done many background subtraction before. Here is an example of removal of background by comparison.

    if (newBitmap.Width == backgroundBitmap.Width && newBitmap.Height == backgroundBitmap.Height)
    {
        for(int x = 0; x < newBitmap.Width; x++)
        {
            for (int y = 0; y < newBitmap.Height; y++)
            {
                Color firstPixel = newBitmap.GetPixel(x, y);
                Color SecondPixel = backgroundBitmap.GetPixel(x, y);

                int Totalfirst = (int)((firstPixel.R * .3) + (firstPixel.G * .59) + (firstPixel.B * .11));
                int Totalsecond = (int)((SecondPixel.R * .3) + (SecondPixel.G * .59) + (SecondPixel.B * .11));
                if (Totalfirst < Totalsecond )
                {

                    Color WhiteColor = Color.FromArgb(0, 0, 0);
                    newBitmap.SetPixel(x, y, WhiteColor);
                }
                else
                {
                    Color SameColor = Color.FromArgb(firstPixel.R, firstPixel.G, firstPixel.B);
                    newBitmap.SetPixel(x, y, SameColor);
                }

            }
        }
    }
    FinalImage = newBitmap; // after substaction of back ground.
}

This code is just a copy and paste from one of the program I did back then. Note that newBitmap is the image of the object and the background, your second picture. While the backgroundBitmap is well, the background.

Use this as your reference on how your code should look like for this segment. Should be good enough.

Also, I am new to this website(my first post here), but I think for image processing questions, you can post it to https://stackoverflow.com/questions next time, there will be more people to help you.

If this method doesn't work, let me know, cause I still got some other methods. But from you question, this should be sufficient.


Python - time frequency spectrogram I have some 64 channel EEG data sampled at 256Hz and I'm trying to conduct a time frequency analysis for each channel and plot a spectrogram.

The data is stored in a numpy 3d array, where one of the dimensions has length 256, each element containing a microvolt reading over all sampled time points (total length is 1 second for each channel of data)

To be clear: my 3D array is 64*256*913 (electrode * voltages * trial). Trial is just a single trial of an experiment. So what I want to do is take a single electrode, from a single trial, and the entire 1D voltage vector and creating a time-frequency spectrogram. So I want to create a spectrogram plot from data[0,:,0] for example.

For each electrode, I want a plot where the y axis is frequency, x axis is time, and colour/intensity is power

I have tried using this in python:

from matplotlib.pyplot import specgram
#data = np.random.rand(256)
specgram(data, NFFT=256, Fs=256)

This gives me something that looks like this:

Right off the bat this looks incorrect to me because the axis ranges are incorrect

Furthermore, when I run the same code for all EEG channels, over all of my data, I end up with the exact same plot (even though I have verified that the data is different for each)

I'm pretty new to signal processing, is there somewhere that I went wrong in either how my data is laid out or how I used my function?

The idea of a spectogram is to split your signal into a number of blocks or frames, which are potentially overlapping. After windowing, an FFT is calculated per frame. The output of these FFTs are collected as column vectors in your graph. Thus, the x-axis is related to time and the y-axis is related to frequency. Since the FFT of a real-valued signal is symmetric, only half the frequencies have to be plotted. Therefore, the y-axis goes from 0 to Fs/2=128.

Since you have chosen the FFT length (NFFT=256) equal to the lenght of your signal, you have only one full frame of data.

Potential solutions are to reduce your FFT length NFFT or increase your data length. For now, the easiest way seems to be in reducing NFFT. Note however, that you also need to adjust noverlap in that case. A widely used amount of overlap is NFFT/2.

Example:
The following command will give you a spectogram.

specgram(data, NFFT=64, Fs=256, noverlap=32)

Alternative implementation:

Nfft=64
specgram(data, NFFT=Nfft, Fs=256, noverlap=Nfft/2)

The FFT length parameter NFFT leads to a tradeoff between time and frequency resolution.

EDIT To answer additional questions in the comments.

From the documentation in
[http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.specgram]:

  Returns the tuple (spectrum, freqs, t, im):

  spectrum: 2-D array

  columns are the periodograms of successive segments

  freqs: 1-D array

  The frequencies corresponding to the rows in spectrum

  t: 1-D array

  The times corresponding to midpoints of segments (i.e the columns in >spectrum)

  im: instance of class AxesImage

  The image created by imshow containing the spectrogram

Thus, the first element in the returned tuple contains the power values.


Implementing a phase-neutral band-pass filter using SciPy Does the following code properly implement a phase-neutral band-pass filter using SciPy?

freq_fraction = FREQ / (SAMPLE_RATE/2)
sharpness = 0.1
filter_length=128
fir = sig.firwin(filter_length, [freq_fraction*(1-sharpness),freq_fraction*(1+sharpness)], pass_zero=False )
filtered = sig.lfilter( fir, 1, data_to_filter )
filtered = sig.lfilter( fir, 1,filtered[::-1] )[::-1]

The code is based on the answer to:

Is it correct to subtract a low-pass filtered signal from the original signal and use the result as a "high-pass"?

which doesn't give a final code listing for how to maintain phase neutrality.  I'm assuming it's just a matter of running the filter forward and backwards but want to double-check.

Also, if there are better/different ways of doing this (other than the remez and firwin2 methods listed in the other question) I'd love to know.

How about scipy.signal.filtfilt and scipy.signal.butter filter? See the example at the bottom of the documentation.

Here's another example with and without phase neutrality: Applying filter in scipy.signal: Use lfilter or filtfilt?


what is the best way of sampling an exponential decay? Some signal is decribed by a formula

$$S(t)=S_0 e^{-tC}\text.$$

I have given $t_\text{min},t_\text{max}$ and number of points. What is the optimal way of sampling such a signal? How to choose the proper values of $t$?

If you know that this is the only signal you're observing, than literally any two sampling points from $[t_\text{min};t_\text{max}]$ will do.

This is a problem involving a bijective function with two unknowns (assuming $C\in \mathbb R$). So, two observation will, as always, do.

To minimize the error, you'd choose sampling times that are as far apart as possible, considering that your function is strictly monotonous.


understanding complex fft results i use  this for complex fft.

Output expected

$fft[3].real= 32$ (peak at 3rd bin)

$fft[61].real= 32$ ((peak at negative frequency pair of 3rd bin))

All other values negligibly small

The input is

$ y.real = \sin (2\pi*3i/64)$ where $i = 0 \to 63$

$ y.imaginary = \sin (0*3i/64)$ where $i = 0 \to 63$ (all zero)

The output i got

$fft[3].imaginary = -32$ (peak at 3rd bin)

$fft[61].imaginary = 32$ ((peak at negative frequency pair of 3rd bin))

All other values negligibly small

This is the first time i am working with complex input fft. Can somebody explain me why i am getting peaks in imaginary part and not in real part of fft. As of my understanding doing Real FFT is nothing but using one half of input as Real part input and other half as Imaginary part input to the Compelx radix 2 FFT algorithm [i remember from John G Proakis textbook]. But i don't understand why this pseudo complex signal produces this kind of output. Also can someone explain me the phase information of this pseudo complex
wave ?

Your expected output can be achieved when the input is as follows:

real part = cos(2*pi*[0:63]*3/64)

imag part = 0

So I think you have a couple of issues:

You are missing a factor 1/64 in your description of the input above.
You should use cos() instead of sin() for the real part of the input.

Try fixing issue 1 in your post, and try fixing issue 2 in your code.


If exponentials are eigenfunctions, where does the transient term come from? Here you can see that the transfer function applied to a cosine input will give you a sinusoid and a transient term:

$$
x(t) = \underbrace{(x(0) + x'(0))(2 e^{-t} - e^{-2t}) + \frac{2}{5} e^{-t}  - \frac{1}{2}e^{-t}}_{{\rm goes\ to\ } 0 {\rm\ as\ } t \rightarrow \infty }\ \ \  + \frac{1}{10} \cos(t) + \frac{3}{10} \sin(t)
$$

However, I don't understand how this can be, aren't exponentials eigenfunctions of LTI systems? so how come it can give an extra (transient) term?

At the same time it makes sense that it has a transient term from the CF of the function. How do I reconcile this?

There is a transient because the input isn't just a cosine, but a causal cosine; that is, the input isn't a cosine from $-\infty \rightarrow \infty$, only from 0 onwards.  If the input were actually just an infinite length cosine, then you would not get any transients.

I notice that the notes you are working from don't show that the input is causal.  One reason I maintain that it actually is, is that otherwise you have an output that is tending to $\infty$ as $t \rightarrow -\infty$.  And also you don't get transients without a change to the system; the change in this case is that the system starts at time $t=0$.


Very basic questions about linear systems Note: This question may seem extremely elementary, but I am not a beginner to signal processing, linear system theory, control theory, etc.
I believe my confusion is over a subtle point, so, before responding, please read the entire question so you understand the source of my confusion first.

Consider this simple autonomous system (autonomous means it has no direct dependence on $t$):

$$\dot{x}(t) = A x(t)$$

Assume the output is the state ($y(t) = x(t)$) and $A \neq 0$. My questions are:

Is this system necessarily linear?
Is this system necessarily time-invariant?
What is the zero-input response of this system?
Would the answers above change if we assumed the system was causal?

Easy, right?

Obviously yes, since it is in the form $\dot{x}(t) = Ax(t)$
Obviously yes, since $A$ does not depend on $t$
The zero-input response is: $x(t) = x(0)\,e^{At}$
No, these are unrelated properties.

But wait!

Linearity means a linear combination of the inputs must give the same linear combination of the outputs. If I scale the input, $u$, then $x$ is not affected at all. So how can it be linear?!
Time-invariance means that if I delay the input, the output is delayed by the same amount of time. But if I delay the input, $u$, then $x$ is not affected at all. So how can it be time-invariant?!
If the initial state is an "input": Doesn't linearity imply $x(0) = 0 \implies x(t) = 0$?
If the initial state isn't  an "input": Isn't the system nonlinear unless $x(0) = 0$? (see #1)
Wouldn't zero-input give zero-output for a causal LTI system (the "ZIZO" property)?

Hopefully it's clear why I'm confused.
Where am I going wrong?

The state-space model of a causal continuous-time LTI system is

$$\begin{align}\dot{x}(t)&=Ax(t)+bu(t)\\y(t)&=cx(t)+du(t)\tag{1}\end{align}$$

where $x(t)$ is the state (vector), $u(t)$ is the input, and $y(t)$ is the output.
If the state $x(t)$ is a length $N$ vector, then $A$ is an $N\times N$ matrix, and $c$ is a $1\times N$ matrix (or row-vector, if you like).

In your example, $b=0$, $c=0$, and $d=0$. So there is basically no input and no output. We can only observe the state, i.e., the "internal behavior" of the system. We can only talk about linearity and time-invariance concerning the states. If $A$ is constant, the system is linear and time-invariant in that sense. But since we can't observe a response to an input signal it is impossible to talk about linearity and time-invariance of the system in terms of its input/output behavior. And, consequently, the corresponding tests become pointless.

And concerning point 4. in your question: yes, an LTI system must have zero output for zero input. It is a non-zero initial state that can produce a non-zero output for zero-input, but with a non-zero initial state, a system is not LTI anymore. If it is LTI otherwise (apart from the non-zero initial state), then it can be modeled as LTI by representing the non-zero initial state as an additional input to the LTI system (which must now have a zero initial state, otherwise it couldn't be LTI).


Frequency specific GCC-PHAT Background
Python GCC-PHAT works for me really well in a quiet environment.
However, in a noisy environment, the results are no good. I want to locate a source of the sound of a known frequency range between x and y in a noisy environment.
The question
How can GCC-PHAT work on specific frequencies?
What I tried
In the gcc_phat I zeroed-out all but the frequencies between x and y. However, this zeroed out the tau result of gcc_phat function.

Since your sound source has a known frequency range, I would recommend a pre-filtering for noise reduction before DOA estimation.
You can use a bandpass filter with the specific frequency range, and then put the filtered signals into GCC-PHAT algorithm.


How to adjust reverse range for Hue in cvInRangeS? cvInRangeS method for HSV colorspace is used like:

cvInRangeS(imgHSV, cvScalar( hMin,sMin,vMin),cvScalar(hMax,sMax, vMax),imgThresh);

For hue value we generally use a range between two numbers, a pair of small and big one like $(10-50)$, which means $10,11,12...48,49,50$ .

Sometimes I need to use a reverse range for red colors. By "reverse range", I mean a range between a big and a small number like $(170-10)$ which loops at the value of $180$: $170,171,172 ... 179,180,0,1,2,3 ... 8,9,10$).

How do I adjust this reverse range to cvInRangeS method? Can anyone explain how to make this work?

You can generate two different threshold images, one with hue in the range from 0 to 10 and the other with hue from 170 to 180. To get an image with hue in range both from 0 to 10 and 170 to 180, just add the two threshold images together.

In C++ code this would look something like:

inRange(imgHSV, Scalar(0, sMin, vMin), Scalar(10, sMax, vMax), imgThresh);
inRange(imgHSV, Scalar(170, sMin, vMin), Scalar(180, sMax, vMax), imgThresh1);

Mat imgThreshCombined = imgThresh + imgThresh1;

You could probably do the same thing by finding the threshold image with hue in the range from 10 to 170 and then just subtracting it from a white image.


Memory for FIR and IIR I have doubt regarding memory of FIR and IIR, I have read that output of FIR depends upon present and previous inputs, while output of IIR depends upon, present and previous inputs as well as previous outputs. Can anyone explain me how FIR consumes more memory than IIR?

Source: https://dspguru.com/dsp/faqs/iir/basics/

Thanks very much.

In order to accomplish a particular specified filtering task, the FIR filter order will almost certainly be much higher than the IIR filter order.  For instance a low-pass filter that requires a particular "sharpness" in cutoff, the transition region between the passband (at low frequencies) and the stopband (at high frequencies).  What a 4th-order IIR filter can do might correspond to what a 20th or 25th-order FIR can do.  The 4th-order IIR requires 10 coefficients and 4 states.  the 20th-order FIR requires 21 coefficients and 20 states.

Feedback can be useful.


Advantages of the Rotation Translation Operation Before Doing FT Smoothing I was reading a relatively old paper from the 1970s on smoothing by FT methods (chemistry applications), where the authors show that if we do rotation translation operation on the signal (y- values) before filtering the data by cutting off higher frequencies, the output is free from various artifacts and oscillations. What that rotational-translation does is that it makes the start point and the end-point equal to zero (see equations in the original text).

I have not seen this type of operation mentioned in modern texts. Here is the explanation by the author of Digital smoothing of electroanalytical data based on the Fourier transformation, 1973.

I did not see any major advantage in MATLAB if we use an example of a noisy Gaussian peak. Was this suggestion an artifact of how FT algorithms worked in the 1970s or there is a theoretical advantage? I assume the author is trying to periodize the signal?

In the figure, please follow the arrows from the bottom to see the rotation-translation operation.  Has anyone seen this equation elsewhere and what could be the fundamental advantage of this operation?  Thanks.

Upon closer inspection, this doesn't even seem to be a rotation:

$R_n$ being the output sample corresponding to the input sample $A_n$, we can understand their equation

  $$R_n = A_n - \Delta_n \tag1$$

as the statement that we add sequence $(-\Delta_n)_{n=1,\ldots,k}$ to the input signal $a$ and get the output signal $r$. I'll call that sequence $s$, $S_n := -\Delta_n$.

Let's look at this not on the individual sample level, but from a signal level:

$$r = a + s\tag{*}\label{additive}$$

From $\eqref{additive}$ we can see that this is just adding another signal to the input signal.

Let's have a closer look at $S$, as they claim it compensates spectral components present in $R$:

\begin{align}
s_n &= \underbrace{-A_1}_{\text{const.}} - \frac{\overbrace{(A_k - A_1)}^{\text{const.}}(n-1)}{\underbrace{k-1}_{\text{const.}}}&
\text{eq. (2) from paper}\\
&= -A_1 + \frac{A_1-A_k}{k-1}(n-1)\\
&=\underbrace{\frac{A_1-A_k}{k-1}}_{c_1}n -\underbrace{\left(\frac{A_1-A_k}{k-1}+A_1
\right)}_{-c_2}\\
&=c_1n+c_2
\end{align}

That is... underwhelming. The magical rotation is just an addition with a ramp of slope $\frac{A_1-A_k}{k-1}$. So let's look at the FFT of that:

\begin{align}\DeclareMathOperator{\dft}{\mathrm{DFT}_k}
\dft\{s\}[l] &= \dft\{c_1n\}[l] + \dft\{-c_2\}[l]\\
&= \begin{cases}
-kc_2+c_1\displaystyle \frac{k(k-1)}{2},&\quad l=0\\
-c_1\displaystyle\frac{k}{1-e^{-j2\pi l/k}},&\quad l\in[1,k-1]
\end{cases} \label{matt}\tag{;}\\
\end{align}

Eq. $\eqref{matt}$ follows according to Matt.

I see absolute no sense in these operation – it  actively introduces an error term, that will skew an estimate the original exponential decay signal model, to the signal – without that actually being there.

So, that "pseudoration" actually fakes an exponential decay  component that depends less on the actual exponential decay of the signal than on the observation length $k$. That transform is bad and you shouldn't use it.

Generally, applying a rectangular window in frequency domain – which seems to be what they're aiming for – is simply a bad idea for circular convolution reasons.
We've got a pretty good answer explaining why that's the case. The fact that this isn't bad here is that there's little signal at the end of the observation anyway – but I don't see how folding that over your beginning actually improves the data! It just makes it look smoother.


Phase ambiguity in QAM modulation I'm currently studying a decision directed carrier recovery scheme. Basically a PLL try to estimate the phase error in the current received symbol, the phase detector compensate the incoming symbol's phase with the estimated one and then the result is multiplied my the complex conjugate of the decision taken.
My notes say:

  Another phenomenon that should be taken into account is the phase ambiguity of $n \pi/2$ that affects the decisions with conventional QAM constellations.

What is he talking about? I don't understand why he tells that there is phase ambiguity in the decision of a QAM symbol.
Also citing Wikipedia:

  Most QAM constellations also have $\pi/2$ phase symmetry

To me there is no phase ambiguity, in fact if I have a 4-QAM with 1 symbol in each quadrant I can just sample the I and Q branch after a matched filter and than perform a simple unambiguous atan like this:

Typically, phase ambiguity arises in the demodulation process when the carrier recovery is done via a loop that tracks the fourth harmonic of the received signal.  So why can't the loop track the carrier frequency directly? Well, the received signal doesn't have a carrier signal -- it has a modulated carrier signal whose phase is going to change at every symbol interval which makes the phase untrackable by a typical loop -- and it is necessary to take the fourth power to eliminate the modulation and produce a signal with unchanging phase. As a result, the local oscillator (LO) in the receiver might be in any of four different phases differing by $\pi/2$, and the local reference signals to the mixers in the I and Q branches of the receiver might end up being not the expected
\begin{align}\big(\cos(\omega t), -\sin(\omega t)\big)& & \text{LO and received signal are in phase}\end{align}
as reference signals (note that the Q branch reference leads the I branch reference by $\pi/2$) but instead
\begin{align}\big(\sin(\omega t),\cos(\omega t)\big)& & \text{LO lags received signal by }\pi/2\end{align}
or
\begin{align}\big(-\cos(\omega t), \sin(\omega t)\big)& & \text{LO lags received signal by }\pi\end{align}
or
\begin{align}\big( -\sin(\omega t),-\cos(\omega t)\big)& & \text{LO lags received signal by }3\pi/2.\end{align}
The result is that instead of the I and Q bits appearing at their respective demodulator outputs (the receiver produces $(I,Q)$), the receiver might be producing $(\bar{Q},I)$, or $(\bar{I},Q)$, or $(\bar{Q},I)$ where the overbars denote bit complements.  Differential demodulation is commonly used to eliminate the phase ambiguity.


Is "Windowed Fourier Transform" a synonym for "STFT"? Sometimes I find the notation WFT (Windowed Fourier Transform) while other times I see STFT (Short Time Fourier Transform).
Are they the same?

In Chapter 2.4 Previous Work of  The Short Time Fourier Transform and Local Signals, S. Okamura, 2011, one reads:

  The STFT is also known under many names such as the windowed Fourier
  Transform, the Gabor transform, and the local Fourier transform.

Later, one discovers that the definition may  slightly vary with different authors, but, at first glance, they are the same. Be careful though. For some persons, a "Windowed Fourier Transform" is a mere Fourier transform of  a windowed signal. Because there is no sliding nor shift involved. For some, "short" implies "finite support".

Leaving theses versions aside, my personal interpretation and use  is that a “Windowed Fourier Transform” is a special case of a slightly more generic  “STFT” acceptation, for two (mild) reasons:

"windowed" somehow conveys the idea that only one window is used. But you can easily use different windows for STFT, long ones for low frequency, shorter ones for high frequencies.
the concept behind  "what a window is"  is not evident. However, for most people, a window  often has a non-zero sum (or unit), and often is either increasing/decreasing (unimodal, flat being a trivial version) or monotonic (like exponential window or weighting) for causal implementations. For short-term, it only needs to be "short", or well-concentrated in time, and could be more complicated.


Inverse Fourier Transform It's well known that convolution in the spatial domain is equivalent to multiplication in the frequency domain. i have multiplied the Fourier transform of an image F with H. where H is the FT of a Gaussian function generated by the following code:

function H = gaussian_kernel(size, sigma)

    half_size = floor(size/2);
    [u,v]=meshgrid(-half_size:half_size,-half_size:half_size);
    H = exp(-(u.^2+v.^2) ./ (2 * sigma^2));
end

and the size of H is equal to the size of the image. my question is how to estimate h which is the blurring in spatial domain given that the FT of the gaussian is gaussian.

Suppose in time domain a Gaussian kernel is expressed as:
f(x,y) = exp (-r.^2/r0.^2);

where r.^2 = x.^2 + y.^2;, r0 is the radius of the exp(-1) point.
Then the 2D FFT should be:
F(u,v) = pi/r0.^2 * exp(-w.^2/w0.^2);

where w.^2 = u.^2 + v.^2; , w0 = 1/ (pi * r0);
As a result, I would change your function to:
function H = gaussian_kernel(size, sigma)

    half_size = floor(size/2);
    [u,v]=meshgrid(-half_size:half_size,-half_size:half_size);
    H = exp(-(u.^2+v.^2) ./ (2 * sigma^2));
    r0 = 2 * sigma^2 * pi^2;
    H = H * r0 * pi;
end

And,
H = gaussian_kernel(16, 2);
subplot(2,1,1),imagesc(H)  % frequency domain
subplot(2,1,2),imagesc(real(fftshift((ifft2(fftshift(H))))))  % time domain

result:

Suppose it is (-N/2+1 : N/2) /N * fs in the frequency axis (N is the sampling point number, and fs is the sampling rate), then it is supposed to be (0:N-1)/(N * fs) in spatial axis.


Formulas of the Fourier transform family It has annoyed me that there doesn't seem to be a source online where the complete complex Fourier transform family is presented with every variable defined. The lack of definitions can be a nuisance when the transforms often use same symbols to denote different variables. So, what are the relevant transforms and the equations?

(Note: I will post the answers myself, but feel free to edit the answer to be as clear and concise as possible).

As a complement to the above answers, we can also classify the transformations based on the orthogonality condition(s) used to derive each pair of transformations.

Fourier Transformation:
$$\frac{1}{2\pi}\int^{+\infty}_{-\infty}e^{ikx}dx = \delta(k)$$

Fourier Series and DTFT:
$$\frac{1}{2\pi}\int^{\pi}_{-\pi}e^{ikx}dx = \delta_{k,0}$$
$$\frac{1}{2\pi}\sum^{+\infty}_{k=-\infty}e^{ikx} = \delta(x),\;\;\;\;k\in Z, x\in (-\pi,\pi)$$
I group Fourier Series and DTFT together, because of the duality. The continuous, periodic function in time domain of Fourier Series, can be seen as the frequency domain function in DTFT. Therefore, there is not much difference between the two, at least on mathematical ground.

DFT:
$$\frac{1}{N}\sum^{N-1}_{n=0}e^{i\frac{2\pi}{N}kn} = \delta_{k,0},\;\;\;\;k=0,1,...N-1$$

Note: Take caution when dealing with the Delta functions.

Kronecker Delta: It can be treated as a discrete elementary function, and commonly used under a summation.
Dirac Delta: It cannot be derived by elementary function, although many introductory texts use $$\delta(x-a) = \lim_{n\to\infty}\sqrt{\frac{n}{\pi}}e^{-n(x-a)^2}$$ to explain it. And treat it as an ordinary function with a value of $\infty$ at $0$, and a value of $0$ otherwise. Delta function is a distribution and has it meaning definite when used under an integral.


On the Measurement Matrix Used for Compressing Sensing Assume we have a matrix $x$ of size $(8,8)$ where each column is considered to be sparse with degree of sparsity equals to $4$. it means that every column can have $4$ zeros and $4$ non-zeros values distributed randomly. The matrix $x$ can be written as follows:

Inverse Fourier transformation $(iFFT)$ is performed for every column in the matrix $x$. It means  the matrix $F$ representing the $FFT$ matrix of size $(8,8)$ is multiplied with every column in matrix $x$, So $X = F^Hx$, where $X$ is the $iFFT$ column-wise of matrix $x$. The resulted matrix $X$ is read in row-wise leading to have a new row $X'$ of size $(1,64)$.(I mean the matrix $X$-transposed is reshaped to have one column).

The vector $X'$ is convoluted with such channel $h$ resulting $y$, so $y = HX'$, where $H$ is the toeplitz matrix gotten based on the channel $h$.

My question is to recover the spare vectors in matric $x$ based on the resulted vector $y$ and matrix $H$, I'm using compressive sensing, i.e OMP algorithm.  My question is how to build the measurement matrix in that way.

$NP$: The measurement matrix can't be built straightforward since we reshaped the matrix $X$ in row-wise way. If not, I think we can build it to be $HF^H$ as used in the paper Beyond Nyquist: Efficient Sampling of Sparse Bandlimited Signals.

Thank you in advance.

It is indeed possible to formulate this setting in terms of matrix-vector products. First, let us re-formulate your $x$ (notice throughout that I use bold letters for vectors and matrices):
$$x = \begin{bmatrix}\mathbf{x}_1 & \mathbf{x}_2 & \ldots & \mathbf{x}_8\end{bmatrix}$$
where $\mathbf x_k$ is the $k$ column of $x$.
I define the vertically stacked vector $\mathbf z$:
$$\mathbf z = \begin{bmatrix}\mathbf{x}_1 \newline
\mathbf{x}_2\newline
\vdots\newline
\mathbf{x}_8\end{bmatrix}$$
Now we can perform your column-wise IDFT of $x$ as:
$$\mathbf Z = \begin{bmatrix}\mathbf F^H \newline
& \mathbf F^H \newline
&& \ddots \newline
&&& \mathbf F^H\end{bmatrix} \mathbf z = (\mathbf I_8 \otimes \mathbf F) \mathbf z$$
$\mathbf Z$ is the column-stacked equivalent of your $X$. Reading $X$ row-wise corresponds to reading $\mathbf Z$ through a permutation:
$$\mathbf Z' = \mathbf P \mathbf Z$$
where
$$\mathbf P = \begin{bmatrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \ldots\newline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & \ldots\newline
\vdots\newline
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &\ldots\newline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & \ldots\newline
\vdots\newline
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & \ldots\newline
\vdots\newline
&&&&&&&&&& \ldots & 0 & 1\end{bmatrix}_{(64 \times 64)}$$
I hope the pattern is understandable: the first 8 rows of $\mathbf P$ pick out the 1st, 9th, ... 57th entry of $\mathbf Z$, the next 8 rows of $\mathbf P$ pick out the 2nd, 10th, ... 58th entry. The last 8 rows of $\mathbf P$ pick out the 8th, 16th, ... 64th entry of $\mathbf Z$.
Now $\mathbf Z' = X'$ and we can write:
$$y = H \mathbf P (\mathbf I_8 \otimes \mathbf F) \operatorname{stack}(x)$$
Notice how I mix your notation with mine, so sorry it looks a bit messy. The "stack" operator stacks your matrix $x$ column-wise into the vector $\mathbf Z$.

It is important to notice here that this is mostly a theoretical exercise to see that the problem can indeed be brought on a standard form, familiar in compressed sensing, that can be plugged into for example OMP. You do not necessarily have to implement it this way. You may be able to still use the operations the way you describe them, and in particular it is more computationally efficient to compute the column-wise IFFT of your $x$ rather than forming the actual DFT Kronecker product matrix $\mathbf I_8 \otimes \mathbf F$ that I describe. However, it may be easier to implement by explicitly using my formulas here.
To put it another way, if you have a sparse reconstruction algorithm that understands the general compressed sensing problem $\mathbf y = \mathbf A \mathbf x$, then you can form $\mathbf A = H \mathbf P (\mathbf I_8 \otimes \mathbf F)$ and $\mathbf x = \operatorname{stack}(x)$ and plug those into the reconstruction algorithm. Whether you we will be able to "deconstruct" $\mathbf A$ and retain a more computationally efficient structure depends on your reconstruction algorithm and its interface.


What is benefit when convert signed 16 bit to unsigned 8 bit in medical image? I have an image with its type is 16 bit signed. I see some paper, the author often convert it to unsigned 8 bit. But they did not explain why they do it. Could you explain what is benefit of this work? And How to implement it by matlab code?

Many image display devices only have an 8-bit DAC.  So if one wants to analyze only what will be displayed on such a device (without added dithering or noise filtering, etc.), one can convert to that 8-bit format and potentially save storage memory.  Many standard image file formats also only support 8 bits of luminance, so the same is true for data planned to be stored or transmitted in such a format.

However, this process is lossy, and will remove information that a display supporting greater luminance resolution (10-bit DAC, etc.) can display.


Problem using pitch shifting with TD-PSOLA and formant preservation I tried out this implementation of TD-PSOLA (https://github.com/sannawag/TD-PSOLA) and even though I know how the algorithm works, I can't understand why I get NO audible differences in the output when using f_ratio = 0.5, given that f_ratio = f_new/f_orig.

I've tried printing the length of 'peaks' (that would be the total number of peaks) arrays from both, the original signal and the new signal, by using find_peaks() again after shift_pitch() returns the output. I got these results:

f_ratio                            #orig_peaks                    #new_peaks
1                                    417                            416
0.9                                  417                            389
0.8                                  417                            403
0.7                                  417                            396
0.6                                  417                            324
0.5                                  417                            404

This looks completely nonsense to me.

I've already tried using another windows (numpy's Hamming, Hanning). Of course, it has to have something to do with the overlap-and-add process...right?

Any help would be greatly appreciated.

I believe you're correct, it's an issue with the overlap-add process. The code doesn't correctly calculate the triangular windows for lower pitches.
In short, the windows are too big.
It produces incorrect results for any pitch ratio below 1, such as 0.5 (one octave down).
The diagram below explains the issue as I understand it:

The big windows introduce so many artifacts (additional peaks) that the result sounds almost exactly the same, apart from the "speaking into fan" effect.
To get around this, it would be easiest to add an if statement which would use the existing windows only when the pitch ratio is equal or greater than 1, otherwise it should use the 'ideal windows'.
However, I'm not sure to implement the 'ideal windows' in Python. The existing algorithm uses newPeaks to calculate windows. For lower pitches it must use the old peaks array. However, this requires new math which I'm not sure how to implement exactly.


Different state-space representations for Auto-Regression and Kalman filter I see that there are different ways to write an AR model into a state-space representation, so that we can apply Kalman filter to estimate the signal. See Example 1, 2 and 3 here.

I wonder what differences are between the different state-space representations on the estimation by Kalman filter?

Thanks!

I unfortunately don't know a whole lot about Kalman filters, but I think I can help you out with the state space stuff.

In Example 1, the AR model is exactly your good old DSP recursive definition of output:

$$ y_t = \alpha + \phi_1y_{t-1} + \phi_2y_{t-2} + \eta_t$$

In this case we write down the state-space model with direct correspondence with the above equation:

$$\begin{pmatrix}y_t \\ y_{t-1}\end{pmatrix} = \begin{pmatrix}\phi_1 & \phi_2 \\ 1 & 0\end{pmatrix}\begin{pmatrix}y_{t-1} \\ y_{t-2} \end{pmatrix} + \begin{pmatrix}\alpha \\ 0 \end{pmatrix} + \begin{pmatrix}1 \\ 0 \end{pmatrix}\eta_t$$

Note that in this case, the states of the system are current and previous values of the output.

In the second example, you're separating your states $c$ from your output values. This means that the states can now be anything, even though they still directly map onto output values. This way we get

$$y_t = \mu + c_t$$

$$ c_t = \phi_1c_{t-1} + \phi_2c_{t-2} + \eta_t$$

And therefore

$$\begin{pmatrix}c_t \\ c_{t-1}\end{pmatrix} = \begin{pmatrix}\phi_1 & \phi_2 \\ 1 & 0\end{pmatrix}\begin{pmatrix}c_{t-1} \\ c_{t-2} \end{pmatrix} +  \begin{pmatrix}1 \\ 0 \end{pmatrix}\eta_t$$

You should also recognize this as the standard state-space representation of a linear system, because you equations for state evolution and state-dependent output are two different equation. This separation is trivial in case of an AR model, but this latter notation is how we think of all linear state-space models in general.

The third example is a curious one. If you multiply out all the coefficients you will realize that it is actually equivalent to the first and the second examples. So why do it? I turns out that example 2 (being the proper state-space representation of the system) is called the Controllable Canonical Form of this system. If you do some reading or simply analyze the system carefully, you will realize that we can put this system into any state we like provided well-behaved values for $\phi_1$ and $\phi_2$ with the single input $\alpha$. Therefore we call such systems controllable, and it's very easy to see from this form of the state-space equations.

You should notice that two linear systems can be identical up to a change of basis. This means that we can pick a different basis to represent the same linear system. You can convince yourself that that's exactly what we've done to go from second to third example. Particularly, we like this linear transformation to transpose the state transition matrix, so that we would get for some unknown state $\boldsymbol{s}$

$$y_t = \begin{pmatrix}1 & 0\end{pmatrix} \boldsymbol{\alpha_t}$$

$$\boldsymbol{\alpha_t} = \begin{pmatrix}s_t \\ s_{t-1}\end{pmatrix} = \begin{pmatrix}\phi_1 & \phi_2 \\ 1 & 0\end{pmatrix}\begin{pmatrix}s_{t-1} \\ s_{t-2} \end{pmatrix} + \begin{pmatrix}\alpha \\ 0 \end{pmatrix} + \begin{pmatrix}1 \\ 0 \end{pmatrix}\eta_t$$

Now we can use the change of basis to find out what this state $\boldsymbol{s}$ has to be with respect to state $\boldsymbol{y}$. And we can calculate it to be

$$\begin{pmatrix}s_t \\ s_{t-1}\end{pmatrix} = \begin{pmatrix}y_t \\ \phi_2 y_{t-1}\end{pmatrix}$$

This form (transpose of Controllability Canonical Form) is called the Observability Canonical Form because if we can put a system in this form, we can easily deduce which states of the system can be observed by simply looking at the output. For some description of the canonical forms, you can read this document, and of course look around on the web. Note that in the document the states are flipped upside down, which does not change anything about the system representation, simply reordering the rows / columns of the matrices.


How to compute number of samples in 1 cycle of signal? I have 3 signal sampled at 12 KHz frequency obtained from a motor rotating at constant 2000, 2020, and 2040 rpm. For 1 second signal recording I have 12,000 time samples. Depending on the frequency I can say (1/f) is the time period of oscillation.

But I am confused how to segment the signal with one time cycle from the whole signal, so that I can choose those many samples out of 12000 samples.

Is there any mathematical way to compute length of 1 cycle of signal when its sampling frequency and rpm are given, without observing the peak to peak in oscilloscope?

Although it is hardly related to DSP but here is a quick calculation:

$R$ Rotations per minute is $\frac{R}{60}$ Rotations per second which gives $12K$  samples per second. Hence, One rotation is proportional to $\frac{1}{\frac{R}{60}}\times 12000 = \frac{720,000}{R}$ samples. For instance, for  $2000$ RPM there are $360$ samples per rotation.


Impulse response for FIR filter design I have solved the problem of finding the filter's order. The order, $N$ is 144 (N= 4/Normalized of BW).

Since the $F_s$(sampling frequency) is 7200Hz, $f_p$(pass-band edge frequency) is 500Hz and $f_s$(stop-band edge frequency) is 700Hz.

Now, I want to:

Create an ideal impulse response i.e, sinc function and then
Apply a Blackman window function to calculate my filter's coefficients.

I am designing the filter in Matlab and I am having difficulty with those two tasks.

Since the blackman window is defined in relation to its length (at least, that is, with the parameter chosen canonically), there's too few parameters open to match the band edges. You could, however:

Determine $\omega_0 \over \pi$ (the $\pi$ being the scale built in to the $\operatorname{sinc}$ function). $\omega_0$ could be picked as the arithmetic mean of the band edges.

> w0=600/7200
Window a sinc scaled accordingly by a blackman window.

 > F=sinc((-71.5:71.5)*w0)).*blackman(144)'

F now contains a sinc windowed by the blackman function, centered about the middle of the filter (giving a constant phase filter).

The resulting filter is quite steep, but it has a rather leaky stopband. If the task at hand was a different one, please excuse and explain in a little more detail.


how to find frequency response of microphone hello I want to find the frequency response of a microphone.  I give the input signal to my speaker and it produces a specific SPL with specific frequency. on the other side, I read the microphone output voltage.
but I wonder which value of SPL is required to find the frequency response of microphone?
should I fix the SPL and just change the frequency and read the microphone output? if so, which SPL value I must choose? should I fix the SPL to 94 dB SPL?
you can assume that I have a really precise reference microphone and I can find out the SPL of speaker by means of this reference microphone. I also don't really care about direction and distance because I use the reference mic so the only thing that matters to me is SPL. also assume that I afford very precise speaker too and don't care about it's response. I appreciate your help
thanks

This is unfortunately a difficult problem.

A microphone doesn't have a single frequency response. It has a different one for each direction of incidence. Especially for second order microphones (dipoles, cardioids, etc.), the frequency response also depends heavily on the distance and radiation impedance of the source. So first you need to decide: what distance and do you want on-axis, frontal averaged or diffuse field response
Measuring the microphone directly requires an EXTREMELY well behaved loudspeaker. Most microphones are very linear and quite flat, but most speakers are definitely not. A measurement grade speaker is hard to come by and very expensive.
If you really just want the microphone response you will also need an anechoic chamber that covers your frequency range of interest. Otherwise the response of your acoustic environment (room, walls, furniture, etc.) will dominates the microphone response. Buying a a few of those will seriously strain your relationship with the CFO of even a mid-sized high tech company. Don't ask me how I know :-)

Your best shot is to find a good reference microphone that has the same directionality type of the microphone you want to measure. Than excite a loudspeaker with pink noise and measure the power spectrum for both the reference mic and the your test mic at EXACT the same location and orientation. Average over many positions, rooms and a few different speakers.
The frequency response of your microphone is the difference of the mic's power spectrum to that of the reference microphone.


Can complex envelope be writen in the form of quadrature components when it has symmetric spectrum? I am reading a chapter on VCO noise in "Design of CMOS phase-locked loops from circuit level to architecture level by Behzad Razavi";I am confused when the upconverted noise is writen as $N_I\cos\omega t - N_Q\sin\omega t$, which means the baseband noise(complex envelope of the upconverted noise) have a form of $N_I+jN_Q$; I have this confusion because I think the complex envelope has symmetric spectrum which implies it is real signal so Nq should be 0; please help me understand this, thanks!

The noise is actually never upconverted, it was already there in the bandpass channel. The textbook is simply modeling the noise as a quadrature signal.
To recap:

The transmitted passband signal is usually assumed to be noise-free.
The passband signal is real, and it occupies frequencies from $-W_2$ to $-W_1$ and from $W_1$ to $W_2$.
The receiver adds noise to the signal, usually in the first few stages of the analog front-end. This noise is passband: it occupies the same spectrum as the received signal.
The passband noise goes through the quadrature receiver and its complex envelope is calculated. In the process, the spectrum from $-W_2$ to $-W_1$ is discarded, and the spectrum from $W_1$ to $W_2$ is shifted to baseband.
As a result, the spectrum of the complex envelope of the noise is not symmetric.
In conclusion, we can think of the passband noise as a quadrature signal.

Here's an illustration:


Evaluation of phase difference between two acquired signals (expected to be sine waves) with the same frequency I would like to calculate the phase difference between two signals $F(t)$ and $R(t)$ which are expected to be:
\begin{align}
F(t)&=A_F\cos\left(\omega t+\phi_F\right)\\
R(t)&=A_R\cos\left(\omega t+\phi_R\right)
\end{align} What I do at the moment is to acquire these two signals, multiply them and then by averaging I can get the phase difference as follow:

$$\phi=\arccos \left(\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R} \right )$$

Now this works fine, but I don't like the fact the I can represent the phase difference only between $0^{\circ}$ and $180^{\circ}$.

Is there an alternative method that I can use to evalutae the phase difference on the full circle $0^{\circ}$ to $360^{\circ}$?

I already tried using FFT and calculate a Cross Power Spectrum, but unfortunately this method is too slow compared to the one mentioned above. I use LabVIEW for the whole acquisition and analysis. The calculation time is crucial for me and I would like to keep it minimal.

EDIT: Here some additional information about the application. I acquire the two $F(t)$ and $R(t)$ with a sampling frequency of 10 MHz. $F(t)$ is a signal (sinusoidal) generated by a signal generator and it drives a piezoelectric element. $R(t)$ is the response of the device that I drive with the piezo. Both signals are acquired simultaneously and I'm interested in the phase difference between the two signals. After acquiring them I follow the following algorithm (in LabVIEW):

Estimation of $A_F$ and $A_R$ from their RMS values.
Multiplication  $X(t)=F(t)R(t)$
Calculation of $2\textrm{Mean}[F(t)R(t)]$ (twice the DC level of $X(t)$)
Calculation of $\cos({\phi})=\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R}$
Calculation of inverse cosine $\phi=\arccos \left(\frac{2\textrm{Mean}\left[F(t)R(t)\right]}{A_F A_R} \right )$

I then repeat the algorithm for a list of different driving frequencies to build a phase spectrum. The phase spectrum looks as expected (Simple Harmonic Oscillator) but it wraps around $180^{\circ}$. The problem is that it doesn't wraps by restarting from $0^{\circ}$ but once it reaches $180^{\circ}$ it goes backward towards $0^{\circ}$. This means that the phase spectrum is continuous but it presents something like a cusp around $180^{\circ}$. I tried some unwrapping, like flipping the portion of the phase spectrum that should be after $180^{\circ}$ but the signals have some noise and they are not perfect sinusoidal signals. This means that I cannot flip by $180^{\circ}$ because the spectrum then doesn't look continuous but with a small step. Here some screenshots of the unwrapping techniques I tried so far:

a. Calculated phase spectrum:

b. Unwrapping around $180^{\circ}$:

c. ZOOM around $180^{\circ}$:

As you can see if I unwrap around $180^{\circ}$ I get the small step. This is a big problem for me because I then fit a simple harmonic oscillator model to my phase spectrum and with that step the fit is not optimal.

EDIT 2: Some additional parameters:

a. Typical frequency ranges of $F(t)$ are $280-500\ \textrm{kHz}$ in particular the frequency range of the images above is $370-420\ \textrm{kHz}$.

b. The frequency around $180^\circ$ is $403\ \textrm{kHz}$.

c. The two signals are acquired at $10\ \textrm{MHz}$ for $1\ \textrm{ms}$ which means there are $10^{3}$ data points for each signal.

d. The amplitude of the reference signal $F(t)$ is $\pm 2\ \textrm{V}$.

EDIT 3: $F(t)$ is generated by a function generator and acquired by a digitizer which acquires simultaneously $F(t)$ and $R(t)$. The letters were chosen this way because for me $F(t)$ is the driving Force of my mechanical system and $R(t)$ is the Response of the mechanical system. The mechanical system is an oscillating cantilever, which is a beam fixed at one end. The beam has a rectangular cross section with length, width and thickness of $500\ \mu \textrm{m}$, $100\ \mu \textrm{m}$ and $1\ \mu \textrm{m}$, respectively.

I really dislike the convention of notation in the question.  "$f$" is for frequency.  Functions having capital letters are for Fourier or Laplace transforms of the time-domain signal having the lower-case letter denoting it.

Let's say our reference signal is

$$ x[n] \triangleq A_x \cos(\omega n + \theta_x) $$

Somehow we can get the Hilbert transform of that signal

$$ \hat{x}[n] \triangleq \mathscr{H} \big \{ x[n] \big \} = A_x \sin(\omega n + \theta_x) $$

If you're generating the reference signal $x[n]$, you should also be able to synchronously generate $\hat{x}[n]$.   If you're measuring $x[n]$, you have to make a hilbert transformer to create $\hat{x}[n]$ (which will be a special FIR filter) and that will introduce a known delay.  You must delay all other signals by the same amount to get this relative phase measurement right.  The general relation of the discrete Hilbert transform is:

$$ \hat{x}[n] = \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] $$

Of course, you will need to window that infinite sum to a finite sum and delay it to make it causal.  Note that all even $m$ terms in the above summation are zero.

$$\begin{align}
 \hat{x}[n] &= \sum\limits_{m=-\infty}^{\infty} \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\
 &\approx \sum\limits_{m=-2M+1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} x[n-m] \\
 &= \sum\limits_{m=1}^{2M-1} w[m] \frac{1 - (-1)^m}{m \, \pi} (x[n-m]-x[n+m]) \\
 &= \sum\limits_{m=0}^{M-1} w[2m+1]\frac{2}{(2m+1) \pi} (x[n-2m-1]-x[n+2m+1]) \\
 &= \sum\limits_{m=0}^{M-1} h[m] \, (x[n-2m-1]-x[n+2m+1]) \\
\end{align}$$

where $w[n]=w[-n]$ is a decent window (like a Kaiser) of width $4M-1$ and the hilbert coefficients are

$$ h[n] \triangleq \frac{2\,w[2n+1]}{(2n+1) \pi} $$

Since this symmetrical, acausal FIR filter looks $2M-1$ samples into the "future", it must be delayed by $2M-1$ samples to be causal and realizable in a real-time context.  So you can not realize in real-time a true Hilbert transformer, but you can realize a very close approximation with the corresponding delay of $2M-1$ samples.  You must delay all other signals by the same amount to make the following calculations valid.  This means you will really be working with $\hat{x}[n-2M+1]$ and you have to explicitly delay $x[n]$ and $y[n]$ by the same amount.

So, while in real-time reality you will be working with $x[n-2M+1]$, $\hat{x}[n-2M+1]$, and $y[n-2M+1]$, for clarity, I am not showing the signals below as so delayed.

Now the measured signal that has the same frequency $\omega$ but different phase is:

$$ y[n] = A_y \cos(\omega n + \theta_y) $$

and the number we are looking for is $\theta_y - \theta_x$ .

Consider these two products:

$$\begin{align}
y[n] \cdot x[n] &= A_y A_x \cos(\omega n + \theta_y) \cos(\omega n + \theta_x) \\
   &= \frac{A_y A_x}{2} \big( \cos(2\omega n + \theta_y+\theta_x) + \cos(\theta_y - \theta_x)  \big) \\
\\
y[n] \cdot \hat{x}[n] &= A_y A_x \cos(\omega n + \theta_y) \sin(\omega n + \theta_x) \\
   &= \frac{A_y A_x}{2} \big( \sin(2\omega n + \theta_y+\theta_x) - \sin(\theta_y - \theta_x)  \big) \\
\end{align}$$

Like you have done before, you want to compute the mean of these two products:

$$\begin{align}
 \overline{y[n] \, x[n]} & = \frac{A_y A_x}{2} \cos(\theta_y - \theta_x) \\
\\
 \overline{y[n] \, \hat{x}[n]} & = -\frac{A_y A_x}{2} \sin(\theta_y - \theta_x) \\
\end{align}$$

We get the mean by use of any decent low-pass filtering.  We don't even need to worry that the DC gain = 1 as you normally do with the mean, since the scaling will come out in the wash.
Then, using complex expression:

$$\begin{align}
 e^{j(\theta_y - \theta_x)} &= \cos(\theta_y - \theta_x) + j \, \sin(\theta_y - \theta_x) \\
  &= \frac{2}{A_y A_x} \left( \overline{y[n] \, x[n]} - j \, \overline{y[n] \, \hat{x}[n]} \right) \\
\end{align}$$

and

$$\begin{align}
\theta_y - \theta_x &= \arg \bigg\{ \overline{y[n] \, x[n]} - j \,\overline{y[n] \, \hat{x}[n]} \bigg\} \\
 &= -\operatorname{atan2}\left(\overline{y[n] \, \hat{x}[n]}, \ \overline{y[n] \, x[n]}\right) \\
\end{align}$$

To do the four-quadrant arctangent:

$$\begin{align}
\operatorname{atan2}(y,x) &= \begin{cases}
\arctan(\frac y x) &\text{if } x > 0, \\
\arctan(\frac y x) + \pi &\text{if } x < 0 \text{ and } y \ge 0, \\
\arctan(\frac y x) - \pi &\text{if } x < 0 \text{ and } y < 0, \\
+\frac{\pi}{2} &\text{if } x = 0 \text{ and } y > 0, \\
-\frac{\pi}{2} &\text{if } x = 0 \text{ and } y < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\\
&= \begin{cases}
\arctan\left(\frac{y}{x}\right) &\text{if } x > 0, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y > 0, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y < 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\end{align}$$


Difference between rceps and cceps in MATLAB I'm trying to construct a cepstrum using MATLAB built in function but I don't get the difference between MATLAB's functions rceps and cceps.
Can anyone explain to me the difference between rceps and cceps?

RTFM? rceps says:

y = real(ifft(log(abs(fft(x)))));

while cecps says:

h = fft(x);
logh = log(abs(h)) + sqrt(-1)*rcunwrap(angle(h));
y = real(ifft(logh));


What is the gain (advantage) of oversampling and noise shaping in D/A conversion? It is clear that oversampling and noise shaping in A/D conversion can help in shaping
quantization noise.  But for D/A conversion, normally there is no quantization
and hence I did not understand its use.
What is the gain (advantage) of  oversampling  in D/A conversion?

If you play 16-bit audio at 48kHz, you need the DAC analog reconstruction filter to pass 20kHz and attenuate 96dB at 24kHz, which is quite steep and requires complex multistage analog filter.
The advantage of using oversampling is moving the sampling rate much higher, for example oversampling by 4x means the DAC runs at 192kHz, and the analog filter only needs to pass 20kHz and block 96dB at 96kHz, which allows for much simpler filter.
The noise shaping is much like dithering after signal processing, but instead of adding white noise uniformly to the whole band, it weights the noise to high frequencies above the audio band, which can then be filtered away by the DAC reconstruction filter, so less noise remains on the audio band, thus increasing the apparent signal-to-noise ratio of the audio band.