What is the minimum number of time points recommended for accurate instantaneous phase estimates using the Hilbert transform I have two time series.  The two series look like they follow the same function, but one appears to lag slightly behind the other.  I would like to determine if there are phase differences between the points in these time series, but each time series is only 4 points long.  Is this enough data to get accurate phase differences using the Hilbert transform?

Forgive me for not having a numerical answer to the question. But here are a few things to think about. The answer depends on how the Hilbert transform (HT) is implemented. If the HT is implemented using an $N$-point tapped-delay line filter-like block diagram the first and last $N$-1 output samples are incorrect. So in this situation four input times samples is definitely not enough samples to produce useful results.

If HT is implemented by way of a forward FFT, zero the neg-frequency spectral components, and inverse FFT (as is done in Matlab) you have a better chance of computing meaning results. However, in this second implementation the input time sequence must be long enough to contain at least one cycle of all the spectral components of your input signal in order to compute somewhat meaningful results. It seems to me that four times samples is too short to reliably estimate any signal parameters of most real-world (information-carrying) signals.


Code Snippet to Convert ADC Samples to I/Q I'd like to pull a narrowband signal out of ADC samples, where the ADC rate is much higher than the bandwidth/central frequency of the narrowband signal.  I'd like to feed the narrowband signal into a software package that want I/Q samples, so I'd like to know if a code snippet is doing the right thing.  There are lots of descriptions online, which seem to be ambiguous about factors of 2, and I've so far failed at finding a code example.
Let's say the ADC operates at sampling frequency $f_{samp}$, and the narrowband signal has central frequency and bandwidth $f_0$ and $bw$.  I think what I'm supposed to do is take my ADC samples, multiply by cos/sin($2\pi f_0 t$), then lowpass/downsample that to get I/Q. What exactly does that lowpass/downsampling look like, though?  I would have naively expected to set the output sample frequency to be the Nyquist rate for the desired bandwidth for both I and Q.  But... that just seems crazy?  I can't imagine people would carry around two Nyquist-sampled timestreams for one signal, so I have to be missing something.  Since nothing is as concrete as code, here's a snippet (ignoring edge effects, windowing, etc.) that does what I think I've read but which I firmly believe must be wrong.  Any suggestions on what it should actually look like would be most appreciated!
import numpy as np

def adc2iq(samples,f_samp,f0,bw):
    #Convert real ADC samples taken at sampling frequency f_samp to I/Q
    #data with full bandwidth bw centered at f0

    n=len(samples)
    #make vectors of sin/cos at frequency f0
    t=np.arange(n)/f_samp
    mycos=np.cos(t*f0*2*np.pi)
    mysin=np.sin(t*f0*2*np.pi)

    #take the Fourier transform of the samples times the sin/cos vectors
    ift=np.fft.rfft(samples*mycos)
    qft=np.fft.rfft(samples*mysin)

    #low-pass/downsample by inverse FFTing a suitable number of samples
    keep_frac=2*bw/f_samp  #is this right?
    kmax=np.int(keep_frac*len(ift))
    i=np.fft.irfft(ift[:kmax])
    q=np.fft.irfft(qft[:kmax])
    return i,q

(In case anyone is curious, my motivation is trying to process GPS signals inside 550 MHz worth of L-band radio astronomy data, and the GPS software seems happiest with I/Q samples)

Consider a real carrier at frequency $f_o$ with phase modulation given by
$$s(t) = 2\cos(2\pi f_o t + \phi(t))$$
Where for simplicity of explanation the amplitude was normalized to be $2$.
This real carrier is at positive and negative frequencies as given by Euler's formula, both of which contain the phase modulation (conjugate symmetric):
$$s(t) = 2\cos(2\pi f_o t + \phi(t)) = e^{j( 2\pi f_o t + \phi(t))}+  e^{-j (2\pi f_o t + \phi(t))}$$
If we multiply the signal by $e^{-j2\pi f_o t}$, which is by Euler's formula $\cos(2\pi f_o t) - j\sin(2\pi f_o t)$, it will move the positive frequency to baseband (DC) and the negative frequency to a negative frequency at twice the rate of the carrier (which we must filter out):
$$e^{-j2\pi f_o t} s(t) = e^{-j2\pi f_o t}e^{j( 2\pi f_o t + \phi(t))}+  e^{-j2\pi f_o t}e^{-j (2\pi f_o t + \phi(t))}$$
$$= e^{j \phi(t)} + e^{-j(4\pi f_ot + \phi(t))}$$
Thus the low pass filter is very wide, passing bandwidth $bw$ at baseband and filtering out the signal up at $-2f_o$.
This is the primary task in recovering the signal at baseband (here in this example the complex signal $e^{j\phi(t)}$ with $I = cos(\phi(t))$ and $Q= sin(\phi(t))$. Any further down-sampling to lower sampling rates after the filtering is just best practice since it is desirable in terms of cost, size, and power to run at the lowest sampling rate possible (minimizes resources needed).
In general, if a modulated signal does not have a symmetric (complex conjugate symmetric) spectrum about its carrier then a complex signal must be used to represent that signal at baseband (using two real signals as I and Q), since any real signal has a complex conjugate symmetric spectrum so one real signal alone would not be able to represent the asymmetric spectrum. The GPS C/A code waveform as the OP is using is however symmetric but if there is any residual small frequency offset (due to Doppler and local clock offsets) then there will be a small positive frequency spectrum offset mostly overlapping with the small negative frequency offset with only one real signal at baseband. By having a complex baseband we can represent that spectrum offset in frequency and therefore the receiver can properly remove that- so the complex baseband processing simplified the receivers ability to resolve frequency and phase offsets between the received signal and the receiver’s local reference.
Please see this additional post for further examples of translating to other frequencies with complex signals:
Frequency shifting of a quadrature mixed signal


Pitch Detection: HPS not working I'm trying to code the HPS (harmonic product spectrum) algorithm and the problem I'm facing is that using my guitar, the fundemental frequency (Here it's 82Hz) just isn't there, so is the 5th Harmonic at 410Hz (Actually it's the lowest local point).
In my algorithm, I did the step of finding the max of all (Mul(harmonic-multiplications))
What pre/post-process should I do to get to 82Hz here.

Missing f0 is a common situation. Pitch could be defined the lowest common multiple of the spacing between spectral peaks, and they seem to be spaced by 82 Hz or 164 Hz here - so there's nothing wrong with your signal and spectrum. Harmonic sum or product should deal with this case. Could you please post plots of your harmonic sum or product? Which result do they give you?

You can also try looking at another pitch detection method like the Average Magnitude Difference Function or YIN.


Measure and simulate non-linear distortion We can capture frequency response, phase response and reverb in recording of impulse response. Then we can use it to simulate audio system and/or space by convolving it with dry audio signal.

Is there similar technique for measuring and simulating non-linear distortion?

I thought of measuring how the amplitude of output signal corresponds to amplitude of input signal. Then I could simulate this system by waveshaping. Will it work? Is there any better method?

Nonlinear systems are very hard to classify and there's no unified theory like for linear systems. In general, you cannot measure/identify non-linear systems in finite time. There are  some specific classes of nonlinear systems which allow for identification or at least approximation.

You already named a trivial one: The memoryless nonlinear system. In this case you can learn the nonlinear transfer function (i.e. your waveshaping function) by just looking at the map from input to output.

Any interesting real system does have memory however, and there we can start with a specific model for the system and adapt the parameters until the model output and the real output match. There are nonlinear optimization techniques that will help to find the proper set of parameters.

If you know little about the internals of the system you want to simulate, then you can use a generic model. The most popular such model is described by the Volterra series and requires that the system is expandable in a series similar to the taylor series but involving convolution terms and that this series converges quickly enough (and in the entire domain you're interested in) to be practically useful. The theory of Volterra kernels delivers the mathematical tool to find the series expansion from input/output observations without the need for non-deterministic optimization techniques.

So summing up, your best chance to really capture the behavior of a nonlinear system is to understand the system in detail and model its components. The black box approach is very difficult and only works under very specific conditions.


28-Band Graphic Equalizer in Matlab I'm trying to design a MATLAB audio player that includes something like a 31-band equalizer. Right now I have the filters working by using the FIR1() function to create an FIR bandpass filter at each frequency specified on the bands. This is great except that it's not exactly what I'm wanting to do. I want to either boost or cut the desired frequency, not isolate the frequency.

How would I go about using these FIR filters in parallel to boost or cut the signal instead of just isolating the signal?

I would suggest the peakingEQ biquad filter by Robert Bristow-Johnson, you can use it to increase or attenuate one frequency or band-frequency, this can be easy to implement in matlab, I used to test and compare some filters, and it seems to work very well to equalizer, you can choose one dB input to do the work !


Dark image separation So I have an image.
I take it into HSV format and v has its intensity values.

now the intensity values of v vary from 0 to 1 (Black to White)

I want to divide v into 3 parts:

1] v1= with intensity values 0<v<0.3

2] v2= with intensity values 0.3<v<0.75

3] v3= with intensity values 0.75<v<0.1

How do I do this?

All you need to do is masking on your v we remember and image is a [rows,cols,layers] matrix, where the layers are Hue=1,saturation=2,intensity value=3

We then create a mask (array of 0s and 1s). For example your first range would be

v1 = (hsv_im(:,:,v_idx) <0.3);

we are taking the all rows, all columns, of layer v, and if the value is <0.3 make it 1, otherwise make it 0. Now to see only those pixels we can multiply the original value array by our mask.

hsv_im(:,:,v_idx) = v0.*v1_mask;

As you recall any value >0.3 was 0 in our mask, so all the other areas will be blacked out.

here is a full example

im = imread('flower.jpg');
hsv_im = rgb2hsv(im);
v_idx =3;

v0 = (hsv_im(:,:,v_idx));    %the original
v1_mask = (hsv_im(:,:,v_idx) <0.3);
v2_mask = (hsv_im(:,:,v_idx) <=0.3 & hsv_im(:,:,v_idx) <0.75);
v3_mask = (hsv_im(:,:,v_idx) >= .75);

%plots the images
figure(1);
subplot(2,2,1);imshow(hsv2rgb(hsv_im));title('original')

hsv_im(:,:,v_idx) = v0.*v1_mask;
subplot(2,2,2);imshow(hsv2rgb(hsv_im));title('0.00< v<0.30')

hsv_im(:,:,v_idx) = v0.*v2_mask;
subplot(2,2,3);imshow(hsv2rgb(hsv_im));title('0.30<=v<0.75')

hsv_im(:,:,v_idx) = v0.*v3_mask;
subplot(2,2,4);imshow(hsv2rgb(hsv_im));title('0.75<=v<1.00')

%plots the masks
figure(2)
subplot(2,2,1);imshow(v0);title('original v')
subplot(2,2,2);imshow(v1);title('0.00< v<0.30 mask')
subplot(2,2,3);imshow(v2);title('0.30<=v<0.75 mask')
subplot(2,2,4);imshow(v3);title('0.75<=v<= 1  mask')


matched filtering in GSM After the estimation of channel impulse response (using the training sequence in burst), to cancel the effect of channel on the received signal, I think, we should convolve the received signal with the inverse of channel impulse response.

But,in GSM, it convolves the received signal with the conjugate of channel impulse response and calls this process as "matched filtering". Could anyone help me with this?

First, please read this answer of mine for a detailed description of matched filters for real-valued signals. In particular, note that what I called the matched filter for a signal $x(t)$ is a(n LTI) filter with
impulse response $h(t) = x(-t)$ which is better described as the time-reversed signal rather than the "inverse" of the impulse response as you call it: to many people, inverse would suggest $[x(t)]^{-1}$ or $-x(t)$
depending on whether they were thinking of multiplicative or additive inverse.. Note also that convolving the channel output with the
time-reversed channel impulse response does not "cancel" the effect of the channel on the signal at the channel input. In general,
matched filtering distorts the signal being filtered tremendously because
the goal of matched filtering is to produce a peak response at the desired sampling time. This is quite different from the audio enthusiast's desire of
reproducing the input signal with the highest possible fidelity, and
thus using an "inverse" filter to remove all traces of the distortion
caused by the channel filter.  See
the last part of my answer cited above for some pictorial illustrations of how matched filtering distorts the signal so as create a peak response
at the sampling instant.

Second, there are no complex signal in real life; they exist only in the diseased imaginations of DSP engineers.  What we have in real life
is a real-valued bandpass signal
$$x(t) = x_i(t)\cos(2\pi f_0 t) - x_q(t)\sin(2\pi f_0 t)
= \Re\left(x_c(t) e^{j2\pi f_0t}\right)$$
where $x_i(t)$ and $x_q(t)$ are slowly varying signals compared to
the carrier signal of frequency $f_0$ Hz, and
$x_c(t) = x_i(t) + jx_q(t)$ is called the complex baseband equivalent
of the RF signal $x(t)$. The filter matched to $x(t)$ thus has
impulse response
\begin{align}
x(-t) &=  x_i(-t)\cos(-2\pi f_0 t) - x_q(-t)\sin(-2\pi f_0 t)\\
&= x_i(-t)\cos(2\pi f_0 t) - [-x_q(-t)]\sin(2\pi f_0 t)\\
&= \Re\left([x_i(-t) -j x_q(-t)]e^{j2\pi f_0t}\right)\\
&= \Re\left([x_c(-t)]^* e^{j2\pi f_0t}\right)
\end{align}
The complex baseband equivalent of the matched filter impulse
response is thus $[x_c(-t)]^*$ and thus when the DSP engineer
proceeds to match-filter the complex baseband equivalent input $x_c(t)$,
she must convolve $x_c(t)$ not with $x_c(-t)$ but rather with
$[x_c(-t)]^*$, cos doing the former instead of the latter
would result in her being in a state of sin.


Value of power spectral density $N_0$ or effect of scaling bandwidth to SNR I want to compare performance over AWGN channel of two systems that are exactly identical but their bandwidth. I will talk only about baseband system for sake of simplicity. For more details about notations, please take a look at this online course MIT DSP course.

Call $P$ the signal power (joules per second).

First, considering a system using bandwidth $W$. According to sampling theorem, the orthogonal basis set contains $2W$ vector. The single side power spectral density of white Gaussian noise $N(t)$ is defined as $N_0$ if its projection $N_k = <N(t), \Phi_k(t)>$ is set of i.i.d Gaussian variables with zero mean and variance $N_0 / 2$.

Sampling interval is $1/2W$, there is $2W$ signal symbol per second, then a symbol has power of $P/2W$. A noise sample has power of $N_0/2$ by definition, the $SNR = (P/2W)/(N_0/2) = P / N_0 W$. Well-known results !!!

Next, I wonder what happens if I increase the bandwidth, say double it, to $2W$ in the same physical environment.

Double bandwidth, half symbol duration. In one second, I have doubled the number of discrete symbol to $4W$. The signal power per discrete sample is now $P/4W$.

Question 1: But what happens to white Gaussian noise sample ?

If I understand $N_0/2$ as noise power per Hertz, does doubling bandwidth to $2W$ doubles noise power per sample ? If yes, SNR should be decreased by factor 4 (2 of reducing signal power per sample and 2 of increasing noise power per sample). I doubt this.
If I use the projection of the same $N(t)$ to the new orthogonal basis set. Number of vector in the set has increased twice from $2W$ to $4W$. Does projecting the same signal to more basis vectors reduce the power of each projected coordinate ? If yes, is it true that the noise sample $N_k$ now have variance $N_0/4$ instead of $N_0/2$ ? And SNR does not change because $SNR = (P/4W) / (N_0/4) = P/N_0 W$.

Question 2: Is $N_0/2$ watts per Hertz ?

The noise power continues to be $N_0/2$, independent of bandwidth. The reason is that the noise variance at the output of a filter with frequency response $H(f)$ is $$\sigma_n^2=\int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,df.$$ (This is a direct application of the Wiener-Khinchin theorem). Since you're assuming an orhonormal basis, the noise power is always $N_0/2$, independent of the bandwidth.

Then, doubling your signaling rate without a corresponding increase in spent energy results in halving the SNR.

Regarding your second question, $N_0/2$ is the PSD of the white Gaussian process assumed to exist at the input of your measuring devices. The PSD is indeed measured in W/Hz.

Note that it is a bit weird to define $P$ as watts per second. One watt is one  joule spent per second; watts per second would be joules, or pure energy. Usually, you assume you have a source providing $P$ watts; this source operates for as long as necessary (even from $t=-\infty$ to $t=\infty$).


Determining Frequency and Period of a wave I'm gathering temperature data from a refrigerator.  The data looks like a wave.  I would like to determine the period and frequency of the wave (so that I can measure if modifications to the refrigerator have any effect).

I'm using R, and I think I need to use an FFT on the data, but I'm not sure where to go from there.  I'm very new to R and signal analysis, so any help would be greatly appreciated!

Here is the wave I'm producing:

Here is my R code so far:

require(graphics)
library(DBI)
library(RSQLite)

drv <- dbDriver("SQLite")
conn <- dbConnect(drv, dbname = "s.sqlite3")

query <- function(con, query) {
  rs <- dbSendQuery(con, query)
  data <- fetch(rs, n = -1)
  dbClearResult(rs)
  data
}

box <- query(conn, "
SELECT id,
       humidity / 10.0 as humidity,
       temp / 10.0 as temp,
       ambient_temp / 10.0 as ambient_temp,
       ambient_humidity / 10.0 as ambient_humidity,
       created_at
FROM measurements ORDER BY id DESC LIMIT 3600
")

box$x <- as.POSIXct(box$created_at, tz = "UTC")

box$x_n <- box$temp - mean(box$temp)
png(filename = "normalized.png", height = 750, width = 1000, bg = "white")
plot(box$x, box$x_n, type="l")

# Pad the de-meaned signal so the length is 10 * 3600
N_fft  <- 3600 * 10
padded <- c(box$x_n, seq(0, 0, length= (N_fft - length(box$x_n))))
X_f    <- fft(padded)
PSD    <- 10 * log10(abs(X_f) ** 2)

png(filename = "PSD.png", height = 750, width = 1000, bg = "white")
plot(PSD, type="line")

zoom <- PSD[1:300]

png(filename = "zoom.png", height = 750, width = 1000, bg = "white")
plot(zoom, type="l")

# Find the index with the highest point on the left half
index <- which(PSD == max(PSD[1:length(PSD) / 2]))

# Mark it in green on the zoomed in graph
abline(v = index, col="green")

f_s     <- 0.5 # sample rate in Hz
wave_hz <- index * (f_s / N_fft)
print(1 / (wave_hz * 60))

I've posted the R code along with the SQLite database here.

Here is a plot of the normalized graph (with the mean removed):

So far, so good.  Here is the spectral density plot:

Then we zoom in on the left side of the plot and mark the highest index (which is 70) with a green line:

Finally we calculate the frequency of the wave.  This wave is very slow, so we convert it to minutes per cycle, and print out that value which is 17.14286.

Here is my data in tab delimited format if anyone else wants to try.

Thanks for the help!  This problem was hard (for me) but I had a great time!

Interesting project you have going on there! :-)

From a signal analysis POV, this is actually a simple question - and yes, you are right that you would utilize the FFT for this frequency estimation problem.

I am not familiar with R, but what you essentially want to do is take the FFT of your temperature signal. Since your signal is real, you will get a complex result as the output of your FFT. You can now take the absolute magnitude squared of your complex FFT result, (since we dont care about phase). That is, $real^2 + imag^2$. This is your power spectral density.

Then, very simply, find the max of where your PSD sits. The abscissa of that max will correspond to your frequency.

Caveat Emptor, I am giving you a general outlook, and I suspect the result of the FFT in R will be normalized frequency, in which case you would have to know your sample rate, (which you do), in order to convert it back into Hz. There are many other important details that I am leaving out, such as your frequency resolution, FFT size, and the fact that you probably want to de-mean your signal first, but it will be good to see a plot first.

EDIT:

Let us take your signal into account. After I de-mean it, it looks like this:

You want to de-mean because the DC bias is a frequency of 0 Hz technically, and you do not want it to drown out other frequencies of interest. Let us call this de-meaned signal $x[n]$.

Now, when you take the FFT, you must take heed of some details. I use an FFT length of 10 times the length of the signal, so we can say that $N_{fft} = 10 * 3600 = 36000.$ Any FFT size set greater than the size of the original signal has the effect of interpolating your FFT result in the frequency domain. While this doesnt add any information for you from an information theory POV, it does help you better discern where exactly your true peak lies, especially in cases where it lies between two frequency bins. To get true better frequency resolution, you want to get more data. (i.e, let the sensor run for a longer period of time).

Moving along, from the text file, I see that your sensor is taking temperature readings every 2 seconds. This means that your sampling rate, or $f_s = 0.5 Hz$

Thus, let us take an FFT of $x[n]$, and call that result $X(f)$. This result will be complex, and in fact, it is conjugate symmetric, but thats not important here. (It just means for your purposes, half of it is redundant). Now if we take $10log_{10}(|X(f)|^{2})$, this is the power spectral density (in log scale). The result looks like so:

You can see how it is symmetric. If you ignore the last half, and just look at the first half and zoom in you, can see this:

So you have a peak at index 70! But what is index 70 in real frequency terms? Here is where you want your frequency resolution. To compute that, we simply take $\frac{F_s}{N_{fft}} = 1.3889e-005 Hz$. This simply means that every index represents a frequency that is a multiple of this number. Index 0 is $0 * 1.3889e-005 = 0 Hz$. Index 1 is $1 * 1.3889e-005 = 1.3889e-005 Hz$. Index 70 is $70 * 1.3889e-005 = 9.7222e-004 Hz$.

We are almost done. Now that you have this figure, we can convert it to something more palatable. This figure is just saying that your signal completes 9.7222e-004 cycles every second. So we can ask, how many minutes does it take to compute one cycle? Thus, $\frac{1}{(9.7222e-004)*60} = 17.14$ minutes per cycle.


Initial conditions for the LTI systems described as a difference equations Why do we need the initial conditions to be zero for the LTI systems described as a difference equations?

First question is why do we need it for linearity? I can't think of any example of the non linear system described as a difference equation with constant coefficients.
The definition of TI system is that it does not depend on particular time the input is applied. I can't understand how does this definition relate to the initial conditions? Those initial conditions will be the same whenever we apply the input so why they have to be zero?

Thanx

The problem is that non-zero initial conditions cause a term in the output signal that does not depend on the input signal. This explains why a system with non-zero initial conditions can neither be linear nor time-invariant. A linear system must have a zero output for zero input. With non-zero initial conditions the output will generally be non-zero, even for a zero input signal. Alternatively, think of scaling a given input signal. A linear system will have a response that is scaled in the same way. However, the part of the output signal caused by non-zero initial conditions will not scale accordingly, because it's independent of the input signal.

The same is true concerning time-invariance. For a time-invariant system, a shifted version of the input signal must result in an output signal with the same shift. However, the output term caused by non-zero initial conditions will not shift accordingly, as it is independent of the input signal.

Consequently, a system described by a linear difference equation with constant coefficients plus initial conditions is only linear and time-invariant if the initial conditions are zero.


Adaptive fiiler (LMS) step size bounds I am reading a book on adaptive filters.
For the optimal filter (Wiener), you can analytically obtain the step size for the gradient descent ($0<\mu<\frac{2}{\lambda_{max}}$).
This is not scary, because for the algorithm itself, the correlation matrix is also needed, at the same time it is possible to find the eigenvalues (although in practice it will probably take a long time).
When switching to adaptive filtering (the first considered LMS algorithm), averaging is abandoned (you may not know the statistical characteristics of the process, great!). Now the algorithm uses just instantaneous signal values. But at the same time, the requirements for stability in the book indicate the same. But where can we get them if we have abandoned the knowledge of the correlation matrix, how to correctly implement the algorithm?

You can analyse the input online, and get some rough estimates on its power, energy, pdf, autocorrelation matrix, etc. and also its eigenvalues of (time-varying) autocorrelation matrix. Then you may estimate the associated upper bound for the step-size of LMS.
Nevertheless, finding the optimal step-size in applications within unknown environments is not easy. So many of the LMS variations are based on different approaches for efficiently, and reliably finding the best step-sizes.


In process of image zooming, Sampling rate changes or stays constant? When we zoom an image, our sampling rate will be changed(increased) or sampling rate will not change because of zooming operation?

Yes, it changes the sampling rate. Just like with audio, 1 kHz tone is 1 kHz tone after resampling from 48kHz to 96kHz. So the image has same content, but with more pixels.


HHT and how to plot the Hilbert spectrum NOTE:
The following is based on an unanswered question posted here:
How to plot the Hilbert Spectrum in Hilbert-Huang transform?
, based on the usage of the "plot_hht" from here:
(https://www.mathworks.com/matlabcentral/fileexchange/19681-hilbert-huang-transform)

From what I understand the "plot_hht" does,
after applying EMD to a signal, the IMF's of the signal remain.
The Hilbert Transform is then applied to each IMF, and the resulting phase angles are somehow used to help find the instantaneous frequency (with respect to time).

Is it possible to represent these instantaneous frequencies vs time, with a third dimension, its strength/amplitude?

Thanks in advance.

EDIT P.S.: For clarity, the desired result is a graph that looks like this:
,
as taken from an article titled "Detecting position dependent tremor with the Empirical mode decomposition".

I figured out an answer.

From here:
https://web.cs.dal.ca/~tt/CSCI690611/papers/HHT.pdf

The Hilbert spectrum x(t) is:

where "a" is the instantaneous amplitude of the resulting IMF's
(i.e.: absolute value of the Hilbert Transform applied to every point of one IMF),
and "w" is the instantanous frequency of the resulting IMF', calculated from the inverse tan of the imaginary divided by real components.

So adding the results using all the IMF's should provide the 3D time-frequency-amplitude information.


impulse response cascaded with time reversed of itself Consider a filter with real-valued impulse response $h[n]$. The filter is cascaded with another filter whose impulse response is $h'[n] = h[-n]$, i.e. whose impulse response is the time-reversed version of $h[n]$.

$h[n]$ cascaded with $h[-n]$

i think we have to perform convolution

if so what is phase of new cascaded impulse response

HINT: See what happens in the frequency domain. Convolution corresponds to multiplication. You just have to figure out the Fourier transform of $h[-n]$, expressed in terms of the Fourier transform of $h[n]$. Multiplying the two functions results in the Fourier transform of the convolution of $h[n]$ with $h[-n]$. The result has a special property which makes it very easy to determine its phase.


DC bias from ADC in frequency detection scenarios In real-world applications where a frequency or multiple frequencies need to be detected from a digital signal that comes from an analog-to-digital converter (ADC), what are the typical 0 Hz (direct current, DC) bias levels in the digital signal in comparison to general noise levels and the levels of the frequencies of interest?
I'm asking the question because, if the answer is that DC bias typically has a significantly higher level than the general noise level, then it may be important to give special attention to it when designing frequency detection methods, for example by doing the frequency detection after a digital DC-blocking filter, or by ensuring that frequency detection is done using a transform such as discrete Fourier transform (DFT) that has basis vectors (aside from the 0 Hz basis vector) that are orthogonal to the DC bias, or by employing windowing to reduce spectral leakage.

You're spot on: whereas other frequencies are typically subject to noise that is somewhat benign shaped and result of random processes, DC is usually affected by things like a DC offset.
Physically, that happens pretty easily: Say, you've got an ADC with 14 bits of effective number of bits, and it has a differential sensing range of 0 V to 2 V. Then, but 1/8192 V is a single voltage bin, about 122 µV. Even if on the analog side, your DC level is as accurate (and it's often not), even the tiniest shift will be quite measurable in cumulative methods like a DFT spectrum due to the (usually fortunate) properties of dithering due to the naturally present voltage noise.
In fact, even the most basic opamp has an offset voltage, and introduces some static error, which might, to make matters worse, drift over time. You can buy "zero-drift" amplifiers that, in the end, tackle that by either forming a DC level control loop or frequency-modulating the DC component, but as you can imagine, making your linear amplifier an intentionally time-variant nonlinear system comes with its own bag of signal worms.
If you're building a direct conversion / quadrature RF receiver, then you'll also find the LO leakage at DC, and that usually really large. So, on the analog side, control loops trying to compensate DC are a common sight. These are, in fact, DC blocking filters, but as anything, they're imperfect (or you can't make them overly blocking without losing signal that you really wanted to digitize).
As you can imagine, that's a slight problem for radio systems: If you tune exactly to the carrier frequency, you stand to lose your carrier. Thus,

using a transform such as discrete Fourier transform (DFT) that has basis vectors (aside from the 0 Hz basis vector) that are orthogonal to the DC bias,

is exactly what many modern OFDM systems do: The quadrature mixer is the preferred mixer in highly integrated circuitry, mostly because it doesn't need the insane sampling rate of direct sampling, nor does it require the expensive and adjustable RF bandpass of subsampling, nor even the intermediate frequency filters that a superheterodyne receiver has, which usually can't be implemented in silicon technology and especially not on the same die.
So, these systems have to deal with DC in their digital baseband. And instead of doing so, they just ignore it: Wifi and 4G, for example, simply leave the center DFT bin unused while doing OFDM.


Is it theoretically possible to perfectly quantize a continuous signal? So, I'm completely new to digital signal processing, but while reading a piece this morning about quantization it got me daydreaming: could a machine ever be fast enough to sample the position and amplitude of each particle that makes up a wave?

Depends on assumptions you are willing to make and what type of signals are you trying to sample, but in theory I think that sampling rate equal to the Planck time would be a gold standard for anything...

This translates to sampling frequency of $1.855 \times 10 ^ {43} \mathtt{Hz}$ ($18.55$ tredecillion hertz). Personally I believe that machines will never be so fast. Obviously we can limit our assumptions to more realistic ones.

Obviously there is more constraints, such as: what sensor is capable of measuring displacement/velocity with such accuracy? Obviously this assumption is applicable even for subatomic scale, and most likely we don't need such accuracy in measurements of sound signals. Additionally when frequency is increasing, sound attenuation in media (air in this case) is increasing to enormous amounts and ultra-sonic sound waves are not propagating on large distances. At some point you will start measuring the Brownian motion. You also mentioned about possibility of measuring each particle separately and it is another, very difficult problem to tackle.


How to Resample Audio Using FFT or DFT I'm down sampling voice audio by first performing an FFT, then only taking the parts of the result that I need, and then performing an inverse FFT. However, it's only working properly when I'm using frequencies that are both power of two, say down-sampling from 32768 to 8192. I perform an FFT on the 32k data, discard the top 3/4 of the data and then perform an inverse FFT on the remaining 1/4.

However, whenever I try to do this with data that doesn't line up properly one of two things happen: The math library I'm (Aforge.Math) using throws a fit, because my samples are not a power of two. If I try to zero-pad the samples so they become power of twos, it get gibberish out on the other end. I also tried to use a DFT instead, but it ends up being insanely slow (this needs to be done in real time).

How would I go about to zero pad the FFT data properly, both on the initial FFT and the inverse FFT at the end? Assuming I have a sample at 44.1khz that needs to get to 16khz, I currently try something like this, the sample being 1000 in size.

Pad input data to 1024 at the end
Perform FFT
Read the first 512 items into an array (I only need the first 362, but need ^2)
Perform inverse FFT
Read the first 362 items into the audio play buffer

From this, i get garbage out at the end. Doing the same thing but without having to pad at step 1 and 3 due to the samples already being ^2, gives a correct result.

The first step is to verify that both your starting sample rate and your target sample rate are rational numbers.  Since they are integers they are automatically rational numbers.  If one of them wasn't a rational number it would still be possible to make the sample rate change, but it is a much different process and more difficult.

The next step is to factor the two sample rates.  The starting sample rate, in this case, is 44100, which factors to $2^2*3^2*5^2*7^2$.  The target sample rate, 16000, factors to $2^7*5^3$.  Thus, to convert from the starting sample rate to the target rate we must decimate by $3^2*7^2$ and interpolate by $2^5*5$.

The previous steps have to be done no matter how you want to resample the data.  Now let's talk about how to do it with FFT's.  The trick to resampling with FFT's is to pick FFT lengths that make everything work out nicely.  That means picking an FFT length that is a multiple of the decimation rate (441, in this case).  For the sake of the example, let's pick an FFT length of 441, though we could have picked 882, or 1323, or any other positive multiple of 441.

To understand how this works it helps to visualize it.  You start out with an audio signal that looks, in the frequency domain, something like the figure below.

When you're done with your processing you want to decrease the sample rate to 16 kHz, but you want as little distortion as possible.  In other words, you simply want to keep everything from the picture above from -8 kHz to +8 kHz and drop everything else.  That results in the picture below.

Please note that the sample rates are not to scale, they are just there to illustrate the concepts.

The beauty of picking an FFT length that is a multiple of the decimation factor is that you can resample simply by dropping portions of the FFT result, and then inverse FFT what is left.  In the case of our example you FFT 441 samples of data, which gets you 441 complex samples in the frequency domain.  We want to decimate by 441 and interpolate by 160 ($2^5*5$), so we keep the 160 samples that represent the frequencies from -8 kHz to +8 kHz.  We then inverse FFT those samples and presto!  You have 160 time domain samples that are sampled at 16 kHz.

As you might suspect, there are a couple of potential problems.  I will go through each one and explain how you can overcome them.

What do you do if your data is not a nice multiple of the decimation factor?  You can easily overcome this by padding the end of your data with enough zeros to make it a multiple of the decimation factor.  The data is padded BEFORE it is FFT'ed.
Even though the method I explained is very simple, it is also non-ideal in that it can introduce ringing and other nasty artifacts in the time domain.  You can avoid that by filtering the frequency domain data before dropping the high frequency data.  You do this by FFT'ing your filter of length $l$, padding your data (before FFT'ing it) with at least $l-1$ zeros (please note that the number of data samples and the number of padding samples must BOTH be a positive multiple of the decimation factor- you can increase the padding length to meet this constraint), FFT'ing the padded data, multiplying the frequency domain data and filter, and then aliasing the high frequency (> 8 kHz) results down into the low frequency (< 8 kHz) results before dropping the high frequency results.  Unfortunately, since filtering in the frequency domain is a big topic in its own right, I won't be able to go into more detail in this answer.  I will say, though, that if you filter and are processing the data in more than one chunk, you will need to implement Overlap-and-Add or Overlap-and-Save to make the filtering continuous.

I hope this helps.

EDIT:  The difference between the starting number of frequency domain samples and the target number of frequency domain samples needs to be even so that you can remove the same number of samples from the positive side of the results as the negative side of the results.  In the case of our example, the starting number of samples was the decimation rate, or 441, and the target number of samples was the interpolation rate, or 160.  The difference is 279, which is not even.  The solution is to double the FFT length to 882, which causes the target number of samples to also double to 320.  Now the difference is even, and you can drop the appropriate frequency domain samples with no problems.


Whats the technical term for a filter that allows a non unity stride when sliding it's window? So a 1D mean filter of size two with a sliding window defined as W below could be considered as having a stride of one.

Wn=(Xn + Xn+1) / 2
Vn=W2n

where as if I were to define a 1D 'stepped mean filter' with a stride of two and labeled U below.

Un= (X2n + X2n+1) / 2

where the number of samples in U (and also V) would be half that of W.
My questions are the following

Is U still a filter (I'm guessing it is even though not fitting with the traditional kernel convolution implementation I associate with filters)
What is the technical language one could use to communicate and research questions regarding filter/non-filter such as just given?
All though not directly related but since I have the floor is, if  V and U would produce equivalent results (my intuition says no)?

That's a decimating filter. As such, it serves as a 1/2-rate resampler.


Is it possible to generate depth map from graycode sequence without using printed calibration patterns? In a camera-projector system, if only the camera's intrinsics are known, is it possible to generate the depthmap fully automatically through a graycode sequence without any manual calibration process (like using a printed calibration pattern then move it away etc.) ?

... if the camera's intrinsics are known and calibrated, is it possible to generate the depthmap fully automatically through a graycode sequence without any manual calibration process ?

Yes, provided that you are talking about a stereo-vision system.

Simply extracting the depth information, does not require knowledge of the camera model. It only requires knowledge of the pixel positions that refer to the same point in space and this is what the Gray code pattern does.

Every Gray code pattern that is projected on to the image "tags" a certain point in space. It is also useful to think about Gray code patterns as bit planes. During the "encoding" process, you segment the image pixels and assign them to bit planes. In this way, every point in space is addressed by a unique binary number that is formed by stacking the bit planes.

If you repeat this process for the left and right images, you will obtain two pairs of $x,y$ coordinates that will tell you where a particular point in space appears to be from the viewpoint of each camera. This is enough to work out "depth" by triangulation.

Knowing the intrinsic parameters allows you to relate this "pixel" depth map to physical coordinates.

For more information, please see this (or any other relevant introduction to imaging using structured light) and this.

Hope this helps.

EDIT:

The basic formula for recovering depth is $Z = f \frac{D}{d_p}$ where:

$f$ is the focal distance of the camera lenses.
$D$ is the distance between the cameras providing the stereo vision
$d_p$ is the "disparity" or relative difference of the position of two pixels associated with the same region in the two images.
$Z$ is the distance from the baseline plane which in turn depends on the relative placement of the two cameras.

Now, if the cameras are not fully calibrated, it is impossible to relate $Z$ with reality as the relative distances would not be preserved. BUT, in an extreme case, you can set $D=1$ and obtain a depth map that gives you some information about the depth of the scene. This should "work" for reasonable assumptions of rectilinear lens and aligned cameras and probably return enough information to discriminate between a ball and a cube, if that is the interest of the application here (?)


How to Cluster Image Colors Using K Means Clustering in CIE a*b Domain? I am working on building detection from aerial images using Image Processing in MATLAB.

I am doing an image segmentation based on the local homogeneity of image. I want to apply K Means Clustering in CIE LAB Color Space with some clusters on the anistropic diffused image for color quantization and generate the image for further processing.
My RGB image is of 877 X 454 dimension.

How do I do that?

Thanks

You should use makecform to create a color transfer structure.
You should use srgb2lab to create a structure which will convert you to LAB.
To apply the color transformation use applycform.
Then, as a result you'll have an image in the LAB domain.

Then use kmeans to cluster the image colors on the AB channels.

Good Luck...


Beamforming back lobes question When performing a conventional digital beamforming I see in some textbooks they show beam patterns as

where the arrival angle ranges between -90 to 90 degrees. But then other places show the full 360 degree view sometimes in a polar plot like this

My question is, when doing a simulation is it ever warranted to ignore signals with angles of arrival greater than 90 degrees? For example, if my signal of interest has an angle of arrival of 30 degrees and an interference signal is at -80 degrees (difference between these is |-80 - 30| = 110 degrees), is it really like the first figure where the gain plot only exists between -90 to 90 degrees?

EDIT: My simulation setup.

the plots are likely from a line or what is often called a ULA and there is a mirror image response.  It is sometimes called a left/right ambiguity.

Your question’s answer depends on the purpose of your simulation.  The plots are probably equivalent, the response is over 360 degrees in both plots but people often drop the redundant half circle.

if you have a problem where objects are distributed over 360 degrees, simulate that, otherwise 180 is fine.

I usually plot in $\cos(\theta)$ coordinates which range from -1 to 1 for most of my work. For a ULA, the pattern is the same regardless of the steering direction using this convention. If the array weights are symmetric, 0 to 1 is sufficient.

It all depends on what your application is and what you are investigating and the conventions of your community which isn’t specified in your question.

The only caveat is that in some applications you might have some structure blocking the back half of the array. In SONAR it is called a baffle. In RF there are other caveats like a ground plane on one side.

Actual arrays have more complicated responses because what you showed are probably Array Factor responses where each element is assumed ideally “isotropic”.  Real array elements have some directionality. The Array Factor also ignores element to element coupling.

In summary, it could be warranted but in many applications it might not matter.


Audio steganography using phase encoding technique I referred the below papers to understand the working of Phase Encoding,but the explanation is not proper.Just a few steps have been given in brief.

ijcsi.org

www.ijarcce.com

Can someone please explain the concept of how Audio Steganography using Phase Encoding works?

Here is the algorithm:

Divide an input signal $S$ into $n$ consecutive blocks $S_i$ ($i\in\{1,\ldots,n\}$) of length $l$ each
Compute magnitude $A_i$ and phase $\phi_i$ of each block using FFT where A(i)=abs(fft(S{i})) and Phi(i)=angle(fft(S{i})) in MATLAB.

For example, here are the first 3 elements of $\phi$:

(Recall that FFT of real signal has odd symmetry which is clearly visible on plots)
Compute phase differences $\Delta\phi_i=\phi_i-\phi_{i-1}$ for $i\in\{2,\ldots,n\}$
To encode your binary data $d$ of length $m:m<l$ assign
$$\phi_{data}[i]=
\begin{cases}
\pi/2,& d_i=0\\
-\pi/2,& d_i=1
\end{cases}$$ for $i\in\{1,\ldots,m\}$.
Replace following elements of first phase sequence, $ {\phi_1} $,  with ${\phi _{data}}[i]$ for $i\in\{1,\ldots,m\}$ as:
$${\phi '_1}[ {L/2 - m + i}] = {\phi _{data}}[i]\,$$
To maintain odd symmetry property of DFT, repeat same progress as:
$${\phi '_1}[{L/2 + 1 + i}] =  - {\phi _{data}}[{m + 1 - i}]$$
To maintain phase differences sequentially reassign $\phi '_i=\phi_{i-1}+\Delta\phi_i$ for $i\in\{2,\ldots,n\}$
Here are the first 3 elements of $\phi '$ after this step:

(Note how encoded $d$ has affected next blocks)
Reconstruct the signal using inverse FFT applied to each block $A_i exp(j \phi '_i)$ (where $j$ is imaginary unit) and joining all blocks together.

Explanations to many details like the length of blocks and why particular steps are taken are given in the documents you and endolith linked to.

A minimal IPython Notebook version: github.com/danylo-dubinin/secret_in_wav.


How to find the envelope of a wide-band signal? I would like some feedback on possible techniques that one may use to determine the envelope of a broad-band time domain signal. I have heard anecdotally, that it is not as straight-forward as it seems, but I am not sure. (Some previous discussions about similar subjects exist here and here but I would like an explicit discussion on broad band signals here).

For narrow-band signals, one may compute the envelope via the absolute value of the analytical signal as such:

Let $s(t) = A cos(2\pi ft)$ be the real signal. Then, the analytical signal is:

$$
a(t) = Ae^{j2\pi ft } = Acos(2\pi ft) + jH(s(t)) = Acos(2\pi ft) + jAsin(2\pi ft)
$$

where $H(.)$ denotes the hilbert transform of the real signal. Then the envelope is given by $E = |a(t)|$.

What we are really doing is converting a real signal into phasor form, and taking that phasors absolute magnitude. In other words, since any real signal can be described as two phasors of equal magnitude rotating at both negative and positive angular frequencies, we simply discard of those phasors and work with the remaining one, and take its absolute magnitude.

Now to my understanding, this method works well for narrow-band signals, and not for broad-band signals. Thus my questions are as follows:

EDIT: I have edited and narrowed down my questions based on the feedback I have received:

Question: I am wondering, what ways one may use to compute the energy-envelope of the following broad-band signals, assuming the above analytical approach does not work well. (If it doesnt, why doesnt it?)

i) A chirp signal, at some carrier. (linear, exponential, etc)

ii) A very short duration windowed sinusoid, of length $T_p$. (So bandwidth is $2/T_p$)

iii) An OFDM signal at some carrier. OFDM is composed of multiple orthogonal (I do not care about it as much, but put it for completion).

iv) @Andrey's 'mixture of tones'. (I am not clear on this, but I am guessing it means simple summation of multiple sinusoids around a carrier?)

I cannot think of any other types of broad-band signals, but answering those cases might go far in grounding my understanding here.

Thanks.

The analytical signal approach works very well for chirp signals. It stops working for mixtures of different tones. Taking the absolute value is a nonlinear operation that introduces cross terms and thus distorts you envelope.

I'd say that the envelope of a broad band signal as a property that is not well defined. First you need to define what you actually mean by that. Than you can come up with a method. For instance, the phase of a signal is not well defined for broad-band signals as well, because there exist two or more phases simultaneously. But depending on the actual problem you can come up with a special solution.

Just some notes to the phase example, because I'm an expert there. In physics there is the problem of defining a phase to measure phase synchronization of two or more broad-band systems. One solution is to compute the Lyapunov spectrum which measure phase synchronization implicitly without ever defining a concrete phase.

I just write that to make you think about what you actually want to know/measure/manipulate in order to step a bit back and have a more general look at your problem.


Simulation of IQ imbalance using MATLAB I already have the value of Image Rejection Ratio (IRR) in dB for IQ imbalanced modulator and I want to simulate this effect using MatLab. How can I obtain $\alpha$ and $\theta$ to calculate the imbalanced gain for I & Q branches ??

Bottom line:

You cannot uniquely determine $\alpha$ and $\theta$ from IRR given in dB since it has a dependence on both given by the following approximation:

$$IRR \approx 20Log_{10}\bigg(\bigg|\frac{1-  \alpha - j\theta}{1+  \alpha + j\theta}\bigg|\bigg)$$

Where

$IRR$ is the image rejection ratio in dB

$\alpha$ is amplitude imbalance as a magnitude ratio (typically 0.9 to 1.1)

$\theta$ is phase imbalance in radians (typically -0.1 to 0.1)

This uses the small angle approximation such that $sin(\theta) \approx \theta$ which covers most IRR applications.

For cases where there is no quadature imbalance and only amplitude imbalance the following relationship is exact and you can then solve for $\alpha$ from IRR:

$$IRR = 20Log_{10}\bigg(\bigg|\frac{1-  \alpha}{1+  \alpha}\bigg|\bigg)$$

From which:

$$\alpha = 1\pm \bigg(1-\frac{1-r}{1+r}\bigg)$$

Where
$r = 10^{IRR/20}$

With IRR given as a negative quantity indicating the image rejection ratio in dB.

For example the amplitude imbalance to provide -30  dB image rejection would be as follows assuming perfect quadrature:

$r = 10^{-30/20} = 0.0316$

$\alpha = 0.9387, 1.0613$

(Or +.51 dB, -.54 dB)

A similar approximation can be given for small angle $\theta$ for when there is no amplitude imbalance. Bessel functions can be used to derive an exact result for $\theta$ applicable to larger angles.

Details

If you complex waveform is already generated (as it would be trivial to generate otherwise) you can use the corrections as given in the diagram below (if it corrects a phase and amplitude error it will equivalently induce a phase and amplitude error):

So multiply I by your $\theta$ term and add to Q to add quadrature error ($\theta$ will be approximately your angle in radians using small angle approximations), and scale I by your $\alpha$ term to add amplitude error.

$$I_{out} = \alpha I_{in}$$
$$Q_{out} = \theta I_{in} + Q_{in}$$

Both $\alpha$ and $\theta$ will both effect IRR. If the IRR is given in dB, there is no way to uniquely determine $\alpha$ and $\theta$. Note the following to see how IRR can be determined from knowing $\alpha$ and $\theta$:

Graphically it is easy to derive this since a circle represents perfect I and Q balance and results in image rejection. We can see this simply by using Euler's Identities:

Given Euler's Identities as

$$cos(\omega t) = \frac{e^{j\omega t}}{2} + \frac{e^{-j\omega t}}{2}$$
$$sin(\omega t) = \frac{e^{j\omega t}}{2j} - \frac{e^{-j\omega t}}{2j}$$

And a circle with perfect IQ balance as:
$$y(\theta) = cos(\omega t) + jsin(\omega t)$$

$$= \frac{e^{j\omega t}}{2} + \frac{e^{-j\omega t}}{2} +  j\frac{e^{j\omega t}}{2j} - j\frac{e^{-j\omega t}}{2j}$$

$$= \frac{e^{j\omega t}}{2} + \frac{e^{-j\omega t}}{2} +  \frac{e^{j\omega t}}{2} - \frac{e^{-j\omega t}}{2} = e^{j\omega t}$$

Resulting in the well known Euler's relationship as

$$ e^{j\omega t} = cos(\omega t) + jsin(\omega t)$$

The above is insightful as it demonstrates the sideband rejection through phase summation from the I and Q components being in perfect balance and quadrature, and recognizing that $e^{j\omega t}$ represents a single sideband tone (a positive frequency $\omega$ while a sinusoid has two sidebands or a positive and negative frequency).

So similarly we can evaluate the effect of amplitude and phase imbalance. From the block diagram above showing the relationship of $\alpha$ and $\theta$ for small angle approximation when $sin(\theta) \approx \theta$:

$$y(\theta) = \alpha cos(\omega t) + jsin(\omega t) + j\theta cos(\omega t)$$

$$y(\theta) = \frac{ \alpha e^{j\omega t}}{2} + \frac{  \alpha e^{-j\omega t}}{2} + \frac{je^{j\omega t}}{2j} - \frac{je^{-j\omega t}}{2j} + \frac{j \theta e^{j\omega t}}{2} + \frac{ j \theta e^{-j\omega t}}{2} $$

$$y(\theta) = \frac{ \alpha e^{j\omega t}}{2} + \frac{  \alpha e^{-j\omega t}}{2} + \frac{e^{j\omega t}}{2} - \frac{e^{-j\omega t}}{2} + \frac{j \theta e^{j\omega t}}{2} + \frac{ j \theta e^{-j\omega t}}{2} $$

$$y(\theta) = \frac{(1+ \alpha +j\theta)}{2} e^{j\omega t} - \frac{(1-  \alpha - j\theta)}{2} e^{-j\omega t} $$

The IRR in dB would be $20Log_{10}$ of the ratio of the two exponential coefficients given as:

$$IRR \approx 20Log_{10}\bigg(\bigg|\frac{1-  \alpha - j\theta}{1+  \alpha + j\theta}\bigg|\bigg)$$


How to represent half of equilateral triangle signal? How to represent the signal in above image ?

I know a ramp signal can be represented by

$$x(t) = t\cdot u(t)$$

but here amplitude of that signal is not same with $t$

Try:
$$
x(t) = \left[ u(t-T_{\rm on}) - u(t-T_{\rm off}) \right] \cdot \frac{T_{\rm off} - t}{T_{\rm off} - T_{\rm on} }
$$

So, to explain:

$T_{\rm on}$ is the time at which the ramp turns on.
$T_{\rm off}$ is the time at which the ramp turns off.

The piece:
$$
 \left[ u(t-T_{\rm on}) - u(t-T_{\rm off}) \right]
$$
will be zero before $T_{\rm on}$ and one afterwards, until $T_{\rm off}$ at which time it will be zero.  Note that this assumes that $T_{\rm off} > T_{\rm on}$.

That just looks like a pulse of duration $T_{\rm off} - T_{\rm on}$ with amplitude 1.

As you said in the question, the function $-t$ has the right sign of the gradient, but (possibly) the wrong slope and (possibly) the wrong intercept. For the numbers given, $T_{\rm off} = 5$ and $T_{\rm on} = 4$, the gradient is correct but the intercept is wrong.

So, the rest of the equation just fixes that:

We need our function to be one at $T_{\rm on}$.
We need our function to be zero at $T_{\rm off}$.

So

$$
1 = m T_{\rm on} + c  \tag{1}
$$
$$
0 = m T_{\rm off} + c \tag{2}
$$

So that
$$
(1) - (2) \rightarrow 1 = m \left( T_{\rm on} - T_{\rm off} \right) \rightarrow m = \frac{1}{T_{\rm on} - T_{\rm off} } \\
$$
and substituting that into (2) gives
$$
0 = \frac{1}{T_{\rm on} - T_{\rm off} } T_{\rm off} + c \rightarrow c = \frac{ - T_{\rm off}}{T_{\rm on} - T_{\rm off} }
$$


How to create an echo effect in audio samples? I am feeding the samples of pure sin wave to the audio DAC, I would like to create an echo effect. What I am not sure is how to go about it.
I was thinking about superimposing the the current buffer that I am feeding with the delayed version of it with some decaying factor but that does not seem to help at all. I am using FIFO buffer to fill up with the next sample.
Is there any standard approach in creating a good echo effect? How long the delay should be(how many samples)?

As hotpaw2 said, you are not going to hear an echo or delay effect unless your signal changes, so you are not hearing it because you are sending in pure, unchanging signal. You might want to try sending a burst of sine waves every second or better yet a chirp or an actual recording of music or someone speaking. It will work best if you have periods of silence in your input.

Here is an explanation of how to create a delay or echo effect in c-like pseudo-code.


Different amplitude for different time resolution with FFT In the following Matlab code, the amplitude of the spectral domain signal is printed
Temporal_Pulse_Width=1;
t=linspace(0,20,2000);
dt=t(2)-t(1);
Time_Period=round(Temporal_Pulse_Width/dt);
AA=zeros(length(t),1);
AA(Time_Period/2:3*Time_Period/2)=1;
A=20*log10(abs(fftshift(fft(AA/length(t)))));
frequencyAxis = ((0:length(t)-1) -ceil((length(t)-1)/2))/length(t)/(t(2)-t(1));
plot(frequencyAxis,A,'LineWidth',3);
max(A)

If we change the vector t as follows,
t=linspace(0,10,2000);

The result (max(A)) will also change. How can I normalize the fft without having this problem?

As you can see in MATLAB documentation for FFT
If you have Y = fft(x) then
$$Y_k = \sum_{j=1}^{N} x_j \cdot exp\left(-\frac{2 \pi i}{n} (k - 1)  (j - 1)\right)$$
$$x_j = \frac{1}{n} \sum_{k=1}^{N} X_k \cdot exp\left(+\frac{2 \pi i}{n} (k - 1)  (j - 1)\right)$$
Notice that the FFT is just a sum and the IFFT is divided by $n$, with this definition you know that  fft(ifft(x)) = x, but as the number of terms increase $Y_k$ will increase.
Your signal is a fixed duration pulse and you want a fixed amplitude in the frequency domain, so it is like an aperiodic Fourier transform
$$ \int x(t) exp\left( \frac{2\pi i}{n} u t \right) dt $$
A zero order approximation for this can be obtained by computing
A=20*log10(abs(fftshift(fft(AA*dt))));

For the pulse of width 1, max(A) will be 0dB this way.


Simulating room impulse response between rooms I want to simulate a sound source in one room and a microphone in a second, adjacent room. I am familiar with the image method and Habet's very efficient implementation for room impulse response (RIR). This implementation will be problematic when the source is in another, second room. However, the limitations of the image method for rectangular rooms and frequency constant reflection coefficients are good enough for me as it is a simulation.
AS I wish to simulate a large number of examples (~1000) so measuring an RIR is not an option and an inaccurate, fast simulation method is good enough.
I also found this paper but I could find no implementation available. This implies that it is possible to simulate though I am not sure.
I am looking for an implementation of an RIR generator for coupled rooms or another quick solution. I wish to spend as little time as possible on learning the inside of the generation method, as it is not my goal, so another tool\method with a good tutorial and fast ramp-up will be also very helpful.
Can anyone recommend the shortest way to achieve this?

This is an FDTD based, wave equation solution. It is the best free solution I have found. None of the other answers has helped so far.


Entropy : do we prefer higher or lower entropy? The entropy of a source $H(S)$  gives the average codeword length to encode a given source alphabet. i.e. it is the average number of bits per symbol required to encode the information in the source. Let, $S_1 = 101101...$ be a message of length $N$. The entropy of the system that produces the message depends on the number of unique symbols, $k$ and the formula is $H(S) = ln k$ where $k$ is the number of unique symbols. The entropy $H(S_1) = \log(2) = 0.69$ bits/symbol. For another different message, $S_2$ having $k=8$ unique symbols, $H(S_2) = \log(8) = 2.07$ bits/symbol. I have 2 questions

Is there a way to determine $n$, the minimum word length for a message sequence of $N$? For example, $S_1 = (10)(11)(01)$ here $n$ = 2.
It is desired that the entropy of the system be maximized, so $H(S_1) < H(S_2)$ so should I prefer $k=8$ symbols then? Which entropy do we prefer - message with higher entropy and hence is this the reason why we prefer higher modulation?

You seem to have a number of misunderstandings, which I'll try to clarify while also trying to help with your questions.

  The entropy of a source $H(S)$ gives the average codeword length to encode a given source alphabet. i.e. it is the average number of bits per symbol required to encode the information in the source.

While this is true, I think it's not the more enlightening way to look at the entropy. Think of the entropy as an inherent property of a source. The entropy measures the average information of each message produced by the source.

As an analogy, think of the source as a black box with a button; every time you press the button, the source tells you something. A high-entropy source will consistently (not necessarily always) produce messages that leave you saying "Wow! That's amazing". A low-entropy source will be boring most of the time.

Now, Shannon proved that it's possible to find an enconding from source messages to bits, so that the average number of bits required per message is equal to the entropy. This is why your statement above is technically correct. Note that, typical of Shannon, he didn't actually tell us how to design such an encoder (but he provided a few hints).

  Let $S_1=101101...$ be a message of length $N$ [...] The entropy $H(S_1)=\log(2)=0.69$ bits/symbol.

Now you can see this statement makes no sense. $S_1$ is a message, so it cannot have an entropy. What it does have is information. To find the entropy of the corresponding source, you need to look at its message set (alphabet) and its probability distribution. If the source has memory, you also have to take that into account.

  Is there a way to determine $n$, the minimum word length for a message sequence of $N$? For example, $S_1=(10)(11)(01)$ here $n = 2$.

I'm not sure what you're asking here, but I'll make a guess. One of the hints Shannon gave us about source encoding is that efficient encoding requires us to gather many source messages and encode them together. This is equivalent to defining a new source with a different alphabet (whose elements are concatenations of the original symbols) and a different probability distribution.

You seem to be asking about moving in the other direction: breaking up the source alphabet into simpler symbols. In the case of a source whose messages are strings of 0s and 1s, you can always define it as a source with alphabet $\lbrace 0,1 \rbrace$, but you have to be careful when doing so, as you may need to introduce memory to the source.

For example, say you have a source with alphabet $\lbrace 00, 01, 11 \rbrace$. To redefine it a source with alphabet $\lbrace 0,1 \rbrace$, you'll need to introduce a mechanism to avoid messages such as $0,0,1,0$, which the original source cannot produce.

  It is desired that the entropy of the system be maximized, so $H(S1)<H(S2)$, so should I prefer $k=8$ symbols then? Which entropy do we prefer - message with higher entropy and hence is this the reason why we prefer higher modulation?

This question does not make much sense to me, since a given source has a given entropy (remember that entropy is an inherent property of a source). So, what is there to maximize?

Now, if you're given a choice between two sources, one of which has higher entropy, what you should do is not clear either. What do you want to achieve?

From the point of view of a source coder (compressor), which is trying to compress the messages as much as possible, having a low-entropy source is preferable. For example, a PNG compressor which is fed single-color images will be able to compress them a lot.

From the point of view of the communications engineer, what you want to do is to match the information content of the source to the capacity of the channel.

Maybe what you mean is this: remember that Shannon says that we need to group source messages together to achieve larger encoding efficiency. This process will also result in a source of larger entropy -- but this means more bits per message, so how can this be better?!

To clarify this point, consider a binary source with alphabet $\lbrace 0,1 \rbrace$ and probability distribution $\lbrace 0.9,0.1 \rbrace$. Even though this source is low-entropy (most of its messages are boring), it's very hard to compress: you still need one bit per message!

Let's say that you redefine this source to have 1024 messages of 10 bits each, and calculate the new probability distribution. You'll find that a few of the new 10-bit messages are much more likely than the others. You can then use (for instance) the Huffman algorithm to find an encoding that approaches the entropy.

So, from this perspective, redefining your source to have larger entropy is better, because you'll be able to find better encoders.

About higher-order modulation -- that doesn't play a role in this problem at all.


Is there a reason why with symmetric padding, the inverse wavelet transform is not the adjoint of the wavelet transform? I recently stumbled upon a bothering fact when using the pywavelet library in Python. When we use the default "symmetric" padding, the inverse wavelet transform is not the adjoint of the wavelet transform (and it's not the mathematical inverse of the wavelet transform, "only in one direction").

However, if you use the "zero" padding the adjoint property is verified (not the inverse property). And if you use the "periodization" padding then all the properties are verified.

I wanted to know if there was a theoretical reason behind why the padding has such a big influence on these properties or if it was down to implementation.

You can find a code below highlighting the said issue:

import numpy as np
import pywt

wavelet_type = "db4"
level = 4
mode = "symmetric"

coeffs_tpl = pywt.wavedecn(data=np.zeros((512, 512)), wavelet=wavelet_type, mode=mode, level=level)
coeffs_1d, coeff_slices, coeff_shapes = pywt.ravel_coeffs(coeffs_tpl)
coeffs_tpl_rec = pywt.unravel_coeffs(coeffs_1d, coeff_slices, coeff_shapes)
_ = pywt.waverecn(coeffs_tpl_rec, wavelet=wavelet_type, mode=mode)

def py_W(im):
    alpha = pywt.wavedecn(data=im, wavelet=wavelet_type, mode=mode, level=level)
    alpha, _, _ = pywt.ravel_coeffs(alpha)
    return alpha

def py_Ws(alpha):
    coeffs = pywt.unravel_coeffs(alpha, coeff_slices, coeff_shapes)
    im = pywt.waverecn(coeffs, wavelet=wavelet_type, mode=mode)
    return im

x_example = np.random.rand(*coeffs_1d.shape)
y_example = np.random.rand(512, 512)
print("Adjoint:")
x_Tadj_y = np.dot(x_example, np.conjugate(py_W(y_example)))
T_x_y = np.dot(py_Ws(x_example).flatten(), np.conjugate(y_example.flatten()))
print(np.allclose(x_Tadj_y, T_x_y))
print("\n Inverse from image to image:")
print(np.allclose(py_Ws(py_W(y_example)), y_example))
print("\n Inverse from coefficients to coefficients:")
print(np.allclose(py_W(py_Ws(x_example)), x_example))

EDIT:
When reading the documentation about the different paddings, I could see that the problematic paddings compute a redundant representation of the original image (for example you have more coefficients than pixels in your original image). Therefore, it's normal that the inverse property is not verified from coefficients to coefficients (if taken at random without structure). However, this doesn't explain why the adjoint property is not verified.

Also I don't know whether this redundance is theoretical to ensure perfect reconstruction or is an implementation choice.

Actually the redundancy perfectly explains why the adjoint property is not verified. The reconstruction operator is not orthogonal anymore and therefore the inverse is not necessarily its adjoint.

You can compute the adjoint efficiently according to https://arxiv.org/pdf/1707.02018.pdf.

Regarding the necessity of redundance, according to this person it is theoretically needed, although he couldn't provided a resource or a proof of that.
But since it's a separate issue, I am going to ask another question (you can check it here).


Why are males' pitch lower than females? Hi :) I'm not sure this is the right place to ask this question, but I'll try anyway:
Here are two well-known facts, and I don't understand how they work together:

Males usually have deeper voices than females because they have longer vocal folds.

Short and relaxed vocal folds produce a deep voice, long and tense vocal folds produce a high-pitched voice.

So... given that males have longer vocal folds, isn't that supposed to mean that they produced high-pitched voices? What did I get wrong?
Thanks

Males usually have deeper voices than females because they have longer vocal folds.

Yes

Short and relaxed vocal folds produce a deep voice, long and tense vocal folds produce a high-pitched voice.

Yes on "relaxed", no on "short".
Treating the vocal cords as strings: The resonance frequency is a function of mass (or density), tension and length.
Higher mass -> lower pitch
Higher tension -> higher pitch
Higher length -> lower pitch

Male vocal cords are longer and heavier.


Check for time invariance of a _system_ using MATLAB? I am having difficulty in developing an algorithm for checking if a given system is time invariant or not.

For example, let there be 2 systems with input output relation given as:

$$y(t)=x(t-2)$$ or $$ y(t)=tx(t-1)$$

Now the first system is time invariant whereas the second one is not. But how do I write a function in MATLAB that would take the input output relationship and find out if the system is time invariant or not.

You can't do that. What you can do is write a function that computes the output given the input signal, then shifts the input signal (the shift could be an input parameter), and then checks if the output signal is a shifted version of the original output signal. If this is NOT the case, then you know that the system is time varying. However, if the output signal is just a shifted version of the original output signal then you don't know anything, because you would need to try all possible (i.e. infinitely many) shifts (and all possible input signals) and verify this condition for all these cases. Only after having checked infinitely many shifts (and all possible input signals) would you be able to assert that the system is indeed time-invariant.


How to compute the power bands of an eeg signal using python? so I have an eeg signal (edf format) that has 25 channels and 248832 entries, sampling frequency of 512Hz. I have to compute the frequency bands:
– Delta: 0.1-4Hz
– Theta: 4-8Hz
– Alpha: 8-12Hz
– Sigma: 12-16Hz
– Beta: 16-36Hz
– Gamma: >36Hz
and plot them accordingly. I am using Python for this with scipy, numpy, etc. and I should get to something like this:

Does anyone have any point-outs/ideeas/tutorials that could help me compute the bands and then get such a plot(probably a histogram)? Thanks!

Here is some code that may solve your problem:

from scipy.io import loadmat
import scipy
import numpy as np
from pylab import *
import matplotlib.pyplot as plt

eeg = loadmat("mydata.mat");
eeg1=eeg['eeg1'][0]
fs = eeg['fs'][0][0]
fft1 = scipy.fft(eeg1)
f = np.linspace (0,fs,len(eeg1), endpoint=False)
plt.figure(1)
plt.plot (f, abs (fft1))
plt.title ('Magnitude spectrum of the signal')
plt.xlabel ('Frequency (Hz)')
show()

You can also check this other link:

http://forrestbao.blogspot.pt/2009/10/eeg-signal-processing-in-python-and.html


Difference between Digital signal and Discrete signal I have just started studying Digital Signal Processing. Can someone explain what is the difference between a Discrete Signal and a Digital Signal in simple words?

Thanks in advance !

Suppose you have a continuous time analog signal. It is continuous in both time and amplitude. Now when you sample it ,you get discrete samples every Ts seconds. Now you have discrete samples(discrete in time) each of which can take continuous value( in amplitude). This is normally referred to as discrete signal(discrete in time but continuous in amplitude).

Now further when you take this discrete signal and quantise it i.e. say assign each sample which takes continuous amplitude to one of the N discrete quantisation levels of a quantiser, Then your total signal now is a digital signal. So a digital signal is discrete in time and discrete in amplitude.


Sound glitches using fft segmented convolution I’m trying to build a binaural audio player, I'm using MIT Kemar c library to get the hrtf and a fft segmented convolution. It works on real time but the problem is each time I change the hrtf I hear noises like small cracks or something like vinyl noise. But I'm not sure why. Please help.
If you are interested in the project here is my github repo:
https://github.com/Bk8/Binaural

One possible way to reduce cracking noises between convolutional block filtering parameter changes is to use an overlapping of the blocks and do a cross fading (linear or raised-cosine) of the overlapping segments.


Removing background noise algorithm What is the algorithm for removing background noise from an audio clip of human speech? Say I have a 10 second audio clip of somebody talking, and there is background noise of them tapping on their keyboard. Now, I would like to remove the "keyboard tapping" signal while preserving the speech signal (always assuming a single speaker).

What are the necessary steps to achieve this? Is the best approach to use voice activity detection to try and isolate a representative signal of keyboard tapping, and then "remove" that pattern from the rest of the signal? What is required to "remove" the specific sound from the signal once identified?

In order to remove the noise from the signal, you need to subtract the spectral content of the noise from that of the signal ( filtered_signal = unfiltered_signal - noise ). This technique is called spectral subtraction. This will only work well if the noise is rather uniform across the entire length of the signal.

The steps required to do this are

Identify the noise region(s) in your signal and do an FFT on it. A single noise region is enough provided that it's long enough for your FFT.
Do the FFT on the entire signal. This may be done in multiple segments (frames).
In the frequency domain, subtract the result of step 1 from step 2. If you're using multiple segments in step 2, then simply subtract the result from step 1 from every segment. This will effectively filter your signal.
Do an inverse fft, which will bring your signal back into the time domain.

Note that I said that this technique will only work if the noise is uniform and by uniform I mean that its spectral content is rather constant. But you can also have regions in your audio where the noise is not present at all so obviously trying to remove the noise from these noiseless regions should be avoided. If your noise removal process is manual (at least in part), it would be possible (albeit tedious) to label such regions and then have the algorithm skip them.


How to prove this theorem about the Z transform and final value theorem? Claim:

  If $\lim_{k\rightarrow\infty} x[k]$ exists and is finite then $X(z)$, the Z-transform
  of $x[k]$, has no poles in the region $|z|>1$ and at most  1 pole at $z = 1$.

Attempt:
\begin{align*}
   X(z) &= \sum_{k\ge 0} x[k]z^{-k}\\\
   H(z)/G(z) &= \sum_{k\ge 0} x[k]z^{-k}\\\
\end{align*}

First prove that no pole can be in $|z|>1$
\begin{align*}
   H(z)/G(z) &= \sum_{k\ge 0} x[k]z^{-k}\\\
   H(z)/[(z-a)G'(z)] &=  \sum_{k\ge 0} x[k]z^{-k},  a>1\\\
   H(z)/[(z-a)G'(z)] &= x[0] + x[1]z^{-1} +x[2]z^{-2}+x[3]z^{-3}...\\\
   H(z)/G'(z)     &= (z-a)(x[0] + x[1]z^{-1} +x[2]z^{-2}+x[3]z^{-3}...)\\\
\end{align*}

Can someone help carry this further?

Let me show you a simple way to see this property. Assume $x[k]$ is a causal sequence and let

$$x[\infty]=\lim_{k\rightarrow\infty}x[k]$$

be finite. Then the sequence $x[k]$ can be written as

$$x[k]=x[\infty]u[k]+y[k]\tag{1}$$

where $u[k]$ is the unit step sequence, and $y[k]$ is a causal sequence that decays to zero as $k\rightarrow\infty$. Taking the $\mathcal{Z}$-transform of (1) gives

$$X(z)=x[\infty]\cdot\frac{z}{z-1}+Y(z)\tag{2}$$

The poles of $X(z)$ are determined by the two terms on the right-hand side of (2). Since $y[k]$ is a causal decaying sequence, its $\mathcal{Z}$-transform $Y(z)$ must have all its poles inside the unit circle. Since the first term on the right-hand side of (2) can only contribute a single pole at $z=1$, $X(z)$ cannot have any poles outside the unit circle. If $x[\infty]\neq 0$, $X(z)$ has exactly one pole at $z=1$ from the first term on the right-hand side of (2). If $x[\infty]=0$, i.e. $x[k]$ decays to zero as $k\rightarrow\infty$, then that first term disappears and $X(z)$ has no pole at $z=1$. This verifies the claim: if $x[k]$ has a finite limit as $k\rightarrow\infty$, $X(z)$ cannot have any poles outside the unit circle, and, if $x[\infty]\neq 0$, it has a single pole at $z=1$ in addition to its poles inside the unit circle.

Also note that the final value theorem directly follows from (2). Multiplying both sides of (2) by $z-1$ gives

$$(z-1)X(z)=x[\infty]z+(z-1)Y(z)\tag{3}$$

Taking the limit $z\rightarrow 1$ yields the final value theorem of the $\mathcal{Z}$-transform:

$$\lim_{z\rightarrow 1}(z-1)X(z)=x[\infty]\tag{4}$$

Note that the limit $\lim_{k\rightarrow\infty}x[k]$ must be guaranteed to exist in order for (4) to make any sense.


Reconstruction of Contiuous - Time Signals In terms of analog signals, we can represent digital signal as :

$$
x[n] \triangleq x_{a}(nT) = \int_{-\infty}^{\infty}X_{a}(f) \, e^{j2\pi f nT} \ \mathrm{d}f
$$

While if we focused on the integral on the right side and according to the Digital Signal Processing by John Proakis, chapter 6.1, we can rewrite it into :

$$
\int_{-\infty}^{\infty}X_{a}(f) \, e^{j2\pi f nT} \ \mathrm{d}f = \sum_{k=-\infty}^{\infty} \int_{(k-1/2)F_{s}}^{(k+1/2)F_{s}}X_{a}(f) \, e^{j2\pi nf/F_{s}}  \ \mathrm{d}f
$$

where $ F_{s} \triangleq \frac1T $.
My question is how the second equation comes up ? what does the interval of $(k-1/2)F_{s}$ to $(k+1/2)F_{s}$ means ?

Furthermore, it is stated in the book that "observing the  $X_{a}(f)$ in the interval of  $(k-1/2)F_{s}$ to $(k+1/2)F_{s}$ is identical to $X_{a}(f-kF_{s})$ in the interval of $-F_{s}/2$ to $F_{s}/2$". Is there any explanation or derivation how both things are identical ?

Thank you so much, hope that my question is clear enough

Because sampled signals Spectrum will have copies of original spectrum at multiples of $f=F_s$. Even if aliasing happened due to choice of $F_s$, the spectrum of sampled signal $x(n.T_s)$ will be periodic in frequency always with period $F_s$. And you can take any one such $k^{th}$ copy or period and integrate from $\frac{k-1}{2}F_s$ to $\frac{k+1}{2}F_s$. Now if you vary $k$ from $-\infty$ to $\infty$, then it is equivalent to integration on LHS.

Key: Sampling at $F_s$ will create a spectrum which is periodic in frequency with period $F_s$.


What's the noise in this signal ? (Beginner question) There is this raw ECG signal that I have obtained its frequency content in Matlab. The signal was sampled at 600Hz. The following is signal and the frequency content.  I want to know what the interferences are. I'm very new to DSP, and I feel like I'm lacking some common sense to understand a signal by looking at it.
It looks like to me the isolated peaks are the noise. How can I tell which one the noise is?

From the spectrogram (frequency domain plot) you have a large signal at 60 Hz and harmonics. This will be mains pickup if this was recorded in a area with a 60 Hz mains supply. On the upper time plot this signal has a peak to peak amplitude of about 0.5 units. (Note that there are 30 cycles over 500 ms.)

It looks like to me the isolated peaks are the noise.

As JRE and MBaz suggest, it helps to know what your wanted signal looks like. The Wikipedia article on ECG gives an idealized view of a ECG trace (as well as some actual plots). You should be able to see that the isolated peaks are actually part of your wanted signal. Your patient seems to have a pulse rate of about 75 bpm.
All your wanted signal seems to be below 25 Hz. To get rid of the mains interference one way would be to low pass filter the signal or use a comb notch filter to remove 60 Hz and its harmonics. A notch filter may have less effect on your wanted signal but whatever filter is used you need to consider how its amplitude and phase response will affect your wanted signal.
An alternative approach is subtraction. You generate a local copy of the interfering signal and subtract it from the raw input leaving behind the wanted ECG. As this isn't based on a filter there are no unwanted amplitude or phase effects on the signal.


Output of a discrete-time LTI system different form than input? This question is related to this one. I'm going through old exams for a 2nd year systems and transforms course, and came across this question. I'm posting this question just in case my other question doesn't actually represent this issue properly...

The question states:

  Suppose that we are given the following information about an LTI system: if the input to the system is $x_1[n]=(1/6)^n u[n]$, then the output is $y_1[n]=[a(1/2)^n +10(1/3)^n]u[n]$. If $x_2[n]=(-1)^n$, then the output is $y_2[n]=(7/4)(-1)^n$ where $a$ is a real number. Determine the transfer function of this system and the value of the number $a$.

The solution proceeds by determining the z-transform of $x_1$ and $y_1$ and stating the transfer function as $H(z)=Y_1(z)/X_1(z)$.

However, it seems to me that it's impossible for this system to generate an output of $(1/2)^n+(1/3)^n$ if the input is $(1/6)^n$,and hence the question is somewhat contrived if not outright bogus - see the discussion pertaining to my other question.

What am I missing here?

NOTE: The answer is not complete; however, it seems the problem cannot be solved.

We have the following z-transforms of the input signals:
$$x_1[n]=(1/6)^n\cdot u[n]\quad\longleftrightarrow\quad X_1(z)=\frac{1}{1-(1/6)\cdot z^{-1}}$$
$$x_2[n] = (-1)^n\quad\longleftrightarrow\quad X_2(z)=\frac{1}{1+z^{-1}}$$

For the output signals we have the following z-transforms:
$$y_1[n] = [a(1/2)^2 + 10(1/3)^n]u[n] \quad\longleftrightarrow\quad Y_1(z) = a\frac{1}{1-(1/2)\cdot z^{-1}} + 10\frac{1}{1-(1/3)\cdot z^{-1}}$$
$$y_2[n] = (7/4)(-1)^n\quad\longleftrightarrow\quad Y_2(z) = (7/4)\frac{1}{1+z^{-1}}$$

For both pairs of signals we know that the transfer function equals $H(z)=\frac{Y(z)}{X(z)}$. From $\{X_2(z),Y_2(z)\}$ we know that $H(z) = 7/4$ (Which is trivial due to linearity, i.e., it follows also directly from $\{x_2[n],y_2[n]\}$).

What remains is to solve $Y_1(z) = H(z)\cdot X_1(z)$ for $a$; however, I don't get a proper result.

EDIT: Another way to appraoch this last step is to solve the following:
$$a\frac{2}{2-z^{-1}} + b\frac{3}{3-z^{-1}} = c\frac{6}{6-z^{-1}}$$
where you may choose $b=10$ if desired.

In matrix representation:
$$\mathbf{A}\left[\begin{array}{1}a\\b\\c\end{array}\right] = \mathbf{b}$$
the matrix $\mathbf{A}$ is a full rank size 3x3 matrix, which shows that no scalar solution for $H(z)=c$ with $b=10$ exists.

EDIT2:
The solution seems to be that $H(-1)=7/4$ and $a=-9$. We already found that $H(z)=7/4$ for any $z$, based on the signals $\{X_2(z), Y_2(z)\}$. When evaluating $\{X_1(z),Y_1(z)\}$ for $z=-1$, the solution $a=-9$ follows easily. However, for another value of $z$, say $z=-2$ we obtain $a\approx -9.375$.

Conclusion, the decomposition only holds for $z=-1$.


Nyquist frequency , sampling distance I have few questions I tried to solve regarding nyquist theorem, and I would like to see your opinion if I'm doing it correctly?(one I know the answer second one not sure).

1.Let $f(x)$ and $g(x)$ be functions with max frequencys $B_1$ and $B_2$, what is the max sampling distance allowed according to nyquist to the function $h(x)=f(x)g(x)$.

My answer: according to the convolution theorem $H(x)=F(x)*G(x)$(* for convolution), so $H(x)=\int_{-B_1}^{B_1}F(u)G(x-u)du$ so $u \leq B_1$ now we knowthat $x-u \leq B_2$ so we get $x\leq B_2+B_1$ and this $x$ is our max frequency, now I can say that our sample rate is $f_s\geq 2(B_1+B_2)$ so our sample max distance denote $T_s$ is $T_s \leq \frac{1}{2(B_1+B_2)}$
is that correct?

2.lets mark function $f(x,y)=sinx*cosy$ what is the max sampling rate?

to be honest here I'm really not sure, I think its zero by intuition but not really sure if I'm correct and if so why?

Please note: those are the notation I saw when I found those questions, thanks

Your answer to the first question is right, provided that the signals involved are real and lowpass (baseband) nature.

The second question is tricky, assuming ideal infinite length sinusoidals, your intuition is correct in the sense that their Fourier transforms multiplied would yield zero (unless their frequencies were the same), and has zero bandwidth. Also note that the convolution was divergent in the first place, if that's a problem for you.

However, any practical implementation of that algorithm would include windows and therefore result of Frequency domain multiplication will be nonzero and the result would be non bandlimited, hence cannot be sampled without anti-aliasing filtering.


Confusion on using FFTW 3.3.4 (C version) for real and even data sets I am using FFTW3 to transform the autocorrelation function of a particular function. I expect the transform to be real since the autocorrelation function is symmetric with respect to negative and positive lags and its Fourier transform give the power spectrum. However, I get nonzero complex part and negative real parts after the transformation. What I may be doing wrong?

I am providing the code below:

/* Compilation and execution command:
 * gcc q1part2.c nrutil.c qtrap.c trapzd.c -lfftw3 -lm -o q1part2 && ./q1part2
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>
#include "nr.h"
#include "nrutil.h"
double h = 0.01;
double a = 0.0, b = 10.0, p;

double func(double t)
{
    return  cos(10 * t) + sin( 20 * t);
}
double intgrd(double t)
{
    return func(t) * func(t + p);
}

int main()
{
    int  N = round(((b - a) / h) + 1);
    double y[2*N - 1], t[2*N - 1]; /* t array hold time delays */
    int i;
    fftw_complex *in, *out;
    fftw_plan MYPLAN;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
    MYPLAN = fftw_plan_dft_1d(2*N - 1, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    for(i = 0; i < 2*N - 1; i++)
    {
        /* The definition of autocorrelation from the book "A First Course on
        * Computational Physics" by Paul DeVries was used.
        */
        p = i * h - 10.0;
        t[i] = p;
        if(p >= 0)
        {
            y[i] = sqrt(1 / (2 * M_PI)) * qtrap(intgrd, a, b - p);
            in[i] = y[i];
        }
        else
        {
            y[i] = sqrt(1 / (2 * M_PI)) * qtrap(intgrd, -p, b);
            in[i] = y[i];
        }
        //printf("%f \t %f \n", t[i], y[i]);
    }
    fftw_execute(MYPLAN);
    for(i = 0; i < N; i++)
    {
        printf("%f \t %f  %f \n", i * 2*  M_PI / (2 * (b - a)), creal(out[i]),cimag(out[i]));
    }
    fftw_destroy_plan(MYPLAN);
    fftw_free(in);
    fftw_free(out);
    return 0;
}

The function qtrap is a numerical integration algorithm provided in the code from the second edition (plain C) of the book Numerical Recipes. I am really confused about the usage of the library and would like to get pointed in the right direction.

Thanks in advance

If the imaginary part is very small you are probably not doing anything wrong.  This is a common occurrence due to floating point quantization noise.  As long as you know the result should be real and the imaginary part is small, just drop the imaginary part.


PSD from autocorrelation in MATLAB I am trying to simulate a simple stochastic process defined by the equation:
\begin{equation}
\frac{1}{v}\frac{db}{dt} +\Gamma_0 b= \sqrt{\sigma}R(t),
\end{equation}
where $R(t)$ is a zero-mean white noise (delta-correlated in time), and $v$, $\Gamma_0$ and $\sigma$ are all real positive parameters. I am simulating this using an Euler-mayurama scheme, namely
\begin{equation}
b(t+\Delta t) = b(t) -v \Gamma_0 b(t) \Delta t + v\sqrt{\sigma} \sqrt{\Delta t} N(0,1),
\end{equation}
with $N(0,1)$ being a standard normal random number, for a finite interval $t\in [0,T]$. I am trying to compute the power spectral density using the Wiener-Khintchine theorem
\begin{equation}
S[f] = \int_{-\infty}^{\infty} e^{-2\pi f\tau} \left\langle b(t)b^*(t+\tau)\right\rangle d\tau,
\end{equation}
from the autocorrelation function $\left\langle b(t)b^*(t+\tau)\right\rangle$. Analytically I can compute this and it gives me
\begin{equation}
S[f] = \frac{v\sigma}{2\Gamma_0}\frac{2v\Gamma_0}{(v\Gamma_0)^2 + (2\pi f)^2},
\end{equation}
But when I try to compute the same numerically, my function looks completely different. Am I implementing this correctly? (see Matlab code below)

%% Autocorrelation test
close all; clear all; clc;

% Inputs
dt = 1e-3;
T = 2;
tv = 0:dt:T;
Lt = length(tv);

v = 3.5;
Gamma0 = 20.3;
sigma = 0.75;
a = v*Gamma0;

% Frequency spectrum properties
Fs = 1/dt; % (Hz) sampling frequency, based on Ben's paper

sims = 1; % number of simulations
U_Walks = zeros(Lt,sims);
for nn = 1:sims
    b = zeros(Lt,1);
    for ii = 1:Lt-1
        b(ii+1) = b(ii) - v*Gamma0*b(ii)*dt + ...
            v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);
    end
    Rxx(:,nn) = xcorr(b);
    U_Walks(:,nn) = b;
    disp(nn);
end

%% Stats
LR = length(Rxx);
N = 2^nextpow2(LR);
Y = fft(Rxx,N); % Taking only one side of the y-axis symmetric FFT
Y = Y(1:N/2+1);
mY = (Y).^2/N;
F = Fs*(0:(N/2))/N; % calculating the frequency range for the x-axis

S_f = 10*log10(mY);

S_f_analytic = 10*log10((v*sigma/2/Gamma0.*(2*a./(a^2 + (2*pi*F).^2))));

%% PLOTS
close all;

figure('units','normalized','outerposition',[0 0 1 1])
set(gcf,'Color','w');
plot(F,S_f,'r',F,S_f_analytic,'b','LineWidth',3);
legend('Numeric','Analytic');
title(['PSD']);
xlabel('f (Hz)'); ylabel('S[f] (dB)');
set(gca,'FontSize',20);

There are two mistakes in your code/method. The first is the term $\sqrt{\Delta t}$ in your second formula; it should be replaced by $\Delta t$. The second is in the computation of the  power spectrum from the estimated auto-correlation. What you do is square the result of the FFT Y to obtain mY, but that's not correct. First of all, Y is complex-valued, and second, the Fourier transform of the auto-correlation is the power spectrum, so there's no need to square the outcome of the FFT, you just need to compute its magnitude.

In sum, replace the line

v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);

by

v*sqrt(sigma)*dt.*normrnd(0,1);

and

mY = (Y).^2/N;

by

mY = abs(Y)/N;

This should give you the desired result, as shown in the figure below:


Peak Detection on ECG signal I've been working on a ECG signal collected during exercise through my chest mount heart rate monitor. My goal is to accurately detect the location of the R-wave; R-wave is the narrow peak similar to the ones I circled in red.

Can you suggest any method for preprocessing the data so that it will be easier to extract the signal?

I have tried wavelet transform but I am not getting any luck.

You can download the data here. Sampling frequency is 150 Hz.

A typical approach would be the Pan-tompkins algorithm. You can also apply a first order low pass filter or a mean filter to improve results.


Chroma-Subsampling: How to properly calculate the data-rate I am having a hard time understanding how to calculate the data rate when leveraging chroma-subsampling in e.g. a Y'UV image:

I have the following examples for a calculation:

Image resolution: 352*288
Frequency: 25 fps

For (4:4:4) the example calculation goes as follows:

(352px * 288px) * 3 color channels * 25 fps * 8 bit = 60 825 600 bit/s

So far so good.

But now comes (4:2:0):

(352px*288px) * 1.5 color channels * 25 * 8 = 30 412 800 bit/s

Now, trying to transfer this example to e.g. (4:1:1) I realized I am not sure whether I have a correct understanding how the ratio of 1.5 color channels is being calculated.

My first guess for the calculation was in case of (4:2:0):
2/4*3=1.5 color channels
Likewise for (4:1:1) I would calculate the ratio for the color channels as:

1/4*3=0.75 color channels

But I am simply not sure whether this would be the correct way to go.
As an alternative, I was thinking along following lines:

color channels (4:1:1): 1 Y' + 1/4 UV = 1.25 color channels

Now which would be the correct way of doing this properly?

Here is a reference to know the difference. I am adding the essential diagram as well:

what you realize is that 4:1:1 has full vertical resolution but 1/4 of Horizontal resolution, where as 4:2:0 has half of vertical as well half horizontal resolution.

However, overall, 4:2:0 and 4:1:1 will have the same number of samples and hence the same bit rate.