Finding parameters of exponentially decaying sinusoids (Matrix Pencil Method) The matrix pencil method is an algorithm which can be used to find the individual exponential decaying sinusoids' parameters (frequency, amplitude, decay factor and initial phase) in a signal consisting of multiple such signals added. I am trying to implement the algorithm. The algorithm can be found in the paper from this link or this link.

In order to test the algorithm, I created a synthetic signal composed of four exponentially decaying sinusoids generated as follows:

fs=2205;
t=0:1/fs:249/fs;
f(1)=80;
f(2)=120;
f(3)=250;
f(4)=560;
a(1)=.4;
a(2)=1;
a(3)=0.89;
a(4)=.65;
d(1)=70;
d(2)=50;
d(3)=90;
d(4)=80;
for i=1:4
    x(i,:)=a(i)*exp(-d(i)*t).*cos(2*pi*f(i)*t);
end
y=x(1,:)+x(2,:)+x(3,:)+x(4,:);

I then feed this signal to the algorithm described in the paper as follows:

function [f d] = mpencil(y)

%construct hankel matrix
N = size(y,2);
L1 = ceil(1/3 * N);
L2 = floor(2/3 * N);
L = ceil((L1 + L2) / 2);

fs=2205;
for i=1:1:(N-L)
    Y(i,:)=y(i:(i+L));
end

Y1=Y(:,1:L);
Y2=Y(:,2:(L+1));

[U,S,V] = svd(Y);
D=diag(S);
tol=1e-3;
m=0;
l=length(D);
for i=1:l
    if( abs(D(i)/D(1)) >= tol)
        m=m+1;
    end
end
Ss=S(:,1:m);
Vnew=V(:,1:m);
a=size(Vnew,1);
Vs1=Vnew(1:(a-1),:);
Vs2=Vnew(2:end,:);
Y1=U*Ss*(Vs1');
Y2=U*Ss*(Vs2');
D_fil=(pinv(Y1))*Y2;
z = eig(D_fil);

l=length(z);
for i=1:2:l
    f((i+1)/2)= (angle(z(i))*fs)/(2*pi);
    d((i+1)/2)=-real(z(i))*fs;
end

In the output from the above code, I am correctly getting the four constituent frequency components but am not getting their decaying factors. If anybody has prior experience with this algorithm or has some understanding about why this discrepancy might be there, I would be very grateful for your help. I have tried rewriting the code from a scratch multiple times but it has been of no help, giving the same results.

Any help would be highly appreciated.

In his book, Numerical Methods that Work, Forman S. Acton writes, page 253:

  For it is  well known that an exponential equation of this type [$y=Ae^{-at}+Be^{-bt}$] in which all four parameters are to be fitted is extremely ill conditionned. That is, there are many combinations of $(a,b,A,B)$ that will fit most exact data quite well indeed (will you believe four significant figures? and when experimental noise is thrown into into the pot, the entire operation becomes hopeless. But those with Faith in Science...

A four exponential sum seems quite ill-posed to me, yet there have been recent advances using sparsity assumptions or more involved optimization algorithms. Recent examples can be found for instance in Parameter estimation of monomial-exponential sums or Nonlinear Approximation by Sums of Exponentials and Translates, or citations thereof (I have not had the opportunity to test them by myself.


Unit step function present in convolution result We are currently learning about convolution in my signals and systems class, and one of our textbook problems is to compute the convolution of $x(t) = e^{-\alpha t} u(t)$ and $y(t) = e^{-\beta t} u(t)$. After working out the integral, I obtained the following result

$$
x(t) * h(t) = \frac{e^{-\alpha t} - e^{-\beta t}}{\beta - \alpha}.
$$

The solutions manual for my textbook has the same answer, except their result is multiplied by $u(t)$. I am unsure of where this is coming from, since I know that I computed the integral correctly. I used the unit step functions to simplify the bounds on the integral, so I do not see how they are still around in the final answer.

Can anyone explain what is going on here?

A step by step derivation of this convolution would start with the following:

$$
\begin{align}
f(t)
  &= \int_{-\infty}^\infty x(\tau) h(t-\tau) \,d\tau \\
  &= \int_{-\infty}^\infty e^{-\alpha\tau} u(\tau) e^{-\beta(t-\tau)} u(t-\tau) \,d\tau
\end{align}
$$
And this is probably where you probably went on and tried to simplified the bounds on the integral with:
$$
\begin{align}
  f(t) &= \int_0^t e^{-\alpha\tau} e^{-\beta(t-\tau)} \, d\tau \\
\end{align}
$$
forgetting that this is only true if $t \geq 0$.
Indeed for $t < 0$ the product $u(\tau) u(t-\tau)$ is 0 for all values of $\tau$, and the result of the integral is thus 0.

So, the correct expression would be:
$$
\begin{align}
f(t)
  &= u(t) \int_0^t e^{-\alpha\tau} e^{-\beta(t-\tau)} \, d\tau \\
  &= \cdots \\
  &= u(t) \frac{e^{-\alpha t} - e^{-\beta t}}{\beta-\alpha}
\end{align}
$$


Detecting trough widths and locations in 1d signal I have the following function obtained from averaging 2d camera image over one axis, where detection of aligned objects is desired:

I need to detect (roughly speaking) locations and widths of the troughs -- sometimes and object will be missing, thus the trough will be much wider. The approach should have some resilience WRT possible lighting variations and other noise; this, and also a desire to use something new and useful) prevents me from going the naive way of thresholding and clustering or something similar.

I was looking at using wavelets: using ideal trough (Haar?) as the mother wavelet function, and the analysis would give its translation and scale at troughs in the signal.

I have no prior experience with wavelets. My background is numerics (variational analysis/FEM, particle systems programming), c++ and python, but I am quite new to signal processing (and its terminology).

What would be the most suitable thing to do? Suggestions, pointers to articles, books, or, best, code examples are much appreciated.

EDIT:

I finally used continuous wavelet transform (as implemented in scipy) where the result local maxima (in yellow) show both location (x-axis) and width (y-axis) of the trough.

import numpy as np
from scipy import signal
# scale the ricker (mexican hat) function to 10-60px width
widths=np.arange(10,60,.2)
cwtout=signal.cwt(vscan,signal.ricker,widths)
plt.imshow(cwtout,extent=[0,len(vscan),60,10],aspect='auto',vmax=abs(cwtout).max(),vmin=abs(cwtout).min())

A wavelet analysis or wavelet transform will give you scale(frequency) vs time vs amplitude. I am not sure how that will help you with your problem.
Here is an algorithm for wavelet transform :

Choose your wavelet : Haar , Molet etc ..say $ \psi (n) $
Choose Scale : the scale should be a multiple of smallest time interval dt.
Calculate the FFT of wavelet - $\Psi_0 (s\omega _k)$
Normalise FFT of the Wavelet - $ \Psi (s\omega_k)=(2\pi s / \delta t)\Psi_0 (s\omega _k)$
Find FFT of time series x(t) = $ X_k = 1/N  \sum_{k=0}^{N-1} x(n) e^{-2\pi i k n/N } $
Then Wavelet transform is :
  $ W(s)= \sum_{k=0}^{N-1} X_k\Psi^* (s\omega_k) e^ {i\omega_k n \delta t} $    ...where * is complex conjugate
Generally the Fourier transform of the wavelet should be known analytically.

Few good sources :

http://paos.colorado.edu/research/wavelets/bams_79_01_0061.pdf
http://users.rowan.edu/~polikar/WAVELETS/WTpart1.html


Why $x(t) = \cos(\frac{\pi}{3}t)$ and $y(t) = \sin(\frac{\pi}{3}t)$ are considered identical signals? Consider the following question that I toke from Signals and systems: Part I exercises, question 1:

  (b) With $x(t) = \cos(\omega_x(t + \tau_x) + \theta_x)$ and
  $y(t) = sin(\omega_y(t + \tau_y)+ \theta_y)$, determine for which of the
  following combinations $x(t)$ and $y(t)$ are identically equal for all
  $t$.

Given that:
$$
\begin{matrix}
\omega_x & \tau_x & \theta_x & \omega_y & \tau_y & \theta_y \\
\hline
\pi/3    & 0      & 2\pi     & \pi/3    & 1      & -\pi/3
\end{matrix}
$$

The solutions manual shows that $x(t)$ and $y(t)$ are identical with these values:

  $\omega_x = \omega_y$, $\omega_x\tau_x + \theta_x = 2\pi$, $\omega_y\tau_y + \theta_y = \frac{\pi}{3}(1) - \frac{\pi}{3} = 0 + 2\pi k$

  Thus, $x(t) = y(t)$ for all $t$.

It toke me sometime to understand this, but still I'm not sure how these signals are identical.

If we substitute the constants we will get:

$$
\begin{align*}
x(t) &= \cos(\frac{\pi}{3}(t + 0) + 2\pi) \\
     &= \cos(\frac{\pi}{3} t + 2\pi) \\
     &= \cos(\frac{\pi}{3}t) \\
\\
y(t) &= \sin(\frac{\pi}{3}(t + 1) - \frac{\pi}{3}) \\
     &= \sin(\frac{\pi}{3}t + \frac{\pi}{3} - \frac{\pi}{3}) \\
     &= \sin(\frac{\pi}{3}t)
\end{align*}
$$

which is clear that $x(t) = \cos(\frac{\pi}{3}t)$ and $y(t) = \sin(\frac{\pi}{3}t)$ are not identical:

Can someone please explain what I'm missing?

There is obviously an error in the problem formulation. From the given solution, the problem should read

  With $x(t) = \cos(\omega_x(t + \tau_x) + \theta_x)$ and
  $y(t) = \cos(\omega_y(t + \tau_y)+ \theta_y)$, determine for which of the
  following combinations $x(t)$ and $y(t)$ are identically equal for all
  $t$.

In this case the given solution $\omega_x=\omega_y$ and $\omega_x\tau_x+\theta_x=\omega_y\tau_y+\theta_y+2\pi k$ makes sense. The solution to the problem as it is given (with $y(t)=\sin(\omega_y(t + \tau_y)+ \theta_y)$), would be

$$\omega_x\tau_x+\theta_x=\omega_y\tau_y+\theta_y-\frac{\pi}{2}+2\pi k$$


What does the Phase part of a FFT of a received signal represent? Imagine we've receive a signal, $R$, at a distance of $L$ from a transmitter.
Then, we do FFT on the received signal. This is what FFT gives us:

\begin{align}
R(f) &= A e^{j\phi}   \tag1\\
\
\end{align}

where, $A$ is the amplitude, we can take it by abs() function over the FFT. And, $\phi$ is the phase and we can take it by angle() function over FFT.

Now this is my confusion:
If the initial phase, I mean the phase of the transmitted signal, would be $\phi_0$, then,
Is it correct if I say:

$\phi = \phi_0 + \theta$ ?

$\theta$  is a number that is added to the transmitted phase over the link.

or I can also write this:

$\phi_{received} = \phi + \phi_0$ ?

$\phi_{received}$ is the phase of the received signal.

Which one of those is correct?

Please take this into account that I have no information about the transmitted signal. Just assume a transmitted signal with initial phase $\phi_0$.

Also just consider a very simple case, like a line-of-sight signal.
And, I'm using MATLAB for FFT, and angle().

EDIT:
I have two received signals at the distances of $L_1$ and $L_2$, with the phases of $\phi_1$ and $\phi_2$.
I calculated the difference between these two phases. And, plotted them.

And, the phase shift between these two received signal should be

$\phi_2 - \phi_1 = 2\pi f \tau$

So, Can I just go and consider the two signals in the time domain and calculate the time difference, $\tau$, between them?
Any thoughts?
What am I missing here?

The phases of an FFT result are a phase measurement relative to the time that the FFT input sample window was taken.  Each phase measurement is fairly useless unless you know the absolute time that the sample window was taken, or know the offset of the window from another measurement window (another FFT window with a known overlap/offset/delay/etc., or another window sampled at the same time but offset in space/etc., or some combination of the two), or are comparing or using phases across FFT bins (for reconstruction, convolution, etc.)

Whether the phase is positive or negative depends on how the cosine (or evenness or real) and sine (oddness or "imaginary") components of each DFT result bin relate to each other within the FFT aperture or sample window.  Usually a negative phase indicates the opposite oddness component (hi/low versus low/hi asymmetry) than a positive phase.

If you do an FFT-shift before the FFT, then the phase measurement (the evenness/oddness ratio represented by the atan2 function) is relative to the middle of the vector of samples of the window fed to the FFT, independant of whether the signal is exactly integer periodic in aperture or not.


Adding Bass to a signal How can bass/low end be added to a signal?

With a low-pass filtered sub-oscillator running at an octave or two below the signal, there still doesnt seem to be enough bass.

Is there a mathematical function to easily add some warmth or boom to the sub oscillator? Or possibly adding two in unison, slightly detuned?

Let's say the audio is given by $ x $.
Let's define Low Pass Filter by $ h $.

Whet you need to do is:

$$ y = x + \alpha (h \ast x) $$

Where $ \alpha $ is the amplification level of the Low Pass data.
Once you set which frequencies you want to amplify it is easy to design.


How to extract the co-ordinates of an ECG signal from a scanned image I checked the earlier question on "How do I recover the signal from an ECG image"

Currently I am trying to write a program to extract the x and y co-ordinates of an ECG signal embedded within the background grid . I was able to use the image from the earlier post and used some inbuilt functions to crop the image and convert the image to grayscale and then to binary.

This is how I did it:

% Cropping the image and saving the co-ordinates
img=imcrop(imread('ecg_test.jpg'));
% select the boundries in the pop up window , hit on crop image to save the
% cropped image x,y co-ordinates
img1=rgb2gray(img);  % Converting the cropped image to grayscale
Bw=im2bw(img1,0.8);  % converting the grayscale image to binary image
imshow(Bw) ![cropped and binary image of the actual ecg_test.jpg][4]

Now I want to clean the image and then extract the co-ordinates of this image.  How do I implement the MATLAB part of:

img = ColorConvert[
   ImagePad[
    Import["http://i.stack.imgur.com/500Kg.jpg"], {{0, -20}, {0, 0}}],
    "Grayscale"];

Image[
 Transpose[Function[With[{m = Min[#]},
     Map[Function[{v}, If[v == m, 1, 0]], #]]] /@
   Transpose[ImageData[img, "Real"]]
  ]
 ]

I believe this is an optimisation problem. E.g. what is the path from the left hand side of the image to the right hand of the image that is 'most black'.

As such 'rectangular shortest path' method would be ideal. You need to resample the image something like 10x to ensure that the method gets all the peaks properly, but in the end you will get something like this:

Update:

Here is a paper on circular shortest paths, but it is basically the same thing:

http://vision-cdc.csiro.au/changs/doc/sun-pallottino03pr2.pdf

Update 2:

MATLAB script:

IM = imread('500Kg.jpg');
IM = rgb2gray(IM);
IM = double(IM);

% Over sample 8 times in horizontal direction
% Amount to oversample is determined by maximum slope, e.g. if max slope is
% 16, oversample by 16.
IM2 = interp2(IM,linspace(1,size(IM,2),size(IM,2)*8),linspace(1,size(IM,1),size(IM,1))');

% Calculate shortest path using 255 - image intensity as energy function
[spath,senergy,slastindex,im] = shortestPath(255-IM2',255,0,1,100);

% Down sample path to account for original oversampling
p_x = 1:size(IM2,2);
p_x = p_x ./ 8;
p_y = spath(:,2);

% Plot on top of original image
imagesc(IM); colormap gray; hold on; plot(p_x,p_y,'r'); hold off;

Shortest path function:

function [ p, e, l, im ] = shortestPath( m, opt, a, ax, bx  )
%SHORTESTPATH Summary of this function goes here
%   Detailed explanation goes here

%Get size of matrix
[sy, sx] = size(m);
im = zeros(sy,sx);

%Convert to value different from optimal
value = abs(opt - m);

%Preallocate index matrix and energy matrix
lastindex = zeros(sy, sx);
energy = zeros(sy, sx);

%Initialise first rows
energy(1,:) = value(1,:);
lastindex(1,:) = 1:sx;

%Initialise temp matrix
t = zeros(1,3);

%Loop through remaining rows
for row = 2:sy
    for col = ax:bx
        %Get the last energies and current values
        if col == ax
            t = [inf energy(row-1,col:(col+1))];
            cv = [inf value(row,col:(col+1))];
        elseif col == bx
            t = [energy(row-1,(col-1):col) inf];
            cv = [value(row,(col-1):col) inf];
        else
            t = energy(row-1,(col-1):(col+1));
            cv = value(row,(col-1):(col+1));
        end

        %Add energy for moving
        t = t + [a 0 a];

        %Add energy from difference from optimum
        t = t + cv;

        %Find minimum
        [v,i] = min(t);

        %Save new values
        energy(row,col) = v;
        lastindex(row,col) = col + i - 2;

    end
end

[v,li] = min(energy(sy,ax:bx));
li = ax + li - 1;
p(sy,:) = [sy li];
im(sy,li) = 1;
for row = (sy-1):-1:1
    i = lastindex(row,li);
    p(row,:) = [row,i];
    im(row,i) = 1;
    li = i;
end

e = energy;
l = lastindex;

end


Group Delay for Hilbert Transformer and Resulting Dispersion An ideal Hilbert Transformer shifts the phase for all positive frequencies by $-\pi/2$ and all negative frequencies by $+\pi/2$ while maintaining constant magnitude everywhere. Group Delay is the negative derivative of phase with respect to frequency, so apart from the discontinuity at frequency $\omega=0$ the group delay would be 0.

Linear phase filters do not cause distortion to a waveform that has many frequency components as each of these components would remain aligned in time at the output of the filter as they were at the input.

Yet the Hilbert Transformer is clearly dispersive meaning different frequency components of a waveform would experience different delays in order for the phase shift to maintain  a $\pi/2$ relationship between the input and output of the network at all frequencies.

The phase of the Hilbert Transformer is linear, so in this case we observe linear phase but still the system is dispersive. How is this explained as this appears to be an exception that I haven’t really seen in definitions of Group Delay or linear phase?

Further imagine a case that is often used in practice with two phase tracking networks used to split a signal into two over a wide frequency with a quadrature relationship. We could imagine a condition with the first network as a linear phase filter and the second network having the same linear phase with a $\pi/2$ offset at all frequencies of interest (generalized linear phase). Here the derivative of the phase with respect to frequency is not zero but constant yet we would have the same dispersive results in order to maintain the constant $\pi/2$ offset.

A real-valued system that doesn't distort the shape of the input signal must have the following input-output relation:

$$y(t)=Ax(t-t_0)\tag{1}$$

with arbitrary real-valued constants $A>0$ and $t_0$.

In the frequency domain, Eq. $(1)$ corresponds to

$$Y(\omega)=Ae^{-j\omega t_0}X(\omega)\tag{2}$$

Consequently, the corresponding system is an LTI system with frequency response

$$H(\omega)=Ae^{-j\omega t_0}\tag{3}$$

i.e., the magnitude response must be constant and the phase response must be linear:

$$\phi(\omega)=-\omega t_0\tag{4}$$

Here, linear is meant in the narrow sense of the word, i.e., it must not have an additive constant term (hence, it must not be affine). Note that for real-valued systems, the phase must be  an odd function of frequency, so the sign of a possible additive constant in the phase would need to change depending on the sign of the frequency. Consequently, if there's such a constant phase term (as is the case for a Hilbert transformer) then the phase has a discontinuity at $\omega=0$, which implies that the group delay is undefined / infinite at $\omega=0$.

The phase response of an ideal Hilbert transformer with a possibly non-zero delay $t_0$ does not satisfy $(4)$ because it is given by

$$\phi(\omega)=-\frac{\pi}{2}\textrm{sgn}(\omega)-\omega t_0\tag{5}$$

As mentioned above, the phase of an ideal Hilbert transform has a discontinuity at $\omega=0$.


Mistake or not - Fourier Series of x(2t+3) I have a couple of resources I have from my university I had being checking and I found this:

Find Fourier Series coefficients of x(2t+3). x(t) is continuous and periodic by T.

I see this solution:

But the 3/2 is seemingly wrong. Am I right? It supposed to be 3, not 3/2. Or am I missing something? They claim that:

Please let me know if this resource is wrong or I am having the mistake and it shall be 3/2 some why.

Ok, I see now: They did it straight from the computation. The 3/2 is because they divided T by two in $e^{j3k{\frac {2\pi} T}}$.
It’s important, because it tells us that the first intuition we have that we can apply scaling to Fourier Series or Fourier Transform and then apply transformation to the Fourier Transform  is wrong. We might mistakenly think that we can scale time, apply transform, and then apply shift by 3 to the transform to move by 3 units in time. As you can see it does not produce the same result, since the shift is by 3/scaling-factor and not straight 3.


Butterworth Polynomial from Butterworth Filter For the butterworth filter, in Matlab command:

[B,A]=butter(2,0.5,'low')

I find this result:

B = [0.29289  0.58579  0.29289]
A = [1  -3.33067e-16  0.17157]

From here I cannot reach the Normalized Butterworth polynomial:

G = [1  1.4142  1]

Nor the poles:

P = [-0.7071+i0.7071  -0.7071-i0.7071]

Which is the relation or the proper parameter i am missing, if any?

The denominator of the transfer function of an analog Butterworth filter is a (non-normalized) Butterworth polynomial. Also the numerator is just a constant (i.e. the coefficients of higher order are 0).

You are designing a discrete filter here, so this does no longer apply. Your numerator has three non-zero coefficients, and the denominator is not a Butterworth polynomial.

You can check that if you design an analog Butterworth adding the parameter 's' when calling the function, you will indeed get a one non-zero element array B and the coefficients of a non-normalized Butterworth polynomial in A.


Nonuniform sampling points based on the first order derivative of signal amplitude Problem:

I need to measure a experimental data set $y(x)$.

However, the experiment is extremely time-consuming and I don't have the luxury to have linearly spaced sampling point $x$.

I need the density of sample points to be proportional to the first order derivative of the data $y^{'}(x)$.

Progress:

first measure $y(x)$ once with uniform x as a test run, then compute $y^{'}(x)$.
Then fit $y^{'}(x)$ to $pdf(x,\alpha, \beta, ... )$.
Then I use the parameters to calculate new data points with inverse CDF $X=icdf[linspace[pdf(min),pdf(max),N],\alpha, \beta, ...]$

The problem is, this method only apply to data sets that is distributed in common probability distribution functions.
Is there any way to make this work with arbitrary data ?
Thanks in advance !

A method that I use very often to approximate the CDF is based on the observation that given a set of $N$ samples, the CDF at the point x the fraction of the points that are smaller than $x$. So if have a vector of sorted samples $S$ you could approximate the CDF as $cdf(x) = i/N$ if $S_{i-1} < x < S_i$.
A simple way to plot the CDF would be
plot(sort(S), (1:N)/N)
Then you can approximate the derivative by a difference equation. You could use a binary search to find the region of the vector $S$ where $x$ lies, say $idx(S, x)$, then you could compute something like
$i_1 = idx(S, x - \Delta)$, $i_2 = idx(S, s + \Delta)$ and use $(N \cdot (S_{i_2+1} - S_{i_1})) / (i_2 - i_1 + 1)$.
A different way to choose the anchor points would be, let $i = idx(S, x)$ and compute $(N \cdot (S_{i+k} - S_{i-k})) / (2k+1)$
And you can fit use a polynomial fit if you need a better accuracy where the CDF derivatives varies rapidly.


Vestigial Filter problem? I have been stuck on this question for a while now. It has to do with vestigial sideband.

I wasn't sure if I should be dividing $H(\omega)$ graph values by 2 because only the positive side of the filter is given....

The answer that I am getting is:

We know that $VSB(t) = A_c f(t)\cos(\omega_c t) \ast h(t)$

then in frequency form I get:

$$VSB(\omega) = A_c\pi  0.5 (\delta(\omega + (\omega_c + \omega_m)) + \delta(\omega + (\omega_c - \omega_m)) + \delta(\omega - (\omega_c + \omega_m)) + \delta(\omega - (\omega_c - \omega_m))) H(\omega)$$

Now using the graph I am multiplying the appropriate $H(\omega)$ values:

$$VSB(\omega) = A_c \pi 0.5 (0.25 \delta(\omega + (\omega_c + \omega_m)) + 0.75 \delta(\omega + (\omega_c - \omega_m)) + 0.25 \delta(\omega - (\omega_c + \omega_m)) + 0.75 \delta(\omega - (\omega_c - \omega_m)))$$

where carrier wave is given by $A_c \cos(\omega_c t)$

I don't have the solution to this, so I have no idea if my answer or reasoning is correct...
It would be sooooo helpful if someone could confirm/clarify

Your reasoning seems reasonable and correct to me. Assuming that the filter's impulse response $h(t)$ is real, then its magnitude response $\left|H(\omega)\right|$ will be symmetric about $\omega=0$, so the negative-frequency terms will be scaled by the same values as those on the positive side of the frequency axis. I don't see any reason why you would need to multiply or divide by 2 anywhere.

The only thing that might differ if you see a solution somewhere else might be the scaling factor at the front of the expression. The typical transform using angular frequency would be:

$$
\mathcal{F}\{\cos(ax)\} = \pi\left(\delta(\omega - a) + \delta(\omega + a)\right)
$$

which assumes a Fourier transform definition of:

$$
\mathcal{F}\{x(t)\} = X(\omega) = \int_{-\infty}^{\infty}x(t) e^{-j\omega t}dt
$$

It looks like you have an extra factor of $0.5$ at the front of your expression that you might have gotten from the transform. Note that different implementations or expressions of Fourier transforms often have arbitrary scaling factors included; what's most important is that if you apply the forward and inverse transforms in sequence, you get the original signal back:

$$
\mathcal{F}^{-1}\left\{\mathcal{F}\{x(t)\}\right\} = \mathcal{F}^{-1}\left\{X(\omega)\right\} = x(t)
$$

or at least to some known scaling factor.


Floating point Descriptor evaluation I have come up with a new descriptor(floating point). I want to evaluate this descriptor against SIFT, GLOH, etc etc. I am using the K. Mikolajczyk, C. Schmid, approach.  The method is described in the link http://www.robots.ox.ac.uk/~vgg/research/affine/descriptors.html. They have provided the binaries and the data set. I extract the patches(41x41) and apply my descriptor on these patches. But the problem is, these patches are not rotational normalized.

SIFT achieves rotation in-variance by rotating the patch in the direction of the dominant angle of the patch. My query is ,how to rotate the patch or achieve rotation normalization? I want the original patch(41x41) and the rotated patch to have the same dimension (41x41). For example, If i use the matlab function imrotate and rotate the original patch by 60 degrees. The rotated patch has a dimension bigger than (41x41). And there are blank spaces in the border of the rotated image. I think applying my descriptor on these rotated and enlarged patches gives wrong results.

M question is, How do i rotate the patches or rotational normalize the patch(Such that the original patch and the rotated patch have the same size 41x41) to apply my descriptor?

You don't. Only circles can be rotated by an arbitrary degree and still have the same shape. A 41x41 image is an axis-aligned square, and only 90 degree rotations keep the shape.

And so rotating over 60 degrees means that the width and height needed to fit the square have increased by sqrt(3/2)..


How does shift and scaling inside of a function affect its Fourier Transform? The properties aren't entirely clear to me, sorry for the basic question.

I know the Fourier Transform of one function.  Say, $\text{rect}(x,y) \Leftrightarrow \frac{\sin \pi u}{\pi u} \frac{\sin \pi v}{\pi v}$.  Now, I have some related function and I want to easily get the FT based on simple properties.

For example, now I have something like $\frac{1}{3}\text{rect}(8x - 4, 4y - 2)$.  How does this shift and scale inside the function affect its (continuous) Fourier Transform?

For one-dimensional Fourier transforms,
$$\begin{align*}
x(t) &\leftrightarrow X(f)\\
x(t-a) &\leftrightarrow X(f)\exp(-j2\pi fa)\\
x(bt) &\leftrightarrow \frac{1}{|b|}X\left(\frac{f}{b}\right)
\end{align*}$$
all of which are proved by change of variables in the definition of the
Fourier transform
$$X(f) = \int_{-\infty}^{\infty} x(t)\exp(-j2\pi ft) \mathrm dt.$$
It is good exercise to work out the details of each of these results
at least once in one's life, and once one has done so, the two-dimensional case
is as easy as rolling off a log.  If you always look up the answers in a table
of Fourier transform properties, the right result will never come to mind just when
you need it most.

More generally, it is a good idea to remember these results in graphic
terms.  If $x(t)$ is a signal with Fourier transform $X(f)$, and $b > 1$,
say $b = 2$,
then $y(t) = x(2t)$ is a time-compressed form of $x(t)$: while $y(0) = x(0)$,
$y(3) = x(6)$, that is, $6$ seconds worth of $x(t)$ have been squeezed down
into $3$ seconds worth of $y(t)$.  Thus, the frequencies have doubled, and so
$Y(f)$ should be something like $X(f/2)$ (but see below).  More generally,

  Compressing the time axis by a factor of $b$ expands
  the frequency axis by a factor of $b$.

But that is not the complete answer.  One of the more boring results in Fourier
transform that everyone tends to ignore is
$$y(0) = \int_{-\infty}^{\infty} Y(f) \mathrm df.$$
But,
$$y(0) = x(0) = \int_{-\infty}^{\infty} X(f) \mathrm df$$
and since we have stretched the horizontal axis by a factor of $b$
in going from $X(f)$ to $X(f/b)$, we better reduce the magnitude of each
frequency component by a factor of $b$ so that the two integrals (= area
under the curve) work out to have the same value, that is,
$$Y(f) = \frac{1}{b}X\left(\frac{f}{b}\right), ~~ b > 0.$$
The result for $b < 0$ can be worked out similarly.


Why the "20" in calculating db ratio? In wikipedia, I read the formula to convert a power ratio to db. This formula is :

10*Math.log10(power_out / power_in);

But, for the sound pressure, I read the formula is :

20*Math.log10(soundLevel_out / soundLevel_in);

I read in some places a reason : it is because the sound pressure, like the voltage, have to be squared in the ratio and this squared becomes a 20 instead of a 10. Ok, but why ?
When I want to double the sound "feeling", I was told I need to double my power amplifier so, I understand it's a question of power (with no squaring) and not a question of voltage, but it's wrong.

Who can help me to understand where does this square comes from ?

According to this source, the sound power received by an aperture is proportional to the pressure squared, i.e.:
$$
P = \frac{A p^2}{\rho c} \cos \theta,
$$
where:

$P$ is the received power
$A$ is the surface area of the receiving aperture
$p$ is the pressure
$\rho$ is the density of medium
$\theta$ is the angle between the propagation direction and the surface normal


OFDM sampling frequency offset how the 4G LTE Or 5G NR OFDM measures and corrects the sampling frequency offset? if there is an algorithm using matlab please help

Actually in OFDM, it is very easy. You can check Van De Beek algorithm for compensation of carrier frequency offset.


Kalman Filter Covariance I've recently started playing with the Kalman filter for a simple 2D (x,y,dx,dy) tracking toy problem. But I seem to have some misunderstanding on what I can expect from the filter. I'm interested in plotting the uncertainty ellipse from the corrected covariance matrix but noticed a few observations:

The covariance decreases to a steady state regardless of how much error I introduce into the measurement.

The variance for x and y are exactly the same even though I introduce more measurement errors in y.

Staring at the maths for a bit it seems that this is how the vanilla Kalman Filter works. What I'm expecting is the opposite of the two points mentioned above. In the end I want to plot an uncertainty ellipse that reflects the error I'm observing. Is this possible at all? Do I have to do some post processing on the covariance matrix?
ANSWER:
It occurred to me what I needed was a Kalman Filter that has the ability to adapt its covariance. I found this paper that details a few different methods to do this.
Evaluating the performances of adaptive Kalman filter methods in GPS/INS integration by Ali Almagbile, Jinling Wang, and Weidong Ding, 2010

The covariance matrix of a Kalman filter is a function of the $ Q $ and $ R $ matrices of the model.
If you use a model where $ R $ and $ Q $ are time invariant or known in prior then the calculation of the covariance matrix $ P $ can be done offline and isn't a function of the measurements.

In some cases, advanced implementations of Kalman Filter estimate the covariances $ R $ and $ Q $ on line according to some data gathered in the process of calculating the result.


How to identify the beginning and the end of a signal using MATLAB I have multiple signals as shown in this image. I'd like automatically detect time for the earliest detection of the signal (approx 0.7e4 in S1-B3-B (orange)) and last detected occurrence of the signal (approx 4.2e4 in S2-B3-B (purple)).

Is there a smart way to that? I'm quite new to signal processing.

There are many sources on Detection Theory. I think there are two relatively simple approaches to your problem. From the example graph I can clearly see a distinct signal. if the desired detected signal has the same trend/shape a good detection technique would be in the direction of matched filtering.

(1) Matched filtering is basically correlating the test signal with the measured signal. if the signals are very similar to each other, matched filtering would return a high value on which you could detect that with a threshold. A good starting point for signal detection would be this lecture slide.

(2) However, that might already be a bit too difficult since you mentioned to be new to signal processing. Another approach when I look at the signal is seeing that the noise is very low when the signal is not present. I'm not sure if the example signals have these distinct amplitudes when the signal is not present (yellow at ~15, purple at ~-730). But if this amplitude is know, we could simply check when the signal is different from this steady value. Or we could calculate how much the signal changes by calculating the gradient (or simply the difference between samples) for example. If the gradient is high, the signal would be increasing at a higher rate.

Maybe start by calculating the difference between samples, i.e.: $s[i]-s[i-1]$ for each sample. I highly encourage to just play with the signals when starting with signal processing. Try different functions on the samples.


How to Downsample The Chrominance Components So I am doing an image compression project for college and am trying to create a visual demonstration of how downsampling chrominance components can reduce the amount of digital information while being hardly noticeable to the average person. While doing research for this project, I came across this video where the guy gives a demonstration of downsampling that leads to there being 100x less color in the image he uses as an example. I was hoping to do something similar but am unsure as to how he was able to downsample the chrominance components in an image by a factor of 10x in each direction.

Looking at the standard ratios for chroma subsampling(/downsampling) the best I could find seems to be 4:2:0 which would downsample the chrominance by a factor of 2x in each direction. It's possible that might be enough for a demonstration, but I was hoping I might could push it a little further.

As a side question, does anyone have any recommendations as to how best to apply the downsampling color to my raw image? I have some ideas involving either photoshop or GIMP, but I was wondering if anybody had a better suggestion.

For demonstration purposes GIMP is just fine. It can split an image to three separate YUV/YCbCr images. You can then manually resize the UV/CbCr images by any amount to downsample the chroma. Then upsample back to original resolution and recombine the image again.


What is the best way of filtering EEG signal for the specific frequency band? I have an EEG data with a sampling rate is 500 Hz. Actually, there are too many ways to filter EEG signals. However, the abundance of options is not valid for me. I have to use a causal, linear-phase FIR filter with a small group delay as possible. I also want to -3dB attenuation at cut- off frequencies.
The problem arises at this point.  Look at the figure below.

I applied 30th order FIR filter (hamming window 4-30 Hz) on the continuous (non-epoched) EEG. The side-lobe attenuation for higher frequencies is satisfactory but the lower frequency components for this filter failed.
I don't want to use rectangular window because of relatively low side-lobe attenuation. How can I design ideal filter to make possible my intention?

Perhaps you should try some of the filter design tools that are available online.
You can define the required filter behavior (Cutoff frequency, Ripple ..) and you will get a design solution that fits your needs.
This kind of tools will help you know what is possible to achieve in terms of filtering performance given the constraints that you detailed
Checkout the following ones ( this is not an exhaustive list) :
http://t-filter.engineerjs.com/
https://webench.ti.com/filter-design-tool/filter-response
Best regards,
Mourad


Help in understanding from book expression of variance of an estimator : PRBS vs real valued input The question is based on   Book : Fundamentals of Statistical Signal Processing by Steven Kay, Chapter 4 : Eq(4.21). The expression for the variance of the estimated coefficients when the input is PRN as

$$\textrm{Var}\left(\hat{h}_i\right)= \frac{\sigma_w^2}{Nr_{uu}[0]} \tag{1}$$
where $\sigma_w^2$ is the variance of the measurement noise; $\sigma_u^2$ is the variance of the input.

Let the model be:
$$x[n] = \sum_{i=1}^m h[i]u[n-i] + w[n] \tag{2}$$

where $u[\cdot]$ is the input process, and $w[n]$ is a zero mean white gaussian noise. $\sigma_{\rm input}^2$ is the variance of the input, $u[n]$ and   $\sigma_{\rm noise}^2$ is the variance of the measurement noise $w[n]$. The cross-correlation function $R_{x,u}[\cdot]$ between $x[\cdot]$
and $u[\cdot]$ is the periodic convolution between the sequence $h[\cdot]$ and the periodic autocorrelation function $R_{u,u}[\cdot]$

Performance is evaluated using CRLB of channel coefficients.

Can somebody please explain to clear my confusion about what the expression of  the variance of the estimator i.e, Cramer-Rao lower bound for the channel coefficients will be

(a) when $u[n]$ takes real values. Is it $\frac{m\sigma_{\rm noise}^2}{N \sigma_{\rm input}^2}$
(b) when $u[n]$ takes values from a finite symbol set $\{+1,-1\}$. Is it $\displaystyle \frac{m\sigma^2_{\rm noise}}{N}$ where $\sigma^2_u = 1$ and $r_{uu} =1 $

The expression in the book is for PRN sequence, but I don't quite understand what PRN sequence is. is it another term for symbolic values (not real values)? By PRN do we mean the values are ${+1,-1}$?

EDIT of the update based on the answer and offline discussions

Just to confirm if I am on the correct track. The general formula for $var(\mathbf{\Delta h}\mathbf{\Delta h^T}) =\frac{ \text{variance of measurement noise}}{\text{normalizing factor * autocorrelation of input}}$ = $\frac{\sigma^2}{D*\text{autocorre of input}} $ where $D = N * variance = N {\sigma^2_u }$.

Therefore, $var(\mathbf{\Delta h}\mathbf{\Delta h^T}) = \frac{\sigma^2}{N \sigma^2_u \sum u[n]u^*[n]}$ = $\frac{\sigma^2}{N(\sigma^2_u \sigma^{*2}_u)^2}$  [autocorrelation having complex conjugate term] ??

For +1/-1 sequence, this formula evaluates to $var(\mathbf{\Delta h}\mathbf{\Delta h^T}) = \frac{\sigma^2}{N\sigma^2_u} $

For QAM :$var(\mathbf{\Delta h}\mathbf{\Delta h^*}) = \frac{\sigma^2}{D (\sum u[n] u^*[n])}$ = $\frac{\sigma^2} { \big(N \sigma_u \sigma_u^*\big) \big(\sum u[n] u^*[n]\big)}$

A PRN sequence is a Pseudo-Random Noise sequence, often generated by using an Linear Feedback Shift Register (LFSR) with the feedback taps done by using a primitive irreducible polynomial in GF{2}, which is the Golois Field of 2 elements. When a primitive and irreducible polynomial in GF{2} is used, the LFSR will produce a "maximum length sequence", meaning the pattern at its output will be the longest possible for the number of registers in the LFSR.  An example generated from the polynomial $1+x^3+x^4$ is shown in the picture below. This polynomial is primitive and irreducible in GF{2}, and in this case we will use the symbols $-1$ and $+1$ as the two elements in the field), and therefore this LFSR as configured will produce a maximum length sequence; the state of the generator, being the values after each cycle at the output of each register in the generator, will cycle through every possible combination pseudo-randomly before repeating, except for the all zero state. (So in the picture shown the choices are $0001$, $0010$, $0011$, $0100$, $\ldots$, $1111$ in pseudo-random order, so the pattern of the sequence at the output will have $15$ values before repeating) Generally the maximum length is $x^N-1$ where $N$ is the order of the polynomial, the $-1$ is because the all-zero state is not allowed).

For clarity we will refer to each output of the LFSR as a "chip", and the complete sequence is referred to as a "symbol".

When you correlate to a matched PRN sequence + noise (the input), you multiply the input by each coefficient in the sequence and sum the result (correlation = multiply + accumulate).  If you are aligned to the code, the multiplication will cause the $-1$'s in the input to be $+1$'s, and the $+1$'s in the input to be $+1$'s, essentially converting the pseudo-random pattern to be all $1$'s. This output of all ones is then accumulated to achieve the correlator output.  (Below shows the multiplication that would then be fed to an accumulator after each result).

The accumulation on the signal component (which is all $1$'s when aligned) will go up in magnitude by $N$ given $N$ samples in the sequence (if we assume one sample per chip).

Assuming the input noise itself is uncorrelated to the PRN sequence, and white within our sample space (which means that each sample of noise is uncorrelated to all the other samples), then the multiplication by $+1/-1$ will not change the variance of the noise. Therefore the standard deviation of the noise will increase by $\sqrt{N}$ in the accumulator. (When you add N uncorrelated noise samples, the standard deviation goes up by the $\sqrt{N}$).

So the variance of the noise is increasing by $N$, but the variance of the signal is increasing by $N^2$ (since the signal magnitude, not power, went up by a factor of $N$). The end result is the signal-to-noise ratio (SNR) changes by a factor of N in power or $10\log(N)$ in $\textrm{dB}$, and this is referred to as the "processing gain" of the PRN correlator.

If we normalize the signal level to be the same at input and output, the noise standard deviation will have changed by a factor of $1/\sqrt{N}$, and it's variance by $1/N$, matching your equation.

Stating this mathematically, consider our samples at the output of the multiplier as
$$X_i\quad\text{for}\quad i=1, 2, \ldots, N$$

The multiplication by $+1/-1$ does not change the distribution of the noise on the input samples (they will still be independent and identically distributed, IID), with the same variance as the input given as
$$\textrm{Var}(X_i) = \sigma^2$$

Multiplying by $\pm 1$ with a synchronized PRN causes the mean to be 1 as previously described,
$$E(X_i)=1$$

We will denote the output after accumulating over all $N$ samples as $Y$, but we will also divide by $N$ to be consistent with the normalized autocorrelation $r_{uu}=1 = E(Y)$:
$$r_{uu} = \frac{1}{N}\Sigma_1^N X_i = E(Y)$$

Where $X_i$ is the product of the input waveform with the PRN.

This results in the expression for the mean and variance of $Y$ as follows:

\begin{align}
E(Y) &= E\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N}E\left[(X_1+X_2+ \ldots +X_N)\right]\\
\textrm{Var}(Y) &= \textrm{Var}\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N^2}\textrm{Var}\left(X_1+X_2+ \ldots+X_N\right)
\end{align}

\begin{align}
E\left[(X_1+X_2+ \ldots +X_N)\right] &= E[X_1]+E[X_2]+ \ldots E[X_N] = N\\
\textrm{Var}(X_1+X_2+ \ldots +X_N) &= \textrm{Var}(X_1)+\textrm{Var}(X_2) + \ldots \textrm{Var}(X_N) = N\sigma^2
\end{align}

Therefore
\begin{align}
E(Y) &= NE[(X_1+X_2+ \ldots +X_N)] = 1\quad\text{and}\\
\textrm{Var}(Y) &= \frac{1}{N^2}\textrm{Var}(X_1+X_2+ \ldots +X_N) = \frac{1}{N^2}N\sigma^2 = \frac{\sigma^2}{N}
\end{align}


Clarification regarding discrete wavelet transform One of the books on "Conceptual Wavelets" by Fugal explains some major differences between the undecimated discrete wavelet transform (UDWT) vs. discrete wavelet transform (DWT). In UDWT the scale of wavelet is increased continuously just like the continuous wavelet transform, but the scale increases in dyads (powers of 2). In the DWT, which is the most commonly used in MATLAB, the filter size remains the same, but the data is reduced dyadically.

His exact wordings are "As discussed briefly in the preview, instead of dyadically stretching the filters, the conventional (decimated) DWT dyadically shrinks the signal instead"

He is using the example of Haar wavelets on a very small set of data, simply exam scores as a signal of eight exam scores [ 80 80 80 80 0 0 0 0].

So the question is when we downsample by 2, which one is more correct to throw away even samples or odd samples?

Assuming that you have a sufficient number of samples, throwing away odd / even samples does not matter.

The DWT can be thought of as measuring time/frequency content with varying levels of time/frequency resolution.  For non frequency varying signals (like chirps), the odd and even samples both contain the same frequency and time content. (Unless we are talking about some explosive transient that occurs in one sample)  then you would need some a priori knowledge about that.

If the signal is of odd length, taking the odd samples would give you one more sample than the even samples but wouldn't make a difference in terms of information gained.


How to Decompose a Separable Filter? I have done some research on the Internet and I have found that a given 2D mask is separable if it exists only a singular value of that matrix. For example, given the following matrix:

$$A = \begin{pmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{pmatrix}$$

If I do for instance in matlab svd(A), it gives me the vector $[0\; 0\; 6]$, thus, the filter is separable. Nevertheless, I would like to know if it is possible to calculate the vectors that multiplied (or convolved) give raise to $A$. As a matter of fact, I know that:

$$\begin{pmatrix} 1  \\ 2  \\ 1  \end{pmatrix} \begin{pmatrix} 1  & 2  & 1  \end{pmatrix}= \begin{pmatrix} 1 & 2 & 1 \\ 2 & 4 & 2 \\ 1 & 2 & 1 \end{pmatrix} $$

However, I would like to calculate those vector for any separable matrix.

Thank you so much for your responses.

Indeed you can do that.
You may look on my answer to How to Prove a 2D Filter Is Separable?

By the SVD for any filter $ A $:

$$ A = \sum_{i = 1}^{n} {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

Since we're talking about separable filter then:

$$ A = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $$

So the columns filter is $ \sqrt{\sigma}_{1} {u}_{1} $ and the rows filter is $ \sqrt{\sigma}_{1} {v}_{1}^{T} $.

Since what's important is that $ A = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $ you can actually choose any 2 vectors which their outer product is $ {\sigma}_{1} {u}_{1} {v}_{1}^{T} $. The choice isn't unique. For example you could take for columns $ {\sigma}_{1}^{2} {u}_{1} $ and $ \frac{1}{{\sigma}_{1}} {v}_{1}^{T} $ for the rows.


Doubt in state space representation $y$ is scalar observations and so C will be a 1x2 matrix.

I want to represent the following model as a state space representation so as to estimate the hidden states from the noisy observations $y$ using Kalman filter.

The  state space model :

x(t+1) = Ax(t) + w(t)

y(t) = Cx(t) + v(t)

w(t) = N(0,Q)

v(t) = N(0,R)

$w(t)$ is a pseudo-random binary signal that excites/ drives ; and $v(t) = N(0,\sigma^2_v)$ is the measurement noise.

The model is an FIR (MA) filter

$$x(t) = h_1 \epsilon(t-1) + h_2 \epsilon(t-2) + \epsilon(t)$$
$$y(t) = x(t) + v(t)$$
$$ y(t) = h^T \epsilon(t) + v(t)$$

(In vector form)

where $\epsilon(t) = w(t)$.

The other model is an IIR (AR) filter
$$x(t) = ax(t-1) + bx(t-2)+ w(t)$$

The state space representation:

$$x(t+1) = a^Tx(t) + w(t)$$

$$y(t) = h^Tx(t) + v(t)$$

How do I represent these as state space so as to apply Kalman Filter?

There are several ways to represent time series models. This is how I proceeded, but unsure because the output of the log-likelihood is a matrix of 2 by 2 with off diagonal elements being infinity and the diagonal elements are same positive values. So, the dimension and the value of log-likelihood is incorrect, I should get negative instead of positive values.

FIR :

Re-writing the above model as:

$$x(t+1) = h_1 \epsilon(t) + h_2 \epsilon(t-1) + \epsilon(t+1)$$
$$y(t) = Cx(t)+v(t)$$

State Space :

$ \left[ \begin{array}{c}
x(t+1)  \\
x(t)\\
 x(t-1) \end{array} \right] $ = $ \left[ \begin{array}{ccc}
1 & h_1 & h_2 \\
0 & 1 &  h_1 \\
0 & 0 & 1\end{array} \right] $  $\times$ $ \left[ \begin{array}{c}
e(t+1)\\
e(t) \\
e(t-1)\end{array} \right] $

$  y(t)$ = $\left[ \begin{array}{c}
1 \hskip 5 pt 0 \hskip 5 pt 0\end{array} \right] $
 $\times$ $ \left[ \begin{array}{c}
e(t+1)\\
e(t) \\
e(t-1)\end{array} \right] $ + $v(t)$

IIR (AR model)

$ \left[ \begin{array}{c}
x(t+1)  \\
x(t)  \end{array} \right] $ = $ \left[ \begin{array}{cc}
a & b  \\
1 & 0
\end{array} \right] $  $\times$ $ \left[ \begin{array}{c}
x(t)\\
x(t-1)\end{array} \right] $ + $\left[ \begin{array}{c}
1\\
0\end{array} \right] $ $\times$ $\left[ \begin{array}{c}
w(t+1)\\
w(t)\end{array} \right] $

$  y(t)$ = $\left[ \begin{array}{cc}
1 \hskip 5 pt 0\end{array} \right]$ $\times$ $ \left[ \begin{array}{c}
x(t)\\
x(t-1)\end{array} \right] $ + $v(t)$

Your FIR state space representation seems to be doing too much.

The way I would write it is to have the process noise is $\epsilon(t)$ as your input, and your state as two time-delayed copies of it:

$$x(t+1) = \left[ \begin{array}{c}
\epsilon(t+1)\\
\epsilon(t) \\
\epsilon(t-1)
\end{array} \right] = \left[ \begin{array}{ccc}
0 & 0 & 0\\
1 & 0 & 0\\
0 & 1 & 0
\end{array} \right] x(t) + \left[\begin{array}{c}
1\\
0\\
0
\end{array} \right] \epsilon(t + 1)
$$

then your output equation is just:
$$
y(t) = \left[ \begin{array}{ccc}
1 & h_1 & h_2 \end{array} \right] x(t) + v(t)
$$

Your IIR representation is too confused for me to make head or tail out of. Can you clarify that you have the right information there?

Other questions:

What is $w(t)$? It does not appear to take part in any of your equations?
What are the parameters for $u(t)$?


How to calculate the FFT period I have a question about finding the FFT duration for 802.11a preamble.

According to the standards, when the bandwith is $20 \textrm{MHz}$ and for $N=64$ we have $\Delta f = 20\textrm{MHz}/64 =312.5 \textrm{kHz}$ and it says FFT duration should be $1/\Delta f=1/(312.5 \textrm{kHz}) = 3.2 \textrm{us}$.

My question is that why FFT duration should be "1 over subcarrier spacing"?

You must have understood the notion of digital linear modulation or discrete time vs continuos time (see Chapter 2).

Another reference.

OFDM can be thought as FDM with sinc pulse whose delay-$T$-shifted versions form an orthonormal basis. In frequency domain, they are seperated by $1/T$ which is denoted $\Delta f$ , i.e. subcarrier spacing. You fix a $T$, then $\Delta f = 1/T$; you fix $\Delta f$, then $T = 1/\Delta f$.

Look in frequency domain, because data is modulated in seperated subcarriers, the final signal is the sum of these subcarriers. Because the Fourier transform preserves the sum operation, in time domain, the final signal is also the sum of inverse-Fourier-transformed subcarrier signals that have the same duration $T$. Thus the final signal, which is sum of several signals having the same duration $T$, has the same duration $T$ that is called OFDM symbol duration or FFT duration or FFT period.


How to produce uniformly distributed low-frequency noise? I want to design noise for the sake of modulation and it should have a frequency of under 1 Hz. This seems impossible to do for a frequency that low.
Here is my Matlab script. The lowest frequency that gives me something useful is 2 Hz. But also then the noise does not travel over the whole amplitude range from 0 to 1. It jumps to 0.9 and is stuck there.
Is there maybe a filter that is suited better?
fs = 44100;
L = 10 * fs;
fpass = 2;
fstop = fpass*2;
x = rand(L,1);

fs = 44100;
Wp = fpass * (2 * pi)/fs;
Ws = fstop * (2 * pi)/fs;
Rp = 0.01;
Rs = 40;

[n, ~] = ellipord(Wp, Ws, Rp, Rs);
[b, a] = ellip(n, Rp, Rs, Wp);

y = filter(b, a, x);

plot(x)
hold on;
plot(y)
hold off;

Well, this is one way you could do it: with LTspice, use the white() function. There are also rand() ([0,1] V, pulses) and random() ([0,1] V, smooth pulses), but white(2*fs*time) will give you what you want between [-0.5, 0.5] V:

B1 is a behavioural current source (bi, or bi2, as it appears in the component selection dialog, F2), and the .wavecommand exports it as a .wav file named output.wav, with 16 bits resolution (can be anything from 1 to 32) and 44.1 kHz sampling frequency (can be any number from 1 to 4294967295). For this case, the file is 8,820,046 bytes, which can be imported. This is a 2 Hz white noise, band limited, properly sampled signal, 100 sec long. Optional settings can include a tighter timestep, though opt plotwinsize=0, here, disables waveform compression and for a 2 Hz signal, it may not make much sense.


Finite duration vs infinite duration sequences? What are difference between finite duration and infinite duration sequences?
As far as i am able to google and know only one difference that is finite duration sequence ROC exist on
xy plane while infinite duration sequence ROC exist on unit circle plane

I have also attached two snapshots in this regard

There is a slight misunderstanding here. There is no such plane called unit circle plane. Unit circle is in the xy plane (or we call it z plane). The difference is, rather, in the shape of the ROC. For infinite duration sequences, the ROC is in the form of a circular strip in the xy plane (as shown in the second figure).

In the case of finite duration sequences, ROC is the entire xy plane except at z=0 or at z = $\infty $ or at z=0 and $\infty$.

e.g. $x(n) = \{1,\underset{\uparrow}{2},3\} $

$X(z) = \sum \limits_{-\infty}^{+\infty} x(n)z^{-n}$

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

This Z transform will not exist only if $z=\infty$ (for the first term) and $z=0$ (for the last term). This can be generalised for all finite duration signals that their ROC consists of entire z plane (xy plane) except at z=0 or z=$\infty$ or at both.


Determining LTI system response to a scale change We know that for an LTI system, if $y(t)$ is the output for $x(t)$ then the response for $x(t-2)$ will be $y(t-2)$ and so on.

But my question is what will be the system response for the input $x(-2t)$ ?Is it possible to determine the output without knowing $h(t)$?

I think it's instructive to look at this in the time domain. For LTI systems, the output $y(t)$ is given by the convolution of the input $x(t)$ with the impulse response $h(t)$:

$$y(t)=(x\star h)(t)=\int_{-\infty}^{\infty}x(t-\tau)h(\tau)d\tau\tag{1}$$

From $(1)$, the response to $x_a(t)=x(t-a)$, $a\in\mathbb{R}$, is given by $y(t-a)$:

$$\begin{align}(x_a\star h)(t)&=\int_{-\infty}^{\infty}x_a(t-\tau)h(\tau)d\tau\\&=\int_{-\infty}^{\infty}x(t-a-\tau)h(\tau)d\tau\\&=y(t-a)\end{align}\tag{2}$$

which proves time-invariance.

However, the response to $x(at)$, $x\in\mathbb{R}$, is

$$\int_{-\infty}^{\infty}x(a(t-\tau))h(\tau)d\tau\neq y(at)\tag{3}$$

Consequently, in general we cannot deduce from the response $y(t)$ to an input $x(t)$ the response to a time-scaled input $x(at)$.


Why a Convolution Matrix Is Called a Kernel? According to Wikipedia:

  In image processing, a kernel, convolution matrix, or mask is a small matrix.

I am wondering, why the matrix is called a kernel? Does it has anything to do with kernel concepts in mathematics (e.g., as those mentioned in this post)?

In general, a discrete (1D) transform can be defined as:

$$ y[k] = \sum_n \phi[k,n] x[n] $$

where the input function $x[n]$ is mapped into the output function $y[k]$, by the kernel, $\phi[k,n]$, of the transform. For example, the kernel for DFT (discrete Fourier transform) is  $\phi[k,n] = e^{-j \frac{2\pi}{N}kn}$.

In essence the kernel acts on the input and produces the output by modifying it, like in the sum above. Therefore the kernel defines the character of the mapping.

Programmatically, a kernel is the core portion of a function which accepts a sequence as its input and returns a computed output whan called.

The input-output relationship for a discrete-space 2D LTI system (typically an image filter) is expressed as a convolution sum :

$$ y[n_1,n_2] = \sum_{k_1=-\infty}^{\infty} \sum_{k_2=-\infty}^{\infty} h[n_1-k_1, n_2-k_2]x[k_1,k_2] $$

The convolution operation is similar to the transform considered above with its kernel being the impulse response $h[n_1, n_2]$ of the LSI system.

For a discrete-time 2D LSI system whose impulse reponse $h[n_1, n_2]$ has a finite domain of support; i.e, finite size, then programmatically the system can be represented as a matrix as well. And the convolution operation than is performed by calling a function which contains the impulse response matrix as its kernel and producing the output pixels per call. And here is the most typical meaning of a (filter) kernel in image processing.

Furthermore, a matrix called mask operates on a block of pixels (sample wise) on the image being processed. They are more of a programming concept than of signal processing, but also find applications such as windowing or frequency domain filtering. Masks find extensive usage for image processing effects. Their operating matrix can also be called as a kernel.


Replicate MATLAB's `conv2()` in Frequency Domain When conv2d is on valid mode, the image needs no padding, because the result is the same size as the image.
When conv2d is on full mode, the result is (image_width + kernel_width -1) * (image_height + kernel_height -1).
Then how do I pad the image?

I created a MATLAB function which is basically conv2() in Frequency Domain:
function [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mO ] = ImageConvFrequencyDomain( mI, mH, convShape )
% Applies Image Convolution in the Frequency Domain.
% Input:
%   - mI                -   Input Image.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double' (Single Channel).
%                           Range: (-inf, inf).
%   - mH                -   Filtering Kernel.
%                           Structure: Matrix.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - convShape         -   Convolution Shape.
%                           Sets the convolution shape.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mI                -   Output Image.
%                           Structure: Matrix (Single Channel).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'conv2()' - https://www.mathworks.com/help/matlab/ref/conv2.html.
% Remarks:
%   1.  A
% TODO:
%   1.
%   Release Notes:
%   -   1.0.000     29/04/2021  Royi Avital     RoyiAvital@yahoo.com
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

CONV_SHAPE_FULL     = 1;
CONV_SHAPE_SAME     = 2;
CONV_SHAPE_VALID    = 3;

numRows     = size(mI, 1);
numCols     = size(mI, 2);

numRowsKernel = size(mH, 1);
numColsKernel = size(mH, 2);

switch(convShape)
    case(CONV_SHAPE_FULL)
        numRowsFft  = numRows + numRowsKernel - 1;
        numColsFft  = numCols + numColsKernel - 1;
        firstRowIdx = 1;
        firstColIdx = 1;
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
    case(CONV_SHAPE_SAME)
        numRowsFft  = numRows + numRowsKernel;
        numColsFft  = numCols + numColsKernel;
        firstRowIdx = ceil((numRowsKernel + 1) / 2);
        firstColIdx = ceil((numColsKernel + 1) / 2);
        lastRowIdx  = firstRowIdx + numRows - 1;
        lastColdIdx = firstColIdx + numCols - 1;
    case(CONV_SHAPE_VALID)
        numRowsFft = numRows;
        numColsFft = numCols;
        firstRowIdx = numRowsKernel;
        firstColIdx = numColsKernel;
        % The Kernel when transformed is shifted (Namely its (0, 0) is top
        % left not middle).
        lastRowIdx  = numRowsFft;
        lastColdIdx = numColsFft;
end

mO = ifft2(fft2(mI, numRowsFft, numColsFft) .* fft2(mH, numRowsFft, numColsFft), 'symmetric');
mO = mO(firstRowIdx:lastRowIdx, firstColIdx:lastColdIdx);

end

It is fully compatible and validated.
The full code is available on my StackExchange Signal Processing Q74803 GitHub Repository.


Automatically crop out certain region around red objects I have the image seen below and I would like be able to automatically crop out everything around each red section - not enough to overlap but enough that you can see how connected to the green parts they are. Is there an easy way to do this, ideally in python/R? In case it matters, I'll be doing some image recognition on the resulting images.

The top image is an example of what I'd like them to look like, the bottom is what I start with.

Sorry if this is stupidly easy or really general, I've never done any work with images so don't even know where to start.

I can't give you a detailed workflow for it right now, but I'll tell you how I'd go about it:

A simple solution would be to color-threshold the image. Your image only uses red and green so it would work quite nicely. I used a red value of 122-255 and a green value of 0-122. This way you obtain a b/w image with the red regions white and the rest of the image black / or the other way around.

On this segmented image you can perform particle analysis. Particle here refers to one white / black region that is separated from the rest (in your case one red blob). It returns the x/y location as well as the extent etc. of the region. this you can use to automatically crop the image to the region of the particle +/- x pixels to obtain the cropped images.

I guess this method won't spare you of cross checking the automated results and deleting falsely detected noise etc though.

I alwaya use ImageJ/FiJi to explore workflows like this, and once I roughly know how I want to proceed, I look for python solutions.
If oyu want to do this in python I suggest you either look into OpenCV or scikit-image.

Stuff I'd read into for your case:

rag thresholding
particle analysis

Once you got the coordinates of the red regions, cropping it is as easy as slicing a numpy array: img_cropped = img[x_coord - extent : x_coord + extent, y_coord - extent : y_coord + extent]


Input image size in MATLAB Hi I am trying to get the image size of the uploaded image in MATLAB. So far I have got :

read=imread(image.png)

ymax=size(read) (1);
xmax=size(read) (2);

However this gives me an error. The image size is 300px by 256px, and I want my ymax and xmax to automatically get the pixel values. How do I do this?

It's not a MATLAB'ish way of doing that, try:

[xmax, ymax, nchan] = size(read);

Where nchan is the number of channels. Also it can be done by:

xmax = size(read, 1);
ymax = size(read, 2);

Or even:

[xmax, ymax, ~] = size(read);


Multi-statelevel detection I'm trying to develop an automated state-detection algorithm, basically expanding on MatLab's "statelevels" function to multiple (i.e. 3+) states, however I've hit a brick wall in trying to figure out how to define histogram bin edges without prior knowledge of the signal and have not been able to find any literature which isn't locked behind a paywall. Does anyone have suggestions, either for an algorithm to achieve multi-state detection or for publicly-available literature?

You're looking for cluster detection. In your case, it's very basic: there's no more than one dimension along which to detect clusters.
It's probably overkill, but 1-dimensional k-means is probably what you want.
Matlab has good documentation on k-means. (Also, seriously, all this is wikipedia or textbook knowledge, there's very many free resources. You very rarely start reading research papers when you enter a new field.)


All-pass filter implementation in numpy I'm playing with Python & numpy to filter audio. It is relatively straightforward to design and us low-high-bandpass/stop FIR and IIR filters, stuff can be done in a couple of lines.

I'm implementing a small software defined radio software, and I need an all-pass filter that has 90 degrees phase shift everywhere, and is flat (as much as possible).

Unfortunately I was unable to find any quick solution for this.

Please share your thoughts on this.

Thank you,

Tamás.

The Hibert transform will produce a signal with 90 degree phase shift.  Look in Numpy documentation at scipy.signal.hilbert as a start.

I'm not sure what you are actually trying to do, but you mentioned 90 degree phase shift and radio in the same post, so you may find the following article of use:

http://home.comcast.net/~szemengtan/LinearSystems/hilbert.pdf


Separate two sinusoids with very close frequencies (5 Hz difference) in a 44.1 Khz audio signal, with good time-domain resolution Here is an audio signal which is a mix of (WAV file available here):

a 2000 Hz sinusoid, beginning at 1.00 sec,  without fade-in/fade-out

a 2005 Hz sinusoid, beginning at 1.031 sec, with a slow fade-out at the end

background noise

(In reality it's even more complex: the sinusoids can also vary in amplitude...)
How to separate the signal into two signals (the two sinusoids), but also with a good temporal resolution?
Attempt 1
I zero-padded the signal to the next power of 2 (final length: 524288), did a real-FFT.
The size of the real-FFT vector h is 262145 frequency bins to cover the frequency range [0, 22050hz], so each bin has a width 0.084 Hz. Pretty good news, we can distinguish the two sinusoids with this!

Now we can isolate the two sinusoids with:
h1 = h.copy()   # the real-FFT
h1[:23750] = 0  # zeroing bins outside [23750, 23800]
h1[23800:] = 0
x1 = irfft(h1)  # inverse real-FFT

h2 = h.copy()
h2[:23810] = 0  # zeroing bins outside [23810, 23860]
h2[23860:] = 0
x2 = irfft(h1)  # inverse real-FFT

It works, but the time-domain resolution is very bad (it's normal because the frequency resolution is very high!) : the sine is very poorly localized in time-domain. Instead of a fast attack in the separated sinusoid, we have a slow fade-in...
Of course, zeroing bins outside of the frequency range of interest (bins [23750:23800]) is not optimal, I should have used a (non-rectangular) window in the frequency domain instead:
h1[:23750] = 0
h1[:23800] = 0
h1[23750:23800] *= window

But even with such a windowing I doubt I can avoid a slow time-domain resolution after the separation.
Attempt 2
Use a STFT instead of a global FFT of the signal. This helps to localize, but ... in order to have a good frequency resolution to be able to separate the two sinusoids, we have to take a big FFTSIZE, such as 16384. Then each of the
8193 frequency bins (real-FFT) will have a 2.7 Hz width!
Not enough to distinguish or separate the two sinusoids that have only 5 Hz of difference... So this approach will fail.

I know this is probably an example of the time-frequency trade off / uncertainty principle, but in this precise case, is there something else we can do to improve the separation?
Remark: I already read Separating waves of very close wavelengths but it did not help for this case.

I believe there are a couple of reasons why your answere is not coming out as expected.  What you have in your problem is the spectrum of the sum of three signals

$G(\omega) = S_1(\omega) + S_2(\omega) + N(\omega)$

where $G(\omega)$ is the spectrum you showed, $S_1(\omega)$, and $S_2(\omega)$ are the spectrum of the signals of interest and $N(\omega)$ is the noise spectrum.  You are attempting to recover signals $S_1(\omega)$, and $S_2(\omega)$ with high fidelity, which may not be possible in this case.  For simplicity, I won't worry about noise here.

Your solution assumes a couple of things.  First it assumes a very small bandwidth for each signal.  You mention that the recovered signals have a slow fade in, which would indicate a band-limited signal.  As an example, a infinite length sinusoid in the time domain would be delta function in frequency.

Also, the spectrums $S_1(\omega)$, and $S_2(\omega)$ overlap, which is problematic.  Think of it as a system of linear equation.  In your case, you have two unknowns, $S_1(\omega)$, and $S_2(\omega)$, but only one equation.  This leads to an underdetermined system that has 0 or infinite solutions.

To overcome this situation you may attempt to assume additional information about the signals $S_1(\omega)$, and $S_2(\omega)$.  Your previous assumptions of limited bandwidth and non overlapping spectrums lead to your previous answer.  In looking toward a better solution, I would suggest new assumptions.  For example, you could attempt to model your signals $s_1(n)$, and $s_2(n)$ as windowed sinusoids and attempt to estimate window function from the data.  Additionally, you may attempt to use a windowing function on the time domain signal to limit the side lobes of the frequency spectrum to improve isolation as well.


Comparing Images by Defining a Reference Model I have a microscopy system which periodically captures an image to be used for shading correction.  Because these images are captured on tissue slides containing, well, tissue, as well as possibly debris, I would like to define a known, good image as the reference with which to compare subsequent images.

These images are captured over clear glass using a line scan camera, so, optimally, there will be little variation across the Y axis, and a quadratic like curve across the X which models the light drop-off at the edges of the field of view of the camera.  An intensity profile, averaged along each column, looks like this:

As you can see, I'm currently using a quadratic fit to model the data.  What I'm now grappling with is what measures to use to compare this image to a future sample.  My initial approach was/is statistical in nature; mainly, comparing residuals.  The difference in residuals between the reference and an arbitrary sample is currently my strongest measure, but I'm not a statistical whiz, and I'm not convinced it is the correct approach yet.

For reference, here are a few samples which should be rejected.  I can add raw data if needed, but I'm really just looking for ideas from the DSP side of things.  The green line is the plot of the reference image, the red line the fit for the image in question.  That last one was captured over faint tissue, and it's a bugger if I continue to use a quadratic fit model.  I may end up using a different approach for that situation.

When looking at your data I want to say that when you have a good image the intensity along the x axis will be somewhat symetric and be much more random when you have a bad image.
The third order statistical moment, the skewness (1rst and second being mean and variance), is a good measure of assymetry. You may want to compare your image based on the skewness of the intensity along the x axis. Set some thresold either manualy or with basic machine learning.


z-transform causality properties: negative coefficents are zero ($x[-1]z^1=0$, $x[-2]z^2=0$, ...) Let's suppose I have a system:

$$Y(z)=X(z)H(z)$$

If the system is causal, does that mean that all the negative coefficients (example: x[-1]) of the transform for $Y(z)$, $X(z)$, and $H(z)$ are zero?  or does that rule apply only to $H(z)$?  what about $Y(z)$ and $X(z)$?

No. The system response is usually $H(z)$. If $H$ is causal, then it has no non-causal coefficients.  $X$ may very well be a signal that exists for all time (and so existed before time 0) and, therefore, have non-causal coefficients.  That means $Y$ will also have non-causal coefficients.

If $Y$ is taken as the system and $Y$ is causal, then your statement is true. But usually $Y$ is the output, $X$ is the input, and $H$ is the system response.


Number of zeros of a sum of Shah functions by applying Rice's formula? There is a Dirac pulse train following the scheme of the Shah function (or $\delta$-cumb function) with its Fourier series of the form:

$$\varsigma(t,T)=\sum_{n=-\infty}^{\infty}\delta (t-nT)=\frac{1}{T}\sum_{n=-\infty}^{\infty}\exp\left(\frac{i 2\pi n t}{T}\right)$$

here $n\in \Bbb Z$; $t$ time and $T(\in \Bbb N)$ the period between the pulses, which occur exactly at the levels of the integers.

The sum of $\varsigma(t,T)$ over $T$ in a given range $[1,T_{\max}]$ and $n$ in a range $[1,T]$ (instead) is accordingly:

$$S(t,T_{\max})=\sum_{T=1}^{T_{\max}}\varsigma(t,T)=\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t}{T}\right)$$

And the zeros (at the level of some integers) of this function with respect to $t$, hence:

$$\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t_0}{T}\right)=0$$

The Question is to identify the number of such zeros $\mathcal Z$ in a range of $t\in[T_{\max},T_{\max}^2]$ by applying Rice's formula, while assuming that the frequencies $1/T$ are rationally independent?

This is not a complete answer in that it does not apply Rice's formula although it does derive the formula for the counting function $\mathcal{Z}(T_{max})$.

For $t,t_0,T,T_{max}\in \mathbb{N}$:

$$\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t}{T}\right)=\cases{1&$T|t$\cr 0 &otherwise\cr}$$
so:

$$S(t,T_{\max})=\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t}{T}\right)$$
counts the number of divisors of $t$ that are $\le T_{max}$ and $t_0$ is defined to be a zero if it satisfies:
$$\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t_0}{T}\right)=\sum_{T=1}^{T_{\max}}\cases{1&$T|t_0$\cr 0 &otherwise\cr}=0 \tag{1}$$
which holds iff $t_0$ is not divisible by any integer $1\le T\le T_{max}$. However $t_0$ is always divisible by $1$ so the sum in $(1)$ must always evaluate to something $\ge1$ for any $T_{max}$. In which case we are forced to say

there are no integers $t_0,\,T_{max}$ such that $(1)$ holds and thus:
$$\mathcal{Z}(T_{max})=0, \,\forall T_{max}\in \mathbb{N}.\tag{2}$$

Alternatively, if we chose to exclude the trivial case in which $1$ divides all integers by beginning the outer sum in $(1)$ at $T=2$ we may find a less trivial outcome. Define $t_0$ with $2\le T_{max}\le  t_0\le T_{max}^2$ to be a solution to the equation:
$$\sum_{T=2}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t_0}{T}\right)=0$$
then $\mathcal{Z}$ counts the number of integers $t_0 \in \left[T_{max},T_{max}^2\right]$ such that $T\nmid t_0 \forall\, 2\le T \le T_{max}$, in which case:
$$\mathcal{Z}=\pi(T_{max}^2)-\pi(T_{max}) \tag{3}$$
where $\pi(x)$ is the prime counting function.

Proof

If: $$\exists t_0 \,\, \text{such that} \,\,\,T\nmid t_0 \,\forall\, 2\le T \le T_{max},$$ but $t_0$ is not prime, then: $$\exists k>T_{max}\,\, \text{such that} \,\,k|t_0,$$ and:
 $$\exists n=\dfrac{t_0}{k}>T_{max},n\in \mathbb{N},$$
which implies:
$$t_0=nk>T_{max}^2,$$
from which we conclude that $t_0 \in \left[T_{max},T_{max}^2\right]$ is not a solution when $t_0$ is not prime. It is clear that if $t_0$ $\left(2\le T_{max}\le  t_0\le T_{max}^2\right)$ is prime, it is not divisible by $2\le T \le T_{max}$ and is therefore a solution, unless $T_{max}$ is itself prime. From this we conclude that all $t_0$ are primes such that $T_{max}<t_0\le T_{max}^2$ and $(3)$ follows.


Difference equations applied to DSP I've found almost nothing about difference equations on the internet. Can you please recommend me something like books or pdf online that handle deep this topic? I'm searching also some exercises. Thank you!

https://ccrma.stanford.edu/~jos/filters/ and specifically
https://ccrma.stanford.edu/~jos/fp/Difference_Equation_I.html
http://www.dspguide.com/pdfbook.htm
https://www.amazon.com/Digital-Signal-Processing-Alan-Oppenheim/dp/0132146355
    (need to buy the book)


GPS signal acquisition I have a set of GPS satellite signal data, the set is 1 ms long and has the information of two satellites (satellite number, code phase and frequency) buried in noise. This data has a carrier frequency of 4.309 Mhz and the sampling rate is 5.714 Mhz.I've been reading but I don't really have pretty clear what to do:

which are the options to make the acquisition? As far as I've read, there are two options: time and frequency based acquisition. which is the best option to use?
The set of data is noisy and continuous:

I have generated the C/A gold sequencies sampled at 5.714 Mhz in order to perform the correlation with the satellite data. How can I design a narrowband filter to get rid of the noise? The local coded I generated are sequencies of +1's and -1's. Do I need to convert this data also into -1's and +1's to perform the acquisition?

If the GPS signal is centered (carrier) at 4.309 MHz, and this is sampled at
5.714 MHz, then there would be an alias copy at 5.714-4.309 = 1.405 MHz.  To
acquire this signal, you would first "down-convert" the signal to baseband, or
0 frequency by multiplying by $exp(-j\omega t)$ where $\omega$ is
$2\pi(1.405e6)$ Hz. The reason to use an exponential and complex downconverter
such as this process instead of just multiplying with a sine, is that it will
facilitate fine adjustments to carrier frequency and phase as you do your
acquisition (which you can do manually to get a feel for how a carrier tracking
system would work).  After you multiply the signal as such to move it to 0
frequency, then correlate with the code for your specific SV of interest.

When correlating to a known code, the magnitude of the correlation will be
proportional to the residual frequency offset of your carrier as a Sinc
function with the first nulls located at $F= \frac1T$ where $T$ is the length
in time of your correlation. (For example, if you choose to correlate over a
full PRN which is 1 ms, the null to null width will be ± 1kHz, meaning you
should be able to see an SV even if it is still ±500 Hz off. For this reason if
I was doing a manual search with a full length PRN, I would step my frequency
in steps of 500 Hz and do a correlation in each "frequency bin" until I stepped
over the full possible Doppler range for my system. For your case to give a
specific example, I would do a first correlation after multiplying by -1.405
MHz ($e^{-j2\pi1.405e6t}$), and know that I would see and SV if it was within
500 Hz (the null extends to 1kHz offset, but beyond 500 Hz and the correlated
signal would start to compete with noise). If I didn't see anything for an SV
that I knew was present, I would then step my frequency up 1 kHz (to -1.406 MHz)
and try there, and I would continue this over all possible Doppler offsets
(typically up to ±5 kHz). If I knew my signal was strong, I could correlate over
1/2 a sequence to double the "window" in which I would see a signal in Doppler
offset. But this would reduce the processing gain by 3 dB but useful for an
initial acquisition search.

Note, if needed see my code generator available here: https://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator  which includes the ability to resample the code to any arbitrary frequency (in your case 5.7144 MHz).

And if you really want to get fancy (at the expense of resources associated with a 2D FFT), see my implementation for "one-shot" code delay and carrier phase here: https://www.mathworks.com/matlabcentral/fileexchange/26035-joint-frequency-and-delay-correlation   This shows how to determine code delay and frequency offset from one block of data in one matrix operation, using the cross-correlation theorem to relate the fft to a circular correlation;

$\text{xcorr}= \text{ifft}\left(\text{fft}\left(a\right)\text{fft}\left(b\right)^*\right)$

combined with the property that a circularly rotated FFT of a time domain signal is the same as the time domain signal in a different Doppler bin. (So once one FFT is computed for the CDMA PRN, you have the FFT's for the PRN at all other Doppler bins just by rotating that one FFT!).   The code provides a similar plot as shown below which is the correlation for an SV vs code and frequency offset.


Sound card for recording audio I need to record sound using 2 (or 4, eventually) microphones and then calculate the time delay of arrival of sound at the different microphone pairs. But first, I need a microphone array. And since I'm just starting out, I would like to connect just 2 mics first and calculate the time at which sound arrives at them.
Can I do this using a sound card?
And I heard sound cards will have typically one or two input channels. What are these then? http://www.ebay.com/bhp/external-sound-card-for-laptop
The 'channels' mentioned along with the above devices are not exactly input channels?

The number of channels mentioned in the products you linked refer to output channels. A Dolby 5.1 system has 6 speakers (center, rear R/L, front R/L, subwoofer), thus 6 channels.

What you're looking for ist usually called audio interface and mostly used in music production studios. Manufacturers are M-Audio, Roland, Behringer (cheap), Focusrite and many more. You can get an overview here. From your previous question I assume that your application is speaker localization. Therefore you probalby don't need high-class converters and pre-amplifiers. If you're planning to use Linux it will reduce your choice.


DTFT of $(-1)^n \cdot \mathrm{sinc}()$ I'm trying to find the DTFT of $$(-1)^n \cdot \frac{\sin (\pi n/2)}{\pi n}$$ I know the DTFT of $\frac{\sin \pi n/2}{\pi n}$ = a box function of amplitude 1, cutoff $\pi/2$. And I know that multiplication in time domain = convolution in frequency, but I'm still not sure how to piece everything together.

Considering @Stanley Pawlukiewicz's hint, you should take into account that $e^{~j \pi} = -1$, thus $(-1)^n = e^{~j \pi n}$, for any integer $n$.

Therefore, time-domain multiplication by $(-1)^n$ corresponds to a frequency-domain shift by $\pi$ radians in digital frequency.

Recall the Fourier transform's theorem: $\mathscr{F}\{x[n]~e^{j \omega_0 n}\} = X(\omega - \omega_0)$ where $\mathscr{F}\{\}$ is, in this case, the Discrete-Time Fourier Transform (DTFT).

Now, since

$$\mathscr{F}\left\{\frac{\sin (\pi n/2)}{\pi n}\right\} = \sum\limits_{k=-\infty}^{+\infty} \Pi\left(\frac{\omega + 2k\pi}{\pi}\right)$$

then, answering to your question:

$$\mathscr{F}\left\{(-1)^n \cdot \frac{\sin (\pi n/2)}{\pi n}\right\} = \sum\limits_{k=-\infty}^{+\infty}\Pi\left(\frac{\omega - \pi + 2k\pi}{\pi}\right)$$

which is a box function of amplitude $1$ and width $\pi$, centered on $\omega = \pi$.  The box function or rect function is

$$\Pi(u) \triangleq \begin{cases}
0 & \mbox{if } |u| > \frac{1}{2} \\
\frac{1}{2} & \mbox{if } |u| = \frac{1}{2} \\
1 & \mbox{if } |u| < \frac{1}{2}. \\
\end{cases}$$

Since the DTFT is periodic in $\omega$ with period $2\pi$, there is a box centered on $\omega = (2k-1)\pi$, for every integer $k$.


Nyquist rate and oversampling Assuming we have time domain signal $$x(t)\,\,\,\,\,\,\, 0\leq t< T$$ and we decide to upsample the signal by factor $L\geq1$.

Then the samples of this oversampled time domain signal are represent as

$$x(kT/L)$$

I have read in paper that setting $L=1$ corresponding to the case of Nyquist rate sampling, does anyone know why this is case

If $L=4$ how many zeros are padded to obtained the upsampled signal?

I provided the text that made me ask the question in first place.

As stated in the paper citation, $(L-1)N$ zeros are padded in frequency domain to achieve an $L$ times oversampling, where $N$ is the number of subcarriers.

Example: Let $N=256$ and $L=4$. For the PAPR simulations use a 1024-IFFT. The 256 "center" inputs of the IFFT are the original OFDM frequency domain signal, the rest is zero. Thus $(L-1)N = 768\,$ zeros have been inserted.

"Nyquist-rate sampling" here means that all IDFT input bins are used and are in general unequal to zero, as for $L=1$ the number of zero padded subcarriers is 0. Since the IDFT input can be interpreted as a frequency domain signal, this means that the complete available bandwidth from $-f_\mathrm s/2$ up to $f_\mathrm s/2$ carries signal energy. In contrast, some of the available bandwidth is unused if zero padding in the frequency domain is applied. This is referred to as oversampling in analogy to sampling an analog signal at a higher rate than twice its maximum frequency. In the context of this question it is worth noting that the oversampling described in the article is only used for simulation purposes and that the signal will usually be analyzed at the IDFT output. A transmission of the oversampled signal is usually not done as the purpose of oversampling here is only PAPR analysis.


Difference between real fft and complex fft with imaginary part of zero in fftw? I have a real 2d matrix. I am taking its fft using fftw. But the result of using a real to complex fft is different from a complex ( with imaginary part equal to zero) to complex fft.

real matrix

 0     1     2
 3     4     5
 6     7     8

result of real to complex fft

36 -4.5+2.59808i  -13.5+7.79423i
0  -13.5-7.79423i 0
0  0              0

Code:

int r = 3, c = 3;
int sz = r * c;
double *in = (double*) malloc(sizeof(double) * sz);
fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sz);
fftw_plan p = fftw_plan_dft_r2c_2d(r, c, in, out, FFTW_MEASURE);
for ( int i=0; i<r; ++i ){
    for ( int j=0; j<c; ++j ){
        in[i*c+j] = i*c + j;
    }
}
fftw_execute(p);

using a complex matrix with imaginary part of zero

complex matrix

 0+0i     1+0i     2+0i
 3+0i     4+0i     5+0i
 6+0i     7+0i     8+0i

result of complex to complex fft

36               -4.5 + 2.59808i  -4.5 - 2.59808i
-13.5 + 7.79423i 0               0
-13.5 - 7.79423i 0               0

Code:

int r = 3, c = 3;
int sz = r * c;
fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sz);
fftw_complex *inc = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * sz);
p = fftw_plan_dft_2d( r,c, inc, out, FFTW_FORWARD,FFTW_MEASURE);
for ( int i=0; i<r; ++i ){
    for ( int j=0; j<c; ++j ){
        inc[i*c+j][0] = i*c+j;
        inc[i*c+j][1] = 0;
    }
}
fftw_execute(p);

I am after the result of complex to complex fft. But the real to complex fft is much faster and my data is real. Am I making a programming mistake or the result should be different?

I haven't tried this myself but I think it's just a problem of ordering the output data. In the real-to-complex case, your output matrix is not a $3\times 3$ matrix, but - due to the implied symmetry - a $3 \times 2$ matrix. So re-arranging your output gives

 36            -4.5+2.59808i
-13.5+7.79423i  0
-13.5-7.79423i  0

where the "missing" third column would just be the complex conjugate of the second column. This is also explained here.


Practical Least Squares FIR filter design: How to find a realistic desired response I'm using Matt. L.'s lslevin.m function to design long ($N\approx 400$) filters to an arbitrary phase and magnitude specification. The specification is the result of a numeric optimization and does cover frequencies between 200Hz and 19kHz. I found that the designed filters show poor approximation of my desired frequency response $H_1[k]$:

Apparently my desired response it hard to meet. Making the filters longer only amplifies the ripple. So I added a linear phase term $H_2[k] = H_1[k]e^{-j\frac{N-1}{2}\omega[k]}$ in the hope to make the desired response a little more realistic to meet. I am not sure why this is requried or why it makes things better, but it appears to be doing something. (Adding a (predictable) phase term to my filters is okay in my scenario as long as two filters designed with the same length have the same phase term added) Here's the result:

However, outside the range of specified frequencies I see huge spikes in the magnitude response (which is to be expected, because I don't control those regions at all):

So I extended my specification by padding $H_2[k]$ with zeros on both sides. I give a weight of 2 to all points inside my specification and a weight of 1 to all the padded zeros - my hope is that I will get a band-pass behaviour with my desired phase and magnitude response in the passband and some sort of suppression outside the passband. I don#t really care about those regions but the huge magnitude from the previous example makes the impulse response have very large coefficients with varying sign which is impossible to implement. With the padding added, things start to look better.

As you can (maybe) see, there is some ripple above 19kHz but in this case it stays low. If I increase the filter length $N$ to 600 or more - hoping to increase the accuracy of LS approximation in the passband - those ripples become extremely large, with the maxima reaching up to +50dB. Which makes this filter impossible to implement.

So as you can tell, I'm not entirely sure what I'm doing. How can I find a gradient of the linear phase term that yields the best approximation? How can I prevent huge magnitudes outside my important region of 200Hz-19kHz?

In general you seem to be doing all the right things.

The linear phase term is required since an arbitrary magnitude and phase response is typically non-causal. That means it's impossible to fit with a causal filter and your fitting error is very high.

One way to estimate a "good" delay value is to sample your target on an FFT grid that's relatively dense and do an inverse FFT. Ideally this looks like an actual impulse response with a clear peak someplace and trailing off at the edge. Circulate this until you've captured the bulk of the energy.

The other options is to try brute force: just try a lot of delays and pick the one that gives you the smallest error.

Controlling the out of band behavior you can play around with the weights of the
"don't care" bands. It's probably best to reduce the weights until the out of band energy is just about the same as the in-band energy. If you want real bandpass behavior, just apply a regular bandpass on top of your least square design. Trying to get a decent high pass at 250 Hz will cost you a lot of coefficients otherwise.

Your target has a lot of detail at low frequencies which determines the minimum number of FIR coefficients. You could try to prefit the major low frequency wiggles with an IIR filter and then just use an FIR for the remaining difference to the target.


Autocorrelation of noise - negative correlation I am investigating autocorrelation of electrical noise as part of an undergraduate experiment (as detailed in http://physlab.lums.edu.pk/images/a/ab/Correlation.pdf). I sampled noise voltages using an 8-bit AtoD from a noise generator (whose description was not provided). I first sampled the noise directly, then I sampled it after passing it through a low-pass filter. I imported the data into MATLAB and used the autocorr function to autocorrelate it. The no-filter data autocorrelated as expected, being just 1 at t=0 (as in Fig 1 in paper); however, the filtered data decays initially as expected but then dips below 0 before coming back up and staying around 0 as in the image below.

My three issues are:

I am not really sure what this means physically (is the correlation 'direction' changed so-to-speak?),
Why might this be occurring? I know that a description of the noise generator may be important here but I do not have that information myself.
Is there any way to correct this issue, and if not, what would be the best compromise?

EDIT:

Here are snapshots of the noise without and with the filter (RC=0.5ms) respectively:

If the impulse response of the low pass filter goes negative, that means that the filter is adding an inverted or negative image of the input signal to the output at that delay.  Which the autocorrelation is finding it as it is, negative.  The impulse response goes negative as one way of canceling out some of the input spectrum in the filter's stopband.


graph of lowpass filter Check the equation 4.8-1 $H(u,v)=1$ if $D(u,v)\le D_0$

my question is the high frequency is at center.

$D(u,v)$ is smaller if it is close to center.

why $H(u,v)=1$ in that case.
Thanks

There is no high frequency in the center. In the center of your 2D spatial frequency distribution is the 0-frequency (constant DC component of the Fourier-transformation). The text excerpt is in itself consistent. Therefore I suppose you misunderstood the problem or the initial assumptions of its solution.


How to Combine / Fuse 2 Least Squares Estimates? Say you want to compute the least squares estimate of $w$ from a data-set:
$$
\begin{bmatrix}d_1 \\d_2 \\\vdots\\d_N \end{bmatrix} =\begin{bmatrix} x_1 \\x_2 \\ \vdots \\x_N\end{bmatrix}w + \begin{bmatrix} q_1 \\q_2 \\ \vdots \\q_N\end{bmatrix}
$$
Note: $q$ is white noise.
This can be compactly written as
$$
\mathbf{d}= \mathbf{H}w + \mathbf{q}
$$
The least squares estimate of $w$ which I shall call $w^*$ is given by
\begin{align}
w^* =&  (\mathbf{H}^T\mathbf{H})^{-1}\mathbf{H}\mathbf{d}\\
   =  & \Big(\sum_{i =1} ^{N}x_i^2 \Big) ^{-1} \Big(\sum_{i =1} ^{N}x_id_i \Big)
\end{align}

The Problem:

I have $N/2$ of the sample and my friend has the other half of the data set. We compute the least squares estimate separately. How do we now combine our estimates of $w^*$ besides 'averaging' it?

I'm unsure if this is the answer you are looking for, but why not save and share $|H_i|^2$ in addition to the least squares estimates?

If you have $w_1^*=\frac 1 {|H_1|^2} H_1^t d_1$, $w_2^*=\frac 1 {|H_2|^2} H_2^t d_2$, and also know $|H_1|^2$ and $|H_2|^2$, you get the total least squares estimate with:

$w^*=\frac 1 {|H_1|^2+|H_2|^2} (H_1^t d_1 + H_2^t d_2)=\frac{|H_1|^2}{|H_1|^2+|H_2|^2}w_1^*+\frac{|H_2|^2}{|H_1|^2+|H_2|^2}w_2^*$.


What is the difference between `image normalization`, `contrast stretching` and `histogram equalization`? I believe that they all do the same task of distributing the frequency of pixel intensities to a bigger range, but can these terms be used interchangeably?

I believe that they all do the same task of distributing the frequency of pixel intensities to a bigger range

They do.

  but can these terms be used interchangeably?

no, they can't. because the way they increase the dynamic range of the image diffears.

All this methods do a transformation of the pixel intensity $u$ of and image into a new pixel intensity $v$ $$v = f(u) $$ in order to improve the dynamic range of the image.

for constrast stretching (that as far as I know is a more general case of image normalization) the transformation is a piecewise linear function

$$ f(u) = \begin{cases}
      \alpha u & 0 \leq u < a \\
      \beta(u - a) + v_a& a \leq u < b \\
      \gamma(u - b) + v_b & b \leq u < L
   \end{cases}$$

Where L is the total number of gray levels.

On the other hand, Histogram equalization is a more powerful technique for image enhancement that aims to give the histogram of a given image the desired uniform shape. This is achieve by the transformation
$$v_{aux} = F_u(u) = \sum_{i=0}^u p_u(i) $$

$$v = int\left( \frac{v_{aux} - v_{aux_{min}}}{1 - v_{aux_{min}}} (L-1) + 0.5 \right) $$

where $p_u(i) $ is the normalized value of the histogram at the  grey level $i$


How to derive an expression for the optimal importance distribution? I'm trying to answer the exercise 7.6 letter b of this book:

https://users.aalto.fi/~ssarkka/pub/cup_book_online_20131111.pdf

page 133 but I'm having some problems in understanding the question because there is no explanation in the book about how to derive the optimal importance distribution... just references to that ( page 125 ) but the problem is that the references are books that I don't have access to them. If someone can put me on the right direction, it will be very helpful. This is the question:

Update 1: Base on Atul Ingle answer, I believe the answer is something like this:

These images are taken from:

https://www.amazon.com/Beyond-Kalman-Filter-Applications-Hardcover/dp/158053631X

Hint

Well, it looks like the question tells you what the optimal importance distribution is! It's $p(\mathbf x_k | \mathbf x_{k-1},\mathbf y_{1:k})$. So the problem boils down to computing a closed form expression for this conditional density using the given state space model.

More Hints

The transition model tells you $p(\mathbf x_k | \mathbf x_{k-1})$ is a multivariate Gaussian. The observation model tells you $p(y_k|\mathbf x_{k})$ is another multivariate Gaussian.

Can you write these multivariate Gaussians in closed form?
Can you express $p(\mathbf x_k | \mathbf x_{k-1},\mathbf y_{1:k})$ in terms of $p(\mathbf x_k | \mathbf x_{k-1})$ and $p(y_k|\mathbf x_{k})$?

Almost complete solution

\begin{eqnarray}
p(\mathbf x_k | \mathbf x_{k-1}, \mathbf y_{1:k}) &=& p(\mathbf x_k | \mathbf x_{k-1}, \mathbf y_{1:k-1}, y_k) \label{1}\tag{1} \\
&=& p(\mathbf x_k|\mathbf x_{k-1}, y_k) \label{2} \tag{2}\\
&=& \frac{p(y_k|\mathbf x_k,\mathbf x_{k-1}) p(\mathbf x_k | \mathbf x_{k-1})}{p(y_k | \mathbf x_{k-1})} \label{3} \tag{3} \\
&=& \frac{p(y_k|\mathbf x_k)p(\mathbf x_k | \mathbf x_{k-1})}{p(y_k | \mathbf x_{k-1})} \label{4}\tag{4}
\end{eqnarray}
where (\ref{1}) is just another way of writing $\mathbf y_{1:k}$, (\ref{2}) follows from Markov property (Eq. 4.2 of the book), (\ref{3}) follows from Bayes' theorem and (\ref{4}) follows from conditional independence of measurements (Eq. 4.3 of the book).

Using the model given in Eq. (7.49) in the original question we have
$$
p(y_k|\mathbf x_k) \sim  N((1 \;\; 0) \mathbf x_k,10^2) \label{5}\tag{5}
$$
and
$$
p(\mathbf x_{k}|\mathbf{x}_{k-1}) = N\left(\left(\begin{array}{}
1 & 1\\
0 & 1
\end{array} \right)\mathbf x_{k-1}, \left(\begin{array}{}
1/10^2 & 0\\
0 & 1
\end{array} \right)\right). \label{6}\tag{6}
$$
To compute $p(y_k | \mathbf x_{k-1})$ we express $y_k$ in terms of $\mathbf x_{k-1}$ as follows:
\begin{eqnarray}
y_k &=& (1\;\;0)\mathbf x_{k} + r_k \\
  &=& (1\;\;0)
\left(\begin{array}{}
1 & 1\\
0 & 1
\end{array} \right)\mathbf x_{k-1} + (1\;\;0)\mathbf q_{k-1} + r_k \\
&=& (1 \;\; 1) \mathbf x_{k-1} + (1\;\;0) \mathbf q_{k-1} + r_k
\end{eqnarray}
which implies that
$$
p(y_k | \mathbf x_{k-1}) \sim N\left((1 \;\; 1) \mathbf x_{k-1}, \frac{1}{10^2}+10^2\right).\label{7}\tag{7}
$$
Finally you should use Wikipedia or Eq. (A.1) from Appendix A.1 of the book for the expression for a multivariate Gaussian density function.

Time for grunge work. Plug (\ref{5}), (\ref{6}) and (\ref{7}) into (\ref{4}) and finish with some algebra.