What should I do know to know the frequency from complex signal? Now I'm trying to digging the DFT function.
But nobody tell about the what do we need it at least to solve the dft from the complex signal- a signal with many frequencies

If I got the complex signal as the below,

then what should we have to know at the first time to do DFT?
Especially, What kind of things should I have to know to solve a signal with many frequencies?

I've seen some parameter from example as the below when I googling.
sampling rate, N, nyquist rate, ..., sort of things.
In the practically field for engineer, if they got a that kind signal, how do they do at the first time?

as I know the below image is DFT equation.

then how do we decide the "N" and "m"?

When you say "complex" signal, I assume you mean just a signal with many frequencies, and not a "complex-valued" signal, in which the the signal values themselves are complex numbers. Let me know if I am incorrect in this assumption. (Actually, computation of DFT is the same in either case).

I'm also assuming your signal is discrete.

The discrete Fourier Transform, $F(\omega)$, can be found by solving the following equation (where N is the number of samples your signal, $f[n]$, contains)

$$F(\omega) = \sum_{n=0}^{N-1} f[n] e^{\frac{-j 2\pi \omega n}{N}}$$

Alternatively, you could use the function fft() on MATLAB.


Kalman Filter - Non Linear Measurement Model I'm new to Kalman filters and estimation in general. I'm running a simple test of an EKF to check my understanding, but I'm getting some odd results with a particular case.

Given a state vector:
$$ \textbf{x} =
\left(\begin{matrix}
p_x \\
p_y \\
p_z \\
\end{matrix}\right)
$$

If I have a sensor that measures these values directly, i.e.:
$$ \textbf{z} =
\left(\begin{matrix}
p_x \\
p_y \\
p_z \\
\end{matrix}\right)
$$

Then I get a converging solution. However, if I decide that my sensor actually measures the square of each value, i.e.
$$ \textbf{z} =
\left(\begin{matrix}
p_x^2 \\
p_y^2 \\
p_z^2 \\
\end{matrix}\right)
$$

Then following the process on the EKF wiki page causes my solution to diverge very rapidly. This seems to be irrespective of the dynamic model of the system, but entirely dependent on whether what the sensor is measuring is linear. Is there a standard approach to follow in this situation?

EDIT: Turns out it was a silly mistake on my part. I just implemented the filter incorrectly. The EKF was valid in this case.

So let's assume your model is:
$$
 \textbf{x}_{k} = \textbf{x}_{k-1} +  \textbf{n}_k
$$
where
$$
\textbf{n}_k =
\left(\begin{matrix}
n_x \\
n_y \\
n_z \\
\end{matrix}\right)
$$
and each of the $n_x,n_y,n_z$ are independent, zero-mean Gaussian noise terms with variances $\sigma^2_{nx},\sigma^2_{ny},\sigma^2_{nz}$.

The output equation is then:
$$
\textbf{z}_k = \textbf{f}(\textbf{x}_k) +  \textbf{m}_k = \left(\begin{matrix}
p_x^2 \\
p_y^2 \\
p_z^2 \\
\end{matrix}\right)
+  \textbf{m}_k
$$
where
$$
\textbf{m}_k =
\left(\begin{matrix}
m_x \\
m_y \\
m_z \\
\end{matrix}\right)
$$
and each of the $m_x,m_y,m_z$ are independent, zero-mean Gaussian noise terms with variances $\sigma^2_{mx},\sigma^2_{my},\sigma^2_{mz}$.

So following the EKF Wikipedia page, that means:
$$
\textbf{F}_{k-1} = \textbf{I}\\
\textbf{H}_{k} = {\tt diag}(2\textbf{x}_k)
$$
where ${\tt diag}$ creates a square matrix from the vector argument.

Implementing this in R (code below) shows the follow example of the true state and the state estimate.  As you can see, it has some problems because the output equation masks the sign of the state (that's why the top plot shows a mirror image around the zero axis).

R Code Below

#30103

# First, construct the signal model and generate the data.
T <- 1000

sigma_n <- 0.1
n_k <- array(rnorm(T*3,0, sigma_n), c(3,T))

sigma_m <- 0.1
m_k <- array(rnorm(T*3,0, sigma_m), c(3,T))

F <- array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))

x_0 <- array(1,c(3,1))
x <- array(0,c(3,T))
x[,1] <- x_0

z <- array(0,c(3,T))

for (k in 2:T)
{
  x[,k] <- F %*% x[,k-1] + n_k[,k]
  z[,k] <- `^`(x[,k],2) + m_k[,k]
}

# Next, form the EKF
H <- 2* array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))

Q <- sigma_m^2*array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))
R <- sigma_n^2*array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))

library("MASS") # For pseudo inverse ginv()

xkm1km1 <- matrix(rep(0,3*T+3),3,T+1)
xkm1km1[,1] <- x_0

xkkm1 <- matrix(rep(0,3*T),3,T)
K <- array(rep(0,3*3*T),c(3,3,T))
Pkm1km1 <- array(0,c(3,3,T+1))
Pkm1km1[,,1] <- array(c(1000,0,0 ,0,1000,0, 0,0,1000), c(3,3))

zhat <- matrix(rep(0,3*T),c(3,T))
err <-  matrix(rep(0,3*T),c(3,T))

for (k in 2:T)
{
  xkkm1[,k] <- F %*% xkm1km1[,k-1]
  Pkkm1 <- F %*% Pkm1km1[,,k-1] %*% t(F) + Q

  H <- 2*diag(xkkm1[,k])

  K[,,k] <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  err[,k] <- z[,k] - H %*% xkkm1[,k]
  xkm1km1[,k] <- xkkm1[,k] + K[,,k] %*% err[,k]
  Pkm1km1[,,k] <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K[,,k] %*% H) %*% Pkkm1
  zhat[,k] <- as.numeric(H %*% xkkm1[,k])
}

p1 <- x
p2 <- xkm1km1
lims <- c(-5,5)
par(mfrow = c(3,1), pty="m")
plot(p1[1,], col="grey",ylim=lims)
lines(p2[1,], col="red")
title("True (grey) and Estimated (red) State p_x")
plot(p1[2,], col="grey",ylim=lims)
lines(p2[2,], col="red")
title("True (grey) and Estimated (red) State p_y")
plot(p1[3,], col="grey",ylim=lims)
lines(p2[3,], col="red")
title("True (grey) and Estimated (red) State p_z")


Beyond cross-eyed 3D (VR application) Making static 3D image is very simple - we just take image pair with something like this:

then we can see 3D image by putting the left and right eye image side by side and crossing eyes or overlapping the images and using color filters for each eye. So far so good.

VR headsets (e.g. Oculus Rift) now allows seeing the 3D image directly.

However, I would like to incorporate slight movements and tilting of the head. Of course, this would immediately spoil the 3D effect.

First idea for better capture is using more cameras (e.g. 2x2 or 3x3 array) to acquire necessary viewpoints (e.g. top and bottom camera is necessary to reconstruct 3D when viewed with head leaned on one side). But how to interpolate between these views based on shifted and rotated head pose? Is depth/3D reconstruction necessary or are there any computational photography techniques to generate appropriate left/right eye images (for a given viewpoint/pose)?

I found various articles regarding depth reconstruction, 3D reconstruction and light fields, but nothing leads me toward the goal of interpolating the views.

Ultimately, I am planning to create a panorama viewer for VR glasses (something like krpano's web viewer supporting Oculus Rift). However, the viewer can account for limited head movements (pan, tilt) and the 3D impression does not hold for leaning, slight head shifts and correct for different eye-to-eye distance.

Any pointing to appropriate articles will be appreciated, thanks!

I have finally found a great overview video on the topic from Richard Szeliski: "Reflections of Image-Based Modelling and Rendering"

The research started with stereo and depth estimation. The depth map can be used to "warp" image in 3D to give it some perceived depth. This does not create parallax effects though.

Another approach called Layered Depth Images used segmentation to several depth layers each having its own depth map. This allows parallax effect to occur, but extra care must be taken to matte foreground and background layers during rendering.

Fast 3D reconstruction of the sorrounding is possible using primitives (boxes, cylinders) matched to scene points. The 3D reconstruction is tedious and requires many images and extensive walk around.

Then, more image-based (rather than geometry-based) approach are presented, which are Lumigraphs and Light Fields. These allow some level of reconstruction of the plenoptic function and solve for rays crossing given viewing plane.

Finally, Concentric Mosaics is the no. 1 solution I have been looking for. A camera travelling around circle (tangential direction) is able to capture 3D plenoptic function of the scene. Each point thus can be viewed from different angles and the viewer can move inside the circle, see parallax, reflections, transparency and occlussion.

Extension is possible by adding support for rays in vertical direction - this will lead to capturing 4D plenoptic function from a camera moving on sphere (tangetial direction) or by using a spherical rig of many cameras.

Further advances in the field are expected with the introduction of plenoptic cameras (e.g. Lytro) that can capture the 4D plenoptic function directly. Some work on making plenoptic mosaics is already there (e.g. "Panorama Light Field Imaging").


How to detect resonance with DFT In the last paragraph of this article, the author gives an example of the DFT of a song and his window's resonant frequency. As he points out (and it is a bit intuitive), if he listens to this song, the window is not likely to break, despite the song's spectrum being centered around the natural frequency of the window.
The question is, if we are given a range of resonant frequencies (the "dangerous" range around the resonant frequency), can we decide whether "the window would break", so to speak, based on the DFT alone, without a physics model (and therefore regardless of the problem, i.e. a window, a vibrating wire, etc.)?
Also, if we have two different signals and we compute their DFTs, how does the real-world physical effect of those signals (say, on a window) compare to their difference in amplitude for the resonant frequency? That is, if the amplitudes are equal, are the effects on the window's vibrations also equal? If one amplitude is half the other one, are the vibrations proportional the amplitudes' ratio (1/2)?

The question is, if we are given a range of resonant frequencies (the "dangerous" range around the resonant frequency), can we decide whether "the window would break", so to speak, based on the DFT alone, without a physics model (and therefore regardless of the problem, i.e. a window, a vibrating wire, etc.)?

No. There is a whole lot of physics between the spectrum of a sound track and what actual force is excerted on a window. Speaker, amplifier, wave propagation in a room, pressure distribution across the pane (both directioin and magnitude), etc.
The whole example is kind of silly. The energy required to break a window is many orders of magnitude larger than what a residential sound system can produce. Your ear drums will break long before the window shows any type of distress, resonance or not.

That is, if the amplitudes are equal, are the effects on the window's vibrations also equal?

Sometimes yes, sometimes no. In order for this to hold we have to require that

All other things are being equal
The physical systems between signal and effect is reasonably linear.


What are the units of my data after an FFT? Magnetometer measures the derivative of the magnetic field, or dB/dt, with an  output in microvolts (mV). The Sampling rate is 128 Hz, so if we collect data for 2 minutes, $2 \times 60 \times  128=15360$ points (discrete case). When I perform an FFT on these time series, what will the units of amplitude density be  after the transform?

In short: if you directly use the voltage output, your FFT amplitude should be in mV as well. If you have the sensor calibration curve, the FFT amplitude should be in Teslas per second (T/s), as you look at a derivative.
If you look at a power density spectrum (squared), then the units above should be squared as well (mV$^2$ or T$^2$/s$^2$).
In general, for a given  unit U, albeit in volt (V), tesla (T), whatever. The FFT sums samples in the original units (U) multiplied by unitless complex values (due to discretization). The units after FFT remain the same as for the original signal, i.e. U.
If you take the absolute value, the same again. For instance, the $0$-frequency index gives you the DC, or average value of your signal  (or at least it is proportional, with a factor of the number of samples, or its square root depending on the normalization).
The continuous Fourier transform, in contrast, "sums" the samples times complex cisoids times a time differential (the $\mathrm{d}t$). Its units would be U.s (s for seconds). Unless you normalize it by something with a unit.


Low Pass Filter - pixel correlation I am trying to understand low pass filters for image processing, I came across an article that stated that low pass filters effect noise more than the real data due to that the image

  Noise always changes rapidly from pixel to pixel because each pixel
  generates its own independent noise. The image from the telescope
  isn't "uncorrelated" in this fashion because real images are spread
  over many pixels.

Does this mean that the filter effects noise more due to the noise being correlated from pixel to pixel?

Most images are inherently low-pass: the object that is displayed in the image generally doesn't change from pixel to pixel.

Most image pixel noise is generated by each individual pixel capturing the image. As a result this noise is generally assumed to be uncorrelated from pixel to pixel.  That means the noise is  (bandlimited) white noise, which has uniform components across all frequencies low to high.

So, when you apply a low pass filter to an image, most of the information about the objects imaged will pass straight through because they are in the low pass filter's pass band.

Conversely, noise will have its high frequency components attenuated by the low pass filter.

That is why the statement low pass filters effect noise more than the real data is generally true.

To answer your explicit question:

  Does this mean that the filter effects noise more due to the noise being correlated from pixel to pixel?

The filter affects the noise more due to the noise being uncorrelated from pixel to pixel.


What methods can be used to remove ringing artifacts in the output of a digital filter? I'm using digital filters to apply spectral mangling-type special effects to audio.

When using a digital filter (vsts/standalone DSP programs/outboard digital filter, etc.), especially when using narrow transition bands/brickwall filters, are there any effective ways to remove ringing artifacts introduced by the filter?  Please bear in mind I am new to DSP and my level of understanding is very basic.

Sticking with linear systems, removing the ringing is nearly the same as adding back some of the spectral content that your really steep transition filters removed.  Why use some crazy scheme to add back the stuff in the "softer" transitions that your hard-edged filters cut out?  Just use a more reasonable total filter response in the first place.

Going to non-linear systems, you could use some sort of AI pattern matching to determine what kind of sound waveforms would be perceived by a human a ringing, and just gate those waveforms out.  But that might just add weird sounding artifacts of its own (as well as also spreading spectral content outside of your really steep filters).


"Dirac Comb" vs "Ones Comb" While learning sampling theory - I noticed that examples of continuous signal sampling always achieved the goal via multiplying the signal with a "Dirac Comb".
I was intrigued by the requirement to use a Dirac Comb instead of a "Ones Comb" - I.E: multiply the signal by the value '1' at the sampling times.

The explanation I got: "That's necessary because multiplying by a Dirac function is equivalent to Convolution in the frequency domain which maintains the spectrum of the original sampled signal".

This made a lot of sense to me...But than, when I proceeded to learn about discrete domain sampling - I was surprised to find out that examples used a "Ones Comb" to sample a signal (not the "Dirac Comb").

Why the difference ? How does discrete time sampling get away with multiplying by a train of discrete ones ?

the integral

$$ \int\limits_{-a}^{a} \delta(x) \, \mathrm{d}x = 1  \qquad \forall a>0$$

but this integral

$$ \int\limits_{-a}^{a} g(x) \, \mathrm{d}x = 0  \qquad \forall a>0$$

where

$$ g(x) = \begin{cases}
1 \qquad & \text{for } x=0 \\
0 \qquad & \text{for } x\ne0 \\
\end{cases} $$

this is why we need this dirac delta thing for the sampling or sifting property:

the integral

$$ \int\limits_{-\infty}^{\infty} f(x) \, \delta(x-x_0) \, \mathrm{d}x = f(x_0) $$

but this integral

$$ \int\limits_{-\infty}^{\infty} f(x) \, g(x-x_0) \, \mathrm{d}x = 0 $$

can you see what's wrong with using $g(x)$ for sampling?


What are the implications of resampling the derivative of a signal with a higher frequency than the original signal? Suppose I have a digital signal measured with sampling time, $T_s=1$ sec. If I take it's derivative, it will, naturally have $T_s=1$ sec. But what are the implications  if I re-sample this derivative to have a lower sampling time, such as $T_s=0.1$ sec?

For example, how is the accuracy or representatives of the derivative influenced by re-sampling?

As Laurent said, no amount of interpolation will restore the information that's not in the original $f_s=1\,\text{Hz}$ signal.

It's been discussed in other posts' comments what a suitable digital differentiator looks like; that's not really easy to answer universally. However, it's certainly a high pass filter. So, I'd say: Try to understand your derivate as a filter operation, if possible. Maybe this gives you a suitable approach on how to implement a "good" (by your application's needs) differentiation.

Also, as Laurent said, when you make assumptions on your original signal, you might find a suitable derivative estimator. The typical approach for DSP (and where Nyquist's theorem mentioned above comes form) is to assume that your signal is a finite sum of complex sinusoids, and hence, every digital signal is equivalently representable by its Discrete Fourier Transform. You might find a derivative according to the differentiation theorem of the Fourier transform, if you want to represent your signal as a sum of orthogonal sinusoids, but you might also represent it according to any other base – and derive a differentiation algorithm that fits your needs best.


Improving Canny Edge Detection and Contours Image Segmentation I am working on an image segmentation project to identify playing cards on varying backgrounds. For my use case I care most that I accurately extract all the cards in the image; false-positives are inconvenient but not fatal. I also care that I can run the segmentation without user interaction (i.e., an algorithm like GrabCut does not work for me).

So far I have the following simple algorithm based on Canny edge detection and contour selection:

# Pre-processing: Convert frame to standard size, 1024x768
gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
gray = cv2.medianBlur(gray, 9)
edges = cv2.Canny(gray, 10, 25)
_, contours, _ = cv2.findContours(edges, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
for c in contours:
    rect = cv2.minAreaRect(c)
    box = cv2.boxPoints(rect)
    area = cv2.contourArea(box)
    ratio = area / size
    if ratio < 0.015: # Any contour large enough is a candidate
        continue
    # Mark this box as possible card

For some (many) inputs it works well enough, or even very well:

However for other inputs the contour detection does not successfully disambiguate cards. The contours are grouped in ways that link cards or result in undesirable bounding boxes:

For the record, here are some other things I have considered:

Using a bilateral filter instead of median blur. Might have improved results a bit, e.g., in the bad example above the first column third row card is now correctly segmented.
Using histogram equalization to enhance contrast. Generally obscured some necessary edges.
Attempting to automatically determine Canny thresholds for each image based on the median. Generally this approach missed many of the noisy card edges that we need to consider.
def auto_canny(image, sigma=0.5):
    v = np.median(image)
    lower = int(max(0, (1.0 - sigma) * v))
    upper = int(min(255, (1.0 + sigma) * v))
    return cv2.Canny(image, lower, upper)
Attempting to use watershed segmentation. The fundamental problem is automatically determining the initial segmentation markers, which turns out to be roughly equivalent to identifying the cards in the first place. Finding "white" is particularly hard because its color varies with the ambient lighting conditions.
Attempting to use mean shift filtering / segmentation. I expected it might denoise the edge detection process by simplifying the image structure but it seemed to have more of a "blurring" effect, especially around reflections.

For reference (or, better yet, if you want to play along at home) I have a set of example inputs and the outputs of my current algorithm: album.

You might try a Hough Transform after your canny-Detection instead of contour detection. (Given that you dont have too much perspective distortion, i.e. the cards are rectangular and not trapezoidal).

Have a look at this related question and the links therein.


Finding Correlation response in fourier domain Lets say I have a system that is trying to find a small image (assume all images are grayscale) within in an image by using correlation. So this system has the baseline image, and I input 5 different images (the five input images are bigger or have more data points than the baseline image) and do correlation of my baseline images against the five input images. From this I would get 5 "tables" of correlation values, and from each tables I select the highest value as the metric to how similar the input image is to the baseline image.

Now lets say I want to upgrade my system to do the correlation in the fourier domain. I pad the smaller image with zeros, do the fft, multiply the data by the conjugate of the baseline to do correlation, and I get the power spectral density at each point through some PSD estimator method, like Bartlett, welch, blackman-tukey, etc.

From my understanding, the more spread out, or the less dense the power spectrum is over a set number of frequency bins, the less correlated the two signals are, and vice versa.

But my question is, in the frequency domain, particularly for a two-dimensional signal, would the metric for measuring how spread the power is simply be the maximum value within the table I would attain after the correlation and PSD estimate, as in the spatial domain case?

Also if you could cite your answers whether it be a book or article, that would be much appreciated

From my understanding, the more spread out, or the less dense the power spectrum is over a set number of frequency bins, the less correlated the two signals are, and vice versa.

I'm afraid you have this backward. Narrow peaks in the original (spatial) domain correspond to strong cross-correlations. Narrow peaks in the spatial domain will be relatively spread out in the Fourier domain. The results of the cross-correlation are usually interpreted in the spatial domain (rather than the Fourier domain).

Transform the image and the (zero-padded) template, multiply in the Fourier domain, and then transform back to the spatial domain. Peaks in the spatial domain correspond to (the possibility of) matches between the template and the image.


Impossible to make a scalogram of a signal with 262144 samples I have a signal with 262144 samples and "sampled" with 200Msps.

I want to use the cwt (Continuous Wavelet Transform), for making a time-frequency scalogram.

Also, I worked in some scripts in python and Matlab, but it's very hard for the machine for make the process and empty the memory of the computer with 262144 samples.

The codes and examples than I run, only work with 1024 samples... 512 samples, and it's very simple for the computer to work with it. But for 262144 SAMPLES !, a 1D signal, it's impossible to work well and output a great scalogram (It's the target).

The signal has frequencies around 10KHz, 20KHz... 200KHz, 1MHz... 2MHz.

I'm really comfortable with Python language, very powerfull and great for data mining, also Matlab, but not eficient in memory use... I think...

Regards and I waiting any suggestions or code in Python, also toolboxes in Python.

Usually this is done by chopping a long signal into shorter windows, doing the transform on each window, then recombining the results.  For finding 10 kHz signals from 200 Msps data, you will need either longer windows than 1k samples, or to downsample the data before doing the transform.


Sub Band DFT response With a DFT lets say I had a signal occupying 0Hz to 2048Hz and i am using 2048 bins.

Can a DFT by design (say FFTW) only calculate the responses in a Sub Band like 768Hz to 1280Hz or -64Hz to -256Hz with the same number of bins?

To put in a context im using FFTW for a project.

Generally if you want to use the same number of bins for a sub-band you should increase the length of the FFT accordingly. That is, say you sub-band width is $BW \;Hz$, then you should increase the FFT length by a factor of $2048/BW$. However, as you can see, this only works if $2048/BW$ is a power of $2$, for example in the first case you mentioned ($768-1260 Hz$). Otherwise, it would required resampling of the signal.


Finding the wobbliest signal Background

I have a dataset of signals that change over time and I'd like to find the "wobbliest" of them. By "wobbliest", I mean a signal that returns consistently to zero, but does not stay at zero. This is illustrated in the figure below by the Green signal:

I do not consider the Blue signal "wobbly", since it stay at zero the whole time. Neither is the Red signal, since it seems to converge to a value.

Current Attempts

I tried getting the frequency information of each signal by doing a Fast-Fourier-Transform on each signal and setting the DC term of each signal to zero to not select for shifted signals. However, this still seemed to favour signals with large unidirectional changes in magnitude (like the Red signal) over the desired wobbly signals (like the Green signal).

Summary

How should I process my signals to find the wobbliest ones? I'm assuming I have to do some sort of operation in the frequency domain, but I don't know exactly what.

By taking the absolute value of the derivative of each signal, I was able to get the index of the wobbly signals I desired.

Here are the top 100 wobbliest signals in my wobbly dataset:

Here are the top 100 wobbliest signals in my non-wobbly dataset:


automatic methods to find the delay and to sync the original pleth signals I have computed the heart rate from two pleth signals (from two fingers -both are acquisition-ed from two different devices) heart rate pattern looks similar with some delay in it i could manually see the delay and add those offsets to the original pleth ...Is there any automated method to find the delay and to add those offsets  to the original  pleth signals in Matlab( both the pleth signals are  downsampled to same 10 Hz)

Compute the cross-correlation of the two signals and look at the lag at which the cross-correlation takes its maximum value - this will give you the time offset providing the best alignment of the two signals.


how to calculate cross correlation of transient signal? there.

Let me have a noisy speech signal $y(k)$ as below.

$$ y(k) = x(k) + v(k) $$

where, $k$ is the discrete-time index, $x(k)$ is clean speech signal, $v(k)$ is a zero-mean random process, uncorrelated with $x(k)$ e.g. white noise.

If we put into this as a vector form:

$$ \textbf{y} (k)  = \textbf{x}(k)+\textbf{v}(k)$$

where,

$$ \textbf{y}(k) = [y(k)  \text{ }\  y(k-1) \text{ }\ ... \text{ }\ y(k-L+1)]^T,$$

$L$ is a length of vector.

If we try to calculate the correlation matrix of size $L \times\ L$, the noisy speech signal can be written as

$$ \textbf{R}_y = E[\textbf{y}(k)\textbf{y}^T(k)] = \textbf{R}_x+\textbf{R}_v.$$

Also, $\textbf{R}_x = E[\textbf{x}(k)\textbf{x}^T(k)]$, $\textbf{R}_v = E[\textbf{v}(k)\textbf{v}^T(k)]$.

However, the speech signal is not the stationary signal but the transient signal.
This correlation matrix also can be used for calculating MVDR(Minimum Variance Distortionless Response) process.

Here is what I want to know.

Because $\textbf{y}(k)$ is not the stationary signal, I think there is a method for calculating the correlation matrix for transient signal.

How can we calculate the cross-correlation matrix for transient signal?

Thank you in advance.

Above article is in the below. If there is any problem, I do delete this post. Thank you.

Jingdong Chen, Jacob Benesty, Yiteng Huang, and Tomas Gaensler(2011), "On single channel noise reduction in the time domain,"  Acoustics, Speech and Signal Processing (ICASSP), 2011 IEEE International Conference on, On page(s): 277 - 280.

Correct speech signal is non-stationary. Hence, you have to stop here as whole dsp is based on stationary assumption, as of now at least to a larger extent. But as we cannot sit idle, we go ahead with quasi-stationary assumption. Take a window (rect/hamming/hanning - these have trade-offs) of 20-30 msec (this is different for male/female/music and is taken as 3-4 pitch periods in speech signal, you can keep it fixed at 30 msec (also referred as narrowband speech processing to get an intuition). Using this window break the speech signal into frames by windowing and shifting the window and extracting a frame for every window. You now have a bag of window length vectors of the given speech signal. Find auto-corr of these for any further processing as within the window the signal is assumed to be stationary (! though not exactly but better than the whole signal).


Principal Component Analysis definition I have just learned about this method, so I am not very familiar with it. As far as I know, Principal Component Anlysis (aka PCA) is used to transform a vector $x$ that belongs to a space of $d$ dimensions to a vector $y$ that belongs to a space of $p$ diamensions, where $p\ll d$ . The new vector $y$ consists of p uncorrelated components. I also know that these $p$ components hold the maximum possible energy of the initial vector $x$ . Last but not least, since PCA is an orthogonal transformation, it should minimize the error norm (the difference between the initial vector and its projection). Correct me if I am wrong.

So let's say we have a random vector $x$ that belong to space $\mathbb{R}^{d}$ and we use PCA to project it to a new space of $p$ dimensions where $p\ll d$. I have come up with three statements. However, I can't find out which one is the correct one:

PCA projects vector $x$ to a space of $p$ dimensions where the projection has minimum energy
PCA projects vector $x$ to a space of $p$ dimensions where the difference between the initial vector and the projection has maximum energy.
PCA projects vector $x$ to a space which is  homeomorphic with $\mathbb{R}^{p}$

So which is the correct statement?

By projecting a vector x using PCA (on the PCs), you maximize the variance in the reduced space. Initially, the space is not optimal in terms of maximizing the variance.

So:

  PCA projects vector  to a space of  dimensions where the
  difference between the initial vector and the projection has maximum
  energy.

(initially the is no maximum variance but after projection the space is optimal in terms of maximizing the variance in the PC defined space)

PS:

Minimizing J1 (the error) corresponds to maximizing the quadratic form e'Se. After demeaning this is e'X'Xe. Set now w = Xe (the projection of the data X). The variance in the projected space is w'w that is equal to the e'X'Xe. Thus, minimizing the error means maximizing varinace


Dominant Frequency Peak Decreases with Increasing Window Size I have a signal that looks like this.

I analyse it using fast Fourier transforms to identify the frequency with the largest peak, which is always close to zero. (There are no other clear peaks.) If I use windows of different sizes (by "window" I mean size of the snippet I use – I don't know if this is the correct use of the terminology, however), the result shows that the dominant peak gets closer and closer to zero as the size of the windows increase. That is to say, if I plot 1/frequency against window size, it is linear.

My initial thought is that if I had an infinitely long window, then the peak would be located at zero implying a DC offset, but the signal is centred on zero in an effort to avoid such an offset. Can someone explain why I see this phenomenon and how I can avoid it?

PS I even tried this using sliding windows so that I can get frequencies based on different parts of the data, but the issue persists.

Since your signal isn't sampled uniformly some strange things might happen when you apply FFT and look at the results.

What you should do is estimate the Uniform DFT of the Non Uniform Time Series.

One easy way to do it is use the reference code and analysis I posted on the question - Frequency Analysis of a Signal Without a Constant Sampling Frequency (Non Uniform Sampling in Time Domain).

Apply it on your data and things should behave more consistently.

One more thing, if you're after a peak, something which will make easier analyzing the data will be centering it. Namely make sure the samples have zero mean.


importance of symmetric windows for STFT for practical applications I am wondering why is the linear phase property of filters with symmetric coefficients important for some practical applications of short time fourier transform, such as image processing? I have read that the linear phase property will ensure that there will be no phase distortion of the signal by the filter? Can anyone explain the reason?

This example may help you to understand this.

The plot below shows the step response for two filters that have the same magnitude response, but the first has non constant group delay (non linear phase) while the second has constant group delay (linear phase).

The first filter is a minimum phase version of the second filter.


How do you map indices to twiddle factors in FFT I am kinda lost on how to map my indices for each stage to my look up array of precomputed twiddle factors.
For an $N=8$ point FFT i have $N/2$ twiddle factors in total:
${
[0] = W_8^0, [1] = W_8^1,  [2] = W_8^2, [3] = W_8^3
}$
With this look up table available how do you calculate, which index matches with each other and which twiddle factor index in the array is needed and whether its negative or positive.
I want to be able to recalculate all of this.
For example:

Lets say as a random choice, I'm on the second stage and on index 3. From this map it shows that:
Stage 2 :: Index 3:: $[1] - W_8^2 * [3]$
So i need to know that I pair with index 1, the twiddle array index is [2] and that the sign for it is negative since i am on the bottom half of the wing.
The butterfly images that every website shows - seems to gloss right over how you actually calculate the mapping of it all.
Does any one know the math that links it up as you iterate the indices for each stage?

Here is C++ code snippet that shows a relatively efficient way of doing it. The idea here is that the step size is always half the previous step and that the index wraps around if it exceeds N/2 which we can implement very efficiently with a bit masking operation (if it's a power of two).
constant int N = 8; // FFT size
constant int moduloMask = (N/2)-1;  // implements modulo N/2 adding

// loop over all stages
int stepSize = N;  // initial step size if FFT length
for (int i = 0; i < numStages; i++)
{
  int twiddleIndex = 0;
  stepSize >>= 1;   // divide by two, it's half of the previous stage
  for (int b = 0; b < numButterflies)
  {
    currentTwiddle = W[twiddleIndex];  // current twiddle factor from table
    // update the index: add step size and warp it around if it exceeds N/2
    twiddleIndex= (twiddleIndex+ stepSize) & moduloMask ;
  }
}


Design discrete controller for zero steady state error I have the following system

where $$G(s)=\frac{0.5}{s+1}+\frac{5}{s+10}$$

How can I design the C(z) controller so that the steady state error for a step input r(t)=1(t) is zero?

I know that this has to do with the system type and in this case we have to deal with a type 0 system which for a step input will give a finite steady state error. Adding an integrator we make the type 1 getting the desired result. Now how would we deal with this in a discrete system? Do I get my constant time controller and convert it to a discrete controller?

You can obtain the steady state error using Finite value theorem (http://www.engr.iupui.edu/~skoskie/ECE595_f05/handouts/fvt_proof.pdf) and then you can compute your controller to make it zero.
or if you have sampling fast enough, you can just discretize your continous version of controller.


DFT and IDFT deriviation problems [![enter image description here][1]][1]

For this problem:
Actually,i think that the question asks us to get the IDFT from the DFT,so I did it..but I am not sure,till now, whether this problems asks us to do transformations or?

What you derive in your answer is that the iDFT of the DFT of a signal is the signal itself. The missing step is that $$ \sum_{k=0}^{N-1}\exp(j2\pi uk/N)=\sum_{m\in\mathbb{Z}}N\delta[u-mN].$$ Then, you get the identity $f[n] = \sum_m f[m]\delta[n-m]$.

However, this does not help for your problem. You need to express the IDFT in terms of DFT operations. So, you need to get

$$
X[m] = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}x[k]\exp(j2\pi km/N)
$$

from a function, which can only calculate

$$
A[m] = \frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}a[k]\exp(-j2\pi km/N).
$$

If you conjugate both sides of the first equation, you get:

$$
X^*[m]=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}x[k]^*\exp(-j2\pi km/N).
$$

This equation is of the structure of the second equation, when you substitute $A[m]=X^*[m]$ and $a[k]=x^*[k]$. Hence, if you perform the DFT the conjugated signal, you get the conjugated of the IDFT of signal:

$$DFT(X^*)=IDFT(X)^*.$$

So, conjugating both sides yields
$$DFT(X^*)^*=IDFT(X)$$
which is exactly what you need to do to solve your problem:

Perform elementwise conjugation of the signal
Perform DFT of the signal
Perform elementwise conjugation of the result

Then, you get back the IDFT of the signal. It's analoguous in the 2D-case (i.e. images). Also, note that you might need to consider a different normalization factor (which would correspond to elementwise division by a scalar), depending on your definition of the DFT. If it is with $\frac{1}{\sqrt{N}}$, then normalization is OK as it is.

Edit: As Royi has pointed out, if your original signal is real, you can omit the 2nd conjugation (since the signal is already real after the DFT, and conjugation won't do anything). As images are most often real-valued, this simplication is very likely to hold here.


Estimating pole radius of formant filter from bandwidth I wish to design an all-pole formant filter by specifying frequency θ (between 0 and π radians for Nyquist) and bandwidth Φ such that the response of the filter at θ-Φ and at θ+Φ should be roughly 50% of the response at θ.

I start by specifying a pole at $p=[r,θ]=re^{jθ}$.

My question is this: How to estimate $r$ as a function of $Φ$?

Note: A fairly crude estimate will suffice for my purposes (vowel synthesis).

This is my approach so far:

Adding a conjugate pole, I construct a transfer function $H(z) = \frac{1}{p-z}\frac{1}{p^*-z}$

$$\frac{1}{H(z)} = (p-z)(p^*-z)$$

$$= pp^* -z(p+p^*) + z^2$$

$$ = r^2 -z(2r cos θ) + z^2 $$

The gain at frequency α will be $\text{Gain}(α) = |H(e^{jα})|$, so the gain at θ is as follows:

$$ \frac{1}{H(e^{jθ})} = r^2 -e^{jθ}(2 rcos θ) + (e^{jθ})^2 $$

$$ = r^2 -re^{jθ}(e^{jθ} + e^{-jθ}) + e^{j\cdot2θ} $$

$$ = r^2 -r(e^{j\cdot2θ} - 1) + e^{j\cdot2θ} $$

$$ = r(r-1) - e^{j2θ}(r-1) $$

$$ = (r-1)(r - e^{j2θ}) $$

So $\text{Gain}(θ) = |H(e^{jθ})| = \frac{1}{(r-1)(r - e^{j2θ})} $ as $0<r<1$. So far so good, nice and tidy!

But now it gets tricky. I require Φ such that:

$$ \text{Gain}(θ) = 2 \text{Gain}(Φ)$$

$$ \frac{1}{|(r-1)(r - e^{j2θ})|} = \frac{2}{|r^2 -z(2r cos θ) + z^2|} \text{where} z=e^{jΦ} $$

$$ 2|(r-1)(r - e^{j2θ})| = |r^2 -e^{jΦ}(2r cos θ) + e^{j2Φ}| $$

Now I should be able to split out real and imaginary components inside each modulus: $2|A+jB| = |C+jD| => 4(A^2+B^2) = C^2+D^2$

But this is starting to look very heavy. How would an engineer approach this?

Is there any way I can do it using iPython/numpy/scipy?

I have here a very pretty answer supplied by @Jazzmaniac, which I will transcribe from our conversation (on IRC Freenode ##dsp).

The key insight is that the gain is inversely proportional to the distance from the pole.

So, double the distance, halve the gain.

(This can be seen by imagining a single pole on the origin with the transfer function $H(z) = 1/z$)

The small circle represents the peak gain (at frequency $e^{jθ}$).

The bigger circle represents a locus of half that gain.

Hence, it is clear that Φ = bandwidth/2


MATLAB code to evaluate audio loundness over time can anyone tell me how to evaluate the loudness/intensity over time using MATLAB? By loudness/intensity over time I mean this:
In the images above, the "black" one on the top is a sound's amplitude over time, and the green line is the corresponding loudness/intensity over time. The green line underneath, which is what I want, is generated using Praat.
Can any one tell me, that given a sound input, how can I generate the corresponding loudness/intensity over time using MATLAB?
For example, given a input, I can generate the upper one using this MATLAB script:
[x,fs] = audioread('soundfile.wav');
t = (0:length(x)-1)/fs;
subplot(211);
plot(t,x);

How can I generate the green line then? Thanks.
(ps. Any other programming languages other than MATLAB is welcome as well)

Any audio recording is generally assumed to have units that are proportional to Pascals (Pa).  This is the instantaneous pressure of the sound wave above (positive) or below (negative) the ambient pressure.

(To actually get Pascals from the WAV file, you need to know the calibration factor for your recording system.)

Loudness of sound (in air) is generally reported as Sound Pressure Level (SPL) in decibels. SPL (in dB) and pascals are related as

Pa_ref = 20e-6;  %reference pressure for SPL
SPL_dB = 10*log10( (Pa / Pa_ref).^2 );

If you use the instantaneous pressure values (ie, the values straight from the WAV file), this equation will give you the instantaneous SPL.  Based on the graph that you showed, you don't want the instantaneous SPL.  You want a time-smoothed SPL.

To get a time-smoothed SPL, the standard approach is to extract the envelope of the pressure signal.  There are lots of ways of doing this.  One way is to rectify and low-pass filter the signal.  In this case, your processing might look like this:

%load and calibrate the data
[wav,fs_Hz]=audioread(myfilename);  %load the WAV file
my_cal_factor = 1.0;  %the value for your system to convert the WAV into Pascals
wav_Pa = wav * my_cal_factor;

%extract the envelope
smooth_sec = 0.125;  %"FAST" SPL is 1/8th of second.  "SLOW" is 1 second;
smooth_Hz = 1/smooth_sec;
[b,a]=butter(1,smooth_Hz/(fs_Hz/2),'low');  %design a Low-pass filter
wav_env_Pa = sqrt(filter(b,a,wav_Pa.^2));  %rectify, by squaring, and low-pass filter

%compute SPL
Pa_ref = 20e-6;  %reference pressure for SPL in Air
SPL_dB = 10.0*log10( (wav_env_Pa ./ Pa_ref).^2 ); % 10*log10 because signal is squared

%plot results
figure;
subplot(2,1,1);
t_sec = ([1:size(wav_Pa)]-1)/fs_Hz;
plot(t_sec,wav_Pa);
xlabel('Time (sec)');
ylabel('Pressure (Pa)');

subplot(2,1,2)
plot(t_sec,SPL_dB);
xlabel('Time (sec)');
ylabel('SPL (dB)');
yl=ylim;ylim(yl(2)+[-80 0]);


Difference(relation?) between filter banks and wavelet decomposition I'm implementing multirate system and met concepts of filter banks and polyphase decomposition of filters.

While trying to understand them, I've found that they look very similar to wavelet decomposition, and now I'd like to understand their relation and differences.

For me, the only thing is kernel types of filters (usually we use sinc kernels in filters, and other, like Haar, in wavelets).

If we stick to the linear version and discrete versions of filter banks and wavelets, filter banks represent  the generic tool, and wavelets can be implemented as a specific instance of iterated $2$-band filter banks satisfying some additional properties, namely that low-pass spaces are embedded dyadically.

In other words: get a single-level $2$-band perfect reconstruction filter-bank, add some regularity so that iterations converge toward something legit, and you have a concrete Discrete Wavelet Transformation.

Kernel types have almost nothing to do with it, from a theoretical perspective. Kernels are chosen with respect to, mainly, implementation (fast or not) or signal properties (regularity).

That was the simple answer. If you want to dig deeper into it, you can look at continuous wavelets (that remain filter banks, as the admissibility condition impose a wavelet to be, somehow, a bandpass), nonlinear wavelets, morphological or hybrid filter banks, etc.

When I was a student, I was teached filter banks by Maurice Bellanger, one (the?) father of polyphase decomposition (Digital filtering by polyphase network: Application to sample-rate alteration and filter banks, 1976). As a young person, I just heard about wavelets, and asked him about it. His answer was: "they are just some special filter banks".


Why could window function help to get more accurate specific frequency amplitude? In order to calculate amplitude of sinusoidal signal at a frequency around 60 Hz (but may be not 60 Hz), window function is adopted while using DFT.

hk is the coefficient of window function and Kdc is the sum of these coefficients.
Np is the number of sample in one period when frequency is 60 Hz.
(XL1, XL2, XL3 are time series sample of three phase electrical quantities.
Sqrt(2) is due to some electrical reason.)
As far as I know, window function is helpful to reduce discontinuity when frequency deviates a little bit. However, I am really confused why this could help to get a more accurate sinusoidal signal amplitude.
I tried to derive, but after getting the steps below, I don't know how to prove V1 is closer to A than V2.
h1 and h2 are two types of window function. h1 corresponds to rectangular window and h2 represents triangle window.

As promised, here is a little derivation to give some insight into the effects of windows on DFT values done in reverse.
Suppose that I have a DFT from some signal.
\begin{equation}
\begin{aligned}
X[k] &= \sum\limits_{n=0}^{N-1} x[n] e^{-i\frac{2\pi}{N}kn}
\end{aligned}
\end{equation}
I look at it, decide it needs a little "smoothing", so I decide to average neighboring DFT values to make a new DFT bin.  Let's call it $Y$.  The weights I'm going to use are just a "twist" of the (rescaled) VonHann values.  If you plug in $\omega=0$ this becomes clear.
Plug and chug.
\begin{equation}
\begin{aligned}
Y[k] &= ( -e^{i\omega}, e^{i\omega}+ e^{-i\omega},-e^{-i\omega} ) \cdot (X[k-1],X[k],X[k+1]) \\
&=  -e^{i\omega} X[k-1] + ( e^{i\omega}+ e^{-i\omega} ) X[k] -e^{-i\omega} X[k+1] \\
&= \sum\limits_{n=0}^{N-1} \left[ -e^{i\omega} e^{i\frac{2\pi}{N}n} + e^{i\omega}+ e^{-i\omega} -e^{-i\omega} e^{-i\frac{2\pi}{N}n} \right] x[n] e^{-i\frac{2\pi}{N}kn} \\
&= \sum\limits_{n=0}^{N-1} 2 \left[ \cos ( \omega) - \cos \left( \frac{2\pi}{N}n + \omega \right) \right] x[n] e^{-i\frac{2\pi}{N}kn} \\
\end{aligned}
\end{equation}
Now, my weighted average has all the sudden become a window function.  Again plug in $\omega=0$ and you can see it is the VonHann.
Time for one of those trig identities.
\begin{equation}
\begin{aligned}
\cos(A) - \cos(B) &= \cos\left( \frac{A+B}{2}  + \frac{A-B}{2} \right) - \cos\left( \frac{A+B}{2} - \frac{A-B}{2} \right) \\
 &= -2 \sin\left( \frac{A+B}{2} \right)  \sin\left( \frac{A-B}{2} \right) \\
\end{aligned}
\end{equation}
Substitute it in:
\begin{equation}
\begin{aligned}
Y[k] &= \sum\limits_{n=0}^{N-1} 4 \sin \left( \frac{n}{N}\pi \right) \sin \left( \frac{n}{N}\pi + \frac{\omega}{2} \right)
x[n] e^{-i\frac{2\pi}{N}kn} \\
Y[k] &= \sum\limits_{n=0}^{N-1} w_{\omega}[n] x[n] e^{-i\frac{2\pi}{N}kn} \\
\end{aligned}
\end{equation}
When you plug in $\omega=0$ you see why the VonHann is also called the sine squared window.
So, here is a family of window functions whose base member is the VonHann.
\begin{equation}
w_{\omega}[n] = 4 \sin \left( \frac{n}{N}\pi \right) \sin \left( \frac{n}{N}\pi + \frac{\omega}{2} \right)
\end{equation}
Like I said, I'm not a window expert as I don't care for them much, but I don't think you'll find these on the list in the references.  Trig equations are slippery though, it could be one of those in different form.
Still it is interesting.
Hope this helps.

For those who say it might be worthless, I say not so fast.
Let $\omega = \frac{2\pi}{N}$
\begin{equation}
w[n] = 4 \sin \left( \frac{n}{N}\pi \right) \sin \left( \frac{n+1}{N}\pi \right)
\end{equation}
Which means is it zeroes out the last sample as well as the first and shifts the center of the window to $(N-1)/2$, the center of the sample set, vs $N/2$ which is the center of the repeat frame.

Extra credit:
Which value of $\omega$ cancels the "The twisting caused by being off bin." which is point #2 in the reference in my first answer.


Why is it easier to extract color in HSV model than in RGB? I am following the open-cv tutorial and they claim that:

  In HSV, it is more easier to represent a color than RGB color-space.

after which they convert the image from BGR color model to HSV, claim some color region as blue and bitmask based on that.

My question is: why is it easier to represent color in HSV model?

Say you have an object that is blue, and you adjust the brightness of a lamp in the room, or you move the object so that it does not get as much light. In RGB, all the three components (red, green, blue) will change but in HSV mostly the V (value) component changes. So to recognize pixels belonging to the object by its color, you only need the H (hue) and S(saturation) components of HSV, in contrast to needing all three components of RGB.

Well actually for almost black pixels the H and S components can be anything due to noise, so you will need to see that V is high enough. An equivalent check needs to be done with RGB too.


What is the difference between a Fourier transform and a cosine transform? In speech recognition, the front end generally does signal processing to allow feature extraction from the audio stream. A discrete Fourier transform (DFT) is applied twice in this process. The first time is after windowing; after this Mel binning is applied and then another Fourier transform.

I've noticed however, that it is common in speech recognizers (the default front end in CMU Sphinx, for example) to use a discrete cosine transform (DCT) instead of a DFT for the second operation. What is the difference between these two operations? Why would you do DFT the first time and then a DCT the second time?

The Discrete Fourier Transform (DFT) and Discrete Cosine Transform (DCT) perform similar functions: they both decompose a finite-length discrete-time vector into a sum of scaled-and-shifted basis functions. The difference between the two is the type of basis function used by each transform; the DFT uses a set of harmonically-related complex exponential functions, while the DCT uses only (real-valued) cosine functions.

The DFT is widely used for general spectral analysis applications that find their way into a range of fields. It is also used as a building block for techniques that take advantage of properties of signals' frequency-domain representation, such as the overlap-save and overlap-add fast convolution algorithms.

The DCT is frequently used in lossy data compression applications, such as the JPEG image format. The property of the DCT that makes it quite suitable for compression is its high degree of "spectral compaction;" at a qualitative level, a signal's DCT representation tends to have more of its energy concentrated in a small number of coefficients when compared to other transforms like the DFT. This is desirable for a compression algorithm; if you can approximately represent the original (time- or spatial-domain) signal using a relatively small set of DCT coefficients, then you can reduce your data storage requirement by only storing the DCT outputs that contain significant amounts of energy.


Differentiating between zero and not sending for OOK How would a trigger works for on-off keying (OOK)? How would it differentiate between not receiving anything and zero?

A simple way to do it would be using a start bit like UART. When idle, don't transmit anything. When you want to transmit something, first transmit a start bit which is a '1' (ON). Then transmit a fixed number of bits known to both the transmitter and receiver. The receiver looks for the start bit and knows to begin demodulating after it receives the start bit.

A more rigorous way to do it would be the following:

Assemble bits into a frame of size X bits.
At the beginning of the frame, transmit a preamble sequence. This is some known sequence of bits that will be the same every time. This is what is done in Wi-Fi, Ethernet, etc.
The receiver looks for that preamble sequence by correlating the received signal with the expected preamble signal. Once you get a high peak in the correlation, the receiver knows that is the start of the frame.
The receiver is now synchronized and knows that it needs to demodulate X bits. After that, the receiver goes back to looking for the start of the next frame.

Raw modulation isn't really useful by itself; you need to add some extra stuff to have a practical communication system.


Filtering data on a known signal I have data that is affected by temperature.  How can I filter out the temperature signal considering I have temperature data collected?

The obvious answer: define a mathematical model of how the temperature affects your data, and use that model to remove the effect of temperature on your data.

Sorry, but with the amount of information you give on your data, the way your temperature affects it or any approach, this is the precisest possible answer.

I'd like to add that what signal processing defines as "filter" is probably not the tool you're looking for, unless you now the spectral properties of your data and your temperature.


Obtain a signal's peak value if it's frequency lies between two bin centers Please suppose the following:

The frequency of a signal's fundamental has been estimated using FFT
and some frequency estimation methods and is lying between two bin
centers
The sampling frequency is fixed
Computational effort is not an issue

Knowing the frequency, what is the most accurate way to estimate the corresponding peak value of the signals fundamental?

One way might be to zero-pad the time signal to increase FFT resolution such that the bin center will be closer to the estimated frequency. In this scenario, one point I am not sure about is if I can zero-pad as much as I want or if there are some drawbacks in doing so. Another one is which bin center I should select after zero padding as the one I am obtaining the peak value from (because one may not hit the frequency of interest exactly, even after zero padding).

However, I am also wondering whether there is another method which may deliver better results, f.e. an estimator which uses the peak values of the surrounding two bin centers to estimate the peak value at the frequency of interest.

The first algorithm that springs to mind is the Goertzel Algorithm. That algorithm usually assumes that the frequency of interest is an integer multiple of the fundamental frequency. However, this paper applies the (generalized) algorithm to the case you are interested in.

Another problem is that the signal model is incorrect.  It uses 2*%pi*(1:siglen)*(Fc/siglen). It should use 2*%pi*(0:siglen-1)*(Fc/siglen) for the phase to come out correctly.

I also think there is a problem with the frequency Fc=21.3 being very low. Low-frequency real-valued signals tend to exhibit bias when it comes to phase/frequency estimation problems.

I also tried a coarse grid search for the phase estimate, and it gives the same answer as the Goertzel algorithm.

Below is a plot that shows the bias in both estimates (Goertzel:blue, Coarse:red) for two different frequencies: Fc=21.3 (solid) and Fc=210.3 (dashed).  As you can see the bias for the higher frequency is much less.

The plot $x$-axis is the initial phase changing from 0 to $2\pi$.


What mathematical tools exist for understanding modulated noise? Suppose we have a signal $n$ that consists of Gaussian white noise.  If we modulate this signal by multiplying it by $\sin 2\omega t$, the resulting signal still has a white power spectrum, but clearly the noise is now "bunched" in time.  This is an example of a cyclostationary process.

$$x(t) = n(t) \sin2\omega t$$

Suppose that we now demodulate this signal at a frequency $\omega$ by mixing with sine and cosine local oscillators, forming I and Q signals:

$$I = x(t) \times \sin\omega t$$
$$Q = x(t) \times \cos\omega t$$

Naively observing that the power spectrum of $x(t)$ (taken over a time interval much greater than $1/f$) is white, we would expect $I$ and $Q$ to both contain white Gaussian noise of the same amplitude.   However, what really happens is that the $I$ quadrature selectively samples the portions of the timeseries $x(t)$ with high variance, while $Q$, ninety degrees out of phase, samples the lower variance portions:

The result is that the noise spectral density in I is $\sqrt{3}$ times that of $Q$.

Clearly there must be something beyond the power spectrum that is useful in describing modulated noise.  The literature of my field has a number of accessible papers describing the above process, but I would like to learn how it is treated more generally by the signal processing / EE communities.

What are some useful mathematical tools for understanding and manipulating cyclostationary noise?   Any references to literature would also be appreciated.

References:

Niebauer et al, "Nonstationary shot noise and its effect on the sensitivity of interferometers". Phys. Rev. A 43, 5022–5029.

I'm not sure specifically what you're looking for here. Noise is typically described via its power spectral density, or equivalently its autocorrelation function; the autocorrelation function of a random process and its PSD are a Fourier transform pair. White noise, for example, has an impulsive autocorrelation; this transforms to a flat power spectrum in the Fourier domain.

Your example (while somewhat impractical) is analogous to a communication receiver that observes carrier-modulated white noise at a carrier frequency of $ 2 \omega $. The example receiver is quite fortunate, as it has its an oscillator that is coherent with that of the transmitter; there is no phase offset between the sinusoids generated at the modulator and demodulator, allowing for the possibility of "perfect" downconversion to baseband. This isn't impractical on its own; there are numerous structures for coherent communications receivers. However, noise is typically modeled as an additive element of the communication channel that is uncorrelated with the modulated signal that the receiver seeks to recover; it would be rare for a transmitter to actually transmit noise as part of its modulated output signal.

With that out of the way, though, a look at the mathematics behind your example can explain your observation. In order to get the results that you describe (at least in the original question), the modulator and demodulator have oscillators that operate at an identical reference frequency and phase. The modulator outputs the following:

$$
\begin{align}
n(t) &\sim \mathcal{N}(0, \sigma^2) \\
x(t) & = n(t) \sin(2\omega t)
\end{align}
$$

The receiver generates the downconverted I and Q signals as follows:

$$
\begin{align}
I(t) &= x(t) \sin(2 \omega t) = n(t) \sin^2(2 \omega t)\\
Q(t) &= x(t) \cos(2 \omega t) = n(t) \sin(2 \omega t) \cos(2 \omega t)
\end{align}
$$

Some trigonometric identities can help flesh out $ I(t) $ and $ Q(t) $ some more:

$$
\begin{align}
\sin^2(2 \omega t) &= \frac{1 - \cos(4 \omega t)}{2}\\
\sin(2 \omega t) \cos(2 \omega t) &= \frac{\sin(4 \omega t) + \sin(0)}{2} = \frac{1}{2} \sin(4 \omega t)
\end{align}
$$

Now we can rewrite the downconverted signal pair as:

$$
\begin{align}
I(t) &= n(t) \frac{1 - \cos(4 \omega t)}{2}\\
Q(t) &= \frac{1}{2} n(t) \sin(4 \omega t)
\end{align}
$$

The input noise is zero-mean, so $ I(t) $ and $ Q(t) $ are also zero-mean. This means that their variances are:

$$
\begin{align}
\sigma^{2}_{I(t)} &= \mathbb{E}(I^2(t)) = \mathbb{E}\left(n^2(t) \left[\frac{1 - \cos(4 \omega t)}{2}\right]^2\right) = \mathbb{E}\left(n^2(t)\right) \mathbb{E}\left(\left[\frac{1 - \cos(4 \omega t)}{2}\right]^2\right) \\
\sigma^{2}_{Q(t)} &= \mathbb{E}(Q^2(t)) = \mathbb{E}\left(n^2(t) \sin^2(4 \omega t)\right) = \mathbb{E}\left(n^2(t)\right) \mathbb{E}\left(\sin^2(4 \omega t)\right)
\end{align}
$$

You noted the ratio between the variances of $ I(t) $ and $ Q(t) $ in your question. It can be simplified to:

$$
\frac{\sigma^{2}_{I(t)}}{\sigma^{2}_{Q(t)}} = \frac{\mathbb{E}\left(\left[\frac{1 - \cos(4 \omega t)}{2}\right]^2\right)}{\mathbb{E}\left(\sin^2(4 \omega t)\right)}
$$

The expectations are taken over the random process $ n(t) $ 's time variable $ t $. Since the functions are deterministic and periodic, this is really just equivalent to the mean-squared value of each sinusoidal function over one period; for the values shown here, you get a ratio of $ \sqrt 3 $, as you noted. The fact that you get more noise power in the I channel is an artifact of noise being modulated coherently (i.e. in phase) with the demodulator's own sinusoidal reference. Based on the underlying mathematics, this result is to be expected. As I stated before, however, this type of situation is not typical.

Although you didn't directly ask about it, I wanted to note that this type of operation (modulation by a sinusoidal carrier followed by demodulation of an identical or nearly-identical reproduction of the carrier) is a fundamental building block in communication systems. A real communication receiver, however, would include an additional step after the carrier demodulation: a lowpass filter to remove the I and Q signal components at frequency $ 4 \omega $. If we eliminate the double-carrier-frequency components, the ratio of I energy to Q energy looks like:

$$
\frac{\sigma^{2}_{I(t)}}{\sigma^{2}_{Q(t)}} = \frac{\mathbb{E}\left((\frac{1}{2})^2\right)}{\mathbb{E}(0)} = \infty
$$

This is the goal of a coherent quadrature modulation receiver: signal that is placed in the in-phase (I) channel is carried into the receiver's I signal with no leakage into the quadrature (Q) signal.

Edit: I wanted to address your comments below. For a quadrature receiver, the carrier frequency would in most cases be at the center of the transmitted signal bandwidth, so instead of being bandlimited to the carrier frequency $ \omega\ $, a typical communications signal would be bandpass over the interval $ [\omega - \frac{B}{2}, \omega + \frac{B}{2}] $, where $ B $ is its modulated bandwidth. A quadrature receiver aims to downconvert the signal to baseband as an initial step; this can be done by treating the I and Q channels as the real and imaginary components of a complex-valued signal for subsequent analysis steps.

With regard to your comment on the second-order statistics of the cyclostationary $ x(t) $, you have an error. The cyclostationary nature of the signal is captured in its autocorrelation function. Let the function be $ R(t, \tau) $:

$$
R(t, \tau) = \mathbb{E}(x(t)x(t - \tau))
$$

$$
R(t, \tau) = \mathbb{E}(n(t)n(t - \tau) \sin(2 \omega t) \sin(2 \omega(t - \tau)))
$$

$$
R(t, \tau) = \mathbb{E}(n(t)n(t - \tau)) \sin(2 \omega t) \sin(2 \omega(t - \tau))
$$

Because of the whiteness of the original noise process $ n(t) $, the expectation (and therefore the entire right-hand side of the equation) is zero for all nonzero values of $ \tau $.

$$
R(t, \tau) = \sigma^2 \delta(\tau) \sin^2(2 \omega t)
$$

The autocorrelation is no longer just a simple impulse at zero lag; instead, it is time-variant and periodic because of the sinusoidal scaling factor. This causes the phenomenon that you originally observed, in that there are periods of "high variance" in $ x(t) $ and other periods where the variance is lower. The "high variance" periods are selected by demodulating by a sinusoid that is coherent with the one used to modulate it, which stands to reason.


What do you mean by radio sensing a channel? I am slowly getting into the field of wireless communications, in particular with WiFi based communications.

802.11 MAC mentions that

  A station wishing to transmit first listens to the medium for a predetermined period and If it senses to be IDLE , then it transmits . If the medium is BUSY , it waits for some random time  and then again senses the channel

I have two questions here:

What exactly do you mean by SENSING here ?
I have a radio A which wants to transmit on 2.45 GHZ. Do I transmit anything dummy data to know if its is buysy or idle ? I am thinking of programming wise
how does radio A knows that channel is busy or idle ? Assuming there is no ACK concept here ?

The station here is continuously monitoring the energy in the medium. Let us say you want to transmit in channel 1 of 2.4GHz band. Your RF hardware is always in receive mode if you don't want to transmit anything. You are continuously sampling the incoming signal in this particular channel. You use ADC to sample and compute the energy $E = \sum |x[n]|^2 $ of these samples corresponding to window of say, 100us. This is a simplified concept of sensing. You don't need to transmit anything to sense the channel. The value $E$ is computed in terms of $dBm$ unit.
After the above step, the $dBm$ value is compared against a threshold specified in standard. For example, let us say the minimum sensitivity needed for receiving QPSK modeulated data is -60dBm. You need to make sure that your transmission is above this. So you sense the channel. If the RSSI (that is $E$ computed above) is more than -60dBm, then you cannot transmit as the channel is BUSY. If it is less than -60dBm, you can transmit. But to be on safer side you can give a margin of say 5dBm. For higher modulation order (like 256QAM) the threshold will be lesser as you need to have very less ambient noise while transmitting 256QAM signals.


On generalizing Gabor-Heisenberg uncertainty in a noisy environment The Gabor-Heisenberg uncertainty represents the fundamental limit of time and frequency resolution one can extract from a signal.

$\sigma_t \sigma_f \geq \frac{1}{4\pi}$

My question is: what happen in a noisy environment?

If a bird sings a perfect constant pitch. Depending of the ambient noise, the frequency estimation of this pitch will be affected by the SNR.

Let say you have a time-varying arbitrary tune (like a bird song) played in a arbitrary colored noise environment. Is there a generalized Gabor uncertainty that takes to account the SNR?

I expected something like this:

$\sigma_t \sigma_f \geq function(SNR)$

Given the downvote(s) and the comments, I'll try to rephrase my question:

Frequency estimation is limited by SNR in a noisy environnent.
Frequency estimation is limited by Gabor-uncertainty.

If statements 1 and 2 are correct, how to reconcile the two statements in one general concept? I suppose, there should be a way to compute the 'frequency-estimation precision' as a function of both SNR and time resolution (maybe with some other assumptions).

This uncertainty thing is often misapplied.  Your statement for #2 is such a misapplication.  The bell curve, aka Gaussian, is an eigenfunction of the FT, (only an approximation in the DFT, but that is a different story).  The uncertainty principle is saying when a shape is taken back and forth from the FT and the narrowest cluster of information is desired, this is accomplished by the Gaussian.  It has nothing to do with parameter estimation, but noise indeed does.

You are conflating different concepts.

My latest blog article gives the actual code for calculating the frequency, amplitude, and phase of a pure tone exactly in the single noiseless pure tone case.  It is also robust in the presence of noise.

A Two Bin Solution

Early on, when trying to show "experts" that I had found exact formula frequencies for tones in the DFT, I was often quoted what you said and they would refuse to even look at my math as if I were proposing a perpetual motion machine.

So many misconceptions exist about the DFT, even among supposed learned folks, and that is scary.

Don't believe me?

Here is an actual quote from an expert who shall remain unnamed.

----------------------------------------------------------
An "exact" computation of frequency based on FTs is not possible for
some very simple reasons:  The Fourier Transform is itself an
estimator, and the time-frequency uncertainty principle means
frequency can't be determined exactly without an infinitely long
observation window.  So essentially all frequency computations of time
series or signals are estimates, it's just that some are more accurate
or more efficient than others.

So claiming that you have an exact solution just isn't going to get
the attention of everyone interested in the field.
----------------------------------------------------------

The only truth in there is the last sentence.  So, I wrote it up in a blog where the math is indisputable and anyone who wants to can verify it for themselves.

No, I've never gotten an apology.


how to understand frequncy response and lowpass filter? I am troubled with understanding freqency response. In lecture, my teacher talked about $H(f)$, which i know means that the response depend on frequency, $f$. But when I read text book, it uses a different notation $H(j\omega)$, which i know that $\omega$ is frequency, but what is the character of $j$ in that case? what does "$j\omega$" mean?
and why is that

?

why $V_L(j\omega)$ could be equal to $V_L e^{j \phi_H}$?

what is meant by $e^{j \phi_H}$?

and how does the first row become the second row?

I know that is a bunch of question, thank you for your effort if you could help me with this.

besides, can anyone explain to my how lowpass filter works? some one tells me that because of frequency response, the unexpected low frequency will cause open loop in RC circuit, so the when the frequency is low, the circuit is open therefore filters out the low frequency, is that right?

$H(j\omega)$ comes from the more general response $H(s)$ which comes from the Laplace transform using the complex unit s which has real and imaginary components given as $$s= \sigma+j\omega$$

This may be too much to swallow if you haven't been exposed to the Laplace Transform yet, but let me make an attempt to simplify it: The Laplace transform decomposes analytic single valued functions into a continuous function of s by correlating the function with all values of $e^{st}$  If you are asking about j you are probably equally confused about $e^{st}$. Note that $Ae^{j\theta}$ is simply a function with magnitude A and angle $\theta$, so if s is only imaginary ($s=j\omega$), then $e^{st}=e^{j\omega}$ and the result is a function whos phase is changing with time. I like to view it as a phasor (vector) of constant magnitude rotating around the complex plane, thus we see that the real part of this is a cosine wave and the imaginary part is a sine wave resulting in Euler's Identity $e^{j\omega}= cos(\omega)+jsin(\omega)$.  When s has real components ($\sigma$) then the time domain function also decays or grows depending on the sign of $\sigma$. See the plot below demonstrating this:

Now to relate this to a low pass filter see this next plot below. This shows the surface plot of H(s) for a low pass filter with two poles in the left half plane:

The frequency response is H(s) when we restrict s to just be imaginary such that $s=e^{j\omega}$, therefore we are seeing the result when we only apply fixed magnitude sine and cosine waves to the input of the filter, omitting all the cases where the input signal may decay or grow- which makes sense when we want to measure the frequency response: apply a fixed frequency at the input and measure how much this one frequency has been attenuated at the output (magnitude response), and how much the signal has shifted in phase (phase response). And repeat this for every frequency (which just means sweep s along the $j\omega$ axis from $-\infty$ to $+\infty$.  This is in fact the Fourier Transform as indicated in the plot above showing an expected response for a low pass filter (in this case one that has peaking near the resonance points based on the pole locations.

Here is another example very close to yours but showing the same relationship of H(s) and $H(j\omega)$

And finally the Bode Plot showing the magnitude and phase versus frequency for this low pass filter, which is specifically the magnitude and phase of H(s) when $s=j\omega$ or otherwise stated more directly as $H(j\omega)$. Electrically what is happening is the capacitor becomes a smaller and smaller impedance as the frequency increases ($X_c = \frac{1}{2\pi fC}$). This capacitor is in parallel with what is typically a resistive load. At low frequencies, DC for example, the capacitor is an infinite impedance or open circuit, as if it is not even there and the voltage at the load will simply be $V_{in}\frac{R_L}{R+R_L}$ but as the frequency increases the output will be $V_{}in\frac{R_L||X_C}{R+R_L||X_C}$ where $||$ represents the parallel combination of the two impedances.  What we see from this is that the voltage at the load will decrease as the frequency increases (low pass filter). At high frequencies (RF and microwave) we deal more with reflection coefficients in that the high frequencies are actually reflected back to the source, but I won't get into that here.


How moving part pixel intensity values of video frames becomes dominant compared to stationary part intensities in reconstructed frames? Hello everyone i want to do dynamic texture video sementation using the Fourier transform in MATLAB. I am applying 3-D fft on dynamic texture video frames (using matlab function 'fftn') and reconstructing video frames back from the phase spectrum only (using matlab function 'ifftn').

I have found moving part pixel intensity values becomes dominant
(means its intensity values are increased so much) compared to
stationary part intensities in reconstructed frames of original
frames .(e.g.in waterfall and traffic on road, water part and moving
car's intensity values are increased respectively compared to the stationary background). how does this
happen?
Also What is the relation between time shift and phase change if we take Fourier transform of video frames?
I want know that what changes will occur in magnitude and phase spectrum amongst the frames and can we derive any relation amongst the interframe magnitude or phase spectrum changes?

Q1. As shown by Oppenheim's experiment, the phase spectrum contains most of the structural information about the image. In 2D this are things like lines and edges. In 3D it is things like lines and edges but also movement. Instead of 2D frames and time, imagine the video as a 3D solid where the z-axis is the frame number. If you took a slice along the z-axis (3rd dimension), movement appears like an edge in the signal.

When you reconstruct using just the phase spectrum, you give all the component sinusoids the same magnitude. This basically normalises the brightness everywhere in the image. Try it in MATLAB:

I = double(imread('YOURIMAGE.TIF'));
f = fft(I);
f = f./abs(f);
I = real(ifft(f));
imagesc(I);

In a 3D image, such as a video, this normalisation is also in the time axis. Hence some smaller amplitude moving parts becoming greater in amplitude (and the opposite for high amplitude parts).

Q2. One video frame in your case is the same as 1 pixel in the x or y directions. The phase change depends on the movement.

Q3. Again that depends on the movement. Slow movement will fall into the lower frequency of the spectrum, high movement into the higher frequency. The spectrum is a global response though. It can only tell you about the entire signal. If you want to do analysis locally (i.e. restricted in location and time) you have to use quadrature filters if you want phase information.


How do I deal with pixel intensity overflow when performing discrete convolution on an image and the result is greater than 255? This is my first post on this forum, and I have been working to come up with a solution for a final project in my hardware design class using system verilog on an Altera Cyclone FPGA.  I've been having trouble finding a solution for overflow when convolving a 3x3 pixel intensity matrix (RGB, values from 0-255) with a specific kernel (blurring, canny edge detection, embossing, etc.)

The question I have is this: What should I do to deal with a result that is greater than 255?  I've seen some references suggesting to "wrap the values around" from 255 back to 0 and up, and also have seen some references suggesting to multiply all pixels by a "shrinking factor" (e.g. .8) to avoid ever having to deal with pixel intensity overflow.

Sadly, I've found few references explaining how to deal with this issue in hardware realization.  Am I missing something here?  Any assistance would be greatly appreciated.

EDIT:

Specifically, what if you convolved a matrix of

[255, 255, 255; 255 255 255; 255, 255, 255] * [1, 1, 1; 1, 1, 1; 1, 1, 1]

?  Wouldn't your result contain values outside of your range of 0-255?

There is no perfect solution to the problem you face. You must make a choice and every choice has pros and cons. For example you can apply saturation thresholding. Or you can use normalization. Saturation threshold clipping is very simple to realize and provide excellent performance provided that intermediate processing stages will not magnify the intensities too much above 255 (or below 0). Otherwise, it will severely distort the image. Normalization guarantees a loss free output but may severely degrade the algorithm quality. Besides it requires you to know in advance the normalization scale, which may not be possible always.


Rectangular vs Hann windowing I'm trying to create a basic digital signal processing pipeline, to manipulate a WAV file.
Nothing fancy, I'm just looking to apply some filtering.

My first attempt was to divide the input signal to N rectangular windows (which do not overlap), convert each window to the frequency domain, apply the filtering of my choice to the spectrum, convert it back to the time domain and insert the resulting window in place of the original one.

(By rectangular window, I mean rectangular in the time domain.)

This worked, however, the only downside was that there was a discontinuity at the edges of each window, severe enough to distort the overall signal and introduce harmonics of significant amplitude.

My attempt to resolve this was to replace the rectangular windows with Hann windows, and apply the same process.
With non-overlapping windows, the obvious result was that the output signal became a set of Hann "pulses", which is as bad as the discontinuities.

I seem to be missing a very basic part of digital signal processing.

I understand the windows should overlap, but how exactly? Should the overlapping portions be averaged? How do we ensure that the amplitude of the resulting signal is unchanged despite Hann window overlapping/averaging?

When you filter a finite length signal, the filtered output will have a larger length, either you implement the filter in time domain or frequency domain. Your desired continuity comes from this excess part which must be added to the result of next segment or window.

For example when you take N samples of a signal and filter it with a filter of length M, the result has a length of N+M-1, where the M-1 excessive samples must be added to the next segment. To do so you have to find your filter response in time domain which has a length of M samples, then zero-pad both your signal and filter to N+M-1 samples then take N+M-1 point FT of both and multiply them.


Maximum likelihood estimation complexity computation I have a basic question about maximum likelihood (ML) estimator and its implementation.

I am trying to simulate a communication system, while using ML at the receiver side to find the transmit sequence. The received signal is given by:

$y=Hx+n$

$y, H, x$ and $n$ indicates the received signal, MIMO channel of $N_t$x$N_r$ transmit and receive antennas, transmitted signal and AWGN, respectively.

I'm trying to detecto $x$ based on ML detector, so the estimated signal $x$ can be written :

where $h_{ln}$ represent the $l^{th}$ column of the channel $H$.  I have read in some books and papers, such as in This paper  that the complexity computation of the above equation is $6N_rN_tM$, where M represent the constellation order of signal, where it was mentioned that $h_{ln}x$ requires $4N_r$ for $Hx$ and $2N_r$ for $|y_n-h_{ln}x|^2$.

My question, in $4N_r$ and $2N_r$ what do the $4$ and $2$ represent?? In other words, the comlexity computation doesn't depend on $x$ too? for example, if we up-sampled the signal $x$ in transmitter or we spread it using Walsh code, the complexity computation won't change as long as we are using the same number of $N_t$ and $N_r$ ??

thank you

Ordinary Least Squares problem
Your
$$\hat x = \arg\min_x \sum_{n=1}^{N_r} \left\lvert y_n - h_n x\right\rvert^2$$
is just a way of saying
$$\hat x = \arg\min_x \left\| y - Hx \right \|$$
and that is called an Ordinary Least Squares estimation problem.
Its objective function is of a quadratic form, so we don't have to actually iteratively optimize or anything. There's an analytic solution to the above equation:
$$ \hat x = \left(H^T H\right)^{-1}H^T y\text.$$
As you can see, this has absolutely nothing to do with the constellation order so far – and it doesn't have to. The constellation order only makes a difference to the decision complexity, and by the orthogonality principle, that's a step to be done after you've estimated the transmit signal.
Don't invert that matrix
Now, it's practically never desirable to calculate a matrix inverse, so the thing that you actually solve is
$$\left(H^T H\right)\hat x = H^T y\text.$$
How complex the solution to that is depends on the properties of your MIMO channel matrix!
The quadratic / full-rank MIMO case
In the (very special) case that $H$ is normal (that implies among other things quadratic, ie. as many receive as transmit antennas $N_r=N_t=:N$), you do that by applying a QR decomposition to $H^T H$:
\begin{align}
\left(H^T H\right)\hat x &= H^T y\\
QR \hat x &= H^T y\\
R\hat x &= Q^T H^T y\tag{1}\\
&= \nu y\tag{2}
\end{align}
Since $R$ is a triangular matrix, $(1)$ s the point where you just backsubstitute. Luckily, $R$, $Q^T$ and $H^T$ only need to be calculated once ever, so we can estimate many $\hat x$ from many $y$ as long as the channel doesn't change at fixed complexity.
Complexity
You do the usual paper trick of ignoring the fact that you need to acquire $H$ first. That can, and quite possibly will, dominate the complexity of this (so don't be a bad scientist – figure out how you get your $H$ and how compelx that is, if you're writing some kind of publication).
So, assuming you already have the $H$, the operations you are doing per channel realization (i.e. only once every time $H$ changes, which probably doesn't happen very often):

calculating one half of the hermitian matrix $H^TH$ (no need to calculate the whole thing, it's redundant).
a QR decomposition of a hermitian matrix ($H^TH$)
calculating $\nu = Q^TH^T$

The operation you do per receive vector (i.e. for every symbol time step, typically, so pretty often):

Backsubstituting in $R\hat x$ to finally get $\hat x$.

The complexities per channel realization are:

$\frac{N^3}2$
depends very much on your QR decomposition algorithm (which depends on matrix size, and need for numerical stability, and your platform), but typically no worse than $N^4$
$N^3$

So that's a total of $C_H = \mathcal O \left(\frac 32 N^3 + N^4\right)$ multiply-accumulates.
The complexity  per receive vector:

$\frac12 N^2$

So, $C_y = \mathcal O \left(\frac12 N^2\right)$.
If we introduce a coherence time (measured in symbol durations) $T_C$ for which $H$ doesn't change, you get an overall complexity of
\begin{align}
C_{\text{overall}} &= \mathcal O \left(\frac {\frac 32 N^3 + N^4}{N T_C} + \frac1{2N} N^2\right)\\
 &= \mathcal O \left(\frac {\frac 32 N^2 + N^3}{T_C} + \frac12 N\right)
\end{align}
per receive symbol (i.e. per element of $\hat x$).
Bringing the constellation order into this
As said above, the constellation order only makes a difference to the complexity of the decider. However, we know under ML rule that the likeliest transmit symbol vector is the one geometrically closest to the likeliest transmit symbol vector – and since all elements of that are to be assumed independent from each other, that minimum can be found individually. The complexity of the decision thus only depends on the (SISO) constellation decoder you pick; nothing has been said about that, and there's immensely different implementations of these, so I will not allow myself to give you a complexity estimate here – I'll just allow myself that I can't see the formula as given by the paper appearing anywhere in my head.
It is a very valid question to ask whether it makes sense to apply the ML detection principle to individual elements of the transmit / receive vector. Unless you're in the very massive MIMO case, it makes sense to assume that there's channel coding on these symbols, and that they are, in fact, not stochastically independent from each other. So, classically, a MLSE-based (typically: Viterbi) soft-input channel decoder would be what you use, or some kind of message-passing based channel decoder.


Transform one-pole cascaded system into a biquadratic section I have two cascaded one-pole sections with the same coefficient b for both sections.
$$
y_A[n]=x_A[n]+b\cdot(y_A[n−1]−x_A[n])\\
y_B[n]=y_A[n]+b\cdot(y_B[n−1]−y_A[n])
$$

The final output of the system is $y_B[n]$.

I would like to re-implement the system as one biquad section of the form:

$$
y[n]=a_0x[n]+a_1x[n-1]+a_2x[n-2]-b_1y[n-1]-b_2y[n-2]
$$

The final output of the system is $y[n]$.

I would like to know what calculations I have to perform in order to make the output of the biquad exactly match the output of the one-pole cascaded system.
$$
y_B[n] = y[n]
$$
Obviously the 5 biquad coefficients $a_0\; a_1\; a_2\; b_1\; b_2$ will have to be determined as a function of the $b$ coefficient as found in the cascaded system above.

Can someone lead me through the math?

Thank you.

A single one-pole filter satisfies the following difference equation:

$$y[n]=x[n]+b\left(y[n-1]-x[n]\right)=(1-b)x[n]+by[n-1]\tag{1}$$

which corresponds to this transfer function:

$$H(z)=\frac{1-b}{1-bz^{-1}}\tag{2}$$

If you cascade two of these systems, the transfer functions are multiplied:

$$H^2(z)=\frac{(1-b)^2}{(1-bz^{-1})^2}=\frac{(1-b)^2}{1-2bz^{-1}+b^2z^{-2}}\tag{3}$$

The transfer function in $(3)$ corresponds to the following difference equation for the cascaded system:

$$y[n]=(1-b)^2x[n]+2by[n-1]-b^2y[n-2]\tag{4}$$

Consequently, implementing $(4)$ is equivalent to implementing $(1)$ twice.


Decimation of signal. How to plot frequency spectrum I have an FDM signal from 0 to 70 MHz. The total signal contains 4 FDM signals from (0-fs/4), (fs/4-fs/2), (fs/2-3fs/4) & (3fs/4-fs).
I used fvtool command in matlab to plot these signals.
Then i downsampled the total signal by 4. Downsampling is time is expansion in freq spectrum.
Thus on paper my new signal should occupy freq (0-fs), (fs-2fs), (2fs-3fs) & (3fs-4fs)
When i plot this downsampled signal using fvtool again, i am just able to see my 1st FDM signal ie b/w (0-fs/4) appearing in the plot.
How to view the entire spectrum from (0-4fs) in fvtool.

The root of your problem is that you are not familiar with the Nyquist Sampling Theorem.

First of all, I question your premise that your sample rate is 70 MHz and that you can represent your signals, which range from 0 to 70 MHz, at this rate.  The Nyquist sampling theorem indicates that you can only represent signals up to half the sample rate, which is 35 MHz in your case.

When you downsample by a factor of four you are reducing your sample rate by a factor of four, which in turn reduces the frequencies that you can represent by a factor of four.  You are thus aliasing your signals, causing them to all be at the same frequency.  Thus, the one signal that you see after you downsample is really all four signals packed into one slot.


DFT analysis of a physical signal what is the reason for the result that has been obtained from the DFT of the signal of these two images. The first image is a DFT analysis of a physical signal u which has a time history shown below sampled at a 1000Hz the second image is the square of the signal u which has a time history shown in the image after the time history of u

The first DFT plot represents a typical case of Gibbs phenomenon.


Time invariance from convolution integral I wanted to prove, time invariance property of system by convolution integral(i.e equal time shift in input result in equal shift in output) but as calculation shown in figure, I got unexpected result, so what I'm missing?

If your delayed input is $\hat{x}(t)=x(t-\alpha)$, then $\hat{x}(t-\tau)=x(t-\tau-\alpha)$, and not $x(t-\tau+\alpha)$. From here everything will work out as expected.


Weird ringing at signal start and end for default Matlab lowpass filter **FIND UPDATES BELOW THE ORIGINAL QUESTION TEXT:
I am trying to lowpass filter electrophysiological data and I see ringing near sharp transients at signal start and end.
Ringing at start:

Ringing at sharp transients:

I am simply using Matlab's default lowpass filter:
fpass=1000;
f_sampling=20000;
fil_X=lowpass(X,fpass,f_sampling);

Two questions:

Why do I see the weird ringing effect at signal start and end?
What strategies can I use for reducing the ringing (both at start/end and at sharp transients) other than reducing the steepness of the filter in the frequency domain? Can some fancy windowing be of help? I tried forcing Matlab to use IIR filter with filtfilt and the ringing stays the same.

Thanks,
Anand
UPDATES:
Padding the signal at start and then removing the padded portion from the filtered trace gets rid of the ringing at start. Thanks!!
For the other problem, I tried an IIR filter without filtfilt. Ignore the unnecessary double filtering. I do it because lowpass automatically applies filtfilt.
[~,ft_obj]=lowpass(X,fpass,f_sampling,'impulseresponse','iir');
fil_X=filter(ft_obj,X);

Here is what I get:

As you can see, it gets rid of the pre-transient ringing, while introducing post-transient ringing. The post-transient ringing should not be a big issue. But, there is also a signal delay/distortion that becomes clear at the same time point across different conditions (red line). That can shift the absolute and relative timing of events in the signal.
Is there a way to compensate for that?

You are using a linear phase filter. These are non causal and hence you get pre-ringing at your transients. Use a minimum phase or IIR filter. Do NOT use filtfilt, since that makes an IIR filter linear phase again and you end up with the same problem
The artifact at the beginning is the step response of the filter. The filter needs access to "past" samples and since it doesn't know what these are, it assumes they are 0. In your case, you have a DC jump and it looks like they are more like -0.05. There multiple work-arounds. In your case it's probable the easiest to pre-pend the signal with a few hundred constant samples that are the average of the first few hundred sample of your signal, filter the whole thing and cut the pre-pended section off again.


How to warp a pixel patch to another camera I have two cameras with rotation and translation starting from the world origin and a patch in the first camera's reference coordinate system. I want to warp this patch into the second camera's reference coordinate system. My idea is to transform every pixel like this:

$$T=T_2-T_1$$

$$R=R_1^{T}R_2$$

$$x' = \left[(x  K_{inv} R) + T\right]\;K$$

Where $R$ and $T$ are the transformation matrices between the cameras and $K$ is the affine transformation matrix that describes the camera properties.

Is this correct? I know that when you go from 3D to 2D you need to project, but that matrix is not invertible, so you cannot really do it here? Am I missing something really big?

$x  K_{inv}$, with $x(3) = 1$, produces the pixels as 3D points "stuck on the screen" of the first camera. After transforming you get the same 3D points, but in the coordinate frame of the second camera, and then you project in that camera. So, overall, what you are doing is looking at the first camera from the point of view of the second one.

I suspect that is not what you want....


Pulse Doppler Radar: time delay in down conversion My understanding of Pulse Doppler Radar is the system measures Time Of Flight for the pulse.

My confusion is that during the down conversion process prior to ADC the lowpass filter causes some time delay in addition to the time of flight. In principle the delay should be known beforehand but I'm sure it's not exact and in general the delay is not completely constant even among the band of interest.

How is this delay dealt with?

The delay is dealt with by knowing it. It is calculated and measured and any residual error is minuscule compared with any other uncertainty.

It’s dealt with design and calibration.


Teager-Kaiser Operator vs. Hilbert Transform Since a couple of months I started working on the extraction (estimation) of signal frequency and amplitude components by means of two different time-frequency approaches, namely the Hilbert transform and the Teager-Kaiser energy operator.

I tested both methods on standard signals, such as chirps, sine, cosine, etc. This preliminary analysis seemed to highlight Teager-Kaiser has better resolution in both frequency and amplitude estimation.

Afterwards I applied the two methods on some acceleration signal derived from a tool simulating the dynamic behavior of a given system. Surprisingly enough, Hilbert transform provided more reliable results in terms of frequency estimation. Teager-Kaiser operator shows some high-frequency content estimation, which is hardly expectable.

I developed the two techniques in Simulink, as visible here:

The results on a synthetic acceleration signal are here:

What did I miss in the design of the Teager-Kaiser operator? Is it any possible to avoid those "false" frequency estimates (those with the arrows)?

This paper may be of interest:

David Vakman, "On the Analytic Signal, the Teager-Kaiser Energy Algorithm, and Other Methods for Defining Amplitude and Frequency." IEEE Trans. Signal Processing. (1996)

Summarising from the paper:

$$\Psi(u) = a^2w^2 = [u'(t)]^2 - u(t)u''(t)$$
$$\Psi'(u') = a^2w^4 = [u''(t)]^2 - u'(t)u'''(t)$$

where $u(t)$ is the signal, and $a$ and $w$ iare the amplitude and frequency estimates to solve for.

$$a(t) = \frac{\Psi(u)}{\sqrt{\Psi(u')}}$$
$$w(t) = \sqrt{\frac{\Psi(u')}{\Psi(u)}}$$

But this means if at some point in the signal $u'(t)=1$ and $u''(t)=u''''(t)=0$, for example, then $\Psi(u) = 1$ and $\Psi'(u) = 0$ and thus $a(t) = \infty$ and $w(t) = 0$.

The paper thus explains this causes spikes in amplitude and zero frequencies with the Teager-Kaiser algorithm.

For your signal, the reverse must be true. If $\Psi(u) = 0$ and $\Psi'(u) = 1$ then you would get spikes in frequency when amplitude is zero. This would occur when $[u'(t)]^2 = u(t)u''(t)$

The Hilbert transform uses the entire signal to compute. It is implemented by doubling the positive frequencies and setting the negative frequencies to zero. Therefore the energy of the signal does not change (Parseval's theorem). Thus, intuitively, there can be no infinite spikes, as this would require infinite energy.