What is definition of independent random variable I wan't to ask that if

E{X}=0

E{Y}=0

and    E{XY}=0

then how can I verify if the two random variables are independent or not.

X , Y are both continuous random variables {I am not able to recall the exact equations for X , Y but  I am sure about the information given above and want some hint of how to proceed}.

What I know about independence {the definition that is causing confusion}

And so E{XY} should be = E{X}E{Y}

Which part of the below derivation doesn't sufficiently satisfies that X and Y are independent {where did it went wrong}

Transition from first line to second line is only possible if two random variables are independent and after it it's just rearrangement of terms.

No, NO and $\mathbf{NO}$.

Knowing that $E[X]=E[Y]=0$ and $E[XY]=0$ (or more generally that $E[XY]=E[X]E[Y]$ or equivalently, $E[XY]-E[X]E[Y]=0$) does not help in the least in proving or deducing that $X$ and $Y$ are independent random variables.  The conditions stated in the previous sentence are necessary for $X$ and $Y$ to be deemed independent random variables, but they are a very far cry from being sufficient for independence to hold. That is, if $E[XY]=E[X]E[Y]$, then we might suspect that $X$ and $Y$ might be independent random variables, but we cannot jump to the conclusion that they are independent random variables based on just the evidence that $E[XY]=E[X]E[Y]$.

  Definition: $X$ and $Y$ are said to be independent random variables (sometimes more prolixly mutually independent random variables) if their joint distribution function (joint CDF) $F_{X,Y}(u,v) = P\{X \leq u, Y\leq v\}$ satisfies
  $$F_{X,Y}(u,v)=F_{X}(u)F_{Y}(v) ~\text{for all real numbers}~ u ~\text{and}~v.\tag{1}$$
  Here, $F_{X}(u)=P\{X \leq u\}$ and $F_{Y}(v)=P\{Y\leq v\}$ are the marginal CDFs of $X$ and $Y$ respectively. It is important to remember the qualifier to the equality $F_{X,Y}(u,v)=F_{X}(u)F_{Y}(v)$: it must hold at all points in the $u$-$v$ plane.. Inanities such as "$X$ and $Y$ are independent when $X=5$ and $Y=3$ but not otherwise" are not permitted.

If $X$ and $Y$ are discrete random variables, then $(1)$ is equivalent to their joint _probability mass function (joint pmf) $p_{X,Y}(u_i, v_j)$ factoring into the product $p_X(u_i)p_Y(v_j)$ of the marginal pdfs
for all $i$ and $j$.

If $X$ and $Y$ are jointly continuous random variables, then $(1)$ is equivalent to their joint _probability density function (joint pdf) $f_{X,Y}(u, v)$ factoring into the product $f_X(u)f_Y(v)$ of the marginal pdfs
_for all real numbers $u$ and $v$.

Now, as the OP correctly points out, if $X$ and $Y$ are independent, then it is true that $E[XY]=E[X]E[Y]$, that is,
$$X, Y ~\text{independent random variables} \implies E[XY]=E[X]E[Y]$$ but

  the reverse implication
  $$E[XY]=E[X]E[Y]\implies X, Y ~\text{independent random variables}$$
  $$\mathbf{DOES~NOT~HOLD}$$

Still don't believe a word of all this malarkey? Consider discrete random variables $X$ and $Y$ such that
$$p_{X,Y}(1,0) = p_{X,Y}(-1,0) = p_{X,Y}(0,1) = p_{X,Y}(0,-1) = \frac 14.$$
Verify that

$X$ and $Y$ take on values $\pm 1$ with equal probability $\frac 14$ and value $0$ with probability $\frac 12$.
$XY=0$ always and so $E[XY]=0$ also, and that $E[X]=E[Y]=0$ also.
The joint pmf does not factor into the product of the marginal pdfs for at any of the four points where the joint pmf is nonzero.


Compare feature vectors of an image in MATLAB I am working on OCR. I have extracted features of 26 alphabets, and I want to compare them with an experimental character feature.

I have tried using Euclidean distance and correlation functions, but I am not getting good efficiency. What are some alternate methods for the same?

If you just need OCR to work, then please try the ocr function in the Computer Vision System Toolbox.

If you need to write your own, take a look at this example, showing how to train an SVM classifier using HOG features for recognizing hand-written digits.


How to remove salt & pepper noise from speech signal? I'm having a lot of problems with my project:

I added salt and pepper noise by inbuilt function
I DWT function with Haar family
Now I want to remove noise from original speech signal.
I don't know how to do this

If you are forced to use the Discrete Haar Transform, then note that the Discrete Haar Mother Wavelet is nothing but a [1 -1], or in other words, a discrete differentiator, as shown below:

As such, it can/will severely accentuate sudden transitions in your waveform, ie, salt and pepper noise.

Note, that the first dyadic scale of the DHT will be composed of coefficients computed from the signal correlated with the Haar Mother wavelet. At this level, salt&pepper noise will correlate very well, and will give you large wavelet co-efficients in the wavelet space. One can then simply null the large values out, and then perform an inverse Discrete Haar Wavelet Transform, (IDHT).

It is strange that you would be required to use the Discrete Haar Wavelet Transform, (over, say, the Stationary Wavelet Transform), because nulling co-efficients in the DWT space and inverse transforming can lead to aliasing. This effect might be small for the problem you describe however.

Another possibility would be to simply perform the DWT, identify the large wavelet co-efficients within the last half of the DWT transform, (corresponding to high correlations between SP noise and the Haar Mother Wavelet), identify those large co-efficients, and simply null them out in the time domain, without doing an inversion. This will get around the problem of aliasing in the inverse transform.


How to find a variance of sample sequence I have a sequence such as

$$r[n] = y[n]v[n]$$

$y[n]$ and $v[n]$ are zero-mean and statistically independent.

I need to find a variance of $r[n]$ and show that it is white and equal to $\sigma ^2_y\sigma^2_v$.

I know that $y[n]$ has variance of $\sigma ^2_y$ and I also know that
$$Cor(v_n, v_{n+m}) = \sigma ^2_v\delta[m]$$

This is what I have so far :

Finding the variance of $r[n]$:
$$\sigma_r^2 = \mathbb{E}([r[n]-m_r)^2]=\mathbb{E}[r^2[n]] - 2\mathbb{E}[r[n]]m_r+m_r^2$$

Then

$$\sigma_r^2 = \mathbb{E}[r^2[n]] - 2m_r m_r+m_r^2 = \mathbb{E}[r^2[n]]-m_r^2 $$

I am stuck at this point. Could you suggest how do I proceed from here?

The process is WSS random process.

$\sigma_r^2=\sigma_y^2\sigma_v^2$ iff

1) $\mathbb{E}\{y\} = \mathbb{E}\{v\}=0$.

AND

2) $y$ and $v$ are independent, OR uncorrelated and their squares ($y^2$ and $v^2$) are also uncorrelated.

Proof: assuming independence,

$\sigma_r^2=\mathbb{E}\{(yv-\mathbb{E}\{yv\})^2\}$ (by definition)

$~~~~=\mathbb{E}\{(yv-\mathbb{E}\{y\}\mathbb{E}\{v\})^2\}$ (by independence)

$~~~~=\mathbb{E}\{y^2v^2+\mathbb{E}^2\{y\}\mathbb{E}^2\{v\}-2yv\mathbb{E}\{y\}\mathbb{E}\{v\}\}$

$~~~~=\mathbb{E}\{y^2v^2\}-\mathbb{E}^2\{y\}\mathbb{E}^2\{v\}$

$~~~~=\mathbb{E}\{y^2\}\mathbb{E}\{v^2\}-\mathbb{E}^2\{y\}\mathbb{E}^2\{v\}$ (by independence, since the square function is a measurable function).\

But, $\mathbb E\{y^2\} = \sigma_y^2 + \mathbb E^2\{y\}$
and $\mathbb E\{v^2\} = \sigma_v^2 + \mathbb E^2\{v\}$ and so,
\begin{align}\sigma_r^2&=(\sigma_y^2 + \mathbb E^2\{y\})(\sigma_v^2 + \mathbb E^2\{v\}) - \mathbb E^2\{y\}\mathbb E^2\{v\}\\
&= \sigma_y^2\sigma_v^2 + \sigma_y^2\mathbb E^2\{v\} +
\sigma_v^2\mathbb E^2\{y\}
\end{align}

which does not equal $\sigma_y^2\sigma_v^2$ unless both $y$ and $v$ have zero mean as the OP has stated: having at least one mean be zero (which does not preclude the other mean from being nonzero) is not enough (as was claimed in an earlier version of this answer).

Proof: assuming $y,v$ uncorrelated and $y^2,v^2$ also uncorrelated,

\begin{align}\sigma_r^2 &=\mathbb{E}\left[(yv-\mathbb{E}[yv])^2\right] &{\scriptstyle{\text{(by definition)}}}\\
&=\mathbb{E}\left[(yv-\mathbb{E}[y]\mathbb{E}[v])^2\right]  & {\scriptstyle{E[yv]=E[y]E[v] ~\text{since}~ y ~\text{and}~v~\text{are uncorrelated}}}\\
&=\mathbb{E}\left[y^2v^2\right] &{\scriptstyle{\mathbb E[y] = \mathbb E[v]  = 0~\text{by assumption}}}\\
&=\mathbb{E}[y^2]\mathbb{E}[v^2] &{\scriptstyle{E[y^2v^2]=E[y^2]E[v^2] ~\text{since}~ y^2 ~\text{and}~v^2~\text{are uncorrelated}}}\\
&= \sigma_y^2\sigma_v^2 &{\scriptstyle{\text{since}~\mathbb E[y]= \mathbb E[v] = 0~\text{by assumption}}}
\end{align}

Note that $y$ and $v$ being independent is a subcase of the above assumptions about $y$ and $v$ being uncorrelated and $y^2$ and $v^2$ also being uncorrelated.  When $y$ and $v$ are independent (and hence uncorrelated), so are $y^2$ and $v^2$ independent (which implies that $y^2$ and $v^2$ are uncorrelated).


Image derivative in the Fourier space It can be shown that the image derivative in the $x$ direction is given by:
$$
\frac{\partial f}{\partial x}=\frac{2\pi i}{N} \mathcal F^{-1}\left(u\cdot \mathcal F(f(x,y)\right)
$$
where $N$ is the width of the image in the $x$ direction and $u$ denotes a variable in the fourier space ($F(u,v)$).

Now, in the image space, derivation is a convolution with this matrix:
$$\partial_x =\left(\begin{array}{cc}
-1 & 1\end{array}\right)$$

this means (from to convolution theorem) that the fourier transform of $\partial_x$ is a matrix $D(u,v)$ such that $D(u,v)=u$.

so I tried to see this in python, but something seems to be wrong here:

d = np.zeros((21,21))
d[10,10] = 1
d[10,9] = -1
D = np.fft.fft2(d)
plt.matshow(np.fft.fftshift(D))

Can anyonw see where I am wrong?

The Fourier transform property in your first equation is for continuous time but the image is a discrete time signal. $\partial_x$ isn't exactly differentiation, it's a discrete time filter which does finite differencing.

Since you are only taking derivatives in one direction, let's think about a 1-D signal for simplicity (or, equivalently, let's assume your image has only one row). Let the signal be $f(x)$ and its continuous time Fourier transform $F(u)$. Then we have the transform pair $\frac{df(x)}{dx} \leftrightarrow j 2\pi u F(u)$. In discrete time, if $x[n] \leftrightarrow X[k]$, you can't really say $\{x[n]-x[n-1]\} \leftrightarrow j2\pi k X[k]$.

Instead, we can think of $\partial_x$ as a high pass filter. It then agrees with the intuitive meaning of the continuous time Fourier transform property too (where multiplying by $u$ attenuates lower frequencies more than higher frequencies). Try this code, to see the high pass response:

d = np.zeros(10000)
d[0] = -1.0
d[1] = 1.0
D = np.fft.fft(d)
plt.plot( np.abs(D) )


Bandwidth of a signal using MATLAB I need to find the bandwidth the signal $x = 50 \operatorname{sinc}(125(t-10))^2$.

(Where $\operatorname{sinc}(x)=\dfrac{\sin (\pi x)}{\pi x}$ [normalized] or $\operatorname{sinc}(x)=\dfrac{\sin x}{x}$ [unnormalized])

How would I go about doing this in MATLAB?

That actually would be a great interview questions. Let's derive it from scratch:

$\text{sinc}(x)$: The Fourier Transform of a $\text{sinc}(x)$ is a rectangle. The normalized $\text{sinc}()$ function results in a bandwidth of 1 Hz (from -0.5Hz to 0.5Hz)
$t-10$: The $(t-10)$ term is irrelevant to the bandwidth: time shift
is equivalent to multiplying with $e^{-j\omega t_0}$, where $t_0$ is time shift, in
the frequency domain. That changes the phase but not the amplitude
or bandwidth
$125*t$: that is a time compression which results in frequency
stretching. Since we compress time by 125 the bandwidth increases by
125, i.e. to 125 Hz (from -62.5 Hz to 62.5 Hz)
$(\cdot)^2$: Squaring is multiplication in the time domain which corresponds to convolution in the frequency domain. The convolution of two equal rectangles is a triangle. The bandwidth doubles to 250 Hz (from -125 Hz to 125Hz). In general the bandwidth of the multiplication of two arbitrary signals is the sum of the individual bandwidths. That's a simple consequence of the multiplication laws for sine and cosine: you get only sum and difference frequencies.
$50*(\cdot)$: that's a simple scale factor that doesn't change the shape of phase or amplitude.

So overall we would expect the spectrum to be of triangular shape (on a double linear plot) going from -125 Hz to +125 Hz with a total bandwidth of 250 Hz.

Writing the Matlab code for this is not as simple as it looks. The code itself is trivial, however you have to pick a sample rate and a proper time axis. In order to do this you need to assess the general shape of the signals first. For example the time domain signal technically goes from $-\infty$ to $+\infty$ and is centered around $t$ = 10s. t = 0:.01:5, really won't do the trick. So in essence you have to do the theoretical analysis before you can write the code (which tends to be a recurring theme in DSP :-))

% 1000 HZ sample rate, since we now the BW is smaller than 500 Hz
fs = 1000;
% pick a good time length. Because of the time shift we need significantly more
% than 10s. Technically the sinc() has inifite length so we need to extend
% it far enough to cover the vast part of the energy AND we need to start
% at negative times, not at t = 0;
n = 2^16; % that's about 60s. Should be plenty
t = (-n/2:n/2-1)'/fs;   % time axis in seconds, symmetric around t = 0;
x = 50*sinc(125*(t-10)).^2;
% frequency axis:
f0 = fs/n;  % frequency resoution in Hz
faxis = (-n/2:n/2-1)'*f0;   % frequency axis in Hz seconds, symmetric around f = 0;
fx = fft(x);
plot(faxis,circshift(abs(fx),n/2));
title('Spectrum of 50*sinc(125*(t-10)).^2');
xlabel('Frequency in Hz');


Only random signal has PSD? I heard from someone online who said: "Only random signal has PSD, determinate signal does not have. For example, period signal does not have PSD."

I am very astonishing about this statement, but it seems he is quite sure about his idea. In my opinion, both random signal and determinate signal have PSD, and it is the Fourier transform of the autocovariance of the signal. However, only stationary random and period signal have stable PSD. Am I right?

Autocovariance is defined on random signals.

$$
C_{XX}(t,s) = E\left[(X_t - \mu_t)(X_s - \mu_s)\right] = E[X_t X_s] - \mu_t \mu_s
$$

The expected value only makes sense for random processes. You could treat a deterministic process as a random process where everything except for the one signal you observed has probability zero.


$\int^{\frac12}_{-\frac12}X(f)e^{j2\pi nf}df=\frac{1}{F_s}\int^{\frac {F_s}{2}}_{-{\frac {F_s}{2}}}X(F)e^{j2\pi nF/F_s}dF$ for $f=\frac{F}{F_s}$? I am using Digital Signal Processing Principles, Algorithms, and Applications 4th edition and by Proakis. Here is what I don't understand

say the signal $x_a(t)$ has its Fourier Transform $X_a(F)$ and digital signal $x(n)=x_a(nT)$ which is $x_a(t)$ after sampling with sampling frequency $F_s$  has its DTFT $X(f)$

since
$$x_a(t)=\int^{\infty}_{-\infty}X_a(F)\,e^{j2\pi Ft}dF\tag{1}$$

because $x(n)=x_a(nT)$ and $nT=\dfrac{n}{F_s}$
$$x(n)=x_a(nT)=\int^{\infty}_{-\infty}X_a(F)e^{j2\pi nF/F_s}dF\tag{2}$$

Using inverse DTFT for $x(n)$ then

$$\int^{\frac12}_{-\frac12}X(f)e^{j2\pi nf}df=\int^{\infty}_{-\infty}X_a(F)e^{j2\pi nF/F_s}dF\tag{3}$$

for $f=\dfrac{F}{F_s}$ then

$$\int^{\frac12}_{-\frac12}X(f)e^{j2\pi nf}df=\frac{1}{F_s}\int^{\frac {F_s}{2}}_{-{\frac {F_s}{2}}}X(F)e^{j2\pi nF/F_s}dF\tag{4}$$

That is what the book says and it use that to explain the reconstruction process in theory. However isn't it will be like this according to substitution in integral

$$\int^{\frac12}_{-\frac12}X(f)e^{j2\pi nf}df = \frac{1}{F_s} \int^{\frac {F_s}{2}}_{-{\frac {F_s}{2}}}X\left(\frac{F}{F_s}\right)e^{j2\pi nF/F_s}dF\tag{5}$$

I think this should be in Math Section but since it is about Digital Signal so I post it here.

You are absolutely right. With $f=F/F_s$ you obtain the equality

$$\int_{\frac12}^{\frac12}X(f)e^{j2\pi nf}df=\frac{1}{F_s}\int_{-F_s/2}^{F_s/2}X\left(\frac{F}{F_s}\right)e^{j2\pi nF/F_s}dF\tag{1}$$

which is also what I find in my third edition of the book (the right-hand side of $(1)$ equals the left-hand side of Eq. $(4.2.81)$ in my edition).


Faster Algorithm to convolve/correlate two sparse 1-D signals in python (or any language) I have two signals which I need to correlate or convolve. Each signal is sampled non-uniformly and the values of the signal I have with me are the timestamp and the magnitude of the signal at that timestamp. The values of the signal at all the other times can be assumed to be zero. The timestamps of the signal have a resolution running in microseconds.

An example of how the signal looks is shown below:

As can be seen, the resolution of the signal is in microseconds and signal is mostly sparse.

If I were to the convolve two signals of this type, I would first have to pad the signals with zeros (since I would have to discretise the signal). While the padding can be done with resolution of microseconds, The number of values to be multiplied becomes too big and the operation becomes increasingly slow. Most of the multiplications in this convolution would be multiplication of zeros(which are pretty much useless). I have therefore chosen a round off value of 2 places (0.xxxxxx becomes 0.xx),since I have to perform 40,000 similar convolutions. I have written my resampling function as shown below.

import numpy as np
import math

def resampled_signal_lists_with_zeros(signal_dict, xlimits):
  '''
      resamples the given signal with precision determined by the round function.
      signal_dict is a dictionary with timestamp as key and signal magnitude is the value of the key.
      xlimits is an array containing the start and stop time of the signal.
  '''
  t_stamps_list = list(signal_dict.keys())
  t_list = list(np.arange(int(math.floor(xlimits[0])), int(math.ceil(xlimits[1])), 0.005))
  t_list = [round(t, 2) for t in t_list]
  s_list = list()
  time_keys = [round(t, 2) for t in t_stamps_list]

  i = 0
  for t in t_list:
    if i < len(t_stamps_list):
      if t==time_keys[i]:
        s_list.append(signal_dict[t_stamps_list[i]])
        i+=1

      else:
        s_list.append(0)
    else:
      s_list.append(0)
  return t_list, s_list

The correlation of two signals padded in the above manner is done using scipy as follows:

from scipy.signal import correlate
output = correlate(s_1, s_2, mode='same')

The output calculated in the above manner is pretty slow .Since the signal is pretty sparse and most of the multiplications in the signal are multiplications of zero, I think there should be a better way to do the same operations. Is there a way to get the result of the convolutions of the two sparse signals faster?

Interpolating discontinuous waveforms is usually not a good idea.

The way I would approach your problem would be to recognize your signals as pulse trains. So I would assume that the convolution a single pulse from one waveform with a single pulse from the other waveform was also a pulse.

so signals A and B could be represented as a train of continuous Dirac functions.
$$
s_A(t)=\sum_{i=1}^{N-1} a_i \delta(t-\tau_A(i))
$$
$$
s_B(t)=\sum_{i=1}^{M-1} a_i \delta(t-\tau_B(i))
$$
where $\tau_A(i)$ and $\tau_B(i)$ are your non uniform arrival times.

There is an identity for delta functions
$$
f(t) \star \delta(t-a)= f(t-a)
$$
where $\star$ denotes convolution, so
$$
s_A(t)\star c_1\delta(t-a)= c_1\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+a))
$$
and
$$
s_A(t)\star \left[c_1\delta(t-a) + c_2\delta(t-b)\right]= c_1\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+a))+c_2\sum_{i=1}^{N-1} a_i \delta(t-(\tau_A(i)+b))
$$
which can be generalized to $s_A(t)\star s_B(t)$

you really just need calculate the time shifts and amplitude products.

I will leave it up to you how to handle the amplitudes but the code below assumes all the impulses are magnitude one.

clear all
close all
clc
%%
M=2 ; % number of spikes in waveform A
N=5; % number of spikes in waveform B
T=10; %  arbitrary duration
%%
waveform_a_times=rand(1,M)*T;
waveform_b_times=rand(1,N)*T;

%%
for i=1:N
    time_diff(i,:)=waveform_a_times+waveform_b_times(i);
end
conv_times=sort(reshape(time_diff,1,M*N));
figure(1)
subplot(3,1,1),stem(waveform_a_times,ones(1,M));
title('waveform A')
subplot(3,1,2),stem(waveform_b_times,ones(1,N));
title('waveform B')
subplot(3,1,3),stem(conv_times,ones(1,N*M));
title('convolution of waveform A with waveform B')

for correlation

$$
f(t)\star \delta(t+a)=f(t+a)$$

you’re also more likely to have coincident times for correlation which you should test for.


How to choose an interpolation method for seismic ground motion signals I am working with seismic ground motion records from the PEER database. They are acceleration measurements that are band limited to 50 Hz and sampled at 200 Hz.

I use these records to drive an explicit finite element analysis of a structure to determine its response under the earthquake. This analysis runs with a 1e-5 time step (100,000 Hz). If I were to input the sampled records directly the FEA program would linearly interpolate the time steps between the samples, which would distort the signal.

Therefore I need to reconstruct the original signal and resample it at 100,000 Hz (or high enough that the interpolation doesn't affect accuracy much).

One property that I would like to preserve is the double integral of the signal, which gives displacement. This keeps the structure from drifting away over time. I do not care about the calculation time, as it only needs to be done once and is trivial compared to the FEA calculation time.

A previous answer discussed several methods for interpolation but not how to choose from among them.

I have experimented with sinc interpolation using a Kaiser window with an artificially generated signal that I sampled and then reconstructed and I see decent results, but there are visible differences in the time domain.

Is there some way to determine what interpolation method (and parameters in the case of a windowed method) would give the most accurate results?

If there are differences then some likely sources might be  that the window on the Sinc wasn't wide enough, that the data wasn't really sufficiently bandlimited, or too much quantization noise in the samples.

The double integration of anything with system offsets or quantization noise is unlikely to be zero. (See random walk) One solution would be to compute the double integral of the first pass interpolated result and then bias the interpolator (noise filtered if needed) until the resulting double integral goes to zero.


Estimating a Low Frequency Signal Corrupted by High Frequency Noise This is a follow up to the question How to Pick a Journal for Publishing a Research Work? question. Suppose we take discrete samples of a low frequency signal corrupted by high frequency noise. We know the signal has no frequencies higher than  20%  of the Nyquist rate. We also know all the noise is at or above 50% of the Nyquist rate, and below the Nyquist rate. I will use the following as a specific example to experiment with.

\begin{align}&\textrm{measurement}(t) = \textrm{signal}(t) + \textrm{noise}(t)\\
&\textrm{signal}(t) = 4.4+\sin(0.06\pi t-0.1) + \sin(0.14\pi t+0.07) + \sin(0.2\pi t+0.4)\\
&\textrm{noise}(t) = \sin(0.06+0.5\pi t) + \sin(0.1-0.68\pi t) - \sin(0.04-0.74\pi t) + \sin(0.03-0.93\pi t)
\end{align}

The measurement is sampled at integer values of ($t$). A plot of the measurement samples and the continuous signal is below.

How well can we estimate $\textrm{signal}(t)$ at integer values of ($t$) using the current and past measurements? I was able to filter it using the Kalman filter below.

A plot of the signal and Kalman filter output is in the next plot. The final plot shows the error in the Kalman filter output which has a maximum magnitude of 1.865 and the RMS error is 0.7895. Can we make a Kalman filter that does better than the Kalman filter above? What known filter method would give the best estimates, and how good would the estimates be?

I learned about discrete lowpass filters 25 years ago, and I always assumed a problem such as the one above was a common application for them. At that time I didn't know this problem is called estimation, and I didn't know about the Kalman filter. Then I recently took this course, and another course Georgia Tech provided that goes into the details of estimation algorithms. However, the Georgia Tech courses mention nothing about using a discrete lowpass filter for estimation. I check several books on estimation, and found nothing about discrete lowpass filters. I work with a computer scientist who has been doing research in data fusion for 20 years, and he knows nothing about z-transforms, transfer functions, etc. Why is it that the estimation community ignores the use of a discrete lowpass filter for estimation? As far as I can tell it's the best approach to estimate the signal above.

To answer your final question:

  Why is it that the estimation community ignores the use of a discrete lowpass filter for estimation? As far as I can tell it's the best approach to estimate the signal above.

That's because you're feeding them the wrong signal model.

TL;DR: Use the right tool for the job!

Gory Details

Any time you start in with the Kalman Filter, you are assuming that the signal (measurement in your equations above) was generated as:
$$
x[k+1] = \mathbf{A} x[k] + u[k]\\
\mbox{measurement}[k] = \mathbf{H} x[k] + v[k]
$$
where
$$
\mathbf{A} = \left [ \begin{array}{cc}
1 & 1 \\
0 & 1
\end{array} \right]\\
\mathbf{H} = \left [ \begin{array}{cc}
1 & 0
\end{array} \right]\\
$$
and $u$ and $v$ are zero mean, Gaussian distributed, independent random variables with covariances
$$
\mathbf{Q} = \left [ \begin{array}{cc}
1.6 & 2.4 \\
2.4 & 4.8
\end{array} \right] \times 10^{-9}\\
\mathbf{R} = 2.35 \times 10^{-6}\\
$$
At least for the equations you cite. More general application of the KF is possible, but usually not done.

So, one thing that might be done to improve your estimates is to ensure that $\mathbf{Q}$ and $\mathbf{R}$ are closer to reality.  From the look of things, your process noise ($u[k]$) has about the same variance as your measurement noise ($v[k]$).  Make them closer.

More importantly, your signal is correlated with your noise because, rather than stochastic noise, you have harmonic noise. Even worse: you have synchronized harmonic noise (meaning the noise moves in lockstep with your signal).

So: the problem is that by using the Kalman Filter your signal model is all wrong.

A better model for the signal is in the frequency domain:  the signal is low frequency. The noise is high frequency.  So, a low pass filter will do a much better job of improving things that the Kalman Filter can hope to do.

If I apply the Kalman filter to your problem (with some modifications to your equations to take account of the difference between $\mathbf{X}[k|k-1]$ and  $\mathbf{X}[k|k]$) then I get the picture below.

The black curve with dots is the measurement.  The red curve is signal.  The blue curve is the output of the Kalman filter.  The green curve is the output of a 5 point moving average filter (a simple low pass filter).

The sum of squared errors for the Kalman filter is 118.63722. The same figure for the low pass filter is 14.18046.  Clearly, the signal model of a low pass filter fits better because the error is smaller.

R code to implement this is below.

# 26489
T <- 128;
t <- 0:(T-1)/T*70

signal <- 4.4 + sin(0.06*pi*t - 0.1) + sin(0.14*pi*t + 0.07) + sin(0.2*pi*t + 0.4)
noise <- sin(0.06 + 0.5*pi*t) + sin(0.1 - 0.68*pi*t) - sin(0.04 - 0.74*pi*t) + sin(0.03 - 0.93*pi*t)

measurement <- signal + noise

xkm1km1 <- matrix(c(4.4, 0),2,1)
Pkm1km1 <- matrix(c(1,0,0,1),2,2)
H <- matrix(c(1,0),1,2)
A <- matrix(c(1,1,0,1),2,2)
Q <- 10^-6 * matrix(c(1.6,2.4,2.4,4.8),2,2)
R <- 10^-6 * matrix(c(2.35),1,1)

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

zhat <- t*0

for (k in 1:T)
{
  xkkm1 <- A %*% xkm1km1
  Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
  K <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  z <- matrix(c(measurement[k]), 1, 1)
  xkm1km1 <- xkkm1 + K %*% (z - H %*% xkkm1)
  Pkm1km1 <- (matrix(c(1,0,0,1),2,2) - K %*% H) %*% Pkkm1
  zhat[k] <- as.numeric(H %*% xkkm1)
}

plot(t, measurement, type="b", pch=19)
lines(t, signal, col="red",  lwd=10)
lines(t, zhat, col="blue", lwd=5)

lpf <-  filter(measurement, c(0.2,0.2,0.2,0.2,0.2), circular = TRUE)
lines(t,lpf, col="green", lwd=4)

errs <- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
print(errs)


Non Uniform FFT / DFT with FFTW? Can FFTW perform the FFT on data which is not uniform in time? I can't seem to find a straight-forward answer to this question anywhere online.

The answer is Non-uniform discrete Fourier transform

I suggest you to take a look in here: NFFT library. Tutorial for that purpose: NFFT 3.0 Tutorial.

You can also find a Python wrapper: pyNFFT.

  The generalisations of the NFFT include

  NNFFT - nonequispaced in time and frequency fast Fourier transform,

  NFCT/NFST - nonequispaced fast (co)sine transform,

  NSFFT - nonequispaced sparse fast Fourier transform,

  FPT - fast polynomial transform,

  NFSFT - nonequispaced fast spherical Fourier transform,

  NFSOFT - nonequispaced fast Fourier transform on the rotation group

  The NFFT is a C subroutine library for computing the nonequispaced
  discrete Fourier transform (NDFT) in one or more dimensions, of
  arbitrary input size, and of complex data. New: A Matlab interface is
  part of the NFFT3. We believe that our library, which is free
  software, and based on  FFTW (FFTW 3.x) should become the NFFT library
  of choice for most applications.


Phase locking value / phase synchronization I have pairs of EEG signals filtered to give me just alpha bands. For a machine learning problem I have found that calculating the phase locking value between the alpha signals of these 2 channels would be a useful classification feature.

I'm using this paper as a reference: http://www.ncbi.nlm.nih.gov/pubmed/22661936

I'm wondering what the correct method/implementation is to calculate PLV/synchrony in Python?

I have not been able to find any modules/code that would help me with this problem. That paper, for example, calculates synchrony between independent components after ICA, but I'm not employing that method here so not sure if their formula translates or if any adjustments are necessary for filtered signals

The closest thing I have found is this: https://www.nbtwiki.net/doku.php?id=tutorial:phase_locking_value#.VchjR6ZVhBc

However, while the code is laid out nearer the bottom of the page, it appears to be MATLAB code. It looks like there are some MATLAB methods being employed and I'm not sure how that would translate over to Python/numpy/scipy

I'm not particularly familiar with this domain and don't trust myself to port the code over without making logic errors, so it would be useful if someone could point me toward a Python implementation of this.

There can be several ways to calculate the Phase locking value (PLV). For relatively mono-component and high SNR (well filtered)-Time domain signal can be converted into analytical signal using Hilbert transform to calculate the phase difference. For the right signal it is a very powerful technique as is shown in the tutorial you have referenced.
Here is a code in python :
y1 and y2 are the input signals. The function gives average phase difference between the two signals even when they are not relatively well behaved. If you run it in a loop, it can give you phase difference in parts or even sample by sample (instantaneous phase).

import numpy as np
import scipy.signal as sig

def hilphase(y1,y2):
    sig1_hill=sig.hilbert(y1)
    sig2_hill=sig.hilbert(y2)
    pdt=(np.inner(sig1_hill,np.conj(sig2_hill))/(np.sqrt(np.inner(sig1_hill,
               np.conj(sig1_hill))*np.inner(sig2_hill,np.conj(sig2_hill)))))
    phase = np.angle(pdt)

return phase

For well behaved signal as shown in the tutorial even a more simple instantaneous phase program can be easily written :

import numpy as np
import scipy.signal as sig

def hilphase(y1,y2):
    sig1_hill=sig.hilbert(y1)
    sig2_hill=sig.hilbert(y2)
    phase_y1=np.unwrap(np.angle(sig1_hill))
    phase_y2=np.unwrap(np.angle(sig2_hill))
    Inst_phase_diff=phase_y1-phase_y2
    avg_phase=np.average(Inst_phase_diff)
return Inst_phase_diff,avg_phase


WAV encoding problem When processing wav files, I encounter this problem: the what I extract PCM from a 8-bit encoded .wav file, I got a sequence of integers. However, when verifying my implementation with wavread in Matlab, the PCM code is represented with so-called double-precision normalized samples. My question is what's the difference between these two format. For example, the hex number is 0x80 0x81, the corresponding double-precision normalized samples are 0, 0.0078. IIRC, PCM code in .wav file records sound pressure value of audio. Is double-precision normalized samples sound pressure level? If so the value should be $10log(P_1/P_2)$ (based 10). But $10log(129/128)$is not 0.078 instead $10ln(129/128)$(based $e$) is. Anything wrong with my intuition?

Thanks and best regards.

The PCM values are proportional to the sound pressure. 8-bit PCM is stored with an offset of 128, so PCM coded values from 0 to 255 code for 'pressure values' from -128 to 127.
These values are often normalized to the range -1.0..+1.0.
The PCM values are linear representation of the sound pressure, and a log operation is only required if you want to express the sound pressure in decibels.
To get actual sound pressure levels in Pascal or dB(A) you would need a calibrated .wav signal.


Why are linear phase filters called so, if they provide Constant delay instead of linear Linear phase filters delay all frequencies by the same amount. Why aren't they called Constant phase filters instead of Linear phase?

As I understand, if there is an input signal with two components $f_1$ and $f_2$, and it is passed through a linear phase system which introduces a delay $t$, $f_1$ is delayed by $t$, and $f_2$ is also delayed by $t$.

To what property does the phase have a linear relationship with? Could someone provide some insight into this? Some Math would be helpful.

Your premise is wrong. Linear phase filters offer linear phase.

Phase is not the same as time delay – at $f_1$, time delay $t$ will lead to a phase of $\varphi_1=-\omega_1\cdot t$, whereas at $f_2$, the phase will be $\varphi_2=-\omega_2\cdot t$ ($\omega$ is $2\pi f$); as you can immediately see, phase is a linear function of frequency.

(It's always good to remind oneself that frequency is just the derivative of phase over time, $\omega=\frac{d\varphi}{dt}$, and thus phase is just an integral of frequency over time.)


Where does this star came from? [EDIT - Detailed Question] First of all, it's a rare topic in google, so I can't find intuitive explanation about step-invariant, so I don't understand actually what it is, except to solve zero order hold transfer function.

In this pdf, page 15 to 18 I've been learning about zero order hold (ZOH).
http://www.engr.usask.ca/classes/EE/480/notes/myee480-p4-sampleHold.pdf

As you can see above, suddenly the laplace transform has become starred, why does it happened?

EDIT

First thought, I think it was like distributive starred property. But why suddenly it has a star while you still have starred curly bracket? Is it a typo?

And also, from the wikipedia
https://en.wikipedia.org/wiki/Starred_transform
It's written that $X^*(s) = L[x(t).dirac_T(t)] = L[x^*(t)]$
The star is only on the big X and small x, and it's not on the L one, so why the creator made it star on the L one?

Furthermore, if you also take a look at page 17 (top right). First, (1-e^-Ts) are converted into laplace transformation form, and then it gets starred, and then it is separated, and then it return to the first original (1-e^-Ts) form. What's the use? It's seems illogical.

If you don't mine, please also explain intuitively what is step-invariant z-transform and why the creator use step-invariant z-transform method to find the transfer function and not the other method.

This seems as a starred transformation. It is a discrete version of the Laplace transformation. Laplace transformation is continuous by nature. Starred transformation is discrete, and the star denotes difference between continuous and analog domain.

This is commonly used when teaching about the Z transform. It tends to be confusing, so the teachers typically just skim over it.

Here are some explanation about it:  https://en.wikipedia.org/wiki/Starred_transform

www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=934


FFT plot in GNU Radio i got similar question as in FFT scaling in GNU Radio
and quite not clear about the answer given so far. I understood that the FFT plot is a relative value in which we have to use measurement instrument to get the actual value.

But, without the measurement instrument, in general, what value will be shown on the plot if the given input is for example complex cosine waveform with amplitude 1?

I've tried the following flowgraph: signal_source - throttle - qt_frequency_sink

and vary the cosine waveform amplitude from 0 to 100. other than amplitude, all blocks use default setting. I expect to get 0dB when the amplitude is 1, but the plot shows ~ -10dB instead,and ~32dB when amplitude is 100 (i'm not sure if this value make sense). SO, how actually the frequency sink calculates these value?

But, without the measurement instrument, in general, what value will be shown on the plot if the given input is for example complex cosine waveform with amplitude 1?

You've done the "right" thing by simply trying.

Realize that the frequency sink simply does a forward FFT and windowing, if you don't set the window to rectangular; the window you're using has a response of ca -9.8dB at DC, so that's what you're seeing here. Use "FFT Window" none, and you'll see your perfect 0dB for a constant 1 input, and 40dB for a constant 100 input.

Since you only have a constant value, so all bins are 0, but the middle one, which will, as defined for the DFT without scaling,

$$\begin{align}
 Y[f] &= \sum_{k=0}^{N_{DFT}-1} y[k] e^{j2\pi f \frac kN},\\
\text{for }f&=0\text{, it follows that }
e^{j2\pi 0 \frac kN}\equiv 1,\text{ so}\\
Y[0] &= \sum_{k=0}^{N_{DFT}-1} y[k] &y\text{ being constant}\\
&= N_{DFT}\cdot y[0]
\end{align}$$

Because it's handy for visualization, GR's approach to scaling of these values is to divide them by the DFT length $N_{DFT}$; otherwise, the same constant signal would have more power in a longer DFT.


Object Tracing with Particular Color I started working on a autonomous robot project in my spare time. Everything is finished except the image proccessing part. I want my robot to follow an object with particular properties. For the start, I decided to use color as the particular property. For example if object's color is yellow, then I want my robot to chase it. I am using OpenCV and tried some of the sample codes but the problem is that in real world lightning is not uniform and the sample codes which I tried did not do a good job because of the changing threshold.

So what can I do? What do you prefer to me? (Are there any better sample codes, algorithm suggestions etc?) I will appreciate if you explain your suggestions as if you were explaining it to someone with nearly none experience with image processing.

Regards,

Amadeus

As you're probably starting to realise this is a very big field with a number of different methods which you can try.

If you want to try to use colour, you can try transforming from RGB to HSV space, which attempts to decouple the colour from the "lightness", this may improve your lighting inconsistencies.

The next problem is how to track your object. If you know the object's colour distribution (histogram) then you can use backprojection to assign a likelihood to each pixel being that colour, then perhaps you could simply track the centre of mass of the back projection.

Alternatively, you could look into template matching, where a template is scanned across the image and areas which match this template closely are scored more highly than non-matches, giving you likely areas for where the object is.

If you have any more specific information about what you are trying to acomplish, it will be easier to help you: do you have an example input image, with targeted highlighted?


Algorithm to find spike I've been doing FFT's (on sound) for a while and usually use the highest intensity or the ones above average to get the spikes in an FFT result.

What would be a good method to get all the 'spikes' in an FFT. I can't use averages or peaks because the frequency response drops off. Visually its quite obvious there is a spike but I can't get a way extract the frequency from the spike.

Low frequencies always tend to have a higher itensity than higher ones but don't necessarily have any spikes in them.

Are there any good reliable techniques for this?

A statistical test might be appropriate.  For instance, if you assume the noise floor is Gaussian, then you could compute the standard deviation in any particular octave (or other appropriate sub-band), and look for maxima that are whatever multiple of standard deviations above the local average that you might consider to be statistically significant for your application.


Linear Difference Equation and Method of Least Squares I'm reading the book "Fault-Diagnosis Systems" by Isermann in the par. 9.2.1a.
The author explains how to estimate the parameter of a linear difference equation using Least Squares.
We start with a standard ordinary linear equation with an error term:

the equation is called 9.15.
Then,

The first thing i don't understand is the syntax used by the author. What's the meaning of the | inside $\hat{y}(k|k-1)$ ?

ps: I'm sorry if what i ask is trivial, but when i have a kind of doubt like this, i don't know how to solve it if not asking it in this site...

The notation $ \hat{y} \left( k \mid k - 1 \right) $ usually means this is an estimated value of $ y \left( k \right) $ given all the available data up to time index $ k - 1 $.
So generally speaking, this is a prediction of one step in time of the data.

The case above also suggests linear estimation.
Namely, $ \hat{y} \left( k \right) $ is built using linear combination of all the given data.


Questions in a Character Degradation Model Paper I have questions in the ‘Noise region definition’ and ‘Noise generation process’ of the paper “A character degradation model for grayscale ancient document images”.

In Noise region definition, g controls the flatness of the regions. What does it exactly mean? How can we say that a noise region is flatter compared to another noise region?

Below is the illustration of the ellipse noise region within the document image. The green ellipse shape is defined to be the noise region.

According to paragraph below, the average value bj of its 8-neighbours (in the initial grayscale image) is used for calculating values of all pixels in the line CiBj.

Does the average value of bj  calculated by averaging the greyscale value of the adjacent pixels in north, northeast, east, south east, south, south west, west and north west?

Does the formula for getting Pk is below.

Please refer to the paragraph below:

I just want to assess , if my comprehension in reading the article is right. Thanks.

All g is doing is saying whether the area is closer to a circle (less flat) or a line (completely flat ellipse).
Yes.
No.  $N(\mu, \sigma^2)$ is a distribution.  You can say
$$P_k \tilde{} N(\mu, \sigma^2)$$
read as $P_k$ is distributed normal, with mean $\mu$ and variance $\sigma^2$.

To find a value for $P_k$ you need to generate a random number with the mean given by equation (1) in the last screenshot:

$$
\mu = \bar{c}_i + (\bar{b}_j - \bar{c}_i) \left( \frac{d_{ik}}{d_{ij}} \right)
$$

In matlab, this would be something like normrnd(mu, sigma), though what sigma value to choose is not clear from your extractions from the paper.


Mel separated filter banks I am currently a bit confused.
I am currently trying to understand the MFCC extraction process, and don't understand how the Mel filter banks only provide a scalar output from each filter bank. What does this scalar value mean? I know the filters aren't infinitely narrow, thereby only allowing one frequency so what about the frequencies around the center frequency?

Scalar Value from the output of each mel filter bank is accumulated scaled magnitudes of frequencies present in mel filters range/bandwidth.


Incorrect dominant frequency from FFT DATA source: Bearing vibration data FEMTO IEEE PHM 2012 challenge.

DATA context:

Sample rate: 25.6 kHz
Sample period: 0.1 Seconds (2560 points)
Motor runs at 1800 rpm

Steps taken:

Low pass butter filter with fc=1000Hz
Apply hann window
Apply FFT (np.fft.fft)

Problem:

The dominant frequency is at ~60Hz which is incorrect it is supposed to be nearer to 30Hz since the motor is ran at 1800RPM.(Correct me if im wrong on this.) What went wrong? Please find the code below.

Fs = 25600  # Hz
sampling_period = 0.1
N = int(Fs * sampling_period)

fc = 1000  # Cut-off frequency of the filter
w = fc / (Fs / 2)  # Normalize the frequency
b, a = signal.butter(10, w)
filter_out = signal.filtfilt(b, a, bearing1_1["horiz_acc"])

f = np.fft.fftfreq(int(N), 1/Fs)
window = np.hanning(int(N))
X = 1/N * np.abs(np.fft.fft(filter_out * window)) * 2  # Multiply by 2 for hann window power correction?

Thanks in advance!

What went wrong?

How do you know it's wrong? Could be that you have way more energy in the first harmonic than at the fundamental.
Debug rule #1: test the code first with a signal where you already know the answer (say a mixture of a few sine waves). Once you have made sure that your code produces the expected result EXACTLY you can move on to real world signals.
Your frequency resolution seems insufficient. You only have 10 Hz of resolution, that will not give a lot of detail at 30 Hz. Keep in mind that the energy at each particular frequency will be "smeared out" across a few neighboring FFT bins, so you want enough resolution to see that.
Your lowpass seems way to steep. Together with filtfilt you get a 20th order filter which will have enormous amount of time domain ringing which filtfilt() unfortunately truncates and the truncation can create significant errors


Is inverse matched filter useful? I am aware of matched filter and its application. Now, wondering if there is any application of inverse matched filter? What I mean my inverse matched filter is that convolution of matched filter and the inverse matched filter would lead to close to delta function.

In your question, you postulate the existence of two filters, let's say $p(t)$ and $g(t)$, such that $p(t) \star g(t)=\delta(t)$. This assumption is problematic because it implies that $G(f)=1/P(f)$, which requires $P(f) \neq 0$ for all $f$. Furthermore, for small values of $P(f)$, $G(f)$ will take arbitrarily large values (see Dilip Sarwate's comments above and below this answer).

Assuming (for the sake of discussion) the existance of such a pair of filters, there may not be any advantage to using them. I'll describe a possible application from digital communications. Let's say you want to transmit a number $A$ over an analog channel. You choose an appropriate analog pulse $p(t)$ and transmit $$s(t) = A\delta(t)\star p(t).$$ (The pulse $p(t)$ could be chosen to fit the channel bandwidth, and/or to have a certain energy). The receiver could perform the "inverse filter" operation on the signal $s(t)$: $$s(t) \star g(t)=A\delta(t),$$ where $g(t)$ is your "inverse matched filter". However, in practice it turns out that detecting $A$ in this way is not optimal, because it does not have the best signal-to-noise ratio. In other words, when you add noise to the transmitted signal, the received signal becomes $$r(t)=s(t)+n(t).$$ In this case, filtering $r(t)$ with the "inverse matched filter" $g(t)$ results in a worse signal-to-noise ratio than filtering with the filter matched to $p(t)$.

If you have access to it, I highly recommend "An Introduction to Matched Filters", by G. Turin, IRE Transactions on Information Theory, June 1960. Your question is answered in page 318.


Filter to speed up convergence of geometric series I want to speed up convergence of geometric series
$\ \  lim + c^n \ \ $ where $|c| < 1$ and $lim$ is not known.
A basic low-pass filter, here Butterworth, is good for $-1 < c < 0$ but not $0 < c < 1$:

butter  N 2  Wn 0.5:  [.293 .586 .293] / [1 0 .172]
|response|: [100 100  99  97  88  71  47  25  10   3   0]
(-.8) ^ n  [100 -80  64 -51  41 -33  26 -21  17 -13]  ->  [29 35 -4 -7  1  1  0  0  0  0]
.8 ^ n     [100  80  64  51  41  33  26  21  17  13]  ->  [29 82 90 62 45 38 31 25 20 16]

Is there a filter that does both ?

Why: many iterative processes that generate points $X_n$ in 1d, 3d, 100d
converge to an $X_{lim}$ slowly,
but roughly geometrically:
either one-sided like $.8^n$, or oscillating like $(-.8)^n$ .
It would be nice to have a simple, robust filter that speeds up both.

Shanks transformation
is a nonlinear transform that works really well on geometric series.
But on linear input 1 2 3 it'll divide by 0.
(Check for that, push on regardless ? practical, but unsatisfying.)
(I believe that
Aitken's delta-squared process
is the same as Shanks -- not sure.)

The problem is that in the case with positive $c$ you have a DC desired signal ($lim$) and a low frequency noise ($c^n$). In the case when $c$ is negative you have a high frequency noise signal, which is of course much easier to separate from the desired DC value.

If you figure out some way to estimate $c$ (sign and magnitude) you can use a simple FIR filter

b=[1,-c]/(1-c); a=1;

which will exactly eliminate the exponentially decaying noise:

$$y[n]=\frac{1}{1-c}(x[n]-cx[n-1])=\frac{1}{1-c}(lim+c^n-c\cdot lim - c\cdot c^{n-1})=lim$$


Preprocessing audio signal for neural network classification I need to identify certain features of the audio signal recorded from microphone in stethoscope. These sounds are only samples i've found, but the final signal will be probably a bit noisier (maybe not, i don't know yet). Below are examples of three signals. Normal one and one of the person sick with atrial spatial defect, third one is of person with late aortic stenosis. These two fluctuations are named S1 and S2 respectively and they are my main subject of interest. Depending on a given heart condition S1 may have doubling and there may be some "noisy" sounds (clicks also) between S1 and S2. The length of the whole cycle (S1, S2) may vary depending on the heart rate.

My question is, what kind of processing should i perform to surely extract and recognize those features, i.e. two large fluctuations normally, with slight delay between them and if person is sick it may be either something before S1 or between S1 and S2 or S1 may double a bit. Also normally S1 has a larger amplitude spike than S2 as i understand.

I have very limited knowledge in this field and not sure about what options are available to me. FFT magnitude graphs are showing somewhat recognizable results, but i think there is too much unneeded information there and they are not so distinct.

The samples themselves can be found here: Washington University: Heart sound samples

Normal heart sound:

Atrial spatial defect heart sound:

Late aortic stenosis:

Some tips that might help:

Your FFT should be shorter than the features you are trying to detect.  An FFT will be calculated for a specific number of samples.  This represents a certain length of time (1024 samples at 1kHz sampling rate is 1024 milliseconds.)  If the time period of your FFT is longer than your signal features then they will all be smeared together - you won't be able to identify the features.  You will have to find a compromise for your FFT length between enough frequency resolution and enough time resolution.  The needed frequency resolution depends on the frequency content of features you need to detect, time resolution depndes on the length of the features (as well as the length of time when there are no features - which is itself perhaps a feature you need to take note of.)
Classify the features as to frequency range.  Make an average of the FFT-points for those ranges - this will give you the intensity of that feature.
Plot the intensity of the features against time.
Locate groupings of feature intensities that are unique for each type of heart defect - these are the things you need to be able to recognize.  For example: Pulmonic stenosis has low frequency (normal heart beat sound) with high frequency (murmur) at the same time, but no murmur type sound between beats - but it also has a click that you will need to detect.
You will need to detect the S1 and S2 by the pause between beats - long pause followed by short means that the last beat is the S2.
The examples on the site you referenced have been cleaned up, but they still have things that might try to sneak in as features.  I hear a 60Hz power hum, and some ringing from the filters that were used to clean things up.  There is also some sone background "white" noise in some of them.

I would try loading the files into Baudline and seeing how to identify the features and their combinations by eye first.  Baudline can show you how the freqeuncy content varies with time so that you can see the relationships between features and relative timing.  Use Audacity to convert the mp3 files to wav and to downsample to 1000Hz sampling rate before loading in Baudline.

The Baudline Color Aperture is handy for making certain parts of the diagram standout - this can make it easier to identify features.

Baudline is for Unix type systems only (Linux, BSD, Mac OSX) but I'm sure there's similar software for Windows.  Audacity is available for Windows as well as Linux and OSX.


STFT amplitude normalization, librosa library So librosa.core.stft returns a complex single sided spectrogram.

My question is:

What normalization of the amplitude values should I perform afterwards?

I believe I have to multiply the amplitude outputs by 2 in order to preserve the energy that was assignated to the negative frequencies. I have not seen in the code that librosa.core.stft does this, although I might have missed it.

And then the other normalization step would be to divide the outputs by the length of the window size right?

However, after doing these two amplitude normalization steps, I still lack a factor of 2 somewhere.
This can be seen in the picture below, where in the coloured picture the spectrogram has been calculated and an amplitude value of 0.5 is observed. And in the second plot, the FT of the whole signal
has been computed and, applying the same normalization, the expected value of 1 is observed.

Picture: The frequency axis are not in Hz, the only important thing here are the amplitude values.

Here is the code:

%Import packages
import scipy.fftpack as fft
import librosa
import numpy as np

%Create sine wave of frequency 1000 Hz and 20000 samples
frequency = 1000
fs = 4 * frequency
time_vector = 1/fs * np.arange(20000)

sine_wave = np.sin(2 * np.pi * frequency * time_vector)

window_size = 1024
hop_length = 512

out  = librosa.core.spectrum.stft(sine_wave, n_fft = window_size, hop_length = hop_length)

%Take magnitude and normalize
out = 2 * np.abs(out) / window_size

fig = plt.figure()
axes = fig.add_subplot(111)
im = axes.pcolormesh(out)
fig.colorbar(im)
plt.show()

%And plot also the fourier transform of the whole signal
transform = 2 * np.abs(fft.fft(sine_wave)) / len(sine_wave)

%Plot only positive frequencies
plt.plot(transform[10000:])

I did answer a similar question a few years back, but can't find it. Basically, you are losing the energy because of the windowing. It's true to say that you should multiply the spectra by 2 in order to retrieve the energy from negative frequencies and divide it by the window size. However, this is only true for the rectangular window. The correct normalization should be following:

out = 2 * np.abs(out) / np.sum(window)

where np.sum(window) in case of rectangular (or no window in fact), is equal to the frame size. It kind of makes sense, right? If window attenuates the signal at its ends, then this energy should be recovered somehow.

Now back to your problem, librosa.stft by default uses hann window which looks like:

If you calculate the np.sum(np.hanning(window_size)) you will get $511.5$. If you read up more, the coherent gain of that window is approximately $0.5$. Which is defined as:

$$CG = \frac{\Sigma w(N)}{N} $$

So in order to solve your problem, I suggest to explicitly pass the window vector to the stft function:

import librosa
import numpy as np
import matplotlib.pyplot as plt

frequency = 1000
fs = 4 * frequency
time_vector = 1/fs * np.arange(20000)

sine_wave = np.sin(2 * np.pi * frequency * time_vector)

window_size = 1024
window = np.hanning(window_size)
hop_length = 512

out  = librosa.core.spectrum.stft(sine_wave, n_fft = window_size, hop_length = hop_length, window=window)
out = 2 * np.abs(out) / np.sum(window)

fig = plt.figure()
axes = fig.add_subplot(111)
im = axes.pcolormesh(out)
fig.colorbar(im)

plt.figure()
slice = out[:, 10]
plt.plot(slice)
plt.grid()

plt.show()

Which will return:


Entropy of enhanced image more than the original image It is said that the entropy of an image is the information that it contains. Let us assume, that the entropy of an image f(x,y) is E1. If we enhance f(x,y), we get an enhanced image g(x,y), whose entropy is let us say E2. Now, is $E2 \le E1$ always? According to me, the above relation should always hold true, as how can an enhanced image contain more information than the original image. Also, it is just the count value at each bin of the histogram that changes when a histogram is equalized using certain technique but overall the total count value of all the bins remains same. So $E1 \ge E2$ should always hold true but I have just seent that E2 will be greater than E1 if we enhance the image using CLAHE technique i.e the entropy of the enhanced image will be greater than the original image. I have searched for the explanation and got some interesting links such as
Entropy of an image but none clears my doubt. Kindly give your inputs and explain me where am I getting it wrong. Thanks.

Image used : https://1drv.ms/u/s!ArfYMunvE-HGgbQ5yDBYTj6UfNdWig

According to me, the above relation should always hold true,

According to the data processing theorem (and here), this is always true.

  how can an enhanced image contain more information than the original image [?]

It cannot. A trivial example is that of interpolation. Suppose that you upsample an image by a factor of 2. This means, inserting a "new" pixel between two known ones. The information remains the same though because no "new light" has been captured. We still do not know what exactly happened BETWEEN the known pixels. All we can do effectively is INTERpolate.

Interpolation works out the values of the unknown pixels based on the values of the KNOWN pixels, therefore re-distributing information across a larger surface. Since all processing takes place over the CAPTURED data it cannot but be inferred that a data processing step can only DECREASE the information content of a dataset.

  I have just seent that E2 will be greater than E1 if we enhance the image using CLAHE technique i.e the entropy of the enhanced image will be greater than the original image.

First of all, you are comparing two entropy measures that appear to have no relationship between them, no shared context.

Having said this, is there any imaginable case in which $E2 > E1$?

Consider the processing step as a black box. It accepts an image ($I$) and it produces an image ($O$), only that in this case, all that this black box does is $O = \alpha \cdot I + (1-\alpha) \cdot N(I)$, where $N(I)$ is a noise process that produces a two dimensional noise image with the same dimensions and overall brightness as $I$. The parameter $\alpha$ determines how much image and how much noise is mixed in the output.

So, at $\alpha = 1$, $O$ is essentially exactly the same as $I$. In this case $E(O) \le E(I)$. BUT, at $\alpha = 0$, $O$ is simply a noise image. An incoherent set of pixels whose entropy is very high. Since there is absolutely no relationship between $I$ and $N$, then it can well be that $E(O) > E(I)$ but this does not mean that anything NEW (AND meaningful) was added to the image.

In the place of $N$ we could also have had extraneous information. For example, image from one camera being cross referenced with an image (of the same scene) from another camera.

But, there is also something else that is at play in your case and this is a "simplistic" view of entropy.

Before we move on to CLAHE, here is a simple experiment. I am using Octave, you can try it in MATLAB or anything else with a bit of re-writing:

n=1500;nBin=16;maxNum=32;
a=round(rand(1,n)*maxNum);
Q = hist(a,nBin)./sum(hist(a,nBin));
H = -sum(Q.*log2(Q))

Where n is the length of a random sequence, nBin is the number of histogram bins, a is a quantised sequence of random numbers up to maxNum, Q is the normalised marginal probability of a and H is the calculated entropy of a.

Now, for a given n,maxNum produce a sequence a and then calculate its Shannon entropy for two different nBin, say for instance nBin=16 and nBin=32.

To recap, SAME sequence of numbers a (absolutely no processing), two DIFFERENT histogram resolutions. What are the two different H?

About 4 for nBin=16 and about 5 for nBin=32. Isn't this a bit suspicious, given that $2^4=16$ and $2^5=32$?

Shannon's entropy gives you the upper bound of the entropy of a source and the way it is derived here, assumes that all the characteristics of the source are captured in Q. Function rand returns numbers from a uniform distribution. If you divide the sequence's histogram at 2 levels, you get 2 equiprobable events and still need 1 bit to describe them (try it, for nBin=2). If you divide the histogram of the SAME sequence at 4 levels, you now get 4 EQUIPROBABLE events for which you need 4 bits to describe them. At 5 levels of the SAME sequence you now get.....and so on.

In this example, a remained the same but the expression of entropy changed, purely as a byproduct of how we calculated it.

This way of calculating the entropy suffers from another problem. That of "empty" bins. If you get a bit more adventurous in the above example and set nBin=500 (for example), you will notice that your H might come back as NaN. Once the histogram contains zero bins, log2(0)*0 produces NaN. Again, we get a result which is purely based on the way entropy is calculated. You might say, "I will omit the zero entries" but then again, you might be omitting them in one image because they did not happen but you would not be able to ignore them in another image because they DID happen.

And THIS is what you see in your case, as a byproduct of how adaptive equalisation works.

Your initial $I$ has a given histogram "shape", a given distribution of grey tones. Most likely this is concentrated around a lobe inside a limited region at the low end. Its entropy will appear low. Yes, you don't need many bits to describe this lobe. AFTER equalisation, the algorithm has now redistributed grey values all over the place, therefore it results in a WIDER histogram (with also differences in its overall shape, with tones that now are not 'negligible') appearing to have a higher entropy. Appearing to require MORE bits to describe the "new" states.

Shannon's entropy is not the only way to calculate entropy. There are other definitions as well, such as those that work over the spectrum of the image or based on data compression.


Designing Rayleigh Fading Wireless Channel with particular properties I have some information on particular wireless channel between an Access Point(AP) and K users as described in the picture below:

How i can create this channel in Matlab?


Does whitening with PCA correspond to power spectrum flattening? I'm finally making some headway in understanding data whitening, but I have a question now relating what I think of as "data-based" whitening to the situation where there are time-varying signals.

So in the "data-based" world, suppose I have a mean-centered column matrix $X$ with covariance $\Sigma = \mathbb{E}[XX^\top]$. This data matrix can be transformed into a whitened dataset $W$ by computing $W=\Lambda^{-1/2}EX$, where $\Lambda$ and $E$ are the eigenvalues and eigenvectors of $\Sigma$.

So far, so good. What I'm curious about is, what if my data vectors are a bunch of windowed snippets of audio (for example) ? I could use the DFT to compute a mean power spectrum for these snippets, and then rescale each snippet by the inverse of the spectrum at each frequency. If I understand correctly, this would whiten the power spectrum of my dataset.

Does this process equate to whitening in the "data-based" world above ?

If so, are the eigenvalues of the data covariance (which I understand are also referred to as a spectrum) directly related to the mean power spectrum of my signals ? How ? Any pointers to good reading on the subject would be very welcome.

Sorry if this is an elementary question, but I haven't had much signal processing background and find this to be really fascinating.

I believe the answer here is no.

Remember what 'whitening' is really doing: It is ensuring that the power/variance that correspond to each of your basis vectors are the same. Which basis vectors? The ones you chose.

In the DFT, you are choosing basis vectors comprised of complex exponential functions. In PCA, you are determining the basis vectors empirically. What is the relationship between complex exponential bases and empirical bases? You cannot say.

Thus, when you say you whitened in the spectral domain, you are really saying "I made sure the power/variance associated with all my complex exponential basis vectors representing my data are the same".

When you say you whitened in the PCA domain, you are really saying "I made sure the power/variance associated with of all my empirically derived basis vectors are the same".

Are your empirically derived basis vectors complex exponentials? There is no guarantee. Thus, whitening in either domain is not equivalent to whitening in the other.

  If so, are the eigenvalues of the data covariance (which I understand
  are also referred to as a spectrum) directly related to the mean power
  spectrum of my signals ? How ?

If I understand your question correctly, then the answer to this is a simple yes, for the simple reason of energy conservation. All the eigenvalues together will give you the total energy/variance in your data. The average of this energy, will have to correspond to your mean power in any domain, including the spectral.

Projecting your data onto either the DFT or PCA (empirical) bases dones not add/remove energy, but simply redistributes it. Thus, your average power has to always be the same.


MATLAB: Generating a random bit stream in a frequency Band for PSK Modulation I want to generate a pseudo-random bit stream in MATLAB in a given frequency band for PSK Modulation. I am using the idinput function to generate the bit stream:

u = idinput(2^17,'prbs',[0.8 1],[0 1]);
figure(1);
[psd2,w] = pwelch(u,blackmanharris(2^17/8),[],[],1e9,'onesided');
semilogx(w,10*log10(psd2/2));grid on;

I want to use minimum a bit rate of 0.8GHz and want the bit steam in the frequency band from 0.8GHz to 1GHz. But when I check the Power Spectral Density, I see that it is flat for the whole band until 1GHz (with a DC component coming).

What is wrong with this code? How can I make the bit stream restrict to a given frequency band?
I want to have this frequency band to be within a modulation bandwidth for the PSK transmitter.

You are attempting to simulate an RF signal directly at it’s carrier frequency which will require much more processing so is not recommended. You can achieve the same results by simulating everything at its equivalent complex baseband signal due to the linear properties of the underlying frequency translation. Equivalently everything in your simulation model (IF bands, channel models etc) can/should all be at complex baseband equivalent representation (the analytic signal).

That said, to do what you are trying to accomplish with either method, you would first create the complex modulation signal mapping your data bits to PSK symbols; for example for 8 PSK 3 bits map to every complex signal such as 000 could become $e^{j\pi/4}$, 111 could become $e^{j5\pi/4}$, etc- any mapping you choose but typically “grey coded” such that adjacent symbol locations only change by one bit.

This sequence which is then sampled only 1 sample per symbol then needs to be up sampled to support implementing pulse shaping (2 samples per symbol would be sufficient but in implementation more samples here can simplify later filtering).

This now represents the complex baseband signal with pulse shaping and if 2 samples per symbol the represented spectrum from a PSD would extend from $-f_S/2$ to $+f_S/2$ where $f_S$ is the symbol rate. (Thus If oversampling more this frequency range can be extended further).

This shows the SIGNIFICANT advantage to doing all simulations at baseband - to now continue to to simulating at a carrier frequency at 900 MHz directly you would need to increase the sampling rate to be greater than the carrier frequency plus half the signal bandwidth plus some margin if completely a complex signal or twice that if simulating a real signal such that you want to model each carrier cycle. There is no reason to do this since everything that carrier is doing can be represented at complex baseband.

But if you still wanted to, against my advice, you would first interpolate the baseband waveform to this much higher rate and then multiply the signal by the carrier frequency ($e^{j\omega_c t}$ if complex or $\cos(\omega t)$ if modeling a real carrier).


Simple software low pass filter I am measuring voltages and currents of a 3-phase electrical machine and I need to calculate the power. Every interrupt (frequency between 30-70kHz) I get values (voltages, currents) from analogue-digital converter and I need to do a simple calculation to determine power. However, I was told to use averaging with a first order filter with a time constant of approximate 1s for the power calculation. This, I believe, means low pass filter with a frequency of approximately 1Hz? What is the simplest solution to this? Moving average?

A first order lowpass filter is usually implemented like this:

$$p[n] = \alpha p[n-1] + (1-\alpha) pi[n]$$

Where $p[n]$ is your filtered power estimation, $p[n-1]$ is the previous result, $pi[n]$ is your new measurement (probably the product of instantaneous voltage and current measurements), and $\alpha$ is a positive parameter just less than 1.

The nearer $\alpha$ gets to 1, the larger the time constant (lower cutoff frequency) of your filter. But beware, especially in embedded systems with limited precision, that getting too near to 1 can make your filter unstable, or at least have problems due to numerical precision.

The cutoff frequency for that filter is around $f_s \frac{1-\alpha}{2\pi\alpha}$, where $f_s$ is your sampling frequency.


Square wave function in FFT spectrum The input function is a 60Hz square wave. From Fourier Transform I know that the frequency on the spectrum is 60,180,300....(2k-1). But there's a series of 120,240,360...(2k) signals appear.
I've googled but all I found is Fourier Transform issues, not Fast Fourier Transform. I'm wondering why there are some 2k terms appearing on the spectrum.
Thanks!

Looks like aliasing.

Your 360 Hz is probably not a real component at 360 Hz but it's the 9th harmonic (540 Hz) that aliases back to 360 Hz at a 900Hz sample rate.

That's easy to test: change the sample rate to, 940 Hz and see of the 360Hz moves to 400Hz or stays put. It moves to 400Hz, it's clearly aliasing.


Frequency dependent phase difference between signals arriving at two microphones I wish to calculate the time difference between a signal (may or may not be periodic), arriving at the right and left channel of a stereo microphone pair kept at a distance of about 30 cm. I am doing so to localize the acoustic source.

I compute the FFT of the left channel and the right channel, and then find the values of the FFT at a particular frequency bin of interest. I then divide the (complex) values to get the phase difference in radians. I then convert the phase difference to a time difference ( time difference = phase difference / 2* pi ) and use this to calculate the angle of arrival of the acoustic source.

What are the ambiguities that could occur in this method of calculating the time difference, as compared to a cross-correlation method of finding time difference?

I am using this method rather than cross-correlation due to the frequency dependence of the phase difference.

For a constant time difference, the unwrapped phase difference measured at the centers of the two windows is linearly related to frequency. (To measure phase at the center of an FFT window, use an fftshift prior to the FFT, or flip the phase of the odd numbered FFT result bins.)  So you can use all of the FFT results, not just one bin.   The low frequency bins may help in determining the potentially ambiguous phase unwrapping of the higher frequency bins, if the signal is sufficiently wideband.


Cepstrum of white gaussian noise What are the statistics of the cepstrum of gaussian white noise?

\begin{align}\newcommand{\Nfft}{ {N_{\mathrm{FFT}} }}\DeclareMathOperator{\FFT}{FFT}\DeclareMathOperator{\IFFT}{IFFT}
x_i &\sim \mathcal{N}(0, \sigma^2)& i&\in\{0,\ldots,\Nfft-1\}\\[1.2em]
%
X   &= \FFT(x)\\
X_i &\sim \mathcal{N}\left(0, \frac{\sigma^2}{\Nfft}\right)\\[1.2em]
%
\left\lvert X_i\right\rvert^2 &\sim \chi^2\\[1.2em]
%
C   &= \IFFT\left(\log{\left\lvert X\right\rvert^2}\right)\\
C_i &\sim ~\ ???
\end{align}

What if I apply a window function?

\begin{align}
X &= \FFT(w \cdot x)\\
  &= \FFT(w)*\FFT(x) \\[1.2em]
%
C &= \IFFT\left(\left\lvert X\right\rvert^2\right)\\
C_i &\sim ???
\end{align}

Even an approximation would be nice. I'm working on a detector and understanding the noise structure is important for my approach.

The logarithm of a $\chi^2$-distributed random variable (RV) is a Log Chi-Squared RV.

Bartlett (of the frequency estimator fame) has a famous result that the sum of independent identically distributed $\ln\chi^2$ variables (and every logarithm is just a scaled version of any other logarithm) converges to a normal distribution much faster than a sum of $\chi^2$-distributed RVs.

Thus, I'm pretty confident that your sum of independent random variables will end up being normal in both real and imaginary part. But I'm too lazy to do a proof: I don't think this will help you that much!

  I'm working on a detector and understanding the noise structure is important for my approach.

The noise is additive to $x$, and since the Fourier transform is linear, it's also additive in $X$.

White noise being uncorrelated to the signal, the signal and noise powers are additive, too.

However, latest with the logarithmic transform, you lose that additivity, and knowing what the cepstrum of white noise in isolation looks like alone will not help you, because the cepstrum of (noise + signal) is not the (cepstrum of noise) + (cepstrum of signal).

That means your noise leads to a strange non-linear rotation thing of your complex cepstrum. It's hard to deal with additive noise in a nonlinear transform, which is why practically all approaches to that from low-SNR communications theory very much try to first invert the nonlinearity, before then dealing with the noise/signal separation. (see: digital predistortion, non-linear DFT for fiber optics, you name it)

Since the logarithm is a invertible operation, I'd strongly propose the first step of your cepstrum-based detector to be that of reverting the logarithm. That's kind of funny, because it just dumps you back to using the power spectrum instead of the cepstrum – but sometimes tools make it harder to work with statistically easy things.

Shot in the blue: As for a detector: Maybe do a noise-variance estimator on the power spectrum (which, again, is easy to calculate from the cepstrum), and do a signal power estimator on the cepstrum, and combine these results cleverly. Since a sensible combination for signal detection of these will be relatively smooth, but inherently pretty nonlinear operation, a reasonable approach to finding an optimum detector is probably stochastic gradient descent; "Neural networks" come to mind.


What is the use of the Welch lower bound when k > 1? The Welch lower bound on the maximum cross-correlation (inner product) of $m$ vectors, each of $n$ elements, is given by $c_{max} = f(m, n, k)$ (eg, see Wikipedia Welch bounds). It is always assumed that $m \gt n$ (otherwise the bound is 0).  Apparently the Welch bound has been widely applied in the restricted case that the vectors are real-valued and $k = 1.\;$ In that case, $c_{max} = \left[ (m-n)/(n(m-1)) \right]^{1/2},\;$ and is always real. The function has practical use in allowing us to compare the $c_{max}\;$ calculated from any given set of vectors with the theoretical minimum possible value.

In contrast, when $k > 1,\;$ then $c_{max}\;$ is frequently a complex quantity, even if the vectors are real-valued.  $\;c_{max}$ becomes real only if $m >> n\;,$ or more specifically, when $m \gt \binom{n+k-1}{k}.\;$
It is unclear to me whether $c_{max}$ has any use when $k \gt 1.$

$\quad$ Q1:  Does $k$ have any physical interpretation?

$\quad$ Q2:  Is the Welch formula with $k \gt 1\;$ useful in any area of physics, computational science,
$\quad\quad$ electrical engineering, or telecommunications? How do you interpret complex returned values?

Given any collection of $m$ unit-norm (or unit-Euclidean-length)
vectors in a real or complex $n$-dimensional vector space $\mathbb R^n$ or
$\mathbb C^n$, let $c_{\max}$ be the maximum magnitude of the
inner product between two vectors in the set. How small can $c_{\max}$ be?
One answer to this question is given by the
Welch bounds which are a collection of lower bounds on the minimum possible
value of $c_{\max}$. These bounds apply to all collections of $m$ unit-norm
vectors in $n$-dimensional space,
and they even hold for the case when the collection is a multiset (and thus
two vectors are identical, and so $c_{\max}$ necessarily equals $1$).

The Welch bounds are of the form
$$c_{\max}^2 \geq \sqrt[k]{C_k(m,n)}, ~k = 1, 2, 3, \ldots \tag{1}$$
or more correctly
$$c_{\max}^2 \geq \sqrt[k]{\max\{C_k(m,n), 0\}},~ k = 1, 2, 3, \ldots \tag{2}$$
because the function $C_k(m,n)$ can have negative value. Note for
example, that for $k = 1$, $(1)$ becomes
$$c_{\max}^2 \geq \frac{m-n}{n(m-1)}$$
where the right side is negative when $m < n$ and $0$ when $m = n$.
This is as it should be because we know that in a $n$-dimensional
space, we can always find $m \leq n$ orthonormal vectors and so
$c_\max^2 = 0$ for this set.
Thus, $(2)$ which becomes
$$c_{\max}^2 \geq \max\left\{\frac{m-n}{n(m-1)}, 0\right\}\tag{3}$$
is a much better way of expressing the bound.
When $m \gg n$ and so $(m-n)/(m-1)$ is approximately $1$, then
(as Welch pointed out as also have you in your question)
the right side of $(3)$ approaches $n^{-1}$ as $m \to \infty$ and so
we have that
$$c_\max  \geq \frac{1}{\sqrt{n}}\tag{4}$$
as a limiting form of the lower bound.

Turning to the question asked about Welch bounds for $k \geq 2$,
these are useful when $m \gg n$ and give larger lower bounds than
the right side of $(4)$. If I remember correctly, the Welch bound
for $k = 2$ gives something like
$$c_\max \geq \frac{2^{1/4}}{\sqrt{n}}$$
when $m \geq \frac{n^3+n^2}{2}$ which is an improvement on $(4)$.

For a discussion of some of these points (as applied to
crosscorrelations and autocorrelations instead of inner products),
I refer you
to my paper "Bounds on crosscorrelation and autocorrelation of sequences"
in IEEE Transactions on Information Theory, November 1979 which is
behind IEEE's paywall unless you have institutional access through your
employer or have paid personally for IEEExplore access.


Could someone please explain how to get frequency response of a given Bode Plot? Having trouble with bode plots. I understand that the  is to be converted to $dB$ but after that I'm stuck. Could someone please show me how this graph gives a frequency response of $10/1+jw10$

If your question is how to use this graph to find the transfer function the answer is simple : you can't.

First of all, this graph give only half of the information: you have a gain plot but not the phase plot, and 2 systems with different TF could have the same gain profile.

Second of all, this is obviously an asymptotic graph, as it presents a discontinuity around the cutoff frequency.

However, you could extract some information from this plot. You can see that it cuts the high frequencies with a slope of -20dB per decade, which means that there is a great chance that it is a 1st order low pass filter. Additionally, you can also read that  the cutoff pulsation is about 0.1 rad/s.

Think of it like a normal function, would you be able to tell what a function is if you any saw an asymptotic graph of one half of it? All you can do is make some educated guesses about the form of the TF.


Estimating the standard deviation of Gaussian filter from smoothed image Firstly, let's say that in order to smooth an image, I convolve it with a Gaussian function having standard deviation $\sigma_x$ and $\sigma_y$. I am now interested in knowing if there exist methods for estimating $\sigma_x$ and $\sigma_y$ from the smoothed image?
Secondly, to potentially make things more complicated, let's say the image pixels instead are drawn from a zero mean unit variance normal distribution (top image), which is then convolved as above (bottom image), could the standard deviations be estimated from the smoothed image?
Here is the MATLAB code used to generate the images. What I would like to achieve is to estimate the standard deviations, which in this case were $\sigma_x=2$ and $\sigma_y=2$.
img = randn(128,128);
figure(1); imagesc(img); colormap(gray); axis off; title('Noise image')
nimg = imgaussfilt(img, [2  2]);
figure(2); imagesc(nimg); colormap(gray); axis off; title('Smoothed noise image')

You can exploit the convolution theorem: Convolving with a Gaussian corresponds to Multiplication with a Gaussian of inverse $\sigma$ in frequency. Let x be the input signal and y be the output signal. Further, Let g be the Gaussian filter.

$$
y = g*x
$$

Or, in frequency domain (using captial letters)
$$
Y=GX
$$

Hence, if both Y and X is known, you can do
$$
G=Y/X
$$
to get your Gaussian Kernel. However, due to boundary conditions etc this will not yield the perfect Kernel. You have to fit a Gaussian function to it to get its standard deviation.

Additionally, if you have no knowledge of the input signal, but know its white noise, you can also fit a Gaussian Curve on the FFT of the output (since the input is white noise and hence flat):

N = 16*128;
img = randn(N,N);
s = 2;

G = fspecial('gaussian',[25 25],s);
Ig = imfilter(img,G,'replicate');
subplot(3,2,1);
imagesc(img);colormap(gray); axis off; title('Noise image')
subplot(3,2,2);
imagesc(Ig);colormap(gray); axis off; title('filtered image')

% Perform FFT
Img = fft2(img);
GImg = fft2(Ig);

subplot(3,2,3);
imagesc(abs(Img));
subplot(3,2,4)
imagesc(abs(GImg));

subplot(3,2,5);
fy = -N/2:(N/2-1);
hold off;
plot(fy, cumsum(fftshift(abs(GImg(:,1)))));
hold on;
plot(fy, normcdf(fy, 0, N/(2*pi*s))*sum(abs(GImg(:,1))), 'r');
hold off;

subplot(3,2,6);
div = abs(GImg ./ Img);
plot(fftshift(abs(div(1,:))));
title('Division of known input with output');


Hilbert transform of symbolic function in Matlab I have a signal of the form $f(\omega,x) = g(\omega) e^{i \omega (x +c)},$ for $g: [\omega_1 , \omega_2] \rightarrow C, \, x\in R$ and some constant $c.$  I want to get the function
$$
F (\omega,x) = \mathcal{H} ( \Re \{f(\cdot,x)\}) (\omega),
$$
where $\mathcal{H}$ denotes the Hilbert transform with respect to frequency $\omega,$ and $\Re$ the real part. I need the function $F$ to be again symbolic on $x$ such that I can differentiate it later.

I am trying to implement it in Matlab using:

L = 1e-2; K = 100; c = 1;
omega = (1/L)*(-K:K-1);

g = exp(-omega.^2/4);

syms x; f(x) = real(g.*exp(1i*omega*(x+c)));

but then the hilbert build-in function is not working. Do you have any idea how to do it?

The real part, the Hilbert transform and the derivation are all operators that commute.  So you can move the derivative just close to the function.


How to simulation different kinds of Noises in speech signal? What are the different kinds of noise in speech signal? How can I simulate the same in Matlab for adding to a clean speech signal?

Types of noises which you could try, that come to my mind at this time are:

White noise,
pink noise,
impulsive noise,
babble(multiple persons speaking in the background)

You can generate white(awgn)/pink noise(pass the white noise through a filter) in matlab.
You could record impulsive noises(pen clicks, knocks on the table, etc), and babble(record yourself and or others multiple times and overlap the recordings)


Why is CDMA so complicated, wouldn't this simpler solution work? I am now reading about Code Division Multiple Access method (CDMA), which is used for multiple hosts to send data at the same time.

The main idea behind it is to use orthogonal vector basis to encode multiple bit streams into one, send the joined stream across the link, then decode the stream differently for every separate message.

What confuses me is that in digital logic negative numbers are introduced, as well as 0 indicating no signal.

When I first heard that such encoding exists, my first idea of how it should work was different from the original and it looks way simpler.

Could someone check and add comments to my method? WHY is it not used instead of ORIGINAL CDMA?

Outline of my solution

I take n orthogonal vectors (where xor is the addition operation in my vector space), each vector corresponding to one chip code. I then encode every bit as a chip code if it is 1, as 0 otherwise. Finally, I xor all bit sequences together and send the joined result. At the other end, recipient can decompose the bit stream into unique components and read the message of every host.

Solution in more details

Assumptions

There are n hosts (numbered from 0 up to n - 1), and all of them want to send b bits at the same time using the same link.
The receiver gets a joined bit stream, but is interested in reading the message of only one host (WLOG the message sent from the 0-th host).

Problem: create a protocol that would satisfy the assumptions.
My solution

For i-th host, assign a chip code to him of the form 2^i (power not xor here) written in binary. Call this chip code c(i). For example:

Chip code 1000..00 would be given to host 0
Chip code 0100..00 would be given to host 1
...
Chip code 0000..01 would be given to host n-1.

Every single bit will be encoded using n bits.

if (host i wants to send 0) {
    encode it as 0000...00  // n zeroes
} else {
    encode it as c(i)
}

For host i wishing to send j-th bit, call the encoded bit as e(i, j). For example:

if (host i wants to send a bit stream, such that j-th bit is 1) {
    e(i, j) = c(i)
} else {
    e(i, j) = 0;
}

For every j, compute

d(j) = e(0, j) ^ e(1, j) ^ ... ^ e(n - 1, j)  // xor operation here

Concatenate d(0), d(1), ... d(n - 1) and send these n * b bits through the link. This will be the joined stream that the recipient at the other end should receive.
The recipient wishing to read the data sent by host 0 can look at (0-th, n-th, 2n-th, ..., b(n-1)-th) bits to see what the first host has sent.
The solution is simpler than original, because it uses only xor operation for encoding, and digits are only binary (i.e. no -2, +3 when doing vector addition). QED.

Does this solution work? If yes, why is it not used? The best reason I could guess is that we want to distinguish between hosts sending something and nothing. Are there other reasons?

Aside question: what are "nearly orthogonal" vectors in the context of CDMA?

Your solution works, and it is known by the name "time division multiple access" (TDMA). Basically, you assign time slots to each transmitter.

CDMA has several advantages compared to TDMA, but it also has its own disadvantages. Both schemes are inferior to more modern techniques such as OFDMA.

Main disadvantages of TDMA:

With data rate constant, the time slices become narrower as the number of transmitters increases. This means that the bandwidth also increases. With CDMA, you can increase the number of transmitters without increasing the bandwidth, by adding orthogonal code words.
TDMA requires very precise synchronization between transmitters. It is very hard to keep a large number of mobile units synchronized so that they won't transmit during the time slot of another unit. This may be solved using guardbands, but that decreases efficiency.

A couple of advantages of CDMA:

More resistance to fading and multipath. CDMA reflections are orthogonal to each other, so they cancel out.
If the codes are private, then there's a measure of inherent privacy.


Generating IF signal for FMCW radar I am trying to generate a signal that represents mixer output of FMCW radar. I am using MATLAB and already tried out this tutorial. There, Phased Array System Toolbox is used to showcase how FMCW radar with sawtooth signal works. The part that is of interest for me right now is up until the range-doppler response. Now what I am trying to accomplish is basically getting the same result (as in pass my signal to plotResponse function and get a range-doppler map with a detection of my target), but I want to get what they call "dechirped" signal from the start - without generating transmitted, reflected and received signals. I've been looking through so many papers in order to determine what formula I need to use to get the correct IF signal, but the more I searched, the more different formulas I found. What I've got now is the formula from this TI webinar, which is
So with that, I am trying to plot the response, here is my code:

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*range2time(range_max,c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp

%% target
R0 = 32; %range in meters
V = 40; %radial velocity, m/s

%%
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirps
    td = 2 * R0 / c; %round trip delay
    phi0 = 4*pi*fc*R0/c; %inital phase

    for j=1:1:samples
        a = (-2*pi*fc*V*i*tm/c ...    %phase shift
             -2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequency
        mix(j,i) = 0.5*cos(a);
        t = t + 1/fs;
    end
end

figure(1)
rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
    'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
    'RangeMethod','FFT','SweepSlope',sweep_slope,...
    'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
    'DopplerFFTLengthSource','Property','DopplerFFTLength',256);

clf;
plotResponse(rngdopresp,mix);
axis([-v_max v_max 0 range_max])

The radar that I'm simulating here does't have quadrature channel, so I have to only form real signal (although I tired doing hilbert(...) and cos(...)+1i*sin(...) just to see how that would look). What I assume I get is the Data Cube similar to the one from the tutorial. But no matter how I fiddle with the formula, I can't get the Velocity right. Here's the plot:

The initial phase I tried to form myself, thinking it would be 2*pi*fc*td, where fc is 77 GHz carrier frequency and td is the time delay of 2*R/c.
So what I think is happening is that I don't shift phase properly, since, from my understanding, that is what velocity is primarily estimated upon. The Doppler shift is too small to have significant impact.
Interestingly, I tried to play with the range and velocity values of the target, and noticed, that changing V doesn't affect the target's position on the heatmap too much, it almost doesn't actually. But changing the R0 in increments of 1 makes the target detection on heatmap shift horizontally and cycle to the other side of the map.
So my question is: how do i properly form the IF signal to be able to display it correctly in MATLAB's range-doppler response map?
Also, I acknowledge that I might be plotting or forming the data cube wrong. Or that my math is plain wrong from the beginning. But I've been so desperately trying to solve this by myself for the past week, that I just want someone to show me the light already :)

The principle behind FMCW is that you transmit a chirped signal and receive a time delayed version of it after reflecting from a target. After mixing and filtering, the resulting signal is a sinusoid at a frequency that is a function of the target's range. This frequency is known as the "beat" frequency $f_b$. Thus, the dechirped signal will have a form of

$$x(t) = e^{j(2{\pi}f_bt + \, \phi)} = e^{j2{\pi}f_bt}\,e^{j{\phi}}$$

Where $\phi$ is a general phase term that we'll ignore for now as it won't effect the determination of the beat frequency, and therefore range. Also, let's not worry about Doppler for now...that's an additional phase term we can easily add later. We'll concentrate on the expression for the homodyned (mixed) signal.

Let the chirped signal we transmit be

$$s_{tx}(t) = e^{j\pi\frac{\beta}{\tau}t^2}$$

Where $\beta$ is the sweep bandwidth of the chirp and $\tau$ is the chirp's length, or pulse width. After reflecting from a target we receive the signal after some delay $t_d$, we have

$$s_{rx}(t) = e^{j\pi\frac{\beta}{\tau}(t - t_d)^2} = e^{j\pi\frac{\beta}{\tau}(t^2 - 2tt_d + t_d^2)}$$

After mixing $s_{rx}(t)$ with $s_{tx}(t)$, which is equivalent to a frequency shift, the higher order term containing $e^{j\pi\frac{\beta}{\tau}t^2}$ drops off and we're left with

$$x(t) = e^{j\pi\frac{\beta}{\tau}(-2tt_d + t_d^2)} = e^{-j\pi\frac{\beta}{\tau}2tt_d}\,e^{j\pi\frac{\beta}{\tau}t_d^2}$$

Now compare this with the first equation, paying attention to the first term, again ignoring the constant phase term. We can then equate the phase functions
$$-\pi\frac{\beta}{\tau}2tt_d = 2{\pi}f_bt$$

So that then we have

$$f_b = -\frac{\beta}{\tau}t_d$$

Since we know our pulse travels at the speed of light $c$, we can rewrite the target's delay in terms of range $R$ and yield the mapping between target range and its beat frequency

$$t_d = \frac{2R}{c} => f_b = -\frac{2R\beta}{c\tau}$$

So generating a dechirped signal is straight forward since it's simply a sinusoid at some beat frequency $f_b$.

Note that these equations apply only to upchirps and downchirps. The negative sign takes care of itself in either case. Triangular and more exotic chirps yield additional frequency terms but this process can be expanded to cover that as well.

To add Doppler, you can add a constant phase term that is updated as you collect pulses to form a range-Doppler map. You can actually start at zero phase for the first pulse and progress from there for the purpose of a simulation. Your additional phase term would look something like

$$e^{j2{\pi}f_dnT_c}$$

Where $n$ is the current pulse number starting at 0 and $T_c$ is your equivalent pulse repetition interval (PRI) which is akin to your sweep time for FMCW radars.

EDIT: After having some time to look at your code directly, I found a few issues.

First, you are missing a factor of two in the Doppler component of the phase.

Second, without getting into theory, your particular system supports a Doppler span wider than what you force the horizontal axis to. This will erroneously change where you perceive the target to be.

Third, the time vector you are using to generate the beat frequency must be reset to 0 after every pulse. This is because the time vector needs to be relative to the time delay of the target $t_d$.

Here is your modified code. I currently do not have the Phased Array toolbox to generate and display the range-Doppler map, so I did it manually.

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*(2*range_max/c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp

%% target
R0 = 20; %range in meters
V = 40; %radial velocity, m/s

%%
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirps
    td = 2 * R0 / c; %round trip delay
    phi0 = 4*pi*fc*R0/c; %inital phase
    t = 0; % Reset

    for j=1:1:samples
        a = (-2*pi*fc*2*V*i*tm/c ...    %phase shift
             -2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequency
        mix(j,i) = 0.5*cos(a);
        t = t + 1/fs;
    end
end

%% Form the range-Doppler map (RDM)

% RDM axes
rangeBinAxis = (0:samples-1).*c/(2*bw);
dopplerBinSize = (1/tm)/chirps;
velocityBinAxis = (-chirps/2:chirps/2-1).*dopplerBinSize*lambda/2;

% 2D FFT to perform range and Doppler compression (i.e. form the RDM)
rdm = fftshift(fft2(mix), 2);

% Plot the RDM for the valid ranges of interest - targets ahead of you
figure;
surf(velocityBinAxis, rangeBinAxis(1:ceil(samples/2)), 20*log10(abs(rdm(1:ceil(samples/2), :))));
% surf(velocityBinAxis, rangeBinAxis, 20*log10(abs(rdm)));  % See the entire spectrum
xlabel("Range (m)");
ylabel("Velocity (m/s)");
axis tight;
shading flat;
view(0, 90);
colorbar;

% figure(1)
% rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
%     'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
%     'RangeMethod','FFT','SweepSlope',sweep_slope,...
%     'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
%     'DopplerFFTLengthSource','Property','DopplerFFTLength',256);
%
% clf;
% plotResponse(rngdopresp,mix);
% axis([-v_max v_max 0 range_max])

Some examples to show the target being positioned appropriately in the RDM (scale is in dB):

R0 = 32 m , v = 40 m/s

R0 = 150 m , v = 40 m/s