Image Restoration using kalman filter I have been trying to restore an image that was blurred with a known Point Spreading Function and corrupted with noise using a kalman filter. I have looked at theory and have a basic understanding of how kalman filters work. But I can't find any good material on how I can use it to restore an image. I have tried a method where I considered that the image was the state to be predicted with no control signal and a observation model equal to the blur function with a known measurement noise i.e,

$$X(k+1) = X(k)$$
$$Z(k) = H*X(k) + R$$
where,

X(k) is the 2D image
Z(k) is the observed noisy blurred image
H is the PSF
R is the noise

here is the code is used,

I = im2double(imread('E:\Wallpapers\2.jpg'));
I = I(:,:,1);

LEN = 2;
THETA = 5;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');

SNR = 20;
noisy = awgn1(blurred,SNR);

noise = noisy - blurred;
R = cov(noise);
H = fft2(PSF,size(I,1),size(I,2));
K = 1e5*eye(size(I));
P = 1e5*eye(size(I));
X = zeros(size(I));
S = zeros(size(I));

for i = 1 : 1000
    S = R + H*P*H';
    K = P*H'/S;
    X = X + K*(noisy - H*X);
    P = (I - K*H)*P*(I - K*H)' + K*R*K';
end

MSE_noise = sum(sum((noisy-I).^2));
MSE1_restored = sum(sum((X-I).^2));

The image is filled with NaNs which is because of K = P*H'/S being badly scaled.
Where am I going wrong? Is there a problem with the code or am I supposed to change the model to restore the image?

I suggest you to investigate Image Restoration using deconvolution since it is more appropriate for your application. Kalman filtering is designed to remove or ignore noisy data in a signal but it is not suited for deblurring (as far as I know).

Noting $\epsilon$ the additive noise, $H$ the convolutional mask (your Point Spreading Function), $x$ the original image you want to retrieve and $y$ the blurred image you have:

$y = H*x + \epsilon$

Deconvolution allows you to retrieve $x$ knowing $y$ and $H$.

Using Matlab you can use the following function:

deconvreg(image, PSF)

where PSF is an estimate of the Point Spread filter and image is the blurred and noisy image.

If you want to know more about the theory, this is called an ill-posed problem that you solve by minimizing a Least Square Error through an iterative process:

$$x_{estimated} = \arg\min_{x}||y - Hx||^2$$

Under some approximations this can be computed in a fast way through FFT.

This does not answer your initial question about Kalman Filtering but I hope it will help you! If you are interested I have documents about image restoration using inverse and Wiener filtering.


change in size of an individual pixel? If we want to change resolution of an  image,we can achieve that by increasing or decreasing the number of pixels

Is there any other way ? Can we change dimensions of an indvidual pixel,rather then changing the number/count of total available pixels?

by changing size of pixel,i mean changing the width and height of array/matrix containg pixel value

A 100x100 pixel image is always a 100x100 pixel image.

You can print it with 100 DPI printer, and it will produce 1x1 inch picture.
A 50 DPI printer will produce 2x2 inch picture because it makes dots that are twice as big.

Same thing when viewing 100x100 image on 15" laptop with 1280x720 screen, it is larger because the pixels are larger than when viewing on a 15" laptop with 1920x1080 screen.

Or same thing when viewing the 100x100 image on 24" 1920x1080 screen, the image is smaller than when viewing the same image on 32" 1920x1080 screen.


Sum of equidistant exponents Consider the next sum
\begin{equation}
\sum_{k = 0}^{N - 1}e^{-j\frac{2\pi}{N}k}
\end{equation}
Its geometric meaning is the sum of uniformly distributed vectors on the unit circle.
Thus, we can say that the sum is zero.
If we try to be more precise, then for the case of even N, we can single out pairs of opposite vectors, each of which in the sum 0. What to do in the case of an odd N?
Is there a way to calculate the sum analytically, i.e. to prove this in general?

For me, only the case of an even number of members of the sum is clear. How to prove this in general terms?

ah, "sum orthogonality proof" is something you'll want to google. You can do it (under a relaxed meaning of "proof") very geometrically. Understand each part in your sum as vector in a plane:

take the first vector: $e^{0}\equiv 1$. It goes from point 0+0j to 1+0j.
attach the second vector to the tip of that vector. Note the inner angle between these two. Also note how the position of the tip is the partial sum of the first two elements from your sum.
attach the third to the tip of the second. Note that the angle between these two is the same.
Repeat until you're out of vectors.

You clearly followed a polygonal chain, where each of the $N$ edges had length 1. Now, what is the condition on the angles for such a polygonal chain to be a closed polygon?


What Is the Derivative of the Function $ s \left( t \right) = \left( 1 - {e}^{\frac{-t}{RC}} \right) u \left( t \right) $? What is the derivative of the following function

  $$ s \left( t \right) =  \left( 1 - {e}^{\frac{-t}{RC}} \right) u \left( t \right) $$

  with respect to $t$, where $u(t)$ is a unit step function?

I am getting $ \delta(t)+\frac{1}{RC}e^{-t/RC}u(t) $ as the answer. Is my answer correct?

By the Product Rule (With abuse of the derivative of Unit Step):

$$\begin{align*}
\frac{d}{dt} \left[ \left( 1 - {e}^{\frac{-t}{R C}} \right) u \left( t \right) \right] & = \frac{d}{dt} \left[ \left( 1 - {e}^{\frac{-t}{R C}} \right) \right] u \left( t \right) + \left( 1 - {e}^{\frac{-t}{R C}} \right) \frac{d}{dt} u \left( t \right) \\
& = \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) + \left( 1 - {e}^{\frac{-t}{R C}} \right) \delta \left( t \right)
\end{align*}$$

By the properties of Dirac's Delta $ f \left( t \right) \delta \left( t \right) = f \left( 0 \right) \delta \left( t \right) $ we get:

$$ \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) + \left( 1 - {e}^{\frac{-t}{R C}} \right) \delta \left( t \right) = \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) + 0 \delta \left( t \right) = \frac{1}{R C} {e}^{\frac{-t}{RC}} u \left( t \right) $$

Remark
The above is with some abuse of Math (Derivative of Unit Step as Delta Function, Zero multiplied by Delta is zero, etc...) which is used (Again, it is abuse but usually gets to the right place).
If I remember correctly, even this intuitive abuse is used in Desor.


Exponential weighted moving average time constant I think it is an easy question but I am stuck.

I want to derive that $$ \alpha = -e^{\frac{-T}{\tau}}\,. $$

Can someone provide me with an answer?

[EDIT] Sorry guys, this was really a senseless question. Here are some more informations. The differential equation of the EWMA:
$$y(n)=\alpha\cdot x(n)+(1-\alpha)y(n-1)$$
and the corresponding frequency magnitude response:
$$H_E(z)=\frac{\alpha}{1-(1-\alpha)z^{-1}}$$
and in the frequency domain for $\omega_s=1 \rightarrow T=1 $ :
$$H_E(\omega)=\frac{\alpha}{\sqrt{1-2(1-\alpha)\cos(\omega)+(1-\alpha)^2}}$$
$$\tau=\frac{1}{2\pi\omega_{3dB}}$$

Hence
$$H_E(\omega_{3dB})=\frac{\alpha}{\sqrt{1-2(1-\alpha)\cos(\omega_{3dB})+(1-\alpha)^2}}\overset!=\frac{1}{\sqrt{2}}\text.$$

I'm looking for $\alpha$, given an $\omega_{3dB}$.

If I understood you correctly, you want to compute the value of $\alpha$ that results in a specified 3dB cut-off frequency for an exponentially weighted moving average filter. If you square your last equation, you get

$$\frac{\alpha^2}{1-2(1-\alpha)\cos(\omega_c)+(1-\alpha)^2}=\frac12\tag{1}$$

which can be rearranged into the following quadratic equation:

$$\alpha^2+2\alpha(1-\cos(\omega_c))-2(1-\cos(\omega_c))=0\tag{2}$$

with the positive solution

$$\alpha = \cos(\omega_c)-1+\sqrt{\cos^2(\omega_c)-4\cos(\omega_c)+3}\tag{3} $$

So, e.g., for a desired cut-off frequency $\omega_c=0.2\pi$, you obtain from $(3)$ a value of $\alpha=0.455886780102867$. The figure below shows the magnitude of the frequency response of the resulting exponentially weighted moving average filter, from which you can see that the desired cut-off frequency is achieved.

EDIT: The formula for $\alpha$ in your question should actually be

$$\alpha=1-e^{-T/\tau},\qquad \tau=1/\Omega_c\tag{4}$$

Note that unlike $\omega_c$ in Eq. $(3)$, $\Omega_c$ in Eq. $(4)$ is not normalized by the sampling frequency. So we have $-T/\tau=-\Omega_cT=-\omega_c$.

Eq. $(4)$ is an approximation, and it comes from applying the impulse invariant transformation to the continuous-time transfer function

$$H(s)=\frac{1}{1+s\tau}\tag{5}$$

which has a 3dB cut-off frequency $\omega_c=1/\tau$. Applying the impulse invariant transformation to $(5)$ gives

$$H(z)=\frac{T}{\tau}\frac{1}{1-e^{-T/\tau}z^{-1}}\tag{6}$$

Comparing the denominator of $(6)$ with the denominator of the discrete-time transfer function of an EWMA filter

$$H(z)=\frac{\alpha}{1-(1-\alpha)z^{-1}}\tag{7}$$

results in the given formula.
Note, however, that this is only an approximation. Especially for cut-off frequencies close to Nyquist, the error of formula $(4)$ becomes relatively large.


What is the correct order of operations for a 2D Haar wavelet decomposition? The source code of iqdb contains a 2D Haar transform implementation. The author claims to have implemented it according to the paper "Fast Multiresolution Image Querying", which is freely available here.

This is the relevant text from the paper:

  A standard two-dimensional Haar wavelet decomposition of an image is very simple to code. It involves a one-dimensional decomposition on each row of the image, followed by a one-dimensional decomposition on each column of the result.

  The following pseudocode performs this one-dimensional decomposition on an array A of h elements, with h a power of two:

proc DecomposeArray (A : array[0..h-1] of color):
    A <- A / sqrt(h)
    while h > 1 do:
        h <- h/2
        for i <- 0 to h-1 do:
            A'[i] <- (A[2i] + A[2i + 1]) / sqrt(2)
            A'[h+i] <- (A[2i] - A[2i + 1]) / sqrt(2)
        end for
        A <- A'
    end while
end proc

  In the pseudocode above, the entries of A are assumed to be 3dimensional color components, each in the range [0,1]. The various arithmetic operations are performed on the separate color components individually.

  An entire r x r image T can thus be decomposed as follows:

proc DecomposeImage(T : array[0..r-1, 0..r-1] of color):
    for row <- 1 to r do:
        DecomposeArray(T[row, 0..r-1])
    end for
    for col <- 1 to r do:
        DecomposeArray(T[0..r-1, col])
    end for
end proc

(end quote)

Implementing it this way does not produce results that match the example images in the majority of the articles I have found on the internet covering this topic, including the Wikipedia article.

Note: The image is divided into 4 large squares, and (only) the top left square is further divided into 4 squares.

However, I have also found counterexamples (i.e. examples that follow the scheme used in the paper above), e.g. here.

The question is whether to loop over all rows and columns, and doing the fully recursive transform in an inner loop for each row or column, - OR - doing one pass of recursion in the outermost loop, and within each single pass only processing the remaining rows and columns.

I have implemented both approaches to demonstrate the difference visually: https://bplu4t2f.github.io/wavelet_toy/

In the approach that wikipedia uses (which I call "pass major" because the pass of recursion is the outermost loop), the emerging pattern shows that each pass divides the image into 4 squares, and only the upper left square is modified in the next pass.

In the approach that iqdb uses (which I call "pass minor"), the emerging pattern shows that only the bottom right of the 4 divisions remains unchanged during subsequent passes.

The pass minor approach feels incorrect to me, because, when looking at it intuitively, it reprocesses parts of the image that have already been transformed during each pass, effectively applying a primitive edge detection scheme recursively on previosly detected edges. It doesn't seem to make much sense to me.

Which of these approaches is correctly being referred to as 2D Haar wavelet decomposition? To both approaches have a name?

[Beginning of the story] Remember the discrete wavelet curse: in 1D, with 2-scale or dyadic wavelets, you cannot have finite support, realness, orthogonality and linear phase (symmetry/antisymmetry) at the same time, except for the Haar wavelet, which lacks of regularity and overlap. You have to lift one constraint to have the other fulfilled. For instance :

if you lift symmetry, you get Daubechies wavelets
if you lift the dyadic, 2-scale, you get FIR $M$-band filter banks, etc.

When extending discrete wavelets to 2D, many options appear. Genuine two-dimensional dyadic wavelets exist (e.g. Design of Regular Nonseparable Bi-dimensional Wavelets Using Gröbner Basis Techniques, 1998), but their lack for separability make them cumbersome (and not often used in practice).

So, most people stick to the large background of known 1D dyadic wavelet designss and apply them on rows an columns of cartesian sampled data, like images. Thus is clearly a lack of imagination, and more genuine 2D, oriented, geometrical wavelets are possible, but the SE margin is too small to wrote about it.

Meanwhile, sticking to 1D wavelet designs, the order of operations on rows and columns matter.

Two main schemes are classics, but depending on the literature, they are more or less known, and more or less used, and often under different names, so:

process all rows then all columns (or reverse): this scheme seems better known in numerical analysis and partial differential equations. It can be found under many names: separable, standard, S-form, rectangular, anisotropic, tensor, hyperbolic, separated.
process rows and columns alternatively: this scheme seems better known in image/video processing. It can be found under many names: non-separable wavelet transform, non-standard form, NS-form, square wavelet transform, isotropic wavelet transform, Mallat decomposition, isotropic wavelet transform, combined wavelet transform.

Some papers are given  here, and it is detailed in "Chapter 3. Oriented and geometrical multiscale representations" of the 2011 review paper on 2D wavelets: A panorama on Multiscale Geometric Representations, intertwining spatial, directional and frequency selectivity.

What is the best? I sure do not know, I I am still working on this. But it's geeting late, and I am attending these days a conference in honor of Alexandre Grossmann and Yves Meyer on wavelets in Paris, and I shall wake up early.

[EDIT: Added references to related questions]

2D DWT Image Issue
What is the correct order of operations for a 2D Haar wavelet decomposition?
2D DWT computation order

Let me backlet... rrr


Playing sounds with only one frequency at a time We have a "speaker", that can make a single frequency at a time (0 - 130k Hz). So we can play a tune by making it play a tone for 1 sec and then the next tone for 0.5 sec etc. We can change the length and frequency of the tones pretty quickly. But we can play any sound as they are typically a superposition of various frequencies (Fourier)... Or can we? Is there some algorithm to convert a sound signal into a sequence of sine waves.

Yes, as long as we are also allowed to vary the amplitude versus time and we're allowed any time period in which to change frequency. Any single valued waveform can also be described in terms of instantaneous frequency versus time, where the instantaneous frequency is the derivative of phase so in this case the rate of change in phase at any given point in time.
If the OP is referring to FM modulation itself (encoding the instantaneous amplitude of a real waveform to the frequency of a constant envelope sinusoid), then yes this is also done where the parameter that can be varied is the frequency deviation (linear mapping of waveform amplitude to frequency) resulting in the modulation index and ultimately the bandwidth required of the modulated signal.  This modulated signal would need to be demodulated with a frequency discriminator to recover the original waveform.


What is the maximum allowed transition band of the low-pass filter used in the reconstruction kernel for a CD with maximum speech frequency of 4 kHz? My question is if the transition band in the question below would be from 4Khz to the Nyquist frequency of 22.05Khz or from 4Khz until 44.1Khz-4Khz=40.1Khz?

If you don't want DAC noise near 40-44 kHz to alias down into your audio, you want the narrower/lower transition band reconstruction filter.


Regularization for inverse filter design Given a $2 \times 2$ matrix, $C$, suppose I want to compute a filter matrix $H = C^{-1}$ and that I need to add regularization for practical purposes (e.g., for an audio filter, regularization is needed to limit the gain of the filter at certain frequencies where the condition number of $C$ is large.) The linear system of equations being solved is $CH = I$, where $I$ is the $2 \times 2$ identity matrix.
My approach would be to use Tikhonov regularization where the solution (at each frequency) is given by $H = \left(C^* C + \beta I\right)^{-1} C^*I$, where $\beta$ is a scalar regularization parameter and * denotes taking the complex conjugate, since this minimizes the objective function $\lVert CH-I \rVert^2 + \beta\lVert H \rVert^2$. So, with an appropriately chosen $\beta$, we can essentially constrain the filter gain (since $\lVert H \rVert$ effectively corresponds to the max. gain of $H$), albeit at the expense of $H$ satisfying $CH = I$ less optimally compared to using $\beta = 0$.
Another approach might be to compute $H$ as $H = \displaystyle\frac{(\det C)^*}{|\det C|^2+\beta}\text{adj}(C)$, where $\text{adj}(C)$ corresponds to the adjugate of $C$ and $\det C$ to its determinant. Note that, for $\beta = 0$, this corresponds precisely to $C^{-1}$. This is an approach I came across recently but have not seen before.
I'm pretty sure the two solutions are not equivalent, but if they are, then a proof of the same would be appreciated. Assuming they are not, I've been trying to determine the benefits of adopting one approach over the other. For one, the latter can only be done if $C$ is square (and, from a computation standpoint, not too large since it requires explicit computation of a determinant and the adjugate). Also, $\beta$ enters the solution differently and so optimizing the value of $\beta$ might be easier in one method than the other. Apart from that, I don't have much insight regarding the differences between the two since I'm not really familiar with the latter approach.
As a response from the community, I'm looking for confirmation as to whether or not the two approaches are different, and for comments on when one might want to use one approach over the other, including a more intuitive explanation of the latter approach, if possible. Links to papers/articles/books that use the latter approach would count and be appreciated. Thanks!

This alternative approach seems fishy to me. Consider the two-dimensional case with $C$ having eigenvalues $\lambda_0, \lambda_1$ that are real and positive and $\lambda_0 \ll \lambda_1$. You have
$$
\det C = \lambda_0\cdot \lambda_1
$$
What you want to prevent with the regularisation is that your inverse gets a too big eigenvalue $\lambda_0^{-1}$, i.e. you want to cap that to at most $\beta^{-1}$. Meanwhile, $\lambda_1^{-1}$ is fine. But if you put in the corresponding eigenvector $v_1$ into that alternative formula, you get
$$\begin{align}
  Hv_1 =& \frac{\det C}{|\det C|^2 + \beta}\cdot\operatorname{adj}C\: v_1
   \\  =& \frac{\lambda_0^2\cdot\lambda_1^2}{\lambda_0^2\cdot\lambda_1^2 + \beta}\cdot C^{-1}\: v_1
   \\  =& \frac{\lambda_0^2\cdot\lambda_1}{\lambda_0^2\cdot\lambda_1^2 + \beta}\cdot v_1
   \\  \approx&
          \frac{\lambda_0^2\cdot\lambda_1}{\beta}\cdot v_1
\end{align}$$
I.e. the scaling for all inputs depends on $\beta$, not just for the low-eigenvalue ones.
I suppose this formula might still work in some situations (perhaps when the determinant is mostly imaginary?), but it certainly isn't appropriate for the general case.
Use Tikhonov.


IQ modulation: baseband to bandpass transformation can someone explain me these doubts?

it is a baseband signal representation.
A- complex signal, R- real.

If I multiply baseband signal with FT( cos) then the signal will look like this ( on the plot the first on the left /For this case, convolution changes nothing).
[![enter image description here][1]][1]

I  shift the signal in frequency but my passband signal is still complex-valued, so it cannot be transmitted over a single channel ( Is it a reason why I can.t use multiplication with cos to modulate complex baseband signal?). otherwise, I will need two channels, one for the real part and one for the imaginary part.

if I multiply baseband only with FT(cos), why I get the same copy A-R in negative part and A-R in positive? In IQ modulation it is R-A in negative part and A-R in positive

If you have a complex baseband signal

$$x(t)=x_R(t)+jx_I(t)=|x(t)|e^{j\phi(t)}\tag{1}$$

and you multiply it with a (real-valued) sinusoid, the resulting signal is obviously still complex. What happens in IQ-modulation is that you generate a real-valued band-pass signal containing the same information as $x(t)$:

$$s(t)=\text{Re}\{x(t)e^{j\omega_ct}\}=x_R(t)\cos(\omega_ct)-x_I(t)\sin(\omega_ct)\tag{2}$$

Using the magnitude and phase of $x(t)$, the bandpass signal $s(t)$ can also be written in the form

$$s(t)=|x(t)|\cos[\omega_ct+\phi(t)]\tag{3}$$

which shows that in general $s(t)$ is amplitude and phase modulated.

EDIT:

This is what happens in the frequency domain when you use a sinusoid for modulation as opposed to IQ modulation. Let

$$r(t)=x(t)\cos(\omega_ct)=x(t)\frac12\big[e^{j\omega_ct}+e^{-j\omega_ct}\big]\tag{4}$$

In the frequency domain this corresponds to

$$R(j\omega)=\frac12\big[X(j(\omega-\omega_c)+X(j(\omega+\omega_c))\big]\tag{5}$$

For the IQ modulated signal $s(t)$ we get (from $(2)$)

$$s(t)=\text{Re}\{x(t)e^{j\omega_ct}\}=\frac12\big[x(t)e^{j\omega_ct}+x^*(t)e^{-j\omega_ct}\big]\tag{6}$$

The Fourier transform of $(6)$ is

$$S(j\omega)=\frac12\big[X(j(\omega-\omega_c)+X^*(-j(\omega+\omega_c))\big]\tag{7}$$

Comparing $(5)$ and $(7)$ you see that the part of the spectrum that is centered around $-\omega_c$ is inverted for $s(t)$ (i.e., it is a mirror image of $X(j\omega)$; note the negative sign before $j$), whereas it is just shifted without inversion for the cosine modulated signal $r(t)$.


FIR filter format - Brutefir - Volumio I am trying to compensate for the high frequency droop of a TDA1543 based nos DAC. One straightforward method is described here.

Having played the resulting .wav, which I feel satisfied with, I am trying to replicate the results with Brutefir.

In essence it is a three tap filter:

  1,
  -0.15,
  0.026

The setup I am using is: Volumio + Brutefir plugin.

What I have attempted is to simply save the three numbers as .txt file and load them into Brutefir.

For some reason I am doing something wrong because there is no resulting music.
Is Brutefir expecting an even number of taps? Is the format the problem?

Yeah, I answered my own question.
So tried the same filter adding zeros to a total of 32 taps.
Plays great now and spectrum is better in the highs.
So this is the filter:

1, -0.15, 0.026, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
(Total taps: 32)


How do I find the convolution of the following signals? $h(t)$ is the inputted to convolve the $x(t)$

The signal written in unit step are:

$x(t) = (5-t)u(t-3) - (5-t)u(t-5)$

$h(t) = 2u(t-1) - 2u(t-3)$

So to convolve I first change the function to:

$h(t-\tau$$) = 2u(t-\tau$$-1) - 2u(t-\tau$$-3)$

The problem however is the signal I get after integrating the product of $x(\tau$$)$ and $h(t-\tau$$)$ doesn't even closely resemble the signal I get via $conv$ function on MATLAB.

I'm trying to do this through integration in the time domain.

If anyone can work it out and explain why, would help me understand this.

EDIT: $y(t)$ is the convoluted signal. The graph below shows one using the $conv$ function the other I manually calculated.

yt2 = (t-4).*(14-t).*u4 - (t-4).*(14-t).*u6 + (t+4).*(15-t).*u6 - (t+4).*(15-t).*u7;

Once you understand the intuition of convolution using the link in the comments, you can work on deriving the equations. These problems are a bit tedious at first and can cause confusion if you try to do too much at once. I'm uploading my hand written notes on how to break the problem down. The key is to think about what happens when there is no overlap, partial overlap when the box enters the triangle from the left, full overlap (occurs only when $t=6$), and partial overlap wen the box begins to exit the triangle.

So now we can check our work using MATLAB! We just need to be careful and understand that conv computes the convolution sum, not integral, but we can still use it to approximate. Here is the code to check your work:

%% Do the problem in MATLAB first
% Original time axis for x(t) and h(t)
t = linspace(-10, 10, 1000);

% Make x(t)
x = zeros( size(t) );
x(  t <= 5 & t >= 3  ) = linspace( 2, 0, length( find( t <= 5 & t >= 3 ) ) );

% Plot the triangle
subplot(4, 1, 1)
plot(t, x)

% Make h(t)
h = zeros( size(t) );
h( t >= 1 & t <= 3 ) = 2;

% Plot the box
subplot(4, 1, 2)
plot(t, h)

% Calculate the convolution and new time axis
y = conv(x, h);
ty = linspace(-20, 20, length(x) + length(h) - 1 );
dy = ty(2) - ty(1);

% Need to multiply by dy to get the same result as doing the integral by hand
% MATLAB conv() does discrete time sum for the convolution but we can
% approximate the integral by multiplying by the time step dy
y = y .* dy;

% Plot the convolution result
subplot(4, 1, 3)
plot(ty, y)
ylim([0 5])

%% Verify using "by-hand" equations
% Make the plot using the derived equations
t = linspace(-10, 10, 100);
y = zeros( size(t) );

tt = t( 4 <= t & t <= 6 );
y( 4 <= t & t <= 6 ) = 10*(tt - 1) - (tt - 1).^2 - 21;

tt = t( 6 < t & t <= 8 );
y( 6 <= t & t <= 8 ) = (tt - 8).^2;

subplot(4, 1, 4)
plot(t, y)
xlim([-20 20])
ylim([0 5])
```


Evaluation of incremental sensor Evaluating a volume flow incremental sensor does not deliver expected results.

The pump of which the volume flow is to be measured has a hand-calculated flow rate of

$$Q = \dfrac{n}{ 60 \;\mathrm{s/min}}   e \pi\left(\frac{d}{2}\right)^2  \approx 2 \;\mathrm{cm^3/s}$$

where $n$ is the rotational speed of $2000 \;\mathrm{1/min}$, $e$ the excentricity of the pump of $1.1 \;\mathrm{mm}$, and $\pi\left(\frac{d}{2}\right)^2$ is the area of the piston with $d = 8 \;\mathrm{mm}$.

Now the sensor (VSE, Type: VS 0,1) delivers 10000 impulses per liter.

The calculation counts the "Delta Impulses", call it dI, and divides by the time delta dt in milliseconds.

$$\dfrac{dI}{dt} = \text{impulses/ms}$$

Now I say that $10$ impulses mean $1 \;\mathrm{cm^3}$ as the spec says ($1 \;\mathrm{liter} = 1000 \;\mathrm{cm^3}$).
So
$$1 \;\text{impulse} = 0.1 \;\mathrm{cm^3/ms}$$

Further with $1000 \;\mathrm{ms} = 1\;\mathrm{s}$, it is  $100 \;\mathrm{cm^3/s}$.
I can multiply the $\frac{dI}{dt}$ with $100$ and get the flow in $\mathrm{cm^3/s}$.

$$1 \frac{dI}{dt} = 100 \;\mathrm{cm^3/s}$$

Is that correct?

Now this calculation is done within the program. The factor $100$ is given as a parameter that I have set to this value to translate the flow into $\mathrm{cm^3/s}$. The program always uses milliseconds (PLC-Control with $1\;\mathrm{ms}$ cycle, increments time ticker each time +1).

I don't know the values of the used impuls count, but the result in $\mathrm{cm^3/s}$ is about $10$ times higher than the hand-calculated flow rate. I can't find a logical mistake. Do you?

The result looks right, though you've put "1 impulse = ...", when it should be "1 impulse/ms = ..." in the 3rd equation.

Have you access to the equipment to check how fast the fluid actually flows and/or how many milliseconds there are between pulses from the sensor?


Find the length of the impulse response for the given output and input Homework Question:

Consider a signal $x[n]=\alpha e^{j \omega_{0} n}+\beta e^{j \omega_{1} n}+\gamma e^{j \omega_{2} n} .$ What is the length of impulse
response $h[n]$ of a system (non-trivial) such that $x[n] * h[n]=0$ ?

I have written the form in frequency domain in terms of system function as $$H(\omega)\big(\alpha \delta\left(w-w_{0}\right)+\beta \delta\left(w-w_{1}\right)+\gamma \delta\left(w-w_{2}\right)\big)=0$$
Thus $H(\omega)$ must have zeroes at $\omega$ = $\omega_0,\omega_1,\omega_2$
I am not sure about the meaning of impulse response length and also how to proceed to find the same from here.

You correctly figured out that the system must have $3$ zeros. So you just have to define a corresponding transfer function $H(z)$ that satisfies
$$H(z_i)=0,\qquad z_i=e^{j\omega_i},\quad i\in\{0,1,2\}$$
In other words, what is the minimum number of coefficients of a polynomial with $3$ zeros? Note that you might need to provide two solutions: one general solution, and one with the restriction that the polynomial coefficients be real-valued (as is required for a real-valued system).


Minimum Output Samples needed to crack a "Gold Code" Generator (Dual LFSR) Please preface your answer with spoiler notation by typing the following two characters first ">!"

Below is an implementation of a "Gold Code" Generator formed by adding (in $GF(2)$) the outputs from two linear feedback shift registers (LFSRs). This generates a pseudo-random sequence and is similar to the C/A Code Generator used in GPS.

An LFSR operates by shifting the contents of its shift register right on each cycle and computing the new input on the left side through the $GF(2)$ addition of certain outputs based on its "generator polynomial".

Each LFSR generator, given that is uses a generator polynomial that supports a maximum length sequence (meaning the polynomial is "primitive") produces a pseudo-random sequence which does not repeat for $2^{10}-1 = 1023$ samples, or "chips".  The combined output therefore also does not repeat for 1023 chips.  The feedback taps in each LFSR are set by the specific generator polynomial used (as shown).

Assume a black box with a a similar Gold Code generator, using two 10th order LFSRs but with unknown and primitive generators, and we only have access to the output code in high SNR conditions. What is the minimum number of chips that we would need to observe in order to determine what the generator polynomials are (in other words, to determine what the feedback taps are)?

For a single polynomial generator, the minimum number of samples is 2N+1 where N is the order of the polynomial. If noise free this would result in 10 equations with 10 unknowns, sufficient to solve the linear equations involved. For the case of the Gold Code given above, the generator shown could also be done with a single polynomial which is determined by convolving the two polynomials shown, resuting in a 20th order polynomial (it just would no longer be irreducible and primitive, but would still generate each Gold Code by changing the initial state to 20 consecutive values that are in the sequence desired). Therefore the answer is 21 chips, and the Berleykamp-Massey algorithm can quickly decipher the generator polynomial that could recreate every Gold code.


Differential Equations applications in Signal Processing Does differential and partial differential equations subject has applications in Signal Processing? I have basic DSP background and know Matrix algebra and Fourier transforms very well, but I am not sure how differential equations would fit in.

I am asking this because I have the option of taking 2 graduate level diff. equations courses(linear, non-linear and dynamical systems) and am wondering if it would be insightful or practical.

I don’t think it will hurt or will be a waste of time.

You can think of SP in a number of ways and those ways overlap. It’s a technolgy, a science,a branch of mathematics, an engineering discipline, and your view is largely tied to your self image.

A lot of SP is based on physical models.  A lot is based on hueristics.

SP overlaps. You might benefit from taking music or statistics. A lot of people take the longer path.

Differential equations are a way to describe a lot of phenomena that SP is used to sense and often control.

You should consider what motivates you.

To be honest a lot of what I took had to do with what was offered that semester. Getting that piece of paper on time was a huge practical consideration.


Find Bandwith from Spectral Efficiency and Eb/N0 to give you a short intro of why I am asking here. My Uni is offering a TeleCom Lecture, I love the topic but all they teach is electronics, the exercise sheets are still TeleCom, we get no supporting material and after a few days of Google searching I have reached maximum confusion. Send help ;)

Ok here is the Deal:

We have a BPSK Modulated Signal with spectral efficiency of p = 0.7 bit/s/Hz.
We want to Transmit Rb = 2.048 Mbit/s
In addition, BER should be 10^-6 or Eb/N0 = 10.5 dB

Task: Find the Bandwith B.

I thought that Rb/p = B  will give me the desired result.
But here I never used Eb/N0. So I assume I am missing something.

What am I missing? Or is the simple relation I found already correct and Eb/N0 is just there to confuse me?

Your initial answer is correct -- the bandwidth depends only on the spectral efficiency and the symbol rate. $E_b/N_0$ is a figure of merit largely designed to be independent of the bandwidth.

Let me give you an example of where $E_b/N_0$ might come into play. Say you are required to provide $E_b/N_0 = 10$, and you know that your receiver is rated at $N_0 = 10 ^ {-12}$. Then, $E_b = 10 N_0 = 10^{-11}$. At a rate $R = 2 \times 10^6$, the received power is $(2 \times 10^6)(10^{-11} = 2 \times 10^{-5}$ watts, or $-47$ dB.

If you know the channel attenuation, then you can calculate the transmitted power. From there, you can calculate things like battery life. And so on...


Mapping soft symbols in DQPSK I'm trying to soft decode DQPSK. According to this paper the bits are:

$b_1 = \mathrm{Re}\{s_n s^*_{n-1}\}, b_2 = \mathrm{Im}\{s_n s^*_{n-1}\}$

If I put 0.707 as a 100% probability of 0:

$\begin{eqnarray}
s_n=(0.707 + 0.707i) \nonumber \\
s_{n-1}=(0.707 + 0.707i) \nonumber
\end{eqnarray}$

This gives:

$\begin{eqnarray}
b_1=0.707 * 0.707 - (0.707 * -0.707)=0.999 \nonumber \\
b_2=0.707 * -0.707 + 0.707 * 0.707 = 0.0 \nonumber
\end{eqnarray}$

And the result doesn't make sense:

The real part (0.999) is fine. Original data corresponds to 100% probability of 0 and real part now from -0.999 to 0.999.
The imag part (0.000) is weird. It can span from -0.999 to 0.999. Being 0.000 means bit 0 has 50% probability. Which is wrong, because original $s_n$ and $s_{n-1}$ correspond to 100% probability.

So it looks like the result constellation is no longer the same as original DQPSK: (0,0) bits now mapped to (1,0) point.

Questions are:

Am I right in assumptions?
How can I get soft symbols from the result constellation?

The paper actually deals with $\pi/4$-QPSK, not with plain DQPSK. In $\pi/4$-QPSK, the symbols are always rotated by a minimum of $\pm \pi/4$. Rappaport's Wireless Communications textbook has this table:

\begin{align*}
11 \rightarrow & \, \pi/4 \\
01 \rightarrow & \, 3\pi/4 \\
00 \rightarrow & \, -3\pi/4 \\
10 \rightarrow & \, -\pi/4
\end{align*}

Unfortunately, the paper you point to does not specify what mapping it's using; the only clue is equation (2), where we see that the rotations are $n\pi/4$ for $n=1,3,5,7$.

Let's assume that $\text{00} \rightarrow 7\pi/4$. Then, in your case you would have $s_{n-1} = \sqrt{2}/2 + i\sqrt{2}/2$, and $s_n = e^{i7\pi/4}s_{n-1} = 1$. Then, according to the paper,

\begin{align*}
b_1 = & \text{Re}(s_n s_{n-1}^*) = \sqrt{2}/2 \\
b_2 = & \text{Im}(s_n s_{n-1}^*) = \sqrt{2}/2
\end{align*}

I worked out $\text{00} \rightarrow 7\pi/4$ as the only rotation that produces results that make sense. However, it's not often that I see a paper with so many typos and so poorly typeset, and I wouldn't trust it without working out the details by myself first.


Using fourier coefficients to reconstruct data in matlab When doing a discrete fourier transform on some data using matlab's fft function, its output is a set of fourier coefficients but I was wondering how do I go about converting these into an and bn so I can reconstruct the signal using sines and cosines.

E.g. I'd have a for loop that continually adds up ai cos(ix) + bisin(ix), where i = 1:N

For real signals, it works a little differently than the complex case as there is a symmetry in the DFT that can be exploited.

The DFT and inverse DFT are conventionally defined as:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i n k \frac{2\pi}{N} } $$

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{i n k \frac{2\pi}{N} } $$

Now express your bin value in cartesian coordinates (due to the IMO poor convention of using x for the signal and X for the DFT bins, this is a different x):

$$ X[k] = x_k + i y_k $$

$$ X[N-k] = x_k - i y_k $$

Note:

$$ e^{i n (N - k) \frac{2\pi}{N} } = e^{-i n k \frac{2\pi}{N} }    $$

Evaluate the kth $(k>0)$ term:

$$
\begin{aligned}
T_k &= X[k] e^{i n k \frac{2\pi}{N} } + X[N-k] e^{-i n k \frac{2\pi}{N} } \\
&=  ( x_k + i y_k )( \cos( n k \frac{2\pi}{N} ) + i \sin( n k \frac{2\pi}{N} ) ) \\
&+  ( x_k - i y_k )( \cos( n k \frac{2\pi}{N} ) - i \sin( n k \frac{2\pi}{N} ) ) \\
&=  2 x_k \cos( n k \frac{2\pi}{N} ) - 2 y_k \sin( n k \frac{2\pi}{N} ) )
\end{aligned}
$$

Put that together, with a shift of one for MATLAB indexing fiasco, and you got what you want.

For the continuous reconstruction, then it is typical to rescale the domain from $ 0 \to N $ to $ 0 \to 2\pi $.

$$ t = \frac{n}{N} 2\pi $$

$$ n = t \frac{N}{2\pi} $$

The $k$th term then becomes:

$$ T_k = 2 x_k \cos( k t ) - 2 y_k \sin( k t ) $$

So for a real valued signal, the DFT real and imaginary are half the Fourier Series coefficients.  But you still have to divide them by N if you used the conventional unnormalized DFT.

Now, when we get to the Nyquist frequency when N is even, for real valued signals, to keep the results real:

$$ T_{N/2} = x_{N/2} \cos( \frac{N}{2} t ) $$

For odd N, the last $k$ will be (N-1)/2 and nothing special needs to be done.

Now when you add them up, you still need to divide by N, so in conclusion:

$$ a_k = \frac{2}{N} x_k $$

$$ b_k = -\frac{2}{N} y_k $$

Finally the zeroth bin, $k=0$, aka the DC bin, MATLAB index of one (stupid), for a real signal is the offset you start your series with.

$$ a_0  = \frac{1}{N} x_0 $$
$$ b_0 = 0 $$

For N even:

$$ a_{N/2}  = \frac{1}{N} x_{N/2} $$
$$ b_{N/2} = 0 $$

I know this isn't the clearest or cleanest explanation, but you can figure it out.

Ced


Finding the z-transform of $h[n] = a^n\cos(2\pi \frac{n}{F_s}f_0)$ for $n ≥ 0$ and zero for $n < 0$ So I'm trying to decide whether the cosine part is intended to be plugged in for $z$ or whether it is strictly part of $h[n]$. (the number a lies in the open unit disk)

I mean I was pretty sure it was all part of $h[n]$ but then upon performing the z-transform I get this rational function

$$\frac{1 - a\cos(2\pi\frac{f_0}{F_s})z^{-1}}{1-2a\cos(2\pi\frac{f_0}{F_s})z^{-1} + a^2z^{-2}}$$

The thing is then I'm supposed to evaluate the poles and zeros and if you just ignore the cosine parts you get this really nice rational expression which factors and simplifies down to $\displaystyle\frac{z}{z-a}$.

So that has gotten me thinking that maybe I'm not understanding things correctly and the cosine portion is supposed to be plugged in for $z$ or something.  Can anyone clarify this for me?

The time domain signal (or impulse response)

$$h(n)=a^{n}\cos n\theta_0,\quad
\theta_0=2\pi\frac{f_0}{f_s},\; n\ge 0$$

is very common: it is a damped sinusoidal function (assuming $|a|<1$) which occurs frequently, because it is one possible response of a second order linear time-invariant system. So concerning your doubt, the cosine part is definitely meant to be part of the time domain signal. The poles and zeros can be found by rewriting $H(z)$:

$$\begin{align}H(z)&=\frac{1-az^{-1}\cos \theta_0}{1-2az^{-1}\cos\theta_0+a^2z^{-2}}\\&=
\frac{z(z-a\cos \theta_0)}{z^2-2az\cos\theta_0+a^2}\end{align}\tag{1}$$

From (1) it is easy to determine the zeros of $H(z)$:
$$z_{0,0}=0\quad z_{0,1}=a\cos\theta_0$$

For determing the poles, we write $H(z)$ as a partial fraction expansion:

$$H(z)=\frac{1}{2}\left [ \frac{1}{1-ae^{j\theta_0}z^{-1}}
+ \frac{1}{1-ae^{-j\theta_0}z^{-1}} \right ]\tag{2}$$

From (2) we see that the poles are given by
$$z_{\infty,0}=ae^{j\theta_0}\quad z_{\infty,1}=ae^{-j\theta_0}$$
We have complex conjugate poles because $h(n)$ is real-valued. Assuming that $h(n)$ is an impulse response, we can see from the poles that the system is stable if $|a|<1$ because then the poles are inside the unit circle of the complex plane.


Circular Convolution of h(n)*x(n) So I understand how Convolution works but Circular convolution seems a bit messy to me, all the examples seems to have even series of number like h(n)={1,2,3,4} &  x(n)={1,2,3,4} and they are all starting from n=0.

This is from my last exam and I'm trying to figure out e)
(My exam is tomorrow and I just realized I had missed to study on this, hence I won't have time to respond with Latex etc since I don't remember that and figuring it out would mean that I won't have time to study for my exam.

My attempt so far, but I don't know if I need to do zero padding (as I understand you only need to do it if you want to do the FFT). And I don't know if you are supposed to start from x(0)*h(-1)..

Grateful for all help, I have tried to search on the forum and to google it and searching on youtube.

BR Rickard

The circular convolution is actually mainly an artefact from the Fast Convolution algorithm, that is calculating the convolution via the FFT. I am not aware of any use of the circular convolution at all. Now, to efficiently calculate the convolution between x and h, you will often use the FFT, as a multiplication in the frequency domain equals a convolution in the time domain (and inversely). The algorithm is then:

conv(x,h) = ifft( fft(x) .* fft(h) )

where .* denotes an element-wise multiplication. You will have to use a zero-padding so both vectors have the same length and you can do the multiplication. This way to calculate the convolution is faster than the direct way in the time domain. However, as fft(x), fft(h) and thus also ifft(...) have the same length (here, the length of x) you cannot express the "real", linear convolution, which would be a longer vector. What happens is, that the result of the linear convolution gets "wrapped around" as you indicate in your drawings.

Let me show you that using your example. The linear convolution is (without proof)

conv(x,h) = [2, 5, 8, 8, 5, 4, 4]

Now as in the circular convolution, the vectors are only 5 long, the last 2 entries of the result will be added in the front, so

circ_conv(x,h) = [2+4, 5+4, 8, 8, 5] = [6, 9, 8, 8, 5]

is the circular convolution. The easiest way (imho) is to first calculate the linear convolution and then wrap around that result to achieve the circular convolution.

If you want to calculate the linear convolution with the FFT, you only have to zero-pad both vectors such that the linear convolution fits into these vectors, i.e. to a length of 7 in your example.


How to do DFT of this signal I am trying to get $X[k]$ when $x[n]$ is equal to

$$x[n] = \cos\left(\tfrac{\pi}{4}n-\tfrac{\pi}{4}\right)$$

I'm using this equation:

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j \frac{2 \pi}{N}kn}$$

but I'm really getting stuck.
So far I have done Euler on the cosine and I am now trying to multiply the two exponents I get from this with the exponent on the right hand side of the equation that I have included above but I cannot get it into any form that I can convert using the DFT.

Does anyone know what I can do?

A large share of such questions are addressed via known sum of series. With trigonometrics arguments like $\cos(\alpha n+\beta)$, one solution is to use Euler/De Moivre formulae:

$$\cos(\alpha n+\beta) = \frac{e^{j(\alpha n+\beta)}+e^{-j(\alpha n+\beta)}}{2}$$

and recognize, under the Fourier sum, two geometric series:

$$\sum_{n=0}^{N-1} \frac{e^{j(\alpha n+\beta)}+e^{-j(\alpha n+\beta)}}{2}e^{-j2\pi kn/N}$$
with terms like :

$$\sum_{n=0}^{N-1}\left( e^{j(\pm\alpha -j2\pi k/N)}\right)^n$$
whose sum is well-known. For a given  $k$, you have:

$$
\begin{align}
\sum_{n=0}^{N-1} e^{\pm j(\alpha n+\beta)}e^{-j2\pi kn/N} &= & e^{\pm j\beta}\sum_{n=0}^{N-1} e^{\pm j\alpha n-j2\pi kn/N}\\
&= & e^{\pm j\beta}\sum_{n=0}^{N-1} e^{j(\pm \alpha -2\pi k/N)n}\\
&= & e^{\pm j\beta}\sum_{n=0}^{N-1} \left(e^{j(\pm \alpha -2\pi k/N)}\right)^n\\
\end{align}
$$
and then, noting

$$ \rho = e^{j(\pm \alpha -2\pi k/N)}$$
you can get:

$$\sum_{n=0}^{N-1}\rho^n =\begin{cases}N \textrm{ when } \rho=1, \\\frac{1-\rho^N}{1-\rho} \textrm{ when } \rho\neq 1.\end{cases} $$


The sum part of the formula for a signals power $$P_x=\lim_{M\rightarrow \infty }\frac{1}{2M+1}\sum_{n=-M}^{M}\left | x[n] \right |^2$$

I have used this formula to find out the powers for different signals, but the sum part of the formula still puzzles me.
Lets take one example:
$$x[n] = 4i+1$$
Simply calculating the absolute value, I got this:
$$|x[n]|^2=|4i+1|^2=(\sqrt{4^2+1^2})^2=17$$
17 seems to be the answer I got after going through the rest of the formula, but I am not sure how the sum works. By using Wolframa, I arrived at:
$$\sum_{n=-M}^{M}17 =34M+17$$
Which would complete the limit part of the formula with the final result of 17, as discussed. How does this sum function, as described above?

A more general expression states that for $ M \geq N$:

$$ \sum_{n= N}^{n = M} c =  (M-N+1) \cdot c $$

where the derivation simply relies on fact that the epxression has (M-N+1) terms :

$$ \sum_{n= N}^{n = M} c = \{ c + c + ... + c\} = (M-N+1) \cdot c $$

And when applied for your particular case (with $N = -M$) it becomes:
$$ \sum_{n= -M}^{n = M} c = (M-(-M)+1) \cdot c = (2M+1) \cdot c $$


Why would amplitude modulation result in change of frequency? Question: (single choice correct; source: JEE Main 2015)

  A signal of $5 kHz$ frequency is amplitude modulated on a carrier wave of frequency $2 MHz$. The frequencies of the resultant signal is/are:

  2 MHz only
  2005 kHz, and 1995 kHz
  2005 kHz, 2000 kHz and 1995 kHz
  2000 kHz and 1995 kHz

My attempt:

I understand that amplitude modulation is the superimposition of an audio-frequency on a radio frequency carrier wave, in such a manner such that the frequency of the modulated wave is the same as that of the carrier wave, but its amplitude varies in accordance with the instantaneous amplitude of the modulating wave (message signal).

Note the part I highlighted, because, based on that part, I believe that the answer to the quoted question should be (1). However, the official answer was given as (3), and I have no clue why. It might have been correct in frequency modulation, but I am not sure why it's true in amplitude modulation as well.

Note that this is high-school level question, please answer accordingly. Many thanks!

$$A[1+\mu\cos(2\pi 5000t)]\sin(2\pi 2000000t)= A\sin(2\pi 2000000t) + \dfrac{A\mu}{2}(\sin(2\pi[2000000-5000]t)+\sin(2\pi[2000000+5000]t))$$

If I didn't make an error with my trigonometric identities.

The three distinct frequencies are there: the carrier and the positive and negative frequency version of the tone around the carrier.

The above assumes an AM that is not suppressed carrier, since the problem statement didn't say "suppressed carrier".


Is there any way to do 3D (Voxel Data) feature detection? In the 2D image processing, there are ways to find and match features of the image (like corner of the edges). However, I could not find the equivalent of that in the 3D-voxel data (Here voxel means a 3D grid with a gradation of values allocated to each box).
So, Is there any algorithm for 3D voxel feature matching?

I never personally tried it but using image registration algorithm on Voxel Data you may have some results by trying registration on sub-cluster of voxel.
I would be interested in any feedback about this idea。


What is the Reference tap of an equalizer (FIR Filter) Currently I'm using the built-in Equalizer (EQ) object from Matlab to equalize my signal. For the EQ there is a parameter called the Reference tap. With different value of the reference tap i got also different results and maybe also different delay of the EQ. My question is: what exactly is the reference tap of a FIR filter and how should i set it to achieve the best performance ? Thank you in advance LA

I haven't checked the code, but that parameter refers to the delay of the reference signal when computing the error used for updating the filter coefficients. If $y[n]$ is the output of the linear equalizer, and $r[n]$ is the reference signal (e.g., the training signal when in training mode), the error is computed as

$$e[n]=y[n]-r[n-d]$$

where $d$ is the delay defined by d=RefTap-1. This means that if the channel introduces a pure delay and the parameter RefTap is set to match that delay, the equalizer doesn't need to do anything. Otherwise, the equalizer would need to try to compensate for that delay, which is impossible for a causal filter (because it would try to implement an advance in time). So that parameter is indeed very important in order to achieve the best possible equalization for a given equalizer filter length. Since in practice the channel delay is often unknown, a reasonable choice is to set the parameter RefTap to half the equalizer length. Then, on average, the best use is made of the equalizer's weights. If the performance is not satisfactory, the number of weights needs to be increased (and the parameter RefTap must be changed accordingly).


Relation between Fourier phase and Analytic Signal Phase Do you know if there is a relation between the
phase of the Fourier Transform and the phase of the Analytic signal?
I mean, if i have a signal s(t)  then...
FT_s=FFT(s(t))
phaseFourier=atan(Imag(FT_s),Real(FT_s))

A_s=s(t)+iH(s(t))
phaseAnalytic=atan(Imag(A_s),Real(A_s))

Then...
phaseFourier= ?(phaseAnalytic)

ThankYou

Discussing the relationship between the Fourier transform and the Hilbert transform, let us consider an example function. It can be any square integrable function; for readability of the text, we select a function with both transforms expressible via simple analytical functions.
Let a signal s be a sinc function of time, s(t) = sinc(t). The Fourier transform of the sinc function is a rect function
$$
\mathscr{F}\{s(t)\}(ω)={{1\over{\sqrt{2π}}}·rect\left({ω\over{2π}}\right)}
$$
which is a function of frequency: OP's FT_s is a function of frequency, or, for the discrete case, an array with indices representing frequency bins.
The Hilbert transform of the sinc function
$$
\mathscr{H}\{s(t)\}(t)={{1-cos(t)}\over{t}}
$$
is the function of time, because the Hilbert transform is not a transform of the kind that maps time domain functions to frequency domain functions -- as the Fourier transform does -- but "a specific linear operator that takes a function, u(t) of a real variable and produces another function of a real variable H(u)(t)" (citation from the Wikipedia article on the Hilbert transform).
Comparing these formulas, one sees that the question "if there is a relation between the phase of the Fourier Transform and the phase of the Analytic signal?" is rather pointless.
FT_s, as defined by OP, is an array of values against frequency bins, and atan(Imag(FT_s),Real(FT_s)) is also an array indexed by the same bin enumerables, corresponding to a continuous-time frequency. Therefore, atan(Imag(FT_s),Real(FT_s)) is an array of phases of the signal frequency components.
A_s is a function of time, and atan(Imag(A_s),Real(A_s)) is also the function of time. It is the instantaneous phase of an analytic signal, and it is not the phase of the Fourier transform frequency component, not in any way.
The OP-defined two continuous-variable "phase functions" and corresponding discrete "phase arrays", phaseFourier and phaseAnalytic, are defined on different domains: phaseFourier is defined on a frequency domain, phaseAnalyticis defined on a time domain. Still, you may want to see the graphs of these function for a sinc(t) signal.
For s(t) = sinc(t), the Fourier transform (OP's FT_s) values is a real constant up to a frequency of 1/(2π) and a zero for higher frequencies. Therefore, the array of phases of frequency components (OP's phaseFourier) is zero up to 1/(2π) and undefined in the region where the Fourier transform amplitude is zero. For sinc(t), one can safely assign all phases of FT components to zero in the entire frequency domain. The FT_s graph is trivial, and shown in every textbook and course notes on signal processing.
The instantaneous phase (OP's phaseAnalytic) graph may be shown less frequently, so here you are:

Summing up, one may agree that, talking about the relationship between the Fourier transform and the Hilbert transform, more appropriate is to indicate that the Hilbert transform is a multiplier operator. The multiplier of H is
$$
σ_{\mathscr{H}}(ω)=-i·sgn(ω)
$$
(see https://en.wikipedia.org/wiki/Hilbert_transform#Relationship_with_the_Fourier_transform).


Behavior of tanh IIR filters If we insert a tanh function (or any other activation function) between the feedback summation and the unit delays, how will such an IIR filter behave for values of $|a| < $1 and $a > ±1$ ?

The filter becomes stable for all values of $a$ but what happens to the frequency response of the filter?

A transfer function is only defined for linear time variant systems. Since your proposed structure is non-linear, the concept of a transfer function doesn't apply any more. You could make approximations for certain special cases and assumptions, but overall the answer to "what happens to the frequency response" is "undefined".


Does the Quality of the Pictures Taken by a USB 3.0 Camera Depend on Cable Length? I've got a JAI GO-5000-USB camera and have tried several cable lengths, ranging from 40cm to 5m.

So far, I haven't observed many differences in the quality of the pictures I've taken. Since I want to avoid a later realization that the quality isn't that good after all, I'm hoping that someone with experience with USB 3.0 cameras and longer cable lengths could warn me in advance.

Now my question is:
Has anyone experienced a decreasing quality of taken pictures depending on longer cable lengths? When did the problems start and in which way did they appear? Is it possible to measure or visualize those effects?

I've tried to visualize them with histograms, but that didn't really show anything useful at all.

No, the data in the USB cable is digital.
Namely there is either an error in the data or the data is the same at any place along the cable.


Highest frequency of collected force signal I have collected force at the footrest during paddling at two occasions. During one occasion I accidently collected the force signals at 150 Hz instead of 1500 Hz. During the other time the data was collected at 1500 Hz. I now want to calculate the highest frequency of the signal collected at 1500 Hz to be able to know if the data collected at 150 Hz is within the Nyquist criterion.
How do I calculate the highest frequency in my signal?

In order to determine the spectral content of your signal you can compute and display its Fourier transform. Assuming you have $N$ samples of the signal $x[n]$ sampled at $F_s = 1500$ Hz, then you can simply determine its frequency content by the following Matlab command:

Xk = fft(x ,N);

Then you can plot its spectrum to see if there is significand energy in the frequency region above $150/2 = 75$ Hz by the following code:

figure, plot ( linspace(0,Fs,N) , abs(Xk) ) ;
xlabel('Frequency [Hz]')
title('DFT Spectrum Magnitude |X[k]| of x[n]');

If there is any signifcand energy above $75$ Hz then there will be significand aliasing. However note that the sensor might have applied a $75$ Hz lowpass filter before sampling the signal, so that you can still avoid alising, but your signal will be distorted, nevertheless, as it would have beeen sampled from a lowpass filtered version of the true signal.


Calculate input impedance from S-parameters I wrote this simple matlab code to calculate input impedance of a low pass filter but not sure whether all the steps are correct. * First attempt

z0           = 50;

A            = importdata("sparameters_lowpass.dat");

freq         = A(:,1);
realPart     = A(:,2);
imagPart     = A(:,3);

s11_freq     = realPart + imagPart*i; % In frequency domain

%X           = fft(s11_freq, N); % Not sure whether this is the right way to interpolate

X            = s11_freq;

X_conj       = conj(flip(X));

X_input      = [X_conj(1:end-1)', X(1), X(2:end)']';

M            = length(X_input);

window       = hamming(M);

s11_time     = ifft(X_input.*window, 1024);

z_in         = z0*(1+ (s11_time))./(1- (s11_time));

S11 is in the frequency domain, and the impedance should also be in the frequency domain. If you want to resample the result, simply use the resample command:

z0           = 50;
A            = importdata("sparameters_lowpass.dat");
N            = 5000;      % number of samples in final result desired

freq         = A(:,1);
realPart     = A(:,2);
imagPart     = A(:,3);

s11_freq     = realPart + imagPart*i; % In frequency domain

z_in         = z0*(1+ (s11_freq))./(1- (s11_freq));

freq_resamp = resample(freq, N)
z_in_resamp = resample(z_in, N)

Note that your formula for S11 would be correct when it is given as a complex signal with a magnitude between 0 and 1 and any phase. Just a word of caution in case the result is in dB magnitude and phase (I don't think that is the case given the extraction of the real and imaginary part shown but that should be confirmed).


fundamental frequency contour of a syllable I have an audio signal and I want to calculate the fundamental frequency contour for each syllable. In the related paper is said that only 9-time points are required. (9 values of frequency contour in each syllable).

As I saw in some papers, calculating fundamental frequency contour is equal to pitch contour.

Can someone please explain it to me what is it exactly?

Thanks in advance

Syllables might require a certain duration to be recognizable by humans.  Over part of that time, voiced syllables might have a pitch, or a clear periodicity in waveform.  Pitch can be described as having a fundamental frequency or F0 (the reciprocal of the perceived repetition period). Note that there a lots of other answers ont this Q&A site on various ways to measure pitch, and the potential problems related to such.

Pitch can change over time.  If you measure pitch at several points in time (say 9 points for example) over the duration where a syllable is voiced, you might find that the measured pitch in fact has changed (or not).  With those pitch measurements over time one could interpolate a curve (spline or polynomial regression, etc.) thru those measured points.  You might end up making a bunch of different looking curves from different snippets of audio.  One could also publish nice looking graphs of those interpolated curves and call the shapes of those curves contours.

Feed a lot of those nice looking pitch/F0 frequency contours (plus a large amount of other hopefully related data) to a machine learning algorithm, and it might be able to infer some possibly interesting correlations between things.

Or not.


References on the MP3 encoding algorithm Does anybody have a reference describing the MP3 algorithm in a step-by-step manner and clear to understand?

These references were a little bit confusing:

Audio compression using modified discrete cosine transform,The MP3 Coding Standard, Joebert S. Jacaba, 2001
The Theory Behind Mp3, Rassol Raissi, 2002

I am willing to encode a bitstream to MP3 after reading the reference that you introduce, manually, for a complete understanding of the procedure.

A compression standard is a quite delicate thing, that took years to develop and tune. I suggest Analysis of the MPEG-1 Layer III (MP3) Algorithm Using MATLAB, 2012 by Jayaraman J. Thiagarajan and Andreas Spanias. After a general overview, it leads you step by step though the different blocks of the whole scheme, providing some Matlab code for each of them.


Signal decomposition: getting phase with known amplitude and frequency of main components Suppose you have sampled a signal and you know amplitude $ \alpha_i $ and frequency $ f_i $ of the main components for the decomposition:

$$
y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)}
$$

Which methods would be the fastest and the most accurate to calculate the phase vector $ \phi_i $?

Boundaries:

Signal sample has constant sampling rate
Signal sample length is not guaranteed to be an integer multiple of the signal period, but is longer than one period

EDIT

The simples way I can think of is to minimise the RMS error with some optimisation algorithm as shown in code below. I’m wondering if there is a more efficient and simpler way to achieve the same.

%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as op

# Build a signal as sum of cosines
def f_signal(t, magnitude=[], frequency=[], phase=[]):
    y = 0
    for m, f, p in zip(magnitude, frequency, phase):
        y += m*np.cos(2*np.pi*f*t+p/180.*np.pi)
    return y

# Length [s] and number of samples of the signal
t_length, N = 0.09, 1000

# Equally spaced time axis
t = np.linspace(0, t_length, N)

# This is the reference signal (we know all information)
y = f_signal(t, [8, 4], [80, 240], [-90, 40])

# Now suppose we lost the phase information and want to reconstruct it

# Lets use the root mean squared error as function to be minimised
def f_rms_error(phase):
    return np.sqrt(np.sum((y - f_signal(t, [8, 4], [80, 240], phase))**2)/len(t))

# Initial guess
phase0 = [0, 0]

# Minimize RMS error with fmin
phase1 = op.fmin(f_rms_error, phase0)

# Print result
print(phase1)

# Reconstruct signal
y_r = f(t, phase=phase1)

# Plot results
plt.figure(1, figsize=[15, 5])
plt.subplot(1,1,1)
plt.plot(t,y)
plt.plot(t,y_r, 'r-‘)

Trying to estimate the phase directly from the following equation requires nonlinear optimization.
$$
y(t) = \sum_{i=0}^N{\alpha_i\cos(2\pi f_i t + \phi_i)}
$$
Since the amplitudes and frequencies are known you can transform your problem to
$$
y(t) = \sum_{i=0}^N{\cos(\phi_i) \alpha_i \cos(2\pi f_i t)} - \sum_{i=0}^N{ \sin(\phi_i) \alpha_i  \sin(2\pi f_i t)}
$$
If you have $M$ samples for $t=t_0,... t_{M-1}$ and define
 $$p_i=\alpha_i \cos(\phi_i)$$ $$q_i=- \alpha_i \sin(\phi_i)$$ as your new parameters. Then you have a classic least square problem with $2N$ parameters.
$$\bf{Y} = A X + \text{noise} $$
where
$${\bf{Y}}=[y(t_0),....,y(t_{M-1})]^T$$

$A$ is a $M\times (2N)$ matrix where
$$A(j,i)=\alpha_i \cos(2\pi f_i t_{j}) ~~~\text{for}~~~ i=0,..,N-1$$
$$A(j,i)=\alpha_i \sin(2\pi f_i t_{j}) ~~~\text{for}~~~ i=N,..,2N-1$$
and
$${\bf{X}}=[p_0,....,p_{N-1}, q_0,....,q_{N-1}]^T$$
The unknown vector $X$ (and thus all $p_i$ and $q_i$) can be estimated by normal equation (see for https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)):

$$X=(A^{T}A)^{-1}A^{T}Y$$

After estimating $p_i$ and $q_i$ you can estimate $\phi_i$, with an ambiguity of $\pm 2\pi k$ by taking 4 quadrant inverse $\tan$:
$$\phi_i=\text{atan2}(-q_i,p_i)$$

Also the amplitude $a_i$ can be estimated:
$$a_i =\sqrt{p^2_i+q^2_i}$$


What is the advantage of measuring IQ data non-digitally? So I know In-phase and Quadrature (IQ) data can be used to analyze changes in both phase and amplitude of the signal. I also know that if I know the I component of the signal (which I could get by measuring it with one oscillator), I could digitally get the Q component using Hilbert Transform. However, I know there are devices (like RTL-SDR sticks) that have two oscillators, 90 degrees out of phase of one another, so they can measure both the I and Q data analogically.

Why go through the trouble, engineering wise, of building such a device when one could always use a computer to generate the Q component by Hilbert transform from the I data? What is the advantage of measuring both?

Thank you.

IQ demodulation is normally done with the "I" component being sampled at a rate that is at or near (including below) the carrier frequency of the signal, or the center of the band of interest.  Since this is well below the sample rate required to meet the Nyquist criterion (e.g. sampling at a rate greater than 2X the highest frequency of the signal spectrum), using the "I" component samples alone at that 1X  rate won't work, as you will get aliasing (likely scrambling together both positive and negative frequency spectrum above and below the sample rate into the one channel of strictly real I.F. or baseband data samples).  This will alias (make a mess of) your spectrum before you can even try to do any sort of Hilbert transform filtering.  (Unless you have a really steep filter before sampling that does the split into just USB or LSB for you.)  Sampling at exactly 2X the carrier has a similar problem.  Using the actual sampled Q component data adds information that properly keeps the positive and negative spectrum (the spectrum both above and below your "I" component sample clock) separated in the IQ I.F. or baseband samples.

IQ sampling can be done using a sample clock (oscillator) that isn't at a frequency much higher (e.g. over 2X higher) than the frequency band of interest.  Instead, only a 90 degree phase delay is required at that much lower sample clock rate (Possibly also simplifying the required filtering as well as the sample clock generation.)


What is the relation between the domain of Fractional Fourier transform and time and frequency? This refers to the link where it clearly writes that:

  The FrFT argument $u$ is neither a spatial one $x$ nor a frequency
  $\xi$. We will see why it can be interpreted as linear combination of
  both coordinates $(x,\xi)$

But I did not see the relation between the domain of Fractional Fourier transform ($u$) and space and frequency $(x,\xi)$. I can understand that domain '$u$' is a linear combination of '$(x,\xi)$' but my question what is the relation? Can anyone point out the exact relation?

It think it is more correct to say that the $u$ axis is a linear combination of the $x$ and $\xi$ axes:

$$(x, \xi) = (1, 0)x + (0, 1)\xi = \left((1, 0)\cos(\alpha) + (0, 1)\sin(\alpha)\right)u,$$

that is:

$$ x = u\cos(\alpha)\\
\xi = u\sin(\alpha).$$

The definition of the fractional Fourier transform from that Wikipedia page is:

  $$\mathcal{F}_\alpha[f](u) =  \sqrt{1-i\cot(\alpha)} e^{i \pi
\cot(\alpha) u^2}  \int_{-\infty}^\infty  e^{-i2\pi (\csc(\alpha) u x
- \frac{\cot(\alpha)}{2} x^2)} f(x)\, \mathrm{d}x.$$

If you use another variable, $y$, for the integral, and substitute $u = \sqrt{x^2 + \xi^2}$ and $\alpha = \text{atan}(\xi, x)$, then the fractional Fourier transform can be written as:

$$\mathcal{F}[f](x, \xi) = e^{i \pi \frac{\xi}{x} (x^2+\xi^2)} \sqrt{1-\frac{i \xi}{x}}\int^\infty_{-\infty}e^{i \pi \frac{y}{x} \left(\xi y-2 x^{\frac{1}{\xi}} \sqrt{\xi^2+x^2}\right)} f(y)\, \text{d}y.$$


How to obtain the probability of a detection? I will use a specific example from image processing to illustrate my question, but I'm actually interested in the higher level / abstract procedure. I must lack the specific vocabulary and my hunch is that my question relates to some topic in math / signal processing / computer science that I'm not yet aware of.

I extracted part of an image and I want to check if it contains an edge.
Using some filtering technique (e.g. Sobel) I obtain pixels which are more likely to be part of an edge.

What I would like to obtain is the likelihood of the image to contain an edge.
For example: The probability that this image contains an edge is p=0.82 .

Q1: Which processing steps am I missing here?

My intuitive approach is the following:
Let's assume I want to purchase a magazine. I know that the average price is 4€. By checking how much money I carry, I can predict the likelihood that I can complete the purchase. By example:

With 1€, the likelihood is low, maybe 0.1
With 4€, the likelihood is average, 0.5
With 50€, I can afford almost any magazine, 0.99

Now, I drew the probabilities and created a PDF (probability density function) out of nowhere. Knowing the actual repartition of magazine prices would help me create a more realistic PDF.

Q2: Is this intuitive approach correct?

Q3: How does it apply to the aforementioned problem in image processing?

Q4: Which topic/keywords is my question actually about? Can you recommend some sources / references to develop my knowledge?

You have a fairly common problem in Detection theory. Let's assume you have a $N \times N$ grid of pixel values and $N$ is small enough that the edge is a straight line in the region.  For that given line, you can formulate a hypothesis on whether and edge was on that line or not. You might also have to some other information like the amplitude of the edge and distribution of the non edge pixels.  Unfortunately, you usually don't know where the lines are, so one thing people do is something called a Generalized Likelihood Ratio GLR test.  So, if you can enumerate all the possible lines that can fall in your $N \times N$ grid and test each, you can detect it, assuming you can estimate the other parameters.  Unfortunately, that's a big search space and probably not a good approach. Even if you could do the search, since pixels will be contained in more than one search, the tests will not be independent and a closed form solution for a Detection probability will be difficult if not intractable to calculate. As an alternative, you  can do an $M$ out N test, test if there are $n$ out of $N \times N$ pixels that are part of a line.  We don't know $n$ either but if there is a single line $n$ is bound by $\sqrt{2} N$.  Again we have nonindependent tests and a difficult calculation, and the pixels that are high aren't necessarily connected as in a line.

A more practical approach is to choose $H_0$ to have no edges, and assuming we have a reasonable sampling distribution, we would look for small $p$ values. The original likelihood we mentioned would be more powerful and often people use the ideal case as an upper bound.

You might try a nonparametric test like Wilcoxon rank test on pairs of blocks.

Radon transforms are good at finding lines, so that is a possible approach as well. Picking $N$ essentially a heuristic as well.

Honestly, it all gives me a headache and I try to stick to time series.

Hopefully someone who knows better will give you a better answer


Bandpass filter BPSK I have a question concerning BPSK and the bandpass filter that seems to be sometimes placed before the receiver circuit (the receiver is a mixer and lowpass filter).

Basically, from what I've noticed, there is sometimes a bandpass filter before the receiver of BPSK (and other types of receivers too I think).  Is the only purpose of this bandpass filter to band limit the signal and basically just improve the carrier to noise ratio?  It doesn't actually need to be there does it?  It's optional right?

Thanks.

The filter that you're referring to is called a preselection filter. Its purpose is to filter out everything but the desired signal of interest before mixing to baseband. Unwanted components could include other signals that are nearby in frequency, or just noise that lies outside the desired signal's bandwidth.

Preselection can serve multiple purposes:

It is common to have a low-noise analog amplifier in your receiver chain to reduce the receiver's overall noise figure. Bandlimiting the input of the amplifier using a preselector can result in improved amplifier performance in some cases, especially those where there are strong interfering signals in nearby bands.
In cases where bandpass sampling is used (i.e. the desired signal is sampled directly at its nonzero intermediate frequency), a bandpass filter will typically be required in order to remove other signal or noise components that would otherwise alias into the sampled signal bandwidth.

Whether the filter is optional or not would be application-dependent.


Could a DCT be used for an audio magnitude spectrum rather than DFT? From what I understand, the DCT has half the bin size as a DFT of the same size N.  The DFT also includes phase information, but often this is not needed when only the magnitude spectrum is desired.

Could the DCT be used to provide a magnitude spectrum with twice the density (half the bin spacing) of the DFT or would out of phase information be lost?
How about with a 50% overlap?

Yes, DCT can be used to provide a magnitude spectrum with twice the density. I do not quite understand overlap, but I am assuming that since DCT covers less, you thought there would be an overlap. To provide an eligible answer to the question, let me make a quick review for usage of DCT in mainly image processing.

First, we need to make some assumptions. In order to use DCT, you need to have a real signal. This is by definition. While you are saying, DCT has half the bin size comparing to DFT in size N, you are assuming that the signal is low frequency signal. Otherwise, not so much.

For usage of DCT in compression, since DFT of image will be symmetric, it produces redundant information (one side mirror will be enough to reproduce the signal). Therefore, kernel of DCT is used in order to produce denser information comparing to DFT. This is also true for low frequency audio signals, it can be used in the same way. While it makes it denser, coefficients gets bigger, since kernel of DCT covers both sides(real and imaginary parts) of the signal.

My major is image processing, so I tried to map DCT and DFT concepts and explanations in image processing. One difference between image and audio could be sizes, though. In image processing, you know the sizes(row and columns for FFT and other purpose of processing). I guess that you need to divide the vector of audio data somehow in order to further process. Without knowing the data, this could be troublesome(I am not sure).

Here is an image taken from web, but I did not write it down where I took it, could be wikipedia.;

As you can see, transformed image is represented in DCT by magnitude spectrum with no problem. In a more compact and denser way, and look at the magnitude of coefficients. It is bigger than two times of DFT. DFT is symmetric, you could just divide it into two. One part is redundant. And one more thing, DCT can store the information is not just half of DFT but nearly quarter of DFT. That is generally the case of DCT overcoming to DFT in images.


Line coding for a modulation order that is not power of 2 In case a symbol encodes a non-integer number of bits, how is done the bit-symbol pairing?

For example, we have a 3-levels FSK with bit rate of 1Mbps and a sequence 01011001. What is the symbol period? What is the sequence of symbols?

Note: My guess is that the actual bit period and symbol period are the same. Hence, the use of 3-levels only helps to increase somehow the robustness of the error detection (similar to Hamming coding or 8b/10b encoding). Thus, we apply a bipolar AMI (Alternate Mark Inversion) to the bit sequence to obtain {0,+1,0,-1,+1,0,0,-1}. Finally, in the receiver, we demodulate every symbol 0 as a bit 0, and every symbol +1 and -1 as a bit 1.

Another reason to encode to $3$ level is that the occupied bandwidth is less.   This is often driven by an emission mask constraint  for transmitter certification. To satisfy this constraint,  the tradeoff is made on symbol distance.   It is however a costly tradeoff affecting receiver sensitivity because the multi level scheme is done post detection where every factor of $2$ increase in throughput costs $6$ dB in sensitivity.   Other schemes like QAM cost only $3$ dB in sensitivity when doubling the number of bits per signaling element.  With the multilevel FSK, you do however retain some of the advantages of an FM strategy in  terms of the CIR performance of the system.


Using ICA on EEG signals for feature extraction I am attempting to use ICA (FastICA via scikit-learn) on EEG signals from seven electrodes per subject for feature extraction and identity classification – that is to extract signal which is related to a specific cortical area without noise and feed it to e.g. MLP for classification.

As I currently understand it, I should be able to get at most seven individual components which I seem to be achieving, but I do not understand how to actually choose the one related to a specific area nor how to verify the results. How do I know I am not extracting noise (they are all very similar)?

The data have been pre-processed with a Butterworth bandpass filter (4 – 40 Hz). The mean error for classification of six subjects via MLP is currently around 12% but similar studies achieved far better results.

See for example:

Selecting Relevant EEG Signal Locations for Personal Identiﬁcation using ica and neural network
Personal Identification by EEG Using ICA and Neural Network

Excuse the lameness with this one, please – the research papers and introductions to ICA I have read are not brilliant at explaining and use examples involving ideal components.

Anyone?

How do I know I am not extracting noise (they are all very similar)?

You should first find or record a well understood set of signals with which you can test your feature extraction code to verify your setup is correct and that the extracted components are as expected.

Then you'll have higher confidence that you're not extracting noise due to errors in your code and you can focus on tuning your algorithms or reproducing other's results.


Discrete or continuous Kalman filter? I have position and acceleration measurements and I would like to apply a Kalman filter to estimate the velocity of the system.

I am not sure yet about how to procede, but I will check the already answered questions on this website (like Estimating velocity from known position and acceleration , Kalman filter with accelerometer with DC offset and Applying Kalman filter to a data set).

Once I have understood how to procede, I will implement it using Matlab. There I saw there are 2 different types of Kalman filter: discrete and continuous. What is the difference? I mean, what is the difference of working with one or the other?

EDIT to be more clear, I am referring to function kalman and kalmd

Since I work with a set of data, should I use the discrete one?

I don't have experience with continuous time kalman filters.  However from your description it sounds like you are making measurements of position and acceleration over time.  If these measurements are sampled, meaning you have individual measurements associated with time, you should be using the discrete time kalman filter.  I feel comfortable broadly stating if you are trying to implement the filter on a pc (matlab), microcontroller, fpga, or dsp based on a set of measurements you are implementing the discrete time filter.

EDIT

Based on your edits, I more clearly understand your question.  I don't have experience with Matlab's built-in kalman filter functions but a quick read of the comments in kalmd seem to indicate to me you want to use kalman and not kalmd.  Below is a snippet pasted from kalmd that should make it clear.

%   The LTI system SYS specifies the plant data (A,[B G],C,[D 0]).
%   The continuous plant and covariance matrices (Q,R) are first
%   discretized using the sample time Ts and zero-order hold
%   approximation, and the discrete Kalman estimator for the
%   resulting discrete plant is then calculated with KALMAN.

This text tells me the function is expecting a continuous function which it then evaluates/discretizes, and then calls the kalman function to create a discrete filter.  From the sounds of it you already have discrete data (measurements over time).  Hence kalman is what you want to use.


Energy of compressed signals I have tried a code to compress a signal using Compressed Sensing(CS).
  The input signal is $x$ and the compressed signal  $y$ is given by :
 $y=Φ*x$   where $Φ$ is the sensing matrix.
I have used the following matlab function to compute the energy of $x$ and the energy of $y$ .

  energy_x= sum(x(1:384,:).^2);
  energy_y = sum(y(1:192,:).^2);

The length of $x$ is 384 , and the length of $y$ is 192.
I found that the energy of $x$ is 446910
  and the energy of $y$ is 77651282 .
Is it is reasonable to obtain higher energy for the output than the input?

Without information on  $Φ$, you can obtain almost anything, since $\lambda Φ$ could be a valid CS matrix as well. Generally, one imposes structure contraints, such as unit energy for their rows or columns.

This being said, compressive sensing does not compress data in a strict sense, and energy is a poor measure of compressibility. Entropy and norm ratios could be more interesting measures.


Bias instability in Gyroscopes : AVER / ADEV I am attempting to implement allan variance to determine the amount of bias that is shown in my IMU's gyroscope ( BMI160 ). I am doing this with the following formula, that was captured in the following tutorials: Tutorial 1 and Tutorial 2

Formula for AVAR:

$\sigma_y^2(\tau) = \frac{1}{2(N-1)} \sum_{i=0}^{N-1}(y(\tau)_{i+1} - y(\tau)_i)^2$

From my reading, I am implementing the calculation using the following statements:

$\tau$ - Bin of time ( 5 seconds of data, 10 seconds, 15 sec, etc ) So if my $\tau$ is $5$ sec, then I'd multiply $5$ by my sample rate to get the number of sample points that I can stuff into a bin. ( $5*40=200$, given that my sample rate is $40$ )

$N$ - Total number of bins that I've split my data into. ( Number of sample points per bin divided by total number of data points, without the remainder )

$y(\tau)$ - The average ( mean ) of the data in some bin of time

I then plot the ADEV ( Sqrt of AVAR ), over $\tau$, and see the following result:

What I see in that plot is that different axis have different bias instability levels. I did not collect hours of data to let the graph start rising again, like a true allan variance plot would. But for the purpose of this question, let's assume that the graph starts rising for each axis after the point i've selected on the graph. ( in the image )

Questions:

When a manufacturer releases a datasheet, and they say that the gyro bias instability is "$0.55^{\circ} / hr$ ", are the taking the average of all three axis to come up with one value? Looking at my plot, that would be rather poor, as two of the axis are fairly low as compared to the third.
Does it make sense to say that you'd want to see a low bias instability, and for it to be identified through a very small $\tau$? i.e Manufacturers claim $0.55^{\circ} / hr$ but rarely in a datasheet do i see the result of $\tau$ being shown of where that value was found. A sensor with $0.55^{\circ} / hr$ at a $\tau$ of $5$ should have significantly less bias than the one with the $\tau$ of $25$. No?

According to this app note by Freescale: http://cache.freescale.com/files/sensors/doc/app_note/AN5087.pdf

which came from this recent related question: How to interpret Allan Deviation plot for gyroscope?

with a copied graphic below where they give the bias instability for each axis.  What a specific manufacturer does on their datasheet, you would be best asking them directly.

Yes you would want to see a low bias instability, and I would think that the same number at a lower $\tau$ would actually be worst unless you were only concerned with doing a measurement over that short duration in time (implying the turn-up occurs sooner so that unit would then be drifting sooner and over the longer duration would be a lot worst). You ultimately want the lowest value at the further averaging time indicating the best long term stability, but depending on the process and use, there may be reasons to be concerned with shorter term stability as well such as if you were making difference reading measurements at 1 second intervals (then the 1 sec ADEV would be of interest).


Second moment ergodicity of gaussian random process How can I prove that a WSS Gaussian stochastic process with mean 0 is
mean-square ergodic in the second moment if and only if:

$$\lim_{n \to \infty} \frac{1}{n}\sum_{k=0}^n r_{xx}^2(k) = 0$$

When $r_{xx}(k)=E[x(n)x(n+k)]$ and $x(n)$ is the process.

A WSS Gaussian process is also a strictly stationary Gaussian process, but beyond that, Gaussianity or the lack thereof has nothing to do with the matter; the question might as well have asked for a proof for a general discrete-time stationary process that the process is second-moment ergodic if and only if the following condition holds:
$$\lim_{n\to\infty} \frac{1}{2n+1}\sum_{i=-n}^n |R_X[n]|^2 = 0 \tag{1}$$ (Note that I have changed the condition to be more symmetric).   Regardless of all this, I think that this is a question whose answer (the proof of an ergodic theorem) is way too complicated for anyone on dsp.SE to type up for the edification of the OP and any other readers who happen to stray here.

Random processes don't have to be first-moment or second-moment ergodic, and they fail to be ergodic when the autocorrelation function does not decay away fast enough as $n\to\infty$. To see this, note that the sum is one of nonnegative terms, and so if $R_X[n]$ has finite support, there is no problem at all, but if $R_X[n]$ is nonzero for all $n$, then the sum in $(1)$ is increasing with $n$, and if the increase is not slow enough, that $\frac{1}{2n+1}$ will not enough to damp the sum down enough that the limit will be $0$. $R_X[n] = \exp(-|n|)$ or $\frac{1}{1+n^2}$? No problem at all: the process is ergodic.  But consider a process in which all the random variables $X[n]$ comprising the process are the same random variable $Y$. The process is stationary, but all the sample paths are (different)  constant functions of time. If $Y=a$, the corresponding sample path has value $a$ for all $n$. Thus, while $R_X[n] = E[X[m]X[m+n]] = E[Y^2]$ is a constant, it is not necessarily the same constant as $$\frac{1}{2n+1}\sum_{i=-n}^n y[i]y[i+n] = a^2$$ and thus the process is not second-moment ergodic. It's not mean-ergodic either as noted in this answer of mine. That answer also has an example of a process that is not stationary but is nonetheless second-moment ergodic.


Allpass filter feedforward/feedback design and code I want to design/implement a simple feedforward/feedback all-pass filter, and I am having some troubles and questions with it.

From Shroeder design (based on J. Smith article here):

I came up with the following pseudocode (for first-order all-pass filter, meaning the delay was set to 1 sample, and b0 = aM):

float allPassGain = aM, value from [0.0, 1.0);

float allPassFeedback = 0.0;
float allPassDelayedSample = 0.0;

float processAllPass(float x)
{
    float delayInput = x + allPassFeedback;
    float feedForward = delayInput * allPassGain;

    float delayOutput = allPassDelayedSample;
    allPassDelayedSample = delayInput;

    allPassFeedback = delayOutput * -allPassGain;

    return delayOutput + feedForward;
}

Here x is an input sample, to process some signal with this all-pass I am going through the signal feeding it sample-by-sample into processAllPass, and store the outputs in the parallel output signal.

But unfortunately it seems like this approach doesn't work. For example, if I create a simple Phaser with a set of successful all-pass filters designed like that, I don't get the desired effect, instead the signal becomes unstable and shoots through the roof.

Does the pseudo-code have issues? I am setting aM to an arbitrary value and the delay line length is strictly 1.

Matt L.'s answer helped me to fix the code, and here is the fixed pseudo-code if that would be interesting (please also note that for the Phaser effect I find it better to switch coefficients' signs, it is slightly different version of the all-pass filter, presented in e.g. "Reducing Artificial Reverberation Algorithm Requirements using FDN" by Jasmin Frenette):

float allPassGain = aM, value from [0.0, 1.0);

float allPassDelayedSample = 0.0;

float processAllPass(float x)
{
    float delayOutput = allPassDelayedSample;
    float feedBack = delayOutput * allPassGain;

    float delayInput = inSamp + feedBack;
    float feedForward = delayInput * -allPassGain;
    allPassDelayedSample = delayInput;

    return delayOutput + feedForward;
}

First of all make sure that the magnitude of $a_M$ is strictly less than $1$; this is necessary for the filter's stability. Second, you use several unnecessary variables. This is no serious problem, but it makes it harder to see what's going on. For $M=1$ a simple pseudo-code would be

v = x - aM * v_old;
y = v * b + v_old;
v_old = v;

The important thing is that the state variable v_old must be remembered between function calls. Are you sure that your variables allPassFeedback and allPassDelayedSample are not reset to $0$ every time you call the function?


FFT on a random signal I have applied FFT on a random noisy signal (say for ex, on some stock data).
Below figures shows input and output. Now, my question is, how to extract useful peaks (frequencies) from the output? Can FFT be useful on such random data?

First of all you should take the magnitude of the FFT (use abs function) - what you've plotted is just a real part of FFT.
Secondly, depending on what you want to achieve, I would suggest to detrend your signal in order to remove the linear drift. From that point it's all up to you what you want to infer from this data.

To be honest, I'd suggest you to check out Orange3 with time-series add-on. It's a good starting point for quickly understanding your data. You can either import the CSV file or fetch series directly from Yahoo Finance. Then it's easy to perform any sort of analysis such as periodogram, ARMA, etc. This way you will avoid mistakes using FFT.


Linear Block Codes Digital Communication Problem How to design a $(8,4)$ linear block code with highest possible error correction and error detection capability?

I have studied linear block codes topic but if I do the conventional way, in case of $(n,k)$ block code, I need a generator matrix which is given as $\mathbf G=[\mathbf P|\mathbf I_k]$, where $\mathbf P$ is the parity matrix and $\mathbf I_k$ the identity matrix. But to find parity matrix I need to find the code words. But how to do that? Since third codeword should be the addition of first two code words.
Any kind of help will be appreciated.

In this answer, I will focus on error correction. I will use as little coding theory as possible, and I'll try to approach the problem from the point of view of a novice taking their first steps in code design.

Our objective is to find the best (8,4) linear block code possible, that is, one that corrects as many bit errors as possible. Is it possible to know in advance how good an (8,4) code can be? To answer this question, we have an invaluable tool: the Hamming bound. It says that a $(n,k)$ code that can correct $t$ errors must satisfy

$$ \sum_{i=0}^t \binom{n}{i} \leq 2^{n-k}. $$

It is easy to write a program that takes $n$ and $k$ as arguments, and outputs the largest $t$ that satisfies the bound. For $n=8$ and $k=4$, we find that the largest possible $t$ is 1. Note that, since the minimum distance is given by

$$ d_{min} = \lfloor (t-1)/2 \rfloor, $$

we still don't know whether the best code has $d_{min}=3$ or $d_{min}=4$. Both will have the same error correcting capabilites, but the code with minimum distance 4 has additional error detecting capabilities.

(Note that, even if there exist numbers $(n,k,t)$ that satisfy the Hamming bound, it is not guaranteed that a $(n,k,t)$ code actually exists. However, if the bound is not satisfied, we know for sure that no such code exists.)

Now that we have a goalpost ($d_{min}=4$), we can focus on finding a code that achieves it. We know that the generator matrix $G$ has $k$ rows and $n$ columns. In its systematic form, $G=[I_k|P]$, where $P$ has $k$ rows and $n-k$ columns. For our (8,4) code, we have to find a matrix $P$ with 4 rows and 4 columns.

Since we don't have any other tools at this point, we can try a brute-force approach, since there are "only" $2^{16}$ possible matrices $P$ (actually a bit less, since we don't want to repeat any columns). Let's try this matrix:

$$ P=\begin{bmatrix}1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1\end{bmatrix}. $$

It is not hard to write a program that generates all 16 possible code words using $P$. The final step is to verify if our code has $d_{min}=4$. This is easy, since the code is linear; all we have to do is find the non-zero code word with smallest weigth (with smallest number of ones). For the matrix $P$ given above, it turns out that $d_{min}=2$, so the code is not the best possible.

It is easy to find generator matrices with $d_{min}=3$; extending the distance to 4 takes a bit more work.

In conclusion:

The Hamming bound can be used to determine what values of $t$ and (approximate) $d_{min}$ we can expect from a code with small $n$ and $k$. For large $n$ and $k$, the achievable value of $t$ may be much less than that given by the bound. Surprisingly, determinining $t$ for given $n$ and $k$ is (for the general case) an unsolved problem.
Design by brute force (by trying many $P$ matrices) is feasible only for small codes. Even for an (8,4) code you may have to try more than 60,000 possibilities. The number grows exponentially with the code size. The amount of work required to generate the entire code, and to determine $d_{min}$ by inspection, soon grows so large that it becomes entirely unfeasible.
It is for this reason that so much research has been done into techniques that allow for less complex design of large codes. When you study other codes, such as cyclic and convolutional codes, keep in mind this motivation. You probably shouldn't try to solve your problem using the brute-force technique outlined here except to get a feel for the sheer magnitude of the computational complexity involved in the design of good codes, and also in their encoding and decoding procedures.


Digital Tanlock carrier recovery loop I’m currently trying to understand a digital tanlock carrier recovery loop for QPSK (for example)

I cannot understand how it works.

For example we an ideal QPSK constellation.

phi = atan2(Q/I) = pi/4,3pi/4,5pi/4 or 7pi/4.
After multiplying by M=4 and (mod 2pi) operation the angle is pi.

So after rotating constellation by pi we only gain phase ambiguity and no phase errors.

But how does it work in case of phase/frequency errors?

And one more question. What does "carrier tracking loop BW mean"?

I guess it refers to loop filter?

Is it the bandwidth of the filter? F(z) = Ki + Kp/(1-z^-1).

Jason thanks a lot for the answer.

You understand the basic operation of the phase detector. Reproducing some of the above, the detected phase error is:

$$
e_k = M \arctan\left(\frac{i_k}{q_k}\right) \mod{2\pi}
$$

As you noted, for an ideal QPSK constellation at the input to the phase detector, $e_k = \pi$ for all possible symbol values. However, in the presence of phase error $\theta$, the soft decisions at the phase detector input will have the form:

$$
\tilde{s_k} = s_k + \theta, \ \theta \in \left[-\frac{\pi}{M}, \frac{\pi}{M}\right]
$$

In this model, $s_k$ is the value of the QPSK symbol that is closest in phase to the observation. If you were making a hard decision on the symbol value based on observation $\tilde{s_k}$, this would be the one you would choose. Because the $\arctan$ function directly calculates the phase of the observed symbol, this phase error simply adds to the phase detector output: $e_k = \pi + \theta$. This allows the phase detector to detect the phase error in each symbol, modulo $\frac{2\pi}{M}$.

The job of the loop filter is to turn this stream of error outputs from the phase detector into a control signal that will adjust the frequency of the NCO that is used before the sampler. The integrator in the loop filter that you showed combats any frequency error (which manifests as a component to the phase error that increases linearly each symbol period). The proportional part of the filter will handle any constant component to the phase error as well.

Carrier tracking loop bandwidth can have multiple interpretations: sometimes it can refer to the pull-in bandwidth of the loop. For example, how large of a frequency error can I have and still expect to achieve lock? If the frequency error is too large, then the loop cannot detect and/or correct for it properly.

As a degenerate example, consider a frequency error $\omega = \frac{\pi}{2T}$ radians per second: this means that the carrier would spin exactly one quarter cycle between successive symbols. The phase detector described here would not detect any phase error at all (beyond the constant $e_k = \pi$ component discussed above).


Computer Vision/Image Processing Riddles I know this topic is not clearly relevant to the forum, but I think it's fun and interesting.

I'm looking for computer vision riddles - fun and interesting questions that have a clear and precise solution. For example the following: We have a binary image with paths (contours) that are one pixel wide. Come up with an algorithm to find the ends of the paths. I know of two solutions to this riddle.

Do you know of any other riddles? I'll be happy to hear some!

Thanks,

Gil.

Nice idea :)

A few years ago, the Instagram team proposed the Unshredder challenge. The chalélenge is over now, but the game is still interesting and could be extended to more difficult puzzles.


DSP interview question: use of the identity in development of a significant transform I'm preparing interview and found this question. But I don't really understand what is the question. Does it ask about Fourier transform or Z transform?

How the simple identity

$$xy=\frac{1}{2}x^2 + \frac{1}{2}y^2 - \frac{1}{2}(x - y)^2$$

was used in developing a significant transform in digital signal processing?

This is related to Chirp Z-transform (CZT) (refer to the Bluestein's algorithm). Using this identity, the CZT can be expressed in terms of a convolution. Hence, it can be efficiently implemented using FFT.