estimate the impulse response after equalization I want to estimate the impulse response of the channel at the receiver. Assuming
some arbitrary impulse response: h=[1  0.2 -0.4 0.0 0.6]. Once the equalizer is constructed, I get the equalizer weights in  w1.

Question1) How to get back the channel estimates: Say,

x is the input to the channel,

d is the output of the channel which is the input to the equalizer,

w is the equalizer,

y is the output of the equalizer or is the equalized signal.

Ideally, y and x must be same. Then would the channel estimates be the least square solution i.e., inv(w^Tw)*w*d' or  inv(w)?

Question2) Can I apply the LMS method for any source input --- Gaussian and non-Gaussian?

Below is the code for LMS equalizer

N=1000;
% number of samples

np = 0.01;
% noise power is 0.01
sp = 1;
% signal power is 1 which implies SNR = 20dB
h=[1  0.2 -0.4 0.0 0.6];  %unknown   impulse  response

x = sqrt(sp).*randn(1,N);
d = conv(x,h);
d = d(1:N) + sqrt(np).*randn(1,N);
w0(1) = 0;
% initial filter weights are 0
w1(1) = 0;
mu = 0.005;
% step size is fixed at 0.005
y(1) = w0(1)*x(1);
% iteration at “n=0”
e(1) = d(1) - y(1);
% separate because “x(0)” is not defined
w0(2) = w0(1) + 2*mu*e(1)*x(1);
w1(2) = w1(1);

for  n=2:N

    y(n) = w0(n)*x(n) + w1(n)*x(n-1);
    e(n) = d(n) - y(n);
    w0(n+1) = w0(n) + 2*mu*e(n)*x(n);
    w1(n+1) = w1(n) + 2*mu*e(n)*x(n-1);
end

Note that what you're trying to do is equalization, as opposed to channel estimation. If we ignore the noise for the moment then, ideally, the concatenation of the equalizer and the channel (modeled as a linear system) should be a pure delay:

$$(h\star w)[n]\stackrel{!}{=}\delta[n-K]\tag{1}$$

where $h[n]$ is the channel impulse response, $w[n]$ are the coefficients of the linear equalizer, and $K$ is some delay that depends on the system design. So by using the LMS algorithm to adapt a linear equalizer you do not directly estimate the channel coefficients $h[n]$. As far as I know there is no straightforward way to reliably estimate the channel weights directly from the equalizer weights.

Regarding your second question, the LMS algorithm doesn't require the source signal to be Gaussian.

A final comment on your code: with only $2$ filter taps you won't be able to sufficiently equalize the given channel. And I would recommend to use vectors in the coefficient update, so you can use different filter lengths without changing your code.


how to understand the MIMO zero forcing as a bank of decorrelators? I am studying the zero forcing as a simple linear receiver for MIMO detection. according to following notes:
http://www.ee.ic.ac.uk/bruno.clerckx/All.pdf.

The zero forcing can be understood as the following two ways.

For a single stream $\mathbf{q}$, we can project the received $\mathbf{y}$ onto a subspace orthogonal to channel matrix $\mathbf{H}$ without $\mathbf{h}_q$, then apply the matched filter to the projected signal or,
We can just work on all streams in one shot with the pseudo inverse of the channel matrix, i.e., $(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T$.

I can understand that the second approach certainly cancels the interference as $(\mathbf{H}^T\mathbf{H})^{-1}$ exists if columns of $\mathbf{H}$ are independent, and obviously $(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{H} = \mathbf{I}$ such that $(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{y}=\mathbf{x} + (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T\mathbf{n}$.

It seems to me that according to (2), a decorrelator for the $q^{th}$ stream is just the $q^{th}$ row of the pseudo-inverse of the matrix $\mathbf{H}$.

My question is, can we prove, or conclude (2) from (1)?

The following is the MATLAB script implemented according to the methods described in the above-mentioned notes.

%zero forcing

randn('seed', 0);
m=2; n=2;

x = randn(n,1)+1j*randn(n,1);
nz=0;

H = randn(m,n)+1j*randn(m,n);

y=H*x+nz;

%zeero forcing filter

G0 = inv(H'*H)*H';
x0=G0*y;

%disp(x-x0);

% transform the y into the null space of the

% subspace spanned by H1m (H w/o h1)

h1=H(:,1);
H1m=H(:,2:end);
[U, S, Vc]=svd(H1m);

% Q1 spans null space of the subspace spanned by H1m (H w/o h1)

% (actually it is null space of H1m')

% i.e., Q1'*y projects y into the null space of H1m

% y1 = Q1'*y = Q1'*h1*x1 + Q1'*n + zero (cancelled interference)

Q1=U(:,n:m);

% (Q1'*h1)' is the match filter to the transformed y1

% g1=h1'*(Q1*Q1'), note that Q1*Q1' is not orthogonal but Q1'*Q1=I

g1=(Q1'*h1)'*Q1';

% normalization

g1=g1/(norm(g1)^2);
x01=g1*y;

%x(1)-x01

%snr for stream 1

snr1=(norm(Q1'*h1))^2

nv=inv(H'*H);

%snr1 - 1/nv(1,1)
```

It seems like you understand your question already but are unsure it is correct? Let the decorrelator be $\mathbf{D}=(\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}^T$, then the symbols post-decorrelator are $\mathbf{z}=\mathbf{D}\mathbf{y}=\mathbf{D}(\mathbf{Hx}+\mathbf{n})$. The first stream, $z_1$, is computed by taking the first row of $\mathbf{D}$, call the row vector $\mathbf{d}_1$ (first row of $\mathbf{D}$), and performing the operation: $z_1=\mathbf{d}_1(\mathbf{Hx}+\mathbf{n})$. You get $z_1=x_1 + \text{filtered noise}$ because $\mathbf{d}_1\mathbf{H}=[1, 0, ..., 0]$.


How can I spread a modulated signal over walsh code I need to confirm if spreading signal by using walsh code in MATLAB is right or no,

Suppose I'm using QAM modulation, so I need to spread the signal over walsh code, is the below step is right:

M = 4;
alphabit = modem.qammod('M',  M ) ;             %QAM modulation
mod = transpose(modulate( alphabit , (0:M-1)));
avg = mod'*mod / M ;                % average power of constellation points
mod = mod / sqrt(avg);              % to make E[x^t x ] = 1
H = hadamard(4);   % The walsh code generation
H_use = H(:,2);  %using the second column for spreading.
S_coded = kron(mod,H_use);  % Spread the modulated data mod over the walsh code

My question is about the last step, is that right how we spread the modulated signal? ..

That's right, you can do spread using the kroon function, but be careful what you spread is not the data, it's the whole possibilities of data which is mod.The data is one of them, which is either, 0.7+0.7i or 0.7-0.7i ...etc.

then you should know how you can do despread in the receiver side too.

Good luck


Savitzky-Golay smoothing filter for not equally spaced data I have a signal that is measured at 100Hz and I need to apply the Savitzky-Golay smoothing filter on this signal. However, on closer inspection my signal is not measured at perfectly constant rate, the delta between measurements ranges between 9.7 and 10.3 ms.

Is there a way to use the Savitzky-Golay filter on not equally spaced data?
Are there other methods that I could apply?

Because of the way the Savitzky-Golay filter is derived (i.e as local least-squares polynomial fits), there's a natural generalization to nonuniform sampling--it's just much more computationally expensive.

Savitzky-Golay Filters in General

For the standard filter, the idea is to fit a polynomial to a local set of samples [using least squares], then replace the center sample with the value of the polynomial at the center index (i.e. at 0). That means the standard SG filter coefficients can be generated by inverting a Vandermonde matrix of sample indicies. For example, to generate a local parabolic fit across five samples $y_0\dots y_4$ (with local indicies -2,-1,0,1,2), the system of design equations $Ac = y$ would be as follows:

$$
\begin{bmatrix}-2^0 & -2^1 & -2^2 \\ -1^0 & -1^1 & -1^2 \\ 0^0 & 0^1 & 0^2 \\ 1^0 & 1^1 & 1^2 \\ 2^0 & 2^1 & 2^2 \end{bmatrix} \begin{bmatrix} c_0 \\ c_1 \\ c_2 \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}.
$$

In the above, the $c_0 \dots  c_2$ are the unknown coefficients of the least squares polynomial $c_0 + c_1x + c_2x^2$. Since the value of the polynomial at $x = 0$ is just $c_0$, computing the pseudoinverse of the design matrix (i.e. $c = (A^TA)^{-1}A^T y\space $) will yield the SG filter coefficients in the top row. In this case, they would be

$$
\begin{bmatrix}c_0 \\ c_1 \\ c_2 \end{bmatrix} =
\begin{bmatrix}
-3 & 12 & 17 & 12 & -3 \\
-7 & -4 & 0 & 4 & 7 \\
 5 & -3 & -5 & -3 & 5 \\
\end{bmatrix}
\begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ y_4 \end{bmatrix}.
$$

Note that since the derivative of $c_0 + c_1x + c_2x^2$ is $c_1 + 2c_2x$, the second row of the matrix (which evaluates $c_1$) will be a smoothed derivative filter. The same argument applies for successive rows--they give smoothed higher-order derivatives. Note that I scaled the matrix by 35 so the first row would match the smoothing coefficients given on Wikipedia (above). The derivative filters each differ by other scaling factors.

Nonuniform Sampling

When the samples are evenly spaced, the filter coefficients are translation-invariant, so the result is just an FIR filter. For nonuniform samples, the coefficients will differ based on the local sample spacing, so the design matrix will need to be constructed and inverted at each sample. If the nonuniform sample times are $x_n$, and we construct local coordinates $t_n$ with each center sample time fixed at $0$, i.e.

$$
\begin{align}
t_{-2} & = x_{-2} - x_0 \\
t_{-1} & = x_{-1} - x_0 \\
t_0 & = x_0 - x_0 \\
t_1 & = x_1 - x_0 \\
t_2 & = x_2 - x_0
\end{align}
$$

then each design matrix will be of the following form:

$$
A =
\begin{bmatrix}
 t_{-2}^0 &  t_{-2}^1 & t_{-2}^2 \\
 t_{-1}^0 &  t_{-1}^1 & t_{-1}^2 \\
 t_0^0 &  t_0^1 & t_0^2 \\
 t_1^0 &  t_1^1 & t_1^2 \\
 t_2^0 &  t_2^1 & t_2^2
\end{bmatrix} =
\begin{bmatrix}
 1 &  t_{-2} & t_{-2}^2 \\
 1 &  t_{-1} & t_{-1}^2 \\
 1 &  0 & 0 \\
 1 &  t_1 & t_1^2 \\
 1 &  t_2 & t_2^2
\end{bmatrix}.
$$

The first row of the pseudoinverse of $A$ dotted with the local sample values will yield $c_0$, the smoothed value at that sample.


FIR / IIR and Wavelets I am a self-learner of DSP by reading various books. I have accomplished basic understanding of Signals - CT and DT and a few transforms. I recently started to learn FIR / IIR Filters.

The part that I cannot understand is that they are termed as 'Filters' which for me logically means blocking / allowing from a threshold value - ex. lower values would be passed higher would be filtered or removed. So if we have a low pass filter - it would lead us to remove high values from a sequence and vice versa for high pass filters - is my understanding correct?

Like HPF (High Pass Filter):

x(n)={1,2,3,4,5,6,7,8,9}

Set threshold as >=5 so output sequence would be {5,6,7,8,9}

OK but the document states FIR / IIR about :

  Finite Impulse Response (FIR) : this type of filter gives a finite
  number of nonzero outputs (response) to an impulse function input. It
  does not use feed-back.

While

  Infinite Impulse Response (IIR) : this type of filter uses feed-back,
  so it could have an infinite number of nonzero outputs (response) to
  an impulse function input.

Now I cannot understand what FIR / IIR has to do with my concept of filters - allow / block high / low values. Where does the question of Feedback comes here?

Similarly for wavelets -

  We call an octave a level of resolution, where each octave can be
  envisioned as a pair of FIR filters, at least for the one-dimensional
  case. One filter of the analysis (wavelet transform) pair is a lowpass
  filter (LPF), while the other is a highpass filter (HPF). Each filter
  has a down-sampler after it, to make the transform efficient.  For
  example, a simple lowpass filter may have coefficients {1/2,1/2},
  producing outputs (x[n] + x[n - 1])/2, which is clearly the average of
  two samples. A corresponding simple highpass filter would have
  coefficients {1/2,-1/2}, producing outputs (x[n] - x[n - 1])/2, half
  the difference of the samples.

I am not able to get the concept of how and why here equation : (x[n] + x[n - 1])/2 and (x[n] - x[n - 1])/2 is being referred?

A filter applies a frequency-dependent gain to be applied to a signal.

The process you describe - removing values below a certain threshold -  is not filtering in a DSP sense.

A typical FIR filter is shown below. It consists of delays, z^-1, and gains, the $b$ terms.

Compare with the IIR filter shown below which has part of its output fed back into the delay structure (the gains are missing on this diagram).


Linear prediction (LPC) of Sine wave samples around maximas I think I missed class when this was explained ...

Anyway as part of a bigger project I have to implement a LPC to predict 2-3 future values of a sinusoidal process. I wrote a small Matlab m-file to calculate the predictor coefficients and plot the resulting predicted values.

It does look well so far when moving around $y(x) \cong 0$ but the values around the maxima/minima $y(x) \cong \pm 1$ are completely off (it predicts into the opposite direction):

I was expecting something like this:

I'm using the variables 'phi' to set sort of a operation point and 'unk' to set the starting point from where values should start to get predicted.

Here my m-file:

N = 100;            % samples
fs = 60;            % sample-frequency [Hz]
Ts = 1/fs;          % sample-time [s]
unk = 10;           % future values to predict (mostly used to set x(k) where we start predicting from)
P = 10;             % predictor order
phi = 0;            % default phase-offset
phi = pi/2;

k = 0:1/N:2;
x = sin(2*pi*fs*k*Ts+phi);  % reference values

x_pred = x(1:101);          % intialize with x ...
                            % last 'unk' values will be replaced by the predicted ones
x_kn = x_pred(1:end-unk);   % values we already know (we will try to predict last 'unk' values)

coeffs = aryule(x_kn, P);      % get 'P' filter coeffs using the known input 'x_kn' values

for i = 1:unk
    % x^ = -vec(a) * vec(x_kn)'
    nextValue = -coeffs(2:end) * ...            % P - 1 coefficient vector
                x_pred(length(x_kn)-1+i: ...    % last value we actually know (at pos len(x_kn))
                -1: ...                         % reverse iterate ...
                length(x_kn)-P+i)';             % until P values before len(x_kn)
    x_pred(length(x_kn)+i) = nextValue;         % now we know the value at len(x_kn)+1
end

figure(1);
plot(0:length(x_kn)+unk-1 ,x_pred, 0:length(x)-1, x); grid;
legend('predicted', 'real');

Disclaimer: This is explanation is based on observation of my MatLab plots and my note be 100% textbook correct!

Well after much looking around and experimenting I read something about the 'yule-walker' method for estimating the coefficients assuming the signal to be 'zero' outside of the observation point. This is not the case for my co-/sine-wave, so the estimated coefficients and the resulting predicted samples $\hat{x}(n)$ will trend towards zero. The further away from the $x$-Axis the sharper the slope towards 'zero' when using yule-walker.

For my problem an approach using Burg's method (coeffs = arburg(x_kn, P);) seems to yield the result I'm looking for, since burg does 'not' assume the Signal trending towards 'zero'?? At least this is how I can explain the result of my two plots. The 'star' marks the point where the prediction of future values start.

Yule-Walker:

Burg-Method:


Understanding MATLAB/Mathematica's bilinear interpolation I am trying to get a working understanding of how to resize images using bilinear and bicubic transformations. I coded a resizing function based on things I've read here and Wikipedia, etc., but I get different results than what MATLAB and Mathematica give.

A minimal working example is the image im=[100 50 0] in MATLAB or im={1, .5, 0} in Mathematica. My algorithm (and "common sense") tells me that resizing this image to 5 pixels should give
im=[100 75 50 25 0] and i={1, .75, .5, .25, 0}, respectively. However, I get
im=[100 80 50 20 0] and i={1, .8, .5, .2, 0} from MATLAB and Mathematica when forcing the bilinear method. Specifically, the MATLAB commands

im=[100 50 0];
imresize(im,[1,5],'bilinear')

returns

ans =  100.0000   80.0000   50.0000   20.0000         0

They are obviously doing the same, different thing than my code (and common sense).

Can someone explain their algorithm to me? I am guessing it my mistakes relate to how to "index" the new image. The original image, with only one row, has a constant "y" value of 1, with columns x=1,2,3. I am assuming that the new image could be viewed as having columns indexed as x=1, 1.5, 2, 2.5, 3, and interpolation is needed to compute the values of f(1.5,1) and f(2.5,1).

Thanks.

imresize gets you different results on Matlab and Octave. The one in Octave calls the interp2 function, which with 'bilinear' method will return you [100 75 50 25 0].

However, it seems like there is some extra processing in Matlab.  It performs affine transformation on the image. The interpolation is calculated on each pixel in the output image, then inverse map onto the indices of the pixels in the input image with the contributing weights. So I guess the mapped coordinates are not distributed uniformly (not exactly at 1.5 and 2.5 between [1,2], and [2,3] in your case).

Probably you can edit imresize_old to see more source code content than imresize.


Determine whether a signal is periodic or not and get fundamental period Here is my signal

Cos(n/2)*cos(pi*n/4)

cos(n/2) has period 4pi and cos(pi*n/4) has period 8

Now, the question is will the signal be periodic for fundamental period 32pi ?

$\cos\left(\dfrac{n}{2}\right)*\cos\left(\pi*\dfrac{n}{4}\right)$

Now finding out period of $\cos\left(\dfrac{n}{2}\right)$
$\cos\left(\dfrac{n}{2}\right)=\cos\left(2*\pi*\dfrac{n}{N1}\right) \implies \dfrac{1}{2}=2*\dfrac{\pi}{N1} \implies N1=4*\pi$

Now finding out period of second function $\cos\left(\pi*\dfrac{n}{4}\right)$
$\cos\left(\pi*\dfrac{n}{4}\right) = \cos\left(2*\pi*\dfrac{n}{N2}\right)\implies \dfrac{\pi}{4}=2*\dfrac{\pi}{N2} \implies N2=8$.

Now for total function period is LCM of $(N1,N2)$ i.e, $\mathrm{LCM}(4*\pi,8) \implies32*\pi$.
Therefore, the period of the function is $\boxed{32*\pi}$.


About a filter's rise time I have just a few questions regarding the rise time of filters.

A rise time seems to be measured by using a step input, with an input typically from 0 to 1. However this is what I cannot seem to find out...

-Does the rise time stay the same if input is from say 0 to 0.5 or some other arbitrary input?

-If not, how can I predict the rise time?

-If so, what happens if I input 0 to 0.5, hold at 0.5, and then step to 1? Does another rise time occur?

Thanks for any information!

When a system is linear, effects are proportional to causes. You can multiply the input by any constant, you'll get the same output multiplied by the same constant.

If you look at how the rise time is defined, you'll see that it is independent of any amplification factor.

You can also add two inputs, and the output will be the sum of the separate outputs corresponding to the separate inputs.

So it you apply a second, delayed step, there can be two situations

steady state (constant output) has been reached after the first step: then you can consider two distinct rises with equal times;
steady state has not been reached yet: there is a single, combined rise, with a different rise time.


Relation between Kalman filter and Sequential linear MMSE estimation Are the results of applying  Kalman filter and  recursive linear MMSE estimation process the same ?

The system model for the Kalman filter is:

$
x_k = F_k x_{k-1} + B_k u_{k-1} + w_k \\
z_k = H_k x_k + v_k
$

The system model for sequential linear MMSE estimation might be:

$
x_k = x_{k-1} \\
z_k = H_k x_k + v_k
$

Those are equivalent if $F_k$ is the identity, $B_k=0$ and $w_k=0$ and I believe the estimators turn out to be the same in that case. If $F_k$ is not the identity, or if either $B_k$ or $w_k$ are non-zero, then I see no immediate way to match sequential linear MMSE estimation to the Kalman filter system model.

On the other hand, let's consider the system equations for $k=1$, $k=2$ jointly:

$
z_1 = H_1 F_1 x_0 + F_1 B_1 u_0 + F_1 w_0 + v_1 \\
z_2 = H_2 F_2 (F_1 x_0 + B_1 u_0 + w_0) + F_2 B_2 u_1 + F_2 w_1 + v_2
$

This can be stacked into one equation:

$
z = A x_0 + \tilde{v} \,,
$

where $z$ is a column vector stacking $z_1$ and $z_2$. $A$ contains $H_1, F_1, H_2, F_2$, and $\tilde{v}$ contains $H_2, F_2, B_1, u_0, u_1, w_0, w_1, v_1, v_2$. Both the Kalman filter and the linear MMSE estimator can be applied to this linear model, and I would expect them to produce the same result. The Wikipedia article assumes that the error covariance for $\tilde{v}$ is zero, but we can trivially fix that by subtracting the error mean (which arises from the known inputs $u_0, u_1$).

So if we model the system with joint equations, rather than sequential, iterative equations, then linear MMSE estimation applies directly to the general model, and Kalman filtering gives no improvement for this more general problem.


Simple Image Edge Preserving Filter I am after a simple, yet effective, Image Edge Preserving Filter.
I need something which is faster than the Bilateral Filter with similar quality (The Guided Filter isn't good enough).
Are there such algorithms?

Recently I have seen the paper Hui Yin, Yuanhao Gong, Guoping Qiu - Side Window Filtering.
They suggest a really simple filtering framework for Edge Preserving Filter:

Basically, what they suggest is filtering the image with a set of filter based on the Box Filter.
This filter set is basically composed of 8 filters with different orientations and sub sets of the Box Filter (As seen in the figure above).
Once you apply all filters (Which each of them can be implemented very efficiently by all the efficient implementations available for Box Filtering) you chose, on a per pixel basis, the one most similar to the original pixel.
Applying it by iterations, yields very efficient and very good filter.
I took the Lena Image:

I applied 40 iterations of the filter (Link to 75 Iterations):

The full code is available on my StackExchange Signal Processing Q74674 GitHub Repository.


what is the name of this curve I was drawing a high pass filter response ,in a polar coordinates, the function is (z+1)/z and then it is 2+2*cos(w)
the plot is
What is the name of this plot, is it a frequency response or impulse response or some thing else.

and how to distinguish it from other curves names.

It is called a cardioid pattern. (cardioid because it looks vaguely heart-shaped)

You've shown us a polar frequency magnitude plot with the frequency given in radians. You'll notice by going counter-clockwise around the unit circle that low frequencies are attenuated while the higher frequencies are passed through.


Correlation of independent random processes Suppose $X(t)$ and $Y(t)$ be two independent random processes. Is $E(X(t_1)Y(t_2))$ necessarily zero?

No. Quoting Wikipedia's article Independence (probability theory):

  If $X$ and $Y$ are independent random variables,
  then the expectation operator $\operatorname{E}$ has the property

  $$\operatorname{E}[X Y] = \operatorname{E}[X]\operatorname{E}[Y].$$

Consider your $X(t_1)$ and $Y(t_2)$ as $X$ and $Y$ in this answer. If both $\operatorname{E}[X] \ne 0$ and $\operatorname{E}[X] \ne 0,$ then the product  $\operatorname{E}[X]\operatorname{E}[Y]$ will be non-zero.


Lloyd Max Quantization and Clustering : Part 1 The Lloyd-Max quantizer is a scalar quantizer which can be seen as a special case
of a vector quantizer (VQ) designed with the Linde Buzo
Gray (LBG) algorithm.

In k-means clustering, we are given a set of n data points in d-dimensional space  and an integer k and the problem is to determine a set of k points in $R^d$, called centers, so as to minimize the mean squared distance from each data point to its nearest center. A popular heuristic for k-means clustering is Lloyd's (1982) algorithm.

Confusion :

Is the Lloyd-Max algorithm alternatively in short called the Lloyd's algorithm ? Or are they completely different?

I am very confused, please help. Thank you.

Different.

Lloyd-Max is a special type of scalar quantizer design which is optimized (in terms of MSE) to source pdf. Hence the quantizer is generally non-uniform.
Lloyd's algorithm (and the more generalized LBG algorithm) is a scheme to design vector quantization. More information here and here.

That said, one should note that Lloyd's algorithm for VQ is essentially an extension of the scalar Lloyd-Max quantization to higher dimensions (vectorized sources).


What is proof of Basic Hilbert Transform? Attached is the question. Need to find the output for signal x(t).

Hilbert transform has the definition that its CTFT is:

$$
H(\Omega) =
\begin{cases}
- \frac{\pi}{2}   &,& \text{ for } \Omega > 0 \\
+ \frac{\pi}{2}   &,& \text{ for } \Omega < 0 \\
\end{cases}
$$

from which you define the time-domian impulse response $h(t)$ associated Hilbert transformer from the inverse CTFT as:

$$ h(t) = \mathcal{IFT} \{ H(\Omega) \}  = \frac{1}{2\pi} \int_{-\infty}^{\infty} H(\Omega) e^{j \Omega t} d\Omega $$

You can evaluate the above integral or just use CTFT properties and pairs to conclude that:
$$ \boxed{ h(t) = \frac{1}{\pi t} } $$

Note that the associated impulse response $h(t)$ of the Hilbert transfomer is noncausal, two-sided and unstable and its Fourier transform can be expressed alternately as
$$ H(\Omega) = -j \text{sgn}(\Omega) $$
where $\text{sgn}(\cdot)$ is the sign function.

Hilbert transformer is most useful for defining and generating the analytic signal $x_+(t)$ associated with the signal $x(t)$, with has the property that

$$
X_+(\Omega) =
\begin{cases}
2 X(\Omega)  &,& \text{ for } \Omega > 0 \\
0            &,& \text{ for } \Omega < 0 \\
\end{cases}
$$


Generating smoothed versions of square wave, triangular, etc I'm trying to get functions depicting smoothed out versions of wave forms(triangle,square, sawtooth, reverse sawtooth) which have the same amplitude and frequency as $cos(x)$.

I found a few trigonometric functions depicting such wave forms:

$\delta=0.05$,

$Square_{smooth}=2.1\tan^{-1}(\sin(\pi x/3)/\delta)/\pi$

$Triangle_{smooth}=1.25-5\cos^{-1}((1-\delta)\sin(\pi x/3.2))/(2\pi)$

$h=12.5$,

$ Sawtooth_{smooth}=-1.5 + 3((x/6-0.5) -(\tanh((((x/6-0.5)+0.5)-\lfloor{(x/6-0.5)+0.5}\rfloor-0.5)h)/(2\tanh(0.5h))+\lfloor{(x/6-0.5)+0.5)}\rfloor-0.5))$

but they do not have the same amplitude or frequency unfortunately. I think $-cos(x)$ might work for $Triangle_{smooth}$, but I'm still looking for ways to fix $Square_{smooth}$ and $Sawtooth_{smooth}$, anyone have any ideas?

It seems the amplitude is not scaled properly. Rather than (2*A/pi) using (A/atan(1/delta)) seems more appropriate. In other words I propose:

y = (A/atan(1/delta))*atan(sin(2*pi*t*f)/delta);

Below is a figure illustrating the difference between the two scaling approaches. For low delta values the difference is not clear but for high delta values the amplitude departs from what is desired. The images can be reproduced using the MATLAB code below.

t=linspace(0,2*pi,500);
delta = 0.1;
A = 1;
f = 1/(2*pi);
y1 = A*sin(2*pi*t*f);
y2 = (2*A/pi)*atan(sin(2*pi*t*f)/delta);
y3 = (A/atan(1/delta))*atan(sin(2*pi*t*f)/delta);

figure; hold on;
h1=plot(t,y1,'r-','LineWidth',5);
h2=plot(t,y2,'g-','LineWidth',4);
h3=plot(t,y3,'b-','LineWidth',2);
legend([h1 h2 h3],{'Sine','Smooth square wave','Corrected smooth square wave'},'Location','SouthOutSide');
axis tight; axis square;
set(gca,'FontSize',15); grid on; box on;
drawnow;


Phase tracking in communication receivers Why does Costa's loop PLL bandwidth need to be narrow? Text books mention that and yet give no explanation

The main purpose of a phase synchronization system, whether Costas loop or a PLL, is to track the phase of an input signal. The bandwidth of a PLL determines whether the signal falls within the range it can track or suppress. If the signal can pass through and the bandwidth is narrow, less noise enters the system increasing the loop SNR, while a wider bandwidth allows more noise deteriorating the tracking performance. But do remember that the narrower the bandwidth, the longer the acquisition time due to the smaller variation injected by the PLL.


Find the noise settings to reproduce it I work on images that contain shot noise (Poisson) acquired from a microscop. On few images, I have a "flat" zone that is supposed to have the exact same intensity, but it's not because of the noise. I was thinking that I could use this zone to estimate the quantity of noise, and then estimate the parameters K and Lambda of the Poisson equation.

Is there a method to do it?

What you'd often be looking for would the variance – the expectation of the squared difference to the mean value in that region. If the noise you add to the actual image is zero-mean, then the mean of the flat region is (in expectation) the actual intensity, and subtracting it from all pixels and squaring the result would give you the noise power.

However, shot noise is not zero-mean:

each flat-zone image pixel $v_i$  follows an actual-intensity-"offset" poisson distribution

$$v_i = m + n_i,\quad n_i\sim\text{Pois}(\lambda)$$

with $m$ being the actual intensity of the object, $n_i$ being shot noise samples, and $\lambda$ the intensity of the Poisson variable.

Now, since shot noise and actual image are independent,

$$\mathbb E(v_i) = \mathbb E(m + n_i) = \mathbb E(m) +  \mathbb E(n_i) =   m + \lambda\text.$$

Sadly, since we don't know the actual $m$  a priori, we'll have to dig one moment  deeper:

\begin{align}
\DeclareMathOperator{\Var}{Var} \Var (v_i) &= \mathbb E\left((m + n_i- \mathbb E(v_i))^2\right)\\
& = \mathbb E\left((m + n_i- m - \lambda)^2\right) \\
& = \mathbb E\left(( n_i - \lambda)^2\right) \\
&=: \Var(n_i)\\
&=\lambda
\end{align}

So, with the Variance of your observation sample as estimator for the $\lambda$ of your Poisson variable, you get the sole statistical property of your shot noise for free. If you want to know the actual flat color, you'd also subtract that from the average of the flat region.


Detecting Zero Crossing in Image Filtered by Laplacian Operator This is my first question on this site, complete noob. But here goes.

I have a positive Laplacian operator

[[0,1,0],
 [1,-4,1],
 [0,1,0]]

Now this Laplacian operator is used to find the outward edges of an image , IIRC. The image which I have in binary format is this

[[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Its hard to make out here, but the image represents some staircase type of  image. I have a PY script which does the convolve and gives me back the following result of convolution :

[[ 0.  0.  0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  1.  0.  1.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.  1.  0.]
 [ 0.  0.  0.  1.  0.  0.  0.  0.  1.  0.]
 [ 0.  0.  1.  0.  0.  0.  0.  0.  1.  0.]
 [ 0.  1.  0.  0.  0.  0.  0.  0.  1.  0.]
 [ 1.  0.  0.  0.  0.  0.  0.  0.  1.  0.]
 [ 0.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]]

Drawing straight line through all 1s actually does give me the outward edge of the binary image.

However, the zero crossing algorithm has confused me .
According to this link , the algorithm to detect zero crossing is:

  Do the following for each pixel $I(u, v)$:

  Look at your four neighbors, left, right, up and down
  If they all have the same sign as you, then you are not a zero crossing
  Else, if you have the smallest absolute value compared to your neighbors with opposite sign, then you are a zero crossing

But here, the 0 and 1 have the same signs , right ? Both are positive. I am confused in this.
If someone can take a small sample out of the result image, and mark some zero crossing for me, it will be clear to me.

Thank you in advance, it's a long post. I hope its clear though.


Terminology question on band-shifting I am trying to figure out if there is anything more to the phrase 'complex band shifting' than just, a simple band shift. Disclaimer, this is NOT a communications problem, and I only say this so that some answers do get carried away with bit rates, etc.

Basically, let us say that I am interested in a signal, that exists in the passband, from frequencies $f_1$ through $f_2$. Let us further assume that I band pass filter this signal with a brick wall BPF. Then, I multiply it with a sine wave of frequency $f_1$. This operation shifts my signal 'down to baseband', and at base band, it now has a (single-sided) bandwidth of $f_2$ - $f_1$, (as it did before). (I can now LPF to get rid of the double frequency component). That is, my signal now exists from $0$ to $f_2$-$f_1$.

To me, this is all there is to band shifting. That is, the key component is a multiplication with a real sinusoid.

Why then, and what does, the phrase 'complex band shift' entail exactly?

I understand that a multiplication with a complex sinusoid is done in comm systems because we then try to decipher phase from the I and Q, so I get that, but if I simply want to shift bands of a signal, (even to an arbitrary frequency), what is 'complex band shifting' and is that always necessary? Cant I just use a real sinusoid for general band shifting?

Thanks.

Complex band shifting is a little more flexible than real band shifting and filtering, and may be a little more efficient too, though if you change the sample rates when moving a signal to baseband that may not be true.

You are correct that you can move your signal to 0 to $f_2-f_1$ Hz in the method that you describe.  You could do the same thing in one step by multiplying the pass band signal by a complex sinusoid with a frequency of $-f_1 $.

Also, when you band shift with a real sinusoid you have to watch out for aliasing.  Band shifting with a complex sinusoid never causes aliasing (unless you reduce the sample rate at the same time, of course).


How detect signal from noise? I'm trying to detect when a signal is "on" from samples I have. Those samples come from a software defined radio, if that matters, and the signal i'm trying to detect is a WiFi signal.

The problem is separating "noise" from "signal", as I did by hand on the following image; green parts are when the signal is on, and non colored ones are when there is no signal (x axis: samples, y axis: energy level).

The main problem I have is that I do not now how to detect such signals, automatically. I tried by using threshold, but the noise level is subject to change. I can always determine when a signal is emitted « graphically », but i'd like to do this automatically.

Below is an exemple of "tricky" situation (signal is on on the left of the image):

What method could I use? I'm searching, specifically, to approach the one  in use in Wifi card, to sense signal and avoid interferences.

I have dealt with similar problems before. Probably the first you need to know is how does the noise spectrum look like.

First approach (mostly empirical): I see your signal of interest (SoI) is a train of rectangular pulses, thus the main frequency components are around DC. You can low-pass filter it to remove the main noise components (probably a moving-average filter will do). Then set a detection threshold (empirically), to a value that will be always surpassed by the filtered SoI, and rarely or never by the filtered noise.

Second approach (mostly theoretical): Check the noise components in the frequency domain (apply FFT). If the noise is nearly flat from DC to $f_s/2$, you can apply the optimal solution to this detection problem, which is implementing an energy detector, as follows:

Assuming we're in presence of (zero mean) Gaussian noise, estimate its statistics (noise variance will do).

How to do it?
Since your SoI is a train of rectangular pulses, its main frequency components are around DC. You can high-pass filter it to remove the SoI. If the filter's cut-off frequency is low enough, the filtered noise will keep a similar noise variance after filtering. Otherwise, multiply the after-filter variance by $(f_s/2)/(f_s/2 - W)$, where $W$ is the filter's bandwidth, and $f_s$ is the sampling rate.
Apply the optimal threshold for detecting a signal in AWGN:
$$
\lambda = 2~\sigma^2_\omega ~ \sqrt{N}~ \left(Q^{-1}(P_{fa}) + \sqrt{N}\right)~~~~[1]
$$

$\sigma^2_\omega$ is the noise variance, $N$ is the signal length (number of samples of the sliding buffer, see step 3), $P_{fa}$ is the design false-alarm probability (i.e., how likely you would admit a false alarm to be triggered).
Use a sliding buffer to analyze the whole signal by pieces, and what you need to compare to $\lambda$ is the energy contained in such a buffer: $$E = \sum^{N-1}_{n = 0} |x[n]|^2$$

More details on the second approach can be found, for example, here:

[1]  S. Atapattu, C. Tellambura, and H. Jiang, Energy Detection for Spectrum Sensing in Cognitive Radio, Springer 2014. Chapter 2:
Conventional Energy Detector.


Matching an LPC Magnitude Spectrum to FFT Magnitude Spectrum I'm trying to get an LPC magnitude spectrum to match to an FFT magnitude spectrum.  Basically I want the peak height in the LPC Spectrum to match the peak height in the FFT Spectrum.

Thanks to numerical recipes I have  some code that produces me an LPC spectrum with what looks like the correct shape.

bool CLPC::BuildMagitudeSpectrumFromLPC( float* pOutSpec, unsigned int outSpecSize )
{
    int spec    = 0;
    while( spec < outSpecSize )
    {
        const float theta = (M_PI * 2.0f) * ((1.0f / (outSpecSize * 2)) * spec);
        const float wpr   = cosf( theta );
        const float wpi   = sinf( theta );

        float sumR  = 1.0f;
        float sumI  = 0.0f;

        float wR    = 1.0f;
        float wI    = 0.0f;
        int i   = 0;
        while( i < mNumLPC )
        {
            const float wTemp = wR;
            wR      = wR * wpr - wI     * wpi;
            wI      = wI * wpr + wTemp  * wpi;

            sumR    -= mLPCs[i] * wR;
            sumI    -= mLPCs[i] * wI;
            i++;
        }

        pOutSpec[spec] = mMeanSquareError / (sumR * sumR + sumI * sumI);
        spec++;
    }
}

Its worth noting that I manually re-scale it to a dB scale elsewhere.

My problem appears to be that I can't figure out the scaling on this spectrum.  The above code returns values of greater than 1 (at seemingly the correct peak locations) for values that are not greater than 0dBFS on the FFT magnitude spectrum.  When I overlay the 2 the toughs also appear to be too high as well.

Has anyone any idea how to get the peak heights and troughs to correspond between the 2?

I'm NOT windowing before passing to the LPC.  Should I be?

Aso, Am I right in thinking that this line is a bit odd:

 pOutSpec[spec] = mMeanSquareError / (sumR * sumR + sumI * sumI);

I would have thought I'm taking the absoloute of the complex number represented by sumR and sumI.  Wouldn't this mean calling sqrt on it? ie.

 pOutSpec[spec] = mMeanSquareError / std::sqrt( sumR * sumR + sumI * sumI );

Doing the above doesn't help my problem all it does is squash up the LPC spectrum somewhat ...

I have a similar function to compute the Spectrum from the predicator coefficients, but instead of your code that I do not understand, I have the following :

pOutSpec[spec] = mMeanSquareError / (sumR * sumR + sumI * sumI);

I have something that computes the real and imaginary parts into my buffer :

data[2 * spec] = sumR / (sumR * sumR + sumI * sumI);
data[(2 * spec) + 1] = -sumI / (sumR * sumR + sumI * sumI);

Then I use different functions to display the spectrum according to the chosen format (Real, Magnitude, Square Magnitude or Log).

I also have some problems with scaling, both for the FFT Spectrum and the LPC Spectrum. I have chosen to scale them based on my display area, but as a result, their scale is not exactly the same on the display.


Filter size vs. FFT size and Overlap add First of all I would like to say that :

I'm new to DSP and in a "learning curve" hence the "basicness" of my questions.
I've read as many posts/articles as possible on the internet but still need "guidance".
I'm french hence the bad formulation sometimes.

My mission/data context :
Performing data analysis on real data sets (oil-rig data).
Sometimes data collected on the rig present low frequency component due to rig movements.
Data are collected every 5 seconds so the sampling rate is 0.2 Hz. The size of the data to be analysed will be 1024/2048/4096 points let's call it N is this post.
My mission is to perform a filtering on these data to eliminate the low frequency component.
This is not a real time application. The range of data to be analysed will be chosen and then the analysis performed.
Chosen approach :
Process these data with a FFT, filtering the FFT and perform an IFFT to obtain and clean(er) set of data. Basically the lowest frequency found in the FFT is to be eliminated.
Many users suggested me to use the overlap-add method to obtain the desired result.
First I was thinking of just zeroing low frequency related bins in my FFT but my readings on the subject suggested to rather use a high pass filtering (with spectral inverted windowed-sinc).
Designing this filter and applying it is what I post here for because it turned out to be way more difficult than expected.
Questions :
Do I have to use the overlap-add method considering that my application is not a real time one ?
For what I have understood I need to design a high pass filter kernel with a spectral inverted sinc function. Then I have to perform a FFT on my filter kernel and multiply points by points my filter FFT(1) coefficients with my data FFT(2) coefficient before performing a IFFT on the resulting FFT. Am I correct on this one ?

According to my readings I can obtain a spectral inverted sinc function by inverting (-sin(x)/x) it and adding 1 to the center of symmetry. I am right, do you guys confirm ? Do I add 1 on the one center point only ?

What is the right size to choose for my filter kernel ? Is it the same size N as my data ? As I have to multiply my FFT(1) with my FFT(2) I would think that I have to choose the same size but I'm probably wrong on this one and the overlap add method may be the answer.

The kernel filter designing is tricky, I don't know where to place the central lobe of my spectral inverted sinc, in 0, on a specific offset in the time domain ?

For your non-real-time application, you do not need OS or OA.  Direct form convolution (such as Octave function conv()) is easy and gives the same result.
But since the title of your post is "Filter size vs. FFT size and Overlap add", this graphic directly addresses that issue.


Incremental Online Update of Tf-IDF Indices for Image Retrieval In their paper Scalable Recognition with a Vocabulary Tree, Nister and Stewenius mention that after an offline training stage, new images can be inserted
on-the-ﬂy into the database. However, at the end of the training stage the IDF is computed and normalized. My question is, how would you incrementally update such a database composed of inverted document frequencies?

Ok, apparently, it is not much different from an update of the mean, given the accumulated statistics (mode statistics). Basically, when the mean and the number of elements are known, the idf weights can be re-computed whenever a new element is inserted. This can be achieved in two ways:

Allocate your tree with the maximum number of features. Keep all the elements inserted. Whenever a new insert occurs, re-compute TF-IDF over the entire set.
Only keep the statistics which are required to update the weights. This way, an entire pass over the data is not necessary.

Normalization is again performed over the entire dataset once the idf is updated. Removal is not much different from the insertion, but one needs to use a suitable pointer-wise data structure if it is desired to have the computational complexity both for insertion and deletion.


How to properly filter a BPSK demodulated signal? I am new to this field and yesterday I began a project for BPSK modulation. So this is what i did:

$\ x(t) = 0 $ or$\ 1 $ binary data

$\ Fs = 15 $ sampling frequency kHz

$\ Fc = 10 $ carrier frequency kHz

BPSK modulation

$\ BPSK(t) = A.sin(2.\pi.Fc.t + \pi.x(t) ) $

BPSK demodulation

$\ S(t) = BPSK(t) . (Ac.sin(2.\pi.Fc.t))$

low pass filter

void lowPassFrequency(float input[], float fs, float output[], int points)
{
  float  CUTOFF = 2 * Fc;
  float RC = 1.0/(CUTOFF*2*M_PI);
  float dt = 1.0/fs;
  float alpha = dt/(RC+dt);
  output[0] = input[0];

  for(int i = 1; i < points; ++i)
  {
      output[i] = output[i-1] + (alpha*(input[i] - output[i-1]));
  }
}

it seems to me that the filtering part is not as it should be ... it supposed to be smoothed so i can get values like -1 and 1 ... Any ideas please?

You must use a sampling frequency of at least the Nyquist rate, which is greater than 15 kHz in this case. It should be greater than $2(f_c + W/2) = 2 (10~\text{kHz} + W/2)$, being $W$ the signal's bandwidth.

You did not specify the bit rate, which also must be considered to determine the signal's bandwidth, and thus the Nyquist rate.

You signal is probably having time-aliasing. More details in the answers to the following question: A voice signal is band-limited to 3.3 kHz. What is the Nyquist frequency?


Derivation of Transfer Functions for Kalman Filter Hi All: I'm  somewhat familiar with the kalman filter from a statistical point of view. But lately I've been trying to familarize myself with the linear systems-EE way of looking at it. So, I've been trying to read the Hannan and Deistler 1988 text slowly. I find it difficult and there is
material in there that I don't follow so this  is what my question is about.

I will write down verbatim what the book has and then explain the parts of it that I don't understand. I don't expect the actual detailed answer but rather a book or notes that explain the derivation. I'm totally new to this notation and terminology so maybe my questions can be easily answered but I don't expect it ? In fact, notes or a reference would still be appreciated even if the answers are straightforward.

Quote from Hannan and Deistler: Page 8

  As will be shown, equations ... can be transformed into a state space system

  1.1.22)$ ~~  x(t+1) = Fx(t) + Lz(t) + K \epsilon(t) $

  1.1.23)  $ ~~ y(t) = H x(t) + \epsilon(t) $

  where $F \in R^{n \times n}$,  $L \in R^{n \times m}$,
  $K \in R^{n \times s}$ and $H \in R^{s \times n}$.

Then there about five pages before the authors comes back to these equations above. I won't quote these five pages but they are mostly about how to use the lag operator to and the z-generating function. Although I find these pages confusing, I do understand the concept of a lag operator. I have a feeling that, if someone can answer my question or point to a useful reference or notes, then they won't need these pages. If someone does need them, I can copy them and send them in a pdf.

So, now I continue to where the authors comes back to above and I quote:

Hannan and Deistler, page 14

  For a state space system (1.1.22), (1.1.23) we have det$(I - Fz) = 1$ for
  $z = 0$ and thus there always exists a causal solution

  1.2.21)$ ~~~  x(t) = (I - Fz)^{-1}z[Lz(t) + K e(t)]$

  where  $(I - Fz)^{-1} = \sum_{j = 0}^\infty F^{j}z^{j}$, provided that the sum in 1.2.21) exists. In this sense, causality is built into the definition of state-space systems.

  The transfer functions corresponding to (1.1.22), (1.1.23) are

  $k(z) = H(I_{n} - Fz)^{-1}zK + I_{s} = I_{s} + \sum_{j=1}^{\infty} HF^{j-1}Kz^{j}$

  $l(z) = H(I_{n} - Fz)^{-1} zL = \sum_{j=1}^{\infty} HF^{j-1}Lz^{j}$

  Starting the system with given initial values $x(0)$, $z(0)$, $\epsilon(0)$ at t = 1, we have

  $ x(t) = \sum_{j=1}^{t} F^{j-1}Lz(t-j) + \sum_{j=1}^{t} F^{j-1}K\epsilon(t-1) + F^{t}x(0)$

  $y(t) = \sum_{j=1}^{t} HF^{j-1}Lz(t-j) + \sum_{j=1}^{t} HF^{j-1}K\epsilon(t-1) + HF^{t}x(0) + \epsilon_{t}$

The three parts that I don't understand regarding above are

1) the statement about the determinant of $(I-Fz)$ being 1 when $z$ = zero.
I mean I understand why that's true but what does it have to do with the price of tea in china ? or causality ?

2) the derivation of the transfer functions, $k(z)$ and $l(z)$. I do understand 1.2.21) so I imagine that 1.2.21) is used to derive them ????? I also understand the second equalities in $k(z)$ and $l(z)$. But not the first equality signs.

3) the equation for $x(t)$ at the very end. I do see how $y(t)$ is obtained once one has $x(t)$.

I imagine that the material that I'm confused about above ( which is most of it )  is hidden deeply in the 5 pages that I didn't quote but, like I said, I found those pages to be quite "opaque", to say the least.

My guess is that I need is nice reference or lectures notes that explain why, when it comes to the KF,

A) the determinant of $(I-Fz)$  is important and

B) how the associated transfer response functions are obtained  from a state space representation?

C) How the $x(t)$ equation is obtained given the transfer function ?

Note that I'm definitely opened to fighting with this material on my own so even just a good link is fine.

Finally,  if anyone knows of a book that has similar objectives to Hannan
and Deistler but is less concise and shows more steps to things (i.e. doesn't assume as much mathematical maturity ), it's appreciated.

This is an attempt to answer. It's not complete, and a little fuzzy, but feel free to comment and I can edit / update.

  1) the statement about the determinant of $(I-Fz)$ being 1 when $z$ = zero.
  I mean I understand why that's true but what does it have to do with the price of tea in china ? or causality ?

So let's have a look:
$$
x(t+1) = Fx(t) + Lz(t) + K \epsilon(t)\\
zx(t) - Fx(t) = Lz(t) + K \epsilon(t)\\
(I-Fz^{-1})x(t) = z^{-1}[Lz(t) + K \epsilon(t)]\\
x(t) = (I-Fz^{-1})^{-1}z^{-1}[Lz(t) + K \epsilon(t)]\\
$$

One thing to note about this is that $z$ is used twice: once as the lag operator and once as the "reference" signal.  That confused me for a while. Another thing is that they appear to use $z$ as the lag operator, whereas $z^{-1}$ is usually used. Above I've assumed the conventional approach to $z$ as the advance operator.

My understanding of causality comes from the ROC associated with the $z$ transform.  A $z$-transform is causal if its ROC extends outwards from the outermost pole.  A $z$-transform is BIBO stable if its ROC includes the unit circle.

I believe what they really should have said is:

  Because we assume that $\displaystyle(I - Fz)^{-1} = \sum_{j = 0}^\infty F^{j}z^{j}$, we assume $j\not < 0$ and so assume causality.

which is also sort of implied in their use of state space equations, as they say In this sense, causality is built into the definition of state-space systems.

  2) the derivation of the transfer functions, $k(z)$ and $l(z)$. I do understand 1.2.21) so I imagine that 1.2.21) is used to derive them ????? I also understand the second equalities in $k(z)$ and $l(z)$. But not the first equality signs.

Assuming for now that $L=0$, then the state space form of the $z$ transform of (1.1.22) and (1.1.23) is just:
$$
H(I - z^{-1})K + I
$$
see these notes, particularly equation (10.20):

Though they use $A, B, C,$ and $D$ instead of $F, K, H,$ and $I$.

This is just the standard way to represent transfer functions in state space.

  3) the equation for $x(t)$ at the very end. I do see how $y(t)$ is obtained once one has $x(t)$.

I believe the $x(t)$ equation is just had by the usual solution to a linear, constant coefficient difference equation with non-zero initial conditions (but using their notation). It's not the same notation, but try this.

  My guess is that I need is nice reference or lectures notes that explain why, when it comes to the KF,

  A) the determinant of $(I-Fz)$  is important and

See answer to 2) above.

  B) how the associated transfer response functions are obtained  from a state space representation?

See references above.

  C) How the $x(t)$ equation is obtained given the transfer function ?

See answer to 3) above.


The computation procedure problem of calculating the BER Here is the processing of calculating the BER average

\begin{align}
P_e &=\int^\infty_0Q(\sqrt{2uSNR})e^{-u}du \\
 & =\frac{1}{\sqrt{2\pi}}\int^{\infty}_0\int^\infty_{\sqrt{2uSNR}}e^{-\frac{t^2}{2}}e^{-u}dtdu \\
& =\frac{1}{\sqrt{2\pi}}\int^{\infty}_0\int^\frac{t^2}{2SNR}_{0}e^{-\frac{t^2}{2}}e^{-u}dtdu \\
& =\frac{1}{\sqrt{2\pi}}\int^{\infty}_0(1-e^{-\frac{t^2}{2SNR}})e^{-\frac{t^2}{2}}dt \\
& =\frac{1}{2}(1-\sqrt{\frac{SNR}{1+SNR}})
\end{align}

channel state $h$ is a complex gaussian r.b.$CN(0,1)$,$SNR=E_s/N_0$,the BER of received symbols is $P_{e|h}=Q(\sqrt{2|h|^2SNR})$

I don't really understand these computation procedure step by step,for example.

Why is $\int^\infty_{\sqrt{2uSNR}}$ become $\int^\frac{t^2}{2SNR}_{0}$ ?
Why is $\int^\frac{t^2}{2SNR}_{0}e^{-u}du$ become $(1-e^{-\frac{t^2}{2SNR}})$ ?

Why is $\frac{1}{\sqrt{2\pi}}\int^{\infty}_0(1-e^{-\frac{t^2}{2SNR}})e^{-\frac{t^2}{2}}dt=\frac{1}{2}(1-\sqrt{\frac{SNR}{1+SNR}})$ ?

Can anyone explain them and teach me?

This is pure mathematics. However, these are some hints.

1) In 2-D coordinates $t$ and $u$ the multiple integral scans the region limited by $t = f_u(u) = 0$ and $t = g_u(u) = \sqrt{2u\textrm{SNR}}$ for $t > 0$ and $u >0$. It is the same region limited by two functions $u = f_t(t) = 0$ and $u = g_t(t) = t^2/2\textrm{SNR}$ for $t > 0$ and $u >0$.

2) because $\int_a^b e^{-u} du = e^{-a} - e^{-b}$

3) use 2) with the identity Q(0) = 1/2.


Replace phases of image FFT I have 2 images of the same scene, captured by different cameras. So they a bit displaced. As far as I can understood, phases of FFT contains info about image details. My goal is to take intensity from first image, details from second and produce enhanced highly detailed image with proper intensity.

I tried to solve problem as is - apply FFT for both images, take magnitudes from first image and phases from second and apply iFFT. But I didn't get what I want. Image looked corrupted, see below.

Then I tried to compute Phase correlation , align images and
apply the algorithm described above again. Result was the same.

Is it possible to do what I want? May be I need to do some calculations in frequency domain, like phase shifting? Matlab code is highly appreciated!

What you want to do is impossible with your technique. The phase of an FFT does not provide the details of an image. Even if it did, why would it enhance an image, if you exchanged the "details" of both images without any calculation?

The phase of an FFT of an image mainly holds the structure information, i.e. how the frequencies contribute to the image (in terms of phase offsets). This is valid for fine and coarse structures in the image. The magnitude tells you, how much each frequency contributes to the image. This information is in most cases less important, as can be seen by the general shape of the 2D magnitude spectrum (many strong low frequency contributions and a rapid decay towards higher frequencies).

I assume, what you would like do do is more a quite complex problem to be solved in image space. You might want to look for publications and algorithms on super-resolution for your application.


FFT subset output yields computational efficiency The traditional FFT gives $N=2^k$ outputs for $N=2^k$ inputs (generally complex.)  The time effort is $O(N \cdot \log(N))$.  If all I want is a compact subset of $M \ll N$ output frequency points, is there some method to "depopulate" the FFT late butterflies and achieve more efficiency than the depopulation burden adds?

The question of the FFT-like algorithm for fast DTFT computation? is close, but first misses the mark being really an off-frequency question, then the Goetzel transform efficiency threshold is far too low for my case.

Assuming your question is about bypassing the later stages of the FFT if not all samples are needed, this may give you some insight:

Knowing that the FFT alorithm is a successive decimations into even and odd samples so that smaller FFT's can be computed and thus leading to the great efficiency of the algorithm; here is a view of what the digital spectrums look like prior to combining in the final stages:

Everything between the pciture above showing the FFT decimation and the picture below showing the final FFT combining is covered in this post, so I won't repeat it here, but please review it if the plot below is not immediately clear: Is this signal perfect reconstructable?

From reading the linked post and the plots shown above, it will hopefully be clear that unless the signal is combined in the final stage, the spectrum will be superimposed with other regions of the spectrum.  So if we work back one stage from the final output, as shown in the figure below, we can hopefully associate with the figure above to see that the upper N/2 pt DFT contains ALL the spectrum (folded on itself), and the lower N/2 pt DFT also contains ALL the spectrum, however separable due to the phase relationship described if combined with the appropriate phase rotations as shown.

So that said, we can see the reason for the combining (and that the signal would be irrecoverable if we did not have access to both prior DFT's). However as long as both DFT's are available for a given point of interest, we can use the mapped output from both DFT's for that point without calculating the rest of the final stage. This of course works back in the same pattern through the implementation of the smaller DFTs (since they two are made up of two DFTs, which are each made up of 2 DFT's etc.

We can see from this plot how the final stage works its way back through these successive DFT's to the input. This plot also gives immediate insight into the possible reduction if only a subset of the points were needed (and shows how quickly any savings if diminished- in this simple case, 1 point is touching 21 of the 48 nodes shown!).  I think it is also interesting (perhaps a little mind-blowing) that if you work your way back to the first stage, each 2 pt DFT output for the same reasons described also contains ALL the spectrum, but in this case folded over itself multiple times (depending on how many stages).  Hence every 2 pt DFT at the input is connected to every output in the final DFT.


How to measure noise, consonance, dissonance of a sound? How can I make a graph of the noisiness of some audio durations of 100ms - 10s with 100% being white noise and 0% being a sine wave.

I have FFT at the moment, are there other ways to know the commensurate random nature and perhaps even the dissonance of sound?

If there are spectral peaks and you want to calculate the dissonance of this sound, you can calculate the sensory dissonance. See the Essentia explanation and the references mentioned there.

If you rather want to measure how noisy a sound is, spectral flatness may be the audio feature of your choice. E.g. see this Wikipedia page.


audio for VOIP - filter and resample or viceversa? I am a noob to signal processing, but I need to create some audio files for IVR for Asterisk VOIP.

So my lady friend with a nice voice records her message on a Windows machine and sends me a 44KHz, stereo 32bit wav file.

I currently squash it into a 8KHz, mono, 16bit signed wav file for Linux, and then convert the wav to alaw or ulaw PCM format.

One source of info for this process is this link.

I can do some magic with Audacity:

Stereo to mono
Equalize, using a "voip" or "telephone" filter: discard anything over 3KHz and under 300Hz.
compress volume/levels to -12 to -15 dB
Resample 44 -> 8KHz
Export as 16bit signed PCM wav
Move it to Linux and sox it to alaw or ulaw

The result is acceptable, but I feel it could be improved.
In particular, after resampling, the voice sounds muffled.

Is this order correct?

Do you have suggestions to get better results?

Your order is fine. The biggest question I have is with,

Equalize, using a "voip" or "telephone" filter: discard anything over 3KHz and under 300Hz.

Since you are resampling to 8 kHz in step 4, you already have less than 4 kHz, I don't see the point in the 3 kHz limit. And knowing nothing about the 3 kHz filtering you're using, I wonder how much intelligibility you're losing in that step. This is easy enough for you to test.

I believe you'll get better quality but less dynamic range with ulaw, compared to alaw. Other than this, the compression (step 3) is probably the next most important aspect to experiment with.


FIR Filter design, using different filter length than signal length Given some neurophysiological application I am filtering in real time data of limited length (e.g. in epochs of 100 ms) using a sampling rate of e.g. 100 Hz.
Due to this short time segment I am trying to implement a FIR filter. This filter should be of high order, and therefore have a length of more than 100 samples.
This would lead to the fact that the "to be filtered signal" has a different length than the estimated filter itself. Using the scipy.convolve function I can then specify different methods how to convolve the signal (full, valid, same).
My question is quite trivial, and related to lack of experience in the field, but does it make sense at all to have different filter and signal length?
Is there literature how to interpret then the filtered signal? Since the convolution needs to define some padding for the signal itself the filtered signal should be only in a certain range interpretable right?
EDIT:
Hereby I want to present an example. I plot a random signal in orange, and a FIR filter with a filter length of 10 s (estimated using mne.create_filter):

My understanding problem is visualized in this context: the signal itself is much shorter than the filter. This leads obviously to the fact that the filtered signal has a longer length. How can I now interprete the filtered signal with respect to the timing of the original signal that had fs=128 Hz? Which time points of the filtered signal correspond to the original signal?

There's no requirement that the signal to be filtered has to have the same length as the filter - in fact, that's practically never the case.
So, there's no problem here; convolution doesn't need to "define some padding"; you'll find discrete convolution is just defined as finite sum over the available lengths.
The length difference between signal and filter isn't of relevance to "track back" a phenomen in the input signal: the time is always the time in the output signal, minus the filter's group delay. Assuming your filter is linear phase / time-symmetric, that's simply half the filter length.


What does "ripple limit cutoff" mean? So I'm reading the datasheet for MAX98357A.
It's an audio amplifier that receives digital signal using I²S bus and has integrated digital low-pass filters of the IIR and FIR types.

In addition to the usual -3 dB and -6.02 dB cut-off points, it also has the "Ripple limit cutoff". This is a term I did not see during my education, and I'm not sure exactly what it means.

Here's an example of the specification:

Ripple limit cutoff 0.43 x Fs
-3dB cutoff         0.47 x Fs
-6db cutoff         0.5  x Fs
Stopband cutoff     0.58 x Fs

I would guess that this is the last point in the frequency response with a passband gain that is still inside the ripple spec. Meaning if you use the 3-dB bandwidth, you will have an attenuation of 3 dB (and therefore a "ripple" of at least 3 dB assuming the ripple spec for the digital filter is less than 3 dB). On the other hand, if you use the ripple limit bandwidth, your attenuation at that point will be within whatever the specified ripple is.

I couldn't find the spec for the ripple in the data sheet, but consider the following example. If the ripple were 0.1 dB and the passband gain were 0 dB then the attenuation at the ripple limit point would be 0.05 dB and from there the gain would decrease until reaching the 3 dB point where the attenuation would be 3 dB.

There is a chance I am wrong about this, since it isn't a standard term I'm familiar with, but I'm posting the answer because it seems very likely to be correct. Maybe someone else with more experience can provide a better response or point out my error if there is one.

The reason they probably include both is that a digital filter designer would only be interested in the ripple limit bandwidth (which we would generally refer to as the passband edge when doing an equiripple design) while many analog filters are specified with their 3-dB (or 6-dB) bandwidth. Probably there was some kind of internal discussion about which one to use, and they wound up with both. But that's all speculation on my part.


Why would I want to define a modulation index for each tone (DSB-FC)? So the exercise is basically a signal $f(t)$ that is going to modulate the carrier $A\cos(\omega_ct)$ using a modulation index of $m=1$. I have to find $A$ and the power of the modulated signal:
$$
f(t)=\cos(\omega_mt)+2\cos(2\omega_mt)
$$

The minimum amplitude of $f(t)$ is $-2$. Then $A = 2$. The power of the signal is, assuming that $R = 1 \ \Omega$:
$$
P = P_c+P_s=\frac{A^2}{2}+\frac{\overline{f^2(t)}}{2}
$$

Having in mind that:
$$
\overline{f^2(t)}=\frac{1^2}{2}+\frac{2^2}{2}=\frac{5}{2}
$$

The power is:
$$
P =\frac{A^2}{2}+\frac{\overline{f^2(t)}}{2}=\frac{2^2}{2}+\frac{5}{4}=3.25
$$

In the book the author uses an effective modulation index that is defined as $m_t = \sqrt{m_1^2+m_2^2}$ where $m_1=1/2$ and $m_2 = 2/2$. So the power is:
$$
P = P_c\left(1+\frac{m_t^2}{2}\right)=2\left(1+\frac{1.12^2}{2}\right)=3.25
$$

My question is, why would I want to define a modulation index for each tone? What do I get from that?

Another thing that I don't understand is that according to this guy the condition that ensures that there's no overmodulation, regardless of the tone frequencies, is $m_1+m_2\leq1$. Obviously in this case the condition is not met but I'm pretty sure there's no overmodulation with $A=2$.

I'm no expert on this, but until the big guys give an answer, here's my take:

For $f(t) = \cos(100t)+2\cos(2\cdot 100t)$, so $\omega_m = 100\ \rm Hz$ for example:

The amplitude's peak value$\ =3$ exceeds $A=2$, so there is overmodulation. This answer only addresses your last point of "I'm pretty sure there's no overmodulation".


relation between DFT to CTFT The signal $$x(t)\;\;\;\;0\leq t\leq 0.2s $$

We know that the CTFT of $x(t)=0$ when $|w|>2*\pi*10^4$

We sample $x(t)$ in sample space of $$T=25\mu s$$ or $$F_s=1/T=40000Hz$$and we get a series with length of 8000.

$$x[n]=x(nT)\;\;\;\;0\leq n \leq 7999 $$

We calculate the DFT of $x[n]$ in size 8000 also.

What is the relation between the CTFT to the DFT

$$X^F(\omega) \rightarrow X_N^d[k]$$

The answer is:

$$
\begin{equation}
  X_N^d[k]=\begin{cases}
    \dfrac{N}{T}X^F(10\pi k), & 0\leq k\leq 3999\\ \\
    \dfrac{N}{T}X^F(10\pi (k-8000)), & 4000\leq k\leq 7999
  \end{cases}
\end{equation}
$$

Let's give it a shot.

We know that during sampling, we obtained a discretized version, $x(nT_s)$ of a continuous time signal $x(t)$ sampled every $T_s$ seconds. I will denote their corresponding Fourier Transforms as $X_{DT}(e^{j\omega})$ and $X_{CT}(j\Omega)$, respectively, according to the majority of the bibliography. Indices denote discrete time (DT) and continuous-time (CT) Fourier Transforms.

The Discrete time Fourier Transform (DTFT) of the sampled signal is given by $$X_{DT}(e^{j\omega}) = \frac{1}{T_s}\sum_{l=-\infty}^{+\infty}X_{CT}\Big(j\Big(\frac{\omega}{T_s} - l\frac{2\pi l}{T_s}\Big)\Big)$$

Assuming that there is no aliasing (like in your case), over one period this turns into $$X_{DT}(e^{j\omega}) = \frac{1}{T_s}X_{CT}\Big(j\frac{\omega}{T_s}\Big)$$

Due to the $2\pi$-periodic nature of the DTFT and the fact that we're dealing with real-valued signals, we can rewrite this as $$X_{DT}(e^{j\omega}) = \left\{\begin{array}{ll}
\frac{1}{T_s}X_{CT}\Big(j\frac{\omega}{T_s}\Big), & 0 \leq \omega < \pi \\
\frac{1}{T_s}X_{CT}\Big(j\frac{\omega-2\pi}{T_s}\Big), & \pi \leq \omega < 2\pi \\
\end{array}\right.$$

The Discrete Fourier Transform (DFT) is an $N-$point uniformly sampled version of the DTFT, so $$X_{DFT}[k] = X_{DT}(e^{j\omega})\Big|_{\omega = 2\pi k/N}$$ for $0 \leq k\leq N-1$.

Combining the last two equations we have $$X_{DFT}[k] = \left\{\begin{array}{ll}
\frac{1}{T_s}X_{CT}\Big(j\frac{2\pi k}{NT_s}\Big), & 0 \leq k < N/2 \\
\frac{1}{T_s}X_{CT}\Big(j\frac{2\pi (k-N)}{NT_s}\Big), & N/2 \leq k \leq N-1 \\
\end{array}\right.$$

Except for the constant $N$ in your results - that I have failed to extract - I think what I wrote will help you.


impulse response of a system from its step response;confusion If the step response of a system is given as $$c(t)=(10e^{-5t}-2e^{-10t}+1)u(t)$$
  then its impulse response is _______?

My Approach:
we know: $$h(t)= \frac {d}{dt} c(t)$$
Process(1):
$$h(t)= \{ 10e^{-5t}(-5) - 2e^{-10t}(-10)+0 \} u(t)$$
$$\implies h(t)= \{ -50e^{-5t} + 20e^{-10t} \} u(t)$$
Process(2):
$$h(t)= \{ 10e^{-5t}(-5) - 2e^{-10t}(-10)+ 0 \} u(t) + (10e^{-5t}-2e^{-10t}+1) \delta (t)$$
$$\implies h(t)= \{ -50e^{-5t} + 20e^{-10t} \} u(t) + 9 \delta (t)$$

so which process is correct? please help...

HINT:

For $t\neq 0$ both solutions are identical and correct. So you only need to consider $t=0$ to find out which of the two solutions is correct. Does the step response jump at $t=0$ or is it continuous? What would a jump at $t=0$ imply for the impulse response? Alternatively, you should try to compute the step responses corresponding to the two impulse responses you obtained. Only one of them will be identical to the given step response.


Highly damped IIR/FIR digital filter I want a highly damped highpass filter (damping of at least $2$), with a cutoff somewhere around $1\textrm{ Hz}$ ($51.2\textrm{ Hz}\quad f_s$)

How do i go about designing the filter with adequate control over the damping?

My best guess was to use a standard second order response and frequency transform it, giving

$$
\frac{s^2 \omega_0}{\omega _0^2 s^2 + 2\zeta\omega_0 ^2 + 1}
$$
This gives me the required damping (after using MATLAB to transform it to the discrete domain using Tustin's method with tustin), but the frequency response is terrible (the $\rm dB$ point is much, much higher than it should be).

Is there any way to solve this?
Can I design a FIR filter with control over the damping?
Is there a higher order $s$-domain function that still gives me control over the damping?

Notes: filter order and phase delay are not overly important to my design.

Your second-order high pass transfer function is wrong. You should use

$$H(s)=\frac{s^2}{s^2+2\zeta\omega_0 s+\omega_0^2}$$

But - as mentioned in the comments - a high damping and a sharp cut-off are incompatible requirements.


partial fractions expansion inverse Z-transform, help I have the correct solution from teacher's solution guide, but I was slightly confused by some algebra about the partial fractions expansion evidently

difference equation is as follows

$ y[n] = \frac{1}{4}  \cdot \left( x[n] +x[n-1]+x[n-2]+x[n-3] \right) $

where $x[n]$ will be unit-step sequence $u[n]$

PROBLEM STATEMENT:

Find what will the output sequence y[n] when input x[n] = u[n] = unit step

We can use tactic to get X(z) and multiply X(z)*H(z) so we get Y(z)
Then inverse Z-transform to get y[k]

we know from earlier calculation that the transfer function for this system willbe as follows

$
H(z) = \frac{z^3 + z^2 + z +1}{4z^3} = \frac{1}{4} * \frac{z^3 + z^2 + z +1}{z^3}
$

because unit step input sequence

$
X(z) = \frac{z}{z-1}
$

We can find out  $Y(z) = X(z) * H(z)$

I had some difficulty with the partial fraction expansion because of the $\frac{1}{4}$ multiplier, how does that constant multiplier term become accounted into the partial fraction expansion?

We were tipped by our teacher to find out $Y(z)/z $  , so that we can find good matches in our z-transforms-table. So, I try to make  Y(z)/z into the partial fraction expansion...

$
Y(z)/z = 1/4 * \frac{z^3 + z^2 + z +1}{z^3*(z-1)}
$

So, do we just "ignore" the constant multiplier $\frac{1}{4}$ on the r.h.s of the partial fraction as follows???

$
Y(z)/z = 1/4 * \frac{z^3 + z^2 + z +1}{z^3*(z-1)} = \frac{A}{z} +  \frac{B}{z^2} +  \frac{C}{z^3} +  \frac{D}{z-1}
$

So, after the "equating coefficients technique" I obtain followins system of equations

4A + 4D = 1

-4A + 4B = 1

-4B +4C = 1

1 = -4C

from these I obtianed coefficients as follows

A= -3/2

B= -1/2

C= -1/4

D= 1

With these results we should multiply both sides by z, so we get Y(z) in partial fractions format

$
Y(z) = A + B/z + \frac{C}{z^2} + \frac{D*z}{z-1}
$

From here, we should be able to find inverse-Z-transforms from the table as follows

$
Z^{-1}(Y(z)) =  -3/4 * Z^{-1}(1) - 1/2 * Z^{-1}(1*z^{-1}) - 1/4 * Z^{-1}( 1 * z^{ -2 } ) + Z^{-1}( \frac{ z }{ z-1 } ))
$

$
y[k] = -3/4 * \delta[k] - 1/2 * \delta[k-1] - 1/4 * \delta[k-2] +u[k]
$

where $\delta[k]$, is the unit impulse, aka the deltafunction

I put that into my Texas instruments calculator and the function seemed to hold true for some initial values

expected values are that for negatives indexes output will be zero

then for k=0: output = 1/4

k=1: output = 2/4

k=2: output = 3/4

k>=3: output = 4/4

First convert your input-output relation to Z-domain to see:

$$ Y(z) = 0.25 \left( 1 + z^{-1} + z^{-2} + z^{-3} \right) X(z)$$

From which you find the transfer function of the LTI system as:

$$ H(z) = \frac{Y(z)}{X(z)} = 0.25 \left( 1 + z^{-1} + z^{-2} + z^{-3} \right) $$

Now this is an FIR system with 4 taps. As MattL. tried (without luck) to convince you that for such a system a partial fraction expansion is not the proper way to compute the output. Indeed for any FIR system, the transfer function is of the form of a $H(z) = \sum h[k]z^{-k}$ and for such a system the output is defined through a convolution as the following:

$$ Y(z) = H(z) X(z) =  \left( \sum h[k]z^{-k} \right) X(z) \implies y[n] = \sum h[k]x[n-k] $$

so for your example the proper way to describe the output for the input $x[n]= u[n]$ is

$$ \boxed{ y[n] = 0.25 ~( u[n] + u[n-1] + u[n-2] + u[n-3]) } $$

Nevertheless, assuming you want to insist on a partial fraction expansion to compute the output $y[n]$ then you would do the following:

Given $x[n] = u[n]$ and its Z-transform $X(z) = \frac{1}{1-z^{-1}} $, then the output is:

$$Y(z) = \frac{1}{4} ~\frac{ 1 + z^{-1} + z^{-2} + z^{-3} } {1 -z^{-1} } $$

Now, ignoring the linear scale factor $1/4$ to the end, since the degree of numerator polynomial is larger (in negative) than that of the denominator polynomial, you should first apply a long division of numerator polynomial into denominator polynomial to get :

$$ \frac{ 1 + z^{-1} + z^{-2} + z^{-3} } {1 -z^{-1} }  = -(3 + 2 z^{-1} + z^{-2}) + \frac{ 4} {1 -z^{-1} } $$

hence you have this (including the scale factor):
$$ Y(z) =  -(\frac{3}{4} + \frac{2}{4} z^{-1} +  \frac{1}{4} z^{-2}) + \frac{1} {1 -z^{-1} } $$

And finally you get the output from the PFE method via inverse Z-transformas:
$$ \boxed{ y[n] =  u[n] - \frac{3}{4} \delta[n] - \frac{2}{4} \delta[n-1] -  \frac{1}{4} \delta[n-2] } $$


How to implement path loss in wireless communications? I would like to implement path loss propagation model in a cellular network scenario.

I know that the path loss can be modeled as follows:

$$\mathrm{PL}(d)\propto\left(\dfrac{d_{ij}}{d_0}\right)^\alpha,$$
where $d_0$ is a reference distance and $d_{ij}$ is the distance between transmitter $i$ and receiver $j$.

My question is how to scientifically choose the value of $d_0$? I read a lot of papers that say typical values of $d_0$ are $10, 100$. I do not really get what does $d_0$ mean in real scenario?

Let say I generate randomly the position of my transmitter $i$ and receiver $j$, is there a constraint on $d_0$ that has to be met? What happens if $i$ and $j$ are very close to each others, i.e., $d_{ij} < d_0$?

Can anyone suggest to me a paper or a chapter in a book that I can read to understand the theory and the implementation of path loss models?

Thank you.

A good text is by WIlliam Y. C. Lee. He did alot of the original work on
propagation models. There are many others, just google "path loss model".
This  equation you selected is a very basic one for a cellular network. It doesn't take into account many of the fActors like environment(urban, suburban, rural), clutter or morphology, line of sight propagation or non line of sight, statistical factors and others. There are empirical models for various environments that have been used and modified (Hata, Okumura cost 231 etc) and can be iteritively adjusted with real field data.

With regard to your question about the reference distance. I'm sure you can find the answer in the references I mentioned above and the Goldsmith reference mentioned by ...

There might be a relationship to that distance and the far field or it may be justa standard reference, I'm not sure. In any case you could use a different model if the user is too close if necessary.

The path loss exponent is also a very important parameter(alpha in your equation)

Hope this helps


What would be the proper notation for denoting part of signal that satisfies a condition? I have a signal $E(t)$ which start at $t=0$ and dies out overtime. What would be a proper notation to specify part of which is above a certain threshold (say for e.g. C)?. Is $E_{interested} = E(t)>C$ a proper notation?

What you are looking for is called an Indicator Function which would assign a $1$ for all those elements of $E(t)$ that satisfy the given condition.

The same could also be expressed via the use of Set Builder Notation as $E_{interested} = \{1|e \in E \land e>C\}$ translating to "Elements of $E_{interested}$ take the value of 1 for all elements $e$ in $E$ satisfying the condition $e>C$". You might have to add $E_{interested} = \{0|e \in E\}$ as a pre-initialisation step to show the default value of $E_{interested}$.

The difference between the two is that the indicator function is defined and so you can simply reference it (e.g. "...where $\mathcal{\chi}(E, \Phi(E))$ denotes the Indicator Function using a predicate $\Phi$...", or even drop $\Phi$ if the predicate is the same throughout) while the set builder notation version could be used more as a "one-of" in-line way in writing down an algorithm.

To an extent, it is also a matter of notation already used and what you are trying to convey.

Hope this helps.

EDIT :

Just to clarify, if you are also interested in the actual values of $E(t)$ (interrupted by zeros wherever $E(t) \leq C$), then given $E_{interested}$ you can now $E_{interested} \cdot E$ to produce that sequence that contains the actual values of $E(t)$ satisfying the condition (rather than just the "pulse sequence").


Equation for impulse train as sum of complex exponentials Could someone please break down what's going on in this equation for me? I understand what the left side looks like, but not so much how the right side is the same thing.

Impulse train:

$$\sum_{m=-\infty}^\infty \delta[n-Nm]=\frac 1N\sum_{k=k_0}^{k_0+N-1}e^{ \ j 2\pi kn/N} \\ \\ n, m, N,k_0 \in \mathbb{Z}$$

where $\delta[n]$ is the Kronecker delta:

$$ \delta[n] = \begin{cases}
    1 \quad n = 0 \\
    0 \quad n \ne 0 \\
\end{cases} \\
\\
n \in \mathbb{Z}$$

The summation on the left side of your equation represents a time domain discrete time periodic signal $x[n]$ whose period is N.
$$ x[n] = \sum_{k=-\infty}^{\infty} \delta[n-kN] $$

And the summation on the right side of your equation is just the same (periodic) signal $x[n]$ represented via its Discrete Fourier Series DFS synthesis form which is an alternative representation for periodic signals
$$ x[n] = \sum_{k=0}^{N-1} a_k e^{j \frac{2 \pi}{N} k n} $$

where the DFS coefficients $a_k$ are all $1/N$ found via DFS analysis equation:
$$ a_k = \frac{1}{N} \sum_{n=0}^{N-1} x[n] e^{-j \frac{2 \pi}{N} n k} $$

I hope you can see why all $a_k$ are $1/N$ (it can be easily seen from the properties of impulse inside the summation...)


Extracting amplitudes of components from FFT I have a 4096-point FFT of a continuous pipe organ sound (no windowing).  I wish to extract the relative amplitudes of the various harmonically-related sinusoidal components for voicing analysis.  I can easily find the bins with the local maximum power for each component, but since the frequencies are not exactly in the middle of a bin, some of these peaks are spread between two adjacent bins, and the amplitude in each bin is less than it would be if the frequency were shifted to exactly hit that bin.  For frequency analysis I have used various interpolation formulas with great success.  But now I am doing amplitude analysis, and I want to know if there is a similar interpolation formula to estimate the amplitude of a component in a way that is independent of whether that component is exactly in the center of a frequency bin or is between two bins.  My current best guess is the take the sum of the powers (magnitude squared) of the bin with the peak power together with the bins immediately adjacent to it.  The sum of those three bin powers would then represent the power of that one component of the sound.  Is there a better way?

Once you use parabolic or other interpolation for a frequency peak location estimate, you can use an offset Sinc kernel, centered at that estimated frequency, to interpolate the magnitude of that peak location.  Sum each of the complex components of all the adjacent component values of the FFT result that are above the local noise floor (weighted by the offset Sinc kernel).  Then take the magnitude of that complex result.

If you use any window function (other than rectangular of FFT width) use the transform of that window function instead of a Sinc kernel.


Hidden Markov Model for Speech Rcognition. HMM Number of States This is a question that came to mind as a result of a previous question Hidden Markov Models - Distinct Observation Symbols and subsequent answer from @pichenettes.

One approach to speech recognition is to use Hidden Markov Models (HMM) to identify patterns in speech.  Both discrete HMM models and continuous HMM models have been used, but the continuous HMM approach seems to produce better results.

In continuous approach, the “input” to the HMM models is a series of vectors constructed from frames of processed speech.  So speech as processed in blocks of some time length resulting in vectors of some fixed dimension say N, with the size of N being directly related to the frame length.

My question is this:

Is there a relationship between the frame length (or alternatively the dimension of the vectors, N) and the number of states resulting in the HMM models used for recognition?

Also  is there a relationshiop between the number of HMM models required to predic speech and the size of N?

The dimension of the frame's feature vector N depends on the frequency range and the number of frequency bins in the frame. For example for 16khz audio signal it's common to take 39 feature vectors created from 40 cepstrum values. For 8khz it's enough to have 20 cepstrum values.

Feature dimension is not equivalent to the frame length. Frame length could be rather arbitrary as long as you catch sound variation in time. Frame length is determined by how fast your speech signal changes. It's enough to have 10ms frame length, though 20ms and 5ms work almost equally well. If you have frame of 30ms you will have issues since some speech sounds change faster. Thrilled r for example.

Number of distinct HMM models (tied phone detectors for example) depends on the diversity of the speech signal, on the number of distinct sounds (temporal audio patterns) there. Each distinct sound ideally has to be better modelled with it's own HMM model. The number of states depends on the size of the training database and on the vocabulary of the database, not on the feature extraction parameters. The more words you want to recognize the more HMM models you need.

There is no direct relation between both feature vector size and length of the frame and the number of states in HMM models.


Inverse Fourier transform from plots Given plots of magnitude and phase of a Fourier transform, how can I obtain the analytical formula for the signal?

I know that $$X_1(j\omega) = 3j\omega[H(\omega+3\pi)-H(\omega-3\pi)]$$ but I can't figure out a way from this point to the answer, which is $$x_1(t)=\frac{3}{\pi t^2}[3\pi t\cos(3\pi t)-\sin(3\pi t)]$$

In this case, it is quite simple. In general, it may be more difficult. I'm not sure what you are getting at with $H(j \omega)$ since you don't define it anywhere, but you are trying to get the analytical time domain expression for the signal from the plots. This is possible by noticing that
$$X_1(j \omega) = |X_1(j \omega)| \exp\left( j \angle X_1(j \omega)\right).$$
In this case, due to the plots' simplicity, we can generate a formula by reading directly from the axes. In this case,
$$X_1(j \omega) = 3 j \omega \left[ u(\omega + 3 \pi) - u(\omega - 3 \pi) \right],$$
where $u(t)$ is the heaviside step function.

Now, we can calculate the Inverse Fourier Transform with the formula,
$$ x_1(t) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} X_1(j \omega) \exp\left(j \omega t \right) d \omega. $$
Plugging in our $X_1(j \omega)$ and simplifying leads us to the expression
$$ x_1(t) = \frac{1}{2 \pi} \int_{-3 \pi}^{3 \pi} 3 j \omega \exp\left(j \omega t \right) d \omega.$$

I assume that you can evaluate this integral. If you need more assistance, let me know.


Obtaining a high sample rate digital approximation to an analog filter from lower frequency measurements Consider a set of dense, but generally irregularly-spaced frequency response measurements of some real low-pass analog filter.  Denote the maximum frequency for which a frequency response measurement is available as F_a_max.

I would like to create a digital filter model for this analog filter.  The sample rate of this filter F_s needs to be upsampled with respect to the rate implied by the measurements, i.e., F_s > 2 * F_a_max.  I have used Matlab's invfreqz function to obtain the transfer function of the digital filter from the frequency response measurements.  In principle, I can associate an arbitrarily high sample rate with the frequency response measurements.  However, I'm concerned about what the response of the digital filter will look like for frequencies greater than F_a_max.  While I don't know the exact desired response for such high frequencies (since I don't have frequency measurements for such high frequencies), I don't want the response to go "crazy" in that high frequency band.

An alternative may be using invfreqs to obtain coefficients of an analog transfer function and then use, e.g., a bilinear transform at the desired sample rate to obtain a digital filter.  However, the same concern exists.

Are there standard approaches to this problem?  I have seem some recommend Lagrange interpolation of the invfreqz output where F_s is first set to 2*F_a_max, and then Lagrange interpolation would resample the response at the higher rate.

You mention $2F_{a_{max}}$, which makes me think that you're trying to invoke the Nyquist-Shannon sampling theorem.  But that theorem doesn't deal with any "highest measured values".  It deals with that process of getting a signal from the continuous-time domain (or some higher-rate sampled domain) to a sampled-time domain.

I think if you have an existing filter that you want to measure and duplicate, that your best bet would be to reverse-engineer the analog filter transfer function, and then use any of the numerous filter design methods to design a digital filter that replicates whatever part of the filter behavior you want to replicate.

Note that unless you sample really fast compared to your highest frequency of interest, you're going to have trouble (or not manage to) get both the phase and the amplitude response to match.


How Can Telecommunication System Know Which Part Is Signal and Which Is Noise? When signal is received how can receiver know which part of received signal is signal and which part is noise??

To elaborate on MM's answer just a little bit.  Unless the receiver has some sort of expectation of the nature of the signal there is no way to tell.  With an expectation (for instance a pure tone), only a best estimation can be made with the rest assumed to be noise.  With something like serial communication a voltage past a threshold is considered a "1" and no (or low) voltage is considered a "0".  In noisy channels, parity bits (or error correcting codes) are used to ensure "what was sent".  There are many possibilities, but the receiver has to know the nature of is being sent, which seems to be the answer you are fishing for.


Signal quality metric I have floating point data signal. I apply lossy control to it and I need some kind of metric to determine its quality.

In image processing, there is PSNR (or MSE). I tried to use the same for it, It gives me some kind of result. But is there anything better (or any standard used for this kind or problem) ?

What kind of result ? Are all the floating point data positive ? Anyway, I will list down some of the Quality Metrics and when to use them for which quality of data.

$1$. SNR - If the difference between the original signal and the reconstructed signal can be interpreted as a zero-mean noise process. SNR measurements relating to signal power or energy are less usual, because the result will be biased by the mean signal value. For audio and speech, zero-mean property applies anyway.

$$
SNR[dB] = 10 \log\frac{\mathcal{E}(s(\textbf{n})-m_{s})^{2}}{\mathcal{E}(s(\textbf{n})-\tilde s(\textbf{n}))^{2}}
$$

Sensitive against phase shifts ( whereas linear phase shift are uncritical for high visual quality ). For example, an observer would hardly find any disturbance in an image shifted by half a pixel using a high quality interpolator while the SNR criterion might indicate
a poor reconstruction quality.

$2$. PSNR - For image and video signals which are strictly positive and bounded within a maximum value $A_{max}$. PSNR is widely used for image and motion distortion comparison.

$$
PSNR[dB] = 10 \log\Sigma \frac{A_{max}^{2}}{(s(\textbf{n})-\tilde s(\textbf{n}))^{2}}
$$

$3$. Weighted SNR -  A more general form of SNR which additionally allows individual position-dependent weighting of errors is the weighted SNR.

$$
PSNR[dB] = 10 \log\Sigma \frac{r(\textbf{n})w(\textbf{n})}{w(\textbf{n})(s(\textbf{n})-\tilde s(\textbf{n}))^{2}}
$$

The weight $w(\textbf{n})$ could be used dependent on local structure (e.g. different weights for edge pixels or textured areas). It has value of one for relevant and
zero for irrelevant positions and allows comparison of two images over an arbitrary-shaped area.

Then there is,

$4$ Segmental SNR -

  SNR criteria are average measures, which do not take into account  fluctuations of quality. Often, a signal with only one highly distorted segment or region will be judged to be of bad quality, while the SNR computed over the entire signal may be extremely high. As possible workaround solutions, additional criteria such as SNR measurements over smaller segments, computation of minima, maxima or other criteria which  express fluctuations have
  been proposed. The  segmental SNR, which is the mean or expected value of SNR results from equal-size segments [Jayant, Noll 1984], penalizes fluctuations in SNR.

And also there is Perceptually Weighted Criteria.


Discrete Fourier series for a sum of deltas so Im given a discrete sum $$x[n] = \sum\limits_{r=\infty}^{+\infty}\delta[n-rN]$$
how do I calculate its discrete Fourier series coefficients?

Thank you.

edit: this is what i've come to up to now:

$\sum\limits_{r=-\infty}^{\infty}\delta[n-rN]$ is a periodic Kronecker delta with period $N$, but the sum to get the fourier coefficients only "uses" one period $\sum\limits_{n = 0}^{N-1}$ so you don't have to worry about the periodic extension of $\delta$ and take just $x[n] = \delta[n]$ as your function. This is going to cancel out every exponential in the sum except when $n=0$, so the result will be independent of the particular coefficient and be just $a_k = 1/N$.


Higher order harmonics during sampling I am studying about the sampling theorem in conjunction for ADC. I got little confused while reading about aliased frequencies. I see that as per the Nyquist theorem, the sampling frequency (fs) should be larger atleast twice the signal frequency (fsig=(2*fs)) to avoid aliasing, which would position all the aliases of the form: ((mfs)+fsig) and ((mfs)-fsig) above the nyquist frequency.

Since I see that the sampling is closely related with mixing operation (sort of modulation), like - its just a multiplication of two frequencies - what happens to the alias - ((mfs)+(nfsig)) and ((mfs)-(nfsig))?

Wouldn't that come below the nyquist frequency?

For example, if fs=100 MHz, fsig=50 MHz, (fs-(3*fin))=40 MHz - analogous to Third order intermodulation product.

Normally, when measuring the ADC output spectrum, it has signal power, some harmonics and quantization noise. the Dynamic range is the difference in signal power and highest harmonic content. So, these harmonics come  from ((mfs)+(nfsig)) and ((mfs)-(nfsig))? Is this understanding correct?

The sampling is indeed analogous to mixing as to my understanding. In the sampling process, we multiply the time domain signal with an impulse train  - the impulses in time are represented as impulses in frequency at integer multiples of the sampling rate. So instead of one or two (for a real sine wave) impulses in frequency, we have an infinite number but the process of multiplying in time is otherwise identical. The result you would get (which explains the aliasing well) is the same thing you would get if you had an infinite number of mixers and LO's, one for each harmonic of the sampling clock.

I have explained this with the additional graphics below to help those more familiar with RF mixing to understand sampling and aliasing, and then undersampling as well.

The reference to third order intermodulation product is not what would explain this aliasing, as that is specifically due to non-linearities in the signal chain. This can certainly occur for the same reason in ADC's causing intermodulation distortion but that is not what causes the aliasing.  So when you see other spurious products in your spectrum, these certainly can be due to non-linear distortion (which can be confirmed by modifying the power level of your dominant signal to see if these products change), or can be caused by aliasing from other frequency bands due to insufficient front end filtering prior to ADC conversion (which means they would be present regardless of the presence of your input signal or not). Another source is spurs on the sampling clock itself, which too is explained well as a mixing process.

To understand the aliasing mechanism and how it is identical to mixing, first observe the sampling process for a 3 Hz Cosine Wave as demonstrated in the graphic below. Note that since we multiply in time the cosine wave with the time domain impulses, we convolve in frequency the 2 impulses representing the cosine wave as in the top portion of the graphic below with the impulses in frequency represented by teh middle portion of the graphic, resulting in the digital spectrum as given in the bottom graphic. In this case the signal was sampled at 20 Hz, so the output spectrum repeats every 20 Hz, so really only the spectrum from -10 Hz to +10 Hz needs to be given for the digital spectrum, as shaded in the graphic. (Or from 0 to 20 Hz, basically any 20 Hz will completely represent the output spectrum). However when working with mixed signal or multi-rate systems, I find that it often mentally helps to "un-roll" the digital spectrum and represent it out to $+/- \infty$ as I have done in this graphic.

Note too that the output spectrum can be completely explained as a mixing process: The two impulses in the middle portion of the graphic at +20 Hz and -20 Hz do represent a real sinusoid at 20 Hz, while the two at +/- 40 Hz represent a sinusoid at 40 Hz, etc... Each of the outputs shown in the digital spectrum can be explained using the traditional "sum and difference" frequency outputs that you may be accustomed to when working with an RF mixer (as explained by the multiplication of two real sinusoids).

This explains aliasing quite well, as introduced in the undersampling example shown in the graphic below (Even though I explained the above with real signals and sinusoids, I prefer to work with complex frequencies such as $e^{j\omega t}$ recognizing that a cosine is $e^({j\omega t} + e^{-j\omega t})/2$ and from that we see that each impulse in these frequency domain plots is a single instance of $e^({j\omega t}$. Then instead of dealing with sums and differences, we just add the frequency terms. For example, the graphic below shows that the frequency component that is labeled "IF" results from the time domain multiplication of the input frequency component at "RF" with the specific frequency impulse in the impulse train that is labeled "-LO".

With that view we can easily see how aliasing occurs (as explained by "mixing" operations!), and the next two graphics show the identical result between sampling in the 1st Nyquist Zone vs 3rd Nyquist zone with the first graphic colored to show which frequency components of the digital sampling spectrum are responsible for which outputs in the digital output spectrum.


Converting mel spectrogram to spectrogram I have a set of songs for which I extracted the STFT (Short-Time Fourier Transform) and used the magnitude spectrum $|S|$ to calculate the mel spectrogram by using a mel filterbank matrix $M$, so $X=\log(M\times |S|)$. I want to know is there any method to reverse this process, i.e. convert from the mel spectrogram back to the spectrogram. I performed some dimensionality reduction on the mel spectrogram, and reconstructed the mel spectrogram from lower dimensions. Now I want to regenerate the audio signal from the reconstructed mel spectrogram, so I guess first reconstruct the spectrogram and then the audio signal.

The problem is that the mel filter bank matrix is not a square matrix, since we the reduce the no of frequency bins, so inverse of $M$ cant be used like this : $ \hat{S}=M^{-1}\exp(X)$. So is there any way to generate the inverse mapping, like some inverse transfer function that can convert from $X$ to $S$?

Both taking a magnitude spectrogram and a Mel filter bank are lossy processes.  Important information needed to reconstruct the original will have been lost.  Thus you need to go back and use the original audio samples to do the reconstruction by determining a time or frequency domain filter equivalent to your dimensionality reduction.

You can make assumptions about the lost information, but those assumptions themselves usually sound inaccurate, artificial and/or robotic.  Or you can use only specially synthesized input, where the assumptions will be correct by design of that input.