Defining the observation matrix for a Linear Kalman Filter So I'm working on some content for an exam tomorrow and have hit a logical snag. The question can be seen below:

In particular I am referring to 5)c).

The questions asks for a measurement model to be derived for this problem.
I am aware that our measurement Zt = CtXt + noise, where Zt represents the measurement, Ct the observation matrix at time t and Xt the state vector at time t.

Xt is defined as (Pt, Vt)^T

One 'solution' I have seen for this problem is to define Ct as [1, 0]. However to me this makes little sense. In the question it states we measure both the position of the robot and the velocity. As such should this observation matrix not be [1,1]? However in which case the result of this multiplication would become = Pt + Vt. Which of course is a nonsensical measurement reading.

Instead we want an obersvation matrix defined as:

As this would give us a Zt measurement result of (Pt, Vt)^T Which the question asks for?

A similar question has been asked before, however the corresponding solution only adds to my confusion, defining an observation matrix of [0,0,1] thereby resulting in two of the three variables in the defined state vector (for that question) being multiplied by 0 and ignored?

Another source has an almost identical question, but states we measure solely the robots position and not both the position and velocity.

Question
Is my thinking right, and therefore the 2x2 identity matrix for the observation matrix correct? Or is there something I am mising?

The first question you link to states that it only measures acceleration and the answer suggests to add the acceleration to the state vector as well, $\begin{bmatrix}x^\top&\dot{x}^\top&\ddot{x}^\top\end{bmatrix}^\top$. This is why the observation matrix becomes $\begin{bmatrix}0&0&I\end{bmatrix}$.

So your reasoning for observation matrix is correct. I only have one remark, namely it is given that the drone moves in 3D space, so the observation matrix should actually be

$$
\begin{bmatrix}
I&0\\0&I
\end{bmatrix}
$$

with $I$ a three by three identity matrix. So effectively the observation matrix would be a six by six identity matrix.


Can edge detection be done in the frequency domain? Can we take advantage of the fact that high frequency components in the FFT of an image generally correspond to edges, to implement an edge detection algorithm in the fourier domain?
I did try multiplying a high pass filter with the FFT of an image. Although the resultant image sort of corresponded to edges, it wasn't exactly the edge detection established using convolution matrices.
So is there any way you could do edge detection in the fourier domain, or it's not possible at all?

Since convolution in spatial domain is multiplication in the Fourier (frequency) domain, you can perform edge detection in Fourier domain by multiplying the spectra of image and the edge detection kernel and then perform IFFT on the result.

I think the high-pass filter alone is not appropriate for edge detection since it keeps all high-frequency features (e.g. sharp peaks and corners) which are usually not classified as edges.

More advanced edge detection methods would be tricky in frequency domain since edges are best described in spatial domain (in my opinion).

The question is why to do edge detection using FFT in the first place? Is it because performance considerations? If so, maybe the high-pass filtered image (produced quickly by FFT) can be quickly filtered again to remove non-edge parts.


Why the time step series affect the noise in FFT results? I used Matlab to calculate a sinusoidal function:

Cy2 = 0.0226809015035037 - 0.235834694520537*sin(5.06307730580011 - 2*pi*2*t);

If I choose a set of Matlab generated data:

dt=5.305164999999997e-04;
t=[dt1:dt1:245];

Then conduct FFT on y(t), I have a nice and clean result (black line):

However, if I use my experimental data's time steps 't_exp' (download here)

t=t_exp %t_exp is the time steps corresponding my experimental data

I get much more noisy data (black line):

I just wonder how come the extra noise while the function is exactly the same, and only a different time step series are chosen.

Conversely, I got noisy experimental data with such time step series (t_exp). Is it possible to make it cleaner by modifying the time step series? I almost mis-interpreted my exp data due to the noise.

Appreciate any help.

==========================

Appendix:
my code to do FFT:

function [f, P1]=fftzp2(t,y)
zpf=1;
%zpf=2;
L=length(t);
dt=t(20)-t(19);
Fs = 1/dt;
Y = fft(y,zpf*L);
%Y = fft(y,L);
P2 = abs(Y/L);
P1 = P2(1:ceil(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:ceil(L/2))/(zpf*L);
%% energy conservation check
Ent=sum(y.^2)*dt;
df = Fs/L;
Enf=sum(abs(Y*dt).^2)*df;
if Ent-Enf>1e-9; disp('energy not conserved'); pause;end
end

Just use MATLAB's interp1 function before FFT


inverse FFT when only a few time points are needed I am developing C++ code to do some modeling that will take as input a parametric model of the complex frequency response (amplitude, phase) and produce as output an impulse response, e.g. using the iFFT. I would like to compare this to a "windowed" impulse response that I have measured for a loudspeaker in a room. The measured impulse can be characterized by an initial portion that corresponds to sound traveling directly through the air to the microphone that lasts a few millisecond, after which sounds that reflect off of the surfaces in the room (floor, walls, etc.) contaminates the signal. The uncontaminated portion is only 3% of the total measured impulse. For instance if I record 16k points at 96kHz, its only the first 470 points of the impulse. I want to compare the impulse response that the iFFT generates from the frequency response model to ONLY the uncontaminated portion of the measured impulse.

I could calculate the entire impulse response from the frequency response and then just throw away 97% of the result but this seems to be very inefficient. The iFFT will be calculated many, many times (thousands probably) while my model is being optimized so I want to make sure that I can make it as efficient as possible. At this point my only option seems to be using FFTW and then just throwing away the data that is not needed (for lack of a better idea).

Is there a fast way to calculate the inverse FFT only for those 470 time points of interest, e.g. not for the entire time span that the FFT can access? I am not intimately familiar with the computation of the FFT and iFFT, so I don't have insight on the answer to this question.

What you're looking for is called a pruned DFT. In principle, it is possible to calculate a subset of outputs from a DFT using fewer mathematical operations. In practice, however, existing highly-optimized FFT implementations like FFTW are designed for full-output transforms. You'll find in many cases, unless you're only concerned with a very small proportion of the transform outputs, it is fastest to just calculate the entire thing and throw away what you don't want.

As it turns out, FFTW has a page in their documentation that talks about this specifically.

If you want to put in a lot of work, it is possible that you could obtain some performance gain by rolling your own FFT implementation for your specific case. FFTW is a general-purpose, any-size-transform library that works well in a variety of conditions, but it's possible that you could beat it if you tune something specifically for your application. Before expending this effort, I would think long and hard about whether it is truly a limiting factor in your system's performance.


expected value of two LTI output signals multiplied (cross correlation) I have an input signal x (assumed to be iid Gaussian with $\mu=0$, $\sigma^2$) which is fed into two linear systems:

$y_1 = h_1 * x$
$y_2 = h_2 * x$

Now I would like to calculate $\mathbb{E}[y_1 y_2]$. What is the proper way to do this - preferably in the frequency domain?

Background: I want to calculate something like $\operatorname{var}(y_1 + y_2)$. This expands to $\operatorname{var}(y_1)+\operatorname{var}(y_2)+2\mathbb{E}[y_1 y_2]$. Now I am very familar with the first two components - the variance. I know that the variance of a stochastic signal is given by the autocorrelation sequence at position 0. With the Wiener-Kinchin theorem this translates to:

$$
\operatorname{var}(y_1) = r_{yy}(0) = \int_0^{\infty} \Phi_{yy}(f) df
$$

and then:

$$
\cdots = \sigma_x^2 \int_0^{\infty} |H_1(f)|^2 df
$$

Now for my problem - with $\mathbb{E}[y_1 y_2]$ I arrive at $\mathbb{E}[x^2 \cdots]$ (giving a variance) but the filter is not square magnitude.

Assumimg that the two linear systems are BIBO-stable, the random processes $\{Y_1(t)\}$ and $\{Y_2(t)\}$ are zero-mean WSS Gaussian processes with autocorrelation functions and power spectral densities given by
\begin{align}
R_{Y_1} &= \sigma^2 (h_1 \star \tilde{h}_1)\\
R_{Y_2} &= \sigma^2 (h_2 \star \tilde{h}_2)\\
S_{Y_1} &= \sigma^2 |H_1|^2\\
S_{Y_2} &= \sigma^2 |H_2|^2
\end{align}
In fact, the processes are also jointly Gaussian and jointly WSS processes with cross-correlation function $$R_{Y_1, Y_2}(\tau) = E[Y_1(t), Y_2(t+\tau)] = \sigma^2 (h_1 \star \tilde{h}_2)$$ and cross-power spectral density
$$S_{Y_1,Y_2}(f) = \sigma^2 H_1(f)H_2^*(f).$$
The OP wants to find $E[Y_1(t)Y_2(t)]$ which is given by
\begin{align}E[Y_1(t)Y_2(t)] &= R_{Y_1,Y_2}(0)\\
&= \sigma^2 h_1\star \tilde{h}_2\big|_{\tau=0}\tag{1}\\
&= \sigma^2\int_{-\infty}^{\infty} H_1(f)H_2^*(f) \,\mathrm df\tag{2}
\end{align}
since the OP prefers the frequency-domain calculation. Personally, given only $h_1$ and $h_2$ (and not $H_1$ and $H_2$), I would say that it is easier/cheaper to just grind out $(1)$ rather than use the frequency-domain calculation $(2)$ but the OP might have specific reasons for opting for the frequency-domain calculation. In particular, $\operatorname{var}(Y_1+Y_2)$ which is what the OP seems to really want to find is just
$$\operatorname{var}(Y_1+Y_2) = \sigma^2\begin{cases}\displaystyle\int_{-\infty}^\infty |h_1(t)+h_2(t)|^2 \,\mathrm dt,\\
\displaystyle\sum_{n=-\infty}^\infty |h_1[n]+h_2[n]|^2\end{cases}$$
which seems easier than the frequency-domain version of the same calculation, but ymmv.


Fourier Decomposition Hello everyone have a look at this video of Fourier Decomposition of an image(otherwise you can also refer the image which shows few plots of different extracted waves from an image) . We also know that a Fourier series is given as

$$\frac {a_0} 2 + \sum \limits _{m=1} ^\infty (a_m \cos \frac {2 \pi m t} T + b_m \sin \frac {2 \pi m t} T)$$

1.In the given video at bottom middle there is plot of extracted waves.My doubt is whether these waves means cosine (or sine) functions from the Fourier series formula?

2.What is mean by bases functions(or images)? Whether the plot of Extracted waves are called as bases functions here?

Note:Although above Fourier Formula is for 1D signal and video is of 2D image,please keep in mind Fourier series formula for 2D signal.

1) A Fourier serie is written as $$\frac {a_0} 2 + \sum \limits _{m=1} ^\infty (a_m \cos \frac {2 \pi m t} T + \mathbf{i} \, b_m \sin \frac {2 \pi m t} T)$$
Unless your signal is even, i.e. $f(-t)=f(t)$ the fourier decomposition will have a real (cosine) part and an imaginary (sine) part. The middle bottom plot is more likely to be the power spectral density i.e. the fourrier transform's squared modulus.

2) Indeed, the Fourrier basis is a orthonormal basis, in a mathematical sense.
$$
\zeta _m(t)=\cos \frac {2 \pi m t} T + i \, \sin \frac {2 \pi m t} T
$$
$$
\int_{0}^{T}{\zeta _k(t)\zeta _l(t) \, dt}=\delta _{kl} %\quad \forall k,l\in[0,T-1]
$$

The two conclusions above still holds in the 2D case simply use as basis functions:
$$
\zeta _{u,v}(x,y)=e^{-i2\pi(ux+vy)}
$$
Define an image $f(x,y)$ and it's Fourier transform :
$$
F(u,v)=\int \int f(x,y) \zeta _{u,v} \, d x d y
$$
$F(u,v)$ isn't real unless $f(x,y)$ is even respect to $x$ and even respect to $y$, and so when building the Fourier transform's representatio, it is more likely that the represented value is the Fourier transform's squared modulus.
In the 2D case, Fourier's basis function still form a basis, it is complicated to define a base in a functional space. The two main things to note are :
First, the scalar product between basis functions $\zeta _{u_1,v_1}$ and $\zeta _{u_2,v_2}$ is always equal to 0 unless $u_1=u_2$ and $v1=v_2$ then it is equal to 1 (orthonormality). One must then prove that every function $f(x,y)$ is uniquely writtable as a linear combinaison of $\zeta _{u,v}$'s. It is true for periodic function, but I can't provide you the demonstration.

I let you check the definition of a basis for a vector space, if you understand it it will give you a good grasp of what is a basis in general.


Finding the impulse response, given a differential function involving the input and output I'm starting to teach myself some signal processing, and I'm stuck on this problem where I have to find the impulse response $h(t)$ of a causal LTI system given: $$y''(t)+2y'(t)+y(t)=x(t)$$ So far what I've done is taken the Laplace transform of both sides of this equation to get $$(s^2+2s+1)Y(s)-3sy(0)-y'(0) = X(s)$$ So the problem now is that I have all these inital values that I don't know what to do with. If I understand correctly, for all $ t<0$ in a causal system we have the output as zero, but this doesn't include zero. If anybody can point me in the right direction, I would really appreciate it. Thanks.

If you do not select
$$
y'(0) = 0\\
y(0) = 0
$$
then the system will not be a linear, time-invariant filter, so one usually makes that selection and then deals with the special cases for when the initial conditions are not zero.


Multiple filters or not? I'm processing audio signal sampled at 44.1 kHz. After some other processing, I need to filter the signal between ~5-50 Hz. Achieving this for a signal with such a high sampling frequency appears to be difficult. I would like the magnitude response to be nicely flat within the bass-band. No phase distortion is allowed. I'm using MATLAB toolbox fdatool for the filter design, FIR filter created with least squares approximation. After some testing, I found out I could get better magnitude responses if I first decimate the signal with some small factor, and then to do the final filtering. Is this a common thing to be done (how about in scientific research, I'm trying to do this the right way..)? All help appreciated.

Your goal should be achievable. The key driver for the filter's complexity is the size of the transition band, which you haven't specified. If you can live with a slow rolloff, then you can use a relatively low order filter.

If, as is more likely, you want to isolate the band from 5 to 50 Hz with relatively sharp cutoffs, then as you have started to discover, you will want to use a multi-stage approach. Instead of decimating by a factor of several hundred all at once, you decimate by a series small factors (something like 2, 4, or 8). By gradually decreasing the sample rate, you end up with filter designs that are much less complex.

This previous question has a very similar aim to yours, only with slightly different filter specs and sample rate. You should be able to adapt the information there to your application.


Inverse z transform - help How can I find inverse z transform of $$X(z)=\frac{5}{(z-2)^{2}}$$ ?

It is known that x[n] is causal.

EDIT: Here is what I have done. Since signal x(n) is causal, convergence of z transform of that signal will be outside of circle with radius r:

We have in bracket sum which represents z transform of signal:

But I don't know what to do next.

You know the geometric series

$$\sum_{n=0}^{\infty}r^n=\frac{1}{1-r},\quad |r|<1\tag{1}$$

Due to convergence of (1) (for $|r|<1$) you can take the derivative with respect to $r$ by taking the element-wise derivatives of the left-hand side. Equating the derivatives of both sides of (1) gives

$$\sum_{n=0}^{\infty}nr^{n-1}=\sum_{n=1}^{\infty}nr^{n-1}=\sum_{n=0}^{\infty}(n+1)r^n=\frac{1}{(1-r)^2},\quad |r|<1\tag{2}$$

With $r=az^{-1}$ you get the following $\mathcal{Z}$-transform relation:

$$(n+1)a^nu[n]\Longleftrightarrow \frac{1}{(1-az^{-1})^2}=\frac{z^2}{(z-a)^2}\tag{3}$$

where $u[n]$ is the unit step sequence.
With the time-shifting property of the $\mathcal{Z}$-transform (see here) you get the following pair:

$$(n-1)a^{n-2}u[n-2]\Longleftrightarrow \frac{1}{(z-a)^2}\tag{4}$$

With (4) you immediately get

$$\frac{5}{(z-2)^2}\Longleftrightarrow 5(n-1)2^{n-2}u[n-2]$$

Of course, all of this is only valid inside the region of convergence, i.e. for $|z|>a=2$.


How can I find the actual frequency given a normalized frequency? Suppose that we have  some sample:

y =

   31.2241
   -5.9830
   22.6058
   -2.2309
  -42.8272
   29.2850
   -6.2652
   20.8502
   -2.4639
  -44.7525
   30.1104
   -6.0180
   21.4476
   -1.5089
  -45.1826
   29.2744

and somehow I know that it is obtained via a sinusoidal method which contains two sine functions at different frequencies and amplitudes, so I apply the music method:

pmusic(y,4)

and obtain the plot:

I know the location of peaks at the normalized peaks, but I don't know the actual frequency.  How can I find it?

I know the formula is $\frac{\text{normalized frequency}*f_s}{2}$, but what about $f_s$?  Can I find it?

[Y,FS,NBITS]=wavread(FILE) returns the sample rate (FS) in Hz and the number of bits per sample (NBITS) used to encode the data in the music file.

After that, [S,F] = pmusic(Y,4,F,FS) to compute the pseudospectrum at frequencies F in Hz.


Sampling rate in Wifi 802.11 -> 20 MHz enough with a BW of 20 MHz? In a Wifi 802.11n context, I'm looking to find the sampling rate to perform some processing on the signal and I've something I cannot understand : in the book Next Generation Wireless LAN's: 802.11n and 802.11ac, I would like to know the meaning of this sentence :

  In 802.11a, the fundamental sampling rate is 20 MHz, with a 64-point
  FFT/IFFT. The Fourier transform symbol period, T, is 3.2 μs in
  duration and F is 312.5 kHz.

Why is the sampling rate « 20 MHz » ?  The bandwitch of the Wifi 802.11 is 20 MHz so my first idea would be to double this value to validate the Nyquist-Shannon criterion. In addition, the guard interval is not included ?

The radio-frequency signal of 802.11n can have a bandwidth of 20 MHz or 40 MHz. But even if it is 20 MHz, it can be sampled at 20 MSa/s if an I/Q demodulator is used: at the receiver the radio-frequency signal is downconverted to the baseband. The I/Q demodulator has two lowpass output signals (Inhphase/Quadrature or real/imaginary part) with bandwidth 10 MHz each. Two samplers with 20 MSa/s each are then sufficient to fulfill the sampling criterion.

There are other receiver architectures where the received signal is first downconverted to an intermediate frequency (for example 10 MHz) and is then sampled at a higher rate (for example 40 MSa/s). In this case only one sampler is necessary.

The guard interval has only a negligible influence on the signal bandwidth. It is not included in the "Fourier transform period $T$", because it is only added after the transform. The symbol duration after the transmitter IFFT is therefore
$$
T = N/f_\mathrm s = 64/(20\,\text{MHz}) = 3.2\,\mathrm{\mu s}
$$


Basic difference between Fourier transform and laplace transform? I have read few links about difference between Fourier transform and Laplace transform but still not satisfied

Please correct me if i am wrong
Simply put, the main difference between Fourier transform and Laplace transform is that real part is set to zero in Fourier transform while real part is non zero in laplace transform?

Fourier transform is an intuitive tool that's a bridge between domain of physics and mathematics, as it quantitatively describes the periodic content of the signals and also frequency response characterisation of systems that occur in physical (and engineering) applications. The use of frequencies is quite intuitive and consistent at least for stable systems...

However for unstable systems (and signals) Fourier transforms becomes mostly awkward (if not useless) to deal with. However control engineers unavoidably must make frequent use of unstable systems in their works. And for this purpuse, Fourier transform is either insufficient or awkward, hence a generalisation of the existing Fourier transform is made into the Laplace transform which conveniently yields mathematical (complex algebric) descriptions of stable as well as unstable systems which was not possible with the Fourier. The Laplace transform, therefore, includes a region of convergence parameter into it.

Another difference between the two transforms is in the time-domain transient analysis of output of LTI systems driven under nonzero initial conditions which is successfully captured in the Laplace transform only. In the sense that LCCDE with initial conditions are straightforwardly solvable by (unilateral) Laplace transforms whereas the standard FT can only solve LCCDE with zero initial conditions (initial rest)...

For one sided and two sided differences, I think Stanley has things to say...


Reconstructing band-limited signal's integral from discrete samples I'm thinking about building an e-bike consumed power meter.

I'm planning on sampling a low-pass filtered current sensing and voltage analog signals on the battery to e-bike controller connection, then somehow digitally integrating the product of these 2 (instantaneous power) to get the consumed power.

I could over-sample say by 10 or 20X of the analog bandwidth and then use a simple sum of areas of the trapezoids to get a reasonably accurate integral (consumed power).

But I'm thinking maybe there's a relatively simple DSP technique to get the integral without such significant oversampling that can be implemented in an FPGA? (e.g. to get a 5% accuracy with sampling at 3x of the analog bandwidth)

Some digital low-pass filter and then the same sum of trapezoid areas on the result?

Let me rephrase what you want to do as:

  I want to approximately calculate the integral of the product of two non-negative discrete-time signals over time. Each signal is sampled from a band-limited continuous-time signal at a sampling frequency $2N$ times its highest frequency. A method is preferred that gives about 5 % error at most, is computationally efficient and does not require a large oversampling factor $N.$

Oversampling factor

The bandwidth of the product of two signals is equal to the sum of their bandwidths. Therefore an oversampling factor of $N = 2$ will suffice to represent the product without aliasing.

Integration

Integration is a linear time-invariant operation, basically a continuous-time filter, and applying it does not increase the bandwidth of the signal. What we want to achieve is, with the product signal as "input":

$$\xrightarrow{\text{input}}\boxed{\text{integrate}}\xrightarrow{\text{output}}$$

However, the impulse response of an integrator is of infinite bandwidth and sampling it to form the impulse response of a discrete-time integrator would cause aliasing, corrupting its frequency response. We can insert to the signal path an ideal lowpass filter that has its cutoff at the bandlimit of the input signal. It does nothing to the signal. The output will remain the same as before:

$$\xrightarrow{\text{input}}\boxed{\text{lowpass}}\xrightarrow{\text{input}}\boxed{\text{integrate}}\xrightarrow{\text{output}}$$

We can combine the lowpass filter and integrator to a single filtering operation. The impulse response of this lowpass integrator filter is the convolution of the impulse responses of the lowpass filter and the integrator filter. This shuffling of operations is enabled by the algebraic properties of convolution, namely associativity. We still get the same desired output as before:

$$\xrightarrow{\text{input}}\boxed{\text{lowpass integrate}}\xrightarrow{\text{output}}$$

The lowpass integrator impulse response:

$$\begin{align}
h[k] &= \int_{-\infty}^k\frac{\sin(\pi x)}{\pi x}dx \\
\\
 &= \int_0^k\frac{\sin(\pi x)}{\pi x}dx + \frac{1}{2}, \\
\end{align}$$

is sampled from the integral of the sinc function, without aliasing. Sinc is the impulse response of the ideal continuous-time low-pass filter with a cutoff frequency at $\pi$.

Figure 1. Impulse response $h[k]$ of an integrator. The impulse response continues indefinitely beyond what is shown, approaching values 0 to the left and 1 to the right.

The impulse response almost represents the process of calculating the sum of (optionally the current sample and) the past samples. As this sum becomes larger and larger, the difference to true integration becomes arbitrarily small, proportionally. The values of $h[k]$ approaching 1 for large $k$ means that for very old samples just summing them gives virtually no error in the integral estimate.

Summary

To summarize, oversample the input signals by a factor of two and calculate with a large accumulator a running sum of obtained samples of the product. This will have an arbitrarily small error compared to true integration, as enough data gets accumulated.

Analog filtering

You should also consider what the analog filtering of the current and voltage signals does to the integral of their product. This can be a big source of error. For example let's say that both signals equal $\sin^2(x).$ If this is filtered to remove all oscillation, it becomes a constant signal $\frac{1}{2}$. The average value (power) of the square of the unfiltered signal is $\frac{3}{8}$ and for the filtered signal $\frac{1}{4}$.


Calculation of an autocorrelation function A sample of a random process is given as:

$$
x(t) = A\cos(2\pi f_0t) + Bw(t)
$$

where $w(t)$ is a white noise process with $0$ mean and a power spectral density of $\frac{N_0}{2}$, and $f_0$, $A$ and $B$ are constants.
Find the auto-correlation function.

Here's my attempt at a solution:

Let $a = 2\pi f_0t$, and $b = 2\pi f_0(t+\tau)$

\begin{align}
\text{Autocorrelation of } x(t) & = E\left\{x(t)x(t + \tau)\right\}\\
& = E\left\{\left(A\cos(a) + Bw(t)\right)\left(A\cos(b) + Bw(t+\tau)\right)\right\}\\
& = E\{A^2\cos(a)\cos(b) + AB\cos(a)w(t+\tau) + AB\cos(b)(wt)\\&\quad + B^2w(t)w(t+\tau)\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + E\left\{AB\cos(a)w(t+\tau)\right\} + E\left\{AB\cos(b)(wt)\right\}\\&\quad + E\left\{B^2w(t)w(t+\tau)\right\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + E\left\{B^2w(t)w(t+\tau)\right\}\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(R_w(\tau)\right)\\
& = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
\end{align}

The expectation terms with the noise in them all equal $0$ (the last is just the auto correlation of white noise ... hence the simplification above. Using trigonometric identities:
$$
\cos(a)\cos(b) = \frac 12\left[\cos(a + b) + \cos(a - b)\right]
$$

we have:

\begin{align}
\text{Autocorrelation of } x(t) & = E\left\{A^2\cos(a)\cos(b)\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
& = E\left\{\left(A^2\right)\frac 12\left[\cos(a+b)+\cos(a-b)\right]\right\} + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
& = \left(\frac{A^2}{2}\right)\left[E\{\cos(a+b)\} + E\{\cos(a-b)\}\right] + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))\\
\end{align}

We're dealing with constant terms, so expectation term goes away and subbing in our initial conditions we get:
$$
\frac {A^2}2 \left[\cos(2\pi f_o(2t + \tau) + \cos(2\pi f_o\tau)\right] + B^2\left(\frac{N_0}{2}\right)(\delta(\tau))
$$

For some reason I can't help but feel I did something incorrectly calculating that autocorrelation ... it's supposed to be a function of $\tau$, but has a $t$ is in there ... I would very much appreciate it if someone could point me in the right direction, or explain what I messed up. I don't know whether it matters, but in this class we're dealing with only wide sense stationary processes.

I guess you've done almost everything right, but have a problem at the calculation of the expectation value regarding $t$. You should calculate the expectation value of the cosine function. Sadly, it does not simply "go away" as you wrote.

Have a look at the Wikipedia page. There you can find another, more explicit, formula for the auto-correlation function of a function $f(t)$:

$R_{\textrm{ff}}(\tau) = \int\limits_{-\infty}^\infty f(t+\tau)\bar{f}(t) \, \textrm{d}t$ .

(Note, that compared to the Wikipedia page, I have taken the liberty to use the variable $t$ in the integration instead of $u$, which would be the mathematically more accurate version.)

As you can see from this equation, you "integrate away" the dependency on t, and indeed you should be left with a function that is independent of $t$.

Note that there is also an version that does not go to infinite times, but is constrained to a period $T$. Maybe this version is more appropriate in your case. However, the same holds for this version: $t$ is integrated away and should not be a variable in the resulting formula.


Pole locations of Butterworth filter I am reading Proakis book "DSP using Matlab", 3rd edition.
I am reading chapter 8, section 8.3, p. 402, and I am confused regarding the equation of poles (roots of denominator of system function) eq 8.47.
How this equation has been derived from eq 8.46 especially highlighted line how the right most side of eq 8.47 $Ω_ce^{j(\pi/2N)(2k+N+1)}$ has been derived from left side $(-1)^{1/2N} jΩ_c$.

Assuming that you understand the left-hand side of Eq. $(8.47)$, for understanding the right-hand side you need to know that $-1=e^{j\pi}$, and that $e^{j2k\pi}=1$. So in order to obtain all $2N$ roots of $(-1)^{\frac{1}{2N}}$ you rewrite $-1$ as
$$-1=e^{j\pi}e^{j2\pi k}\tag{1}$$
from which you get
$$(-1)^{\frac{1}{2N}}=e^{j\frac{\pi}{2N}}e^{j\frac{2\pi k}{2N}}=e^{j\frac{\pi}{2N}(1+2k)},\qquad k=0,1,\ldots,2N-1\tag{2}$$
And since $j=e^{j\pi/2}$, the final result is
$$\Omega_cj(-1)^{\frac{1}{2N}}=\Omega_ce^{j\pi/2}e^{j\frac{\pi}{2N}(1+2k)}=\Omega_ce^{j\frac{\pi}{2N}(1+2k+N)},\qquad k=0,1,\ldots,2N-1\tag{3}$$
If this is all new to you, you should really read up on complex numbers. They are fundamental in understanding many concepts in signal processing.

PS: There is a typo in the equations in the book. The numerator on the right-hand of Eq. $(8.46)$ should have $\Omega_c$ instead of $\Omega$. Also, the left-hand side of Eq. $(8.47)$ should have $\Omega_c$ instead of $\Omega$.


What's the difference between using DFT, IDFT or DCT to calculate cepstrum of a power spectrum? I've seen different equations that calculate cepstrum from power spectrum, but the equations are not consistent. Some people use Fourier transform, some use the inverse Fourier transform, and some use the Cosine transform. I.e.:

$$ C(k) = IDFT(\log_{10}P(n)) \\
   C(k) = DFT(\log_{10}P(n)) \\
   C(k) = DCT(\log_{10}P(n)) \\
$$

For example to calculate LPCC (Linear Prediction Cepstral Coefficients), I see either IDFT and DFT being used equally commonly, while MFCC (Mel Frequency Cepstral Coefficients) uses DCT exclusively.

Are they different? If so what's the different? If not why is there an inconsistency?

They are all somewhat different, in a family of potential candidates,  depending on the system you use and the features you are interested in.

Cepstra play on anagrams related to standard Fourier terms:  quefrency alanysis, liftering, cepstrum, saphe (for the phase)  etc. Not surprisingly, one can find some mixing in cepstrum definitions. You can think of it as a spectrum of a spectrum, scrambled.

The original cepstrum was a squared magnitude of the inverse Fourier transform of the logarithm of the squared magnitude of the Fourier transform of a signal (Bogert et al. 1963). But people have developed generalized versions (related to homomorphic systems). This is where the apparent inconsistency arises.

From a very mundane point of view, take a transformation $F$ (often orthogonal) that gives you some kind of spectrum. Separate (or not) a phase component ($\angle\left( F\left( x\right)\right)$, which could be the sign, some official phase or a mere constant) and an amplitude component $|F\left( x\right)|$. Put an exponent $p$ to the latter, add a little constant $k$ to it to avoid it from vanishing. Pour everything in a companding function $L$: the  $\log$, or even a Box-Cox transform: $(x^\lambda-1)/\lambda$. Now, transform the result  with $F^*$, which can either be $F$ or its inverse $F^{-1}$:

$$F^*\left(\angle\left( F\left( x\right)\right).L\left(k+ |F\left( x\right)|^p\right)\right)\,.$$

Now you can look at the absolute value of this, or a power again, or some inverse of the function $L$.

The key features are:

what is a meaningful frequency information for you or your system? If your system natively uses a DCT, then it is natural to compute a spectrum in a DCT domain.
computations: indeed, a few orthogonal transforms have been used for cepstra: sine, cosine, or Hadamard, which can be quite fast.
features looked at: from basic to invertible.


synchronizing two music tracks I would like to synchronize two different music tracks (and maintain synchronization) as a DJ would with two records.

So far, I have the beats mapped for each song, so I am able to adjust their rates and time align them accordingly.  What I've found though, with many sample tracks, is that the two tracks will drift out of sync after a while due to the tempos of the songs not remaining constant. I've tried compensating for this by periodically checking the tempos of the songs and re-adjusting their rates so that they are synchronized again by the start of the following measure.

Is there a common (or better) approach of synchronizing two music tracks?

well, it depends on how similar the audio signal from the two tracks are.  i dunno if the two tracks are identical, to within some degree of error or noise, and to within some degree of time alignment, but  somehow, in a manner that is salient to your model or to your problem, you have to have a measure of something in the two signals to align.

so let

$$ a_x(t) = f\left\{ x(t) \right\} $$

$$ a_y(t) = f\left\{ y(t) \right\} $$

be that something that defines what it is you're trying to synchronize.   so you want to find a delay offset $\tau$, that maximizes the cross-correlation between $a_x(t)$ and $a_y(t)$

$$ R_{xy}(\tau) = \int_{-\infty}^{+\infty} a_x\left(t+\frac{\tau}{2}\right) a_y\left(t-\frac{\tau}{2}\right) w(t) \ dt $$

with whatever window function $w(t) \ge 0$ you like best.  find the $\tau$ that maximizes $R_{xy}(\tau)$

so i'm not telling you how to do this using discrete-time code.  you get to figure that out.


Eigenvalues of FIR Convolution Matrix Are the eigenvalues of the FIR convolution matrix the zeros of the corresponding FIR filter?

Suppose I have an FIR filter $H(z) = h_{0} + h_{1}z^{-1} + h_{2}z^{-2}$. I want to implement it using a matrix, so I have:

$$ H = \begin{bmatrix}
    h_{0} & 0 & 0  \\
    h_{1} & h_{0} & 0  \\
    h_{2} & h_{1} & h_{0}
\end{bmatrix} $$

Can I find the zeros of $H(z)$ by finding the eigenvalues of $H$?

I know that in general, the roots of a polynomial are given by the eigenvalues of its companion matrix. However, can I just find the eigenvalues of $H$ for this example instead?

For any convolution matrix (even truncated), your diagonal entries would be $h_0$; hence, the characteristic polynomial, setting $b_0 = h_0-\lambda$,

\begin{align}|H-\lambda I| = \begin{vmatrix} b_0 & 0  & 0\\h_1 & b_0 & 0\\h_2 & h_1 & b_0 \end{vmatrix} &= b_0\,\begin{vmatrix} b_0 & 0\\h_1 & b_0 \end{vmatrix} - 0\,\begin{vmatrix} h_1&0\\h_2&b_0\end{vmatrix} + 0\,\begin{vmatrix} h_1 & b_0\\h_2&h_1 \end{vmatrix} \\ \\
&=b_0\,\begin{vmatrix} b_0 & 0\\h_1 & b_0 \end{vmatrix}\\
&=b_0\left(b_0\,b_0 - 0\,h_1 \right)\\
&=b_0^3
\end{align}

Little induction is necessary to see that for any square matrix of dimension $N\times N$ with an lower left structure and a constant $d$ in each element of the diagonal, the characteristic polynomial is $(d-\lambda)^N$; the Eigenvalues are roots of that polynomial, and hence there's but one Eigenvalue (with $N$ multiplicity): $d$, which in your case is $h_0$.

So, since not all zeros of $H(z)$ are defined by your first filter tap, no, you can't find all zeros of $H(z)$ through Eigenvalues.


Analyze stability of a closed loop system with Bode How can one know, with just looking at the Bode plots, if the closed-loop system is going to be stable or not?
This closed-loop system is unstable, why?

@DSP Novice
@Jake

It was my fault that the input perturbation was misplaced.

Closed loop stable: Yes .
Thank you.

1) To determine if a system is stable by looking at a bode plot, generally, you just need to look at two things: Gain Margin and Phase Margin.

First find the frequency at which you phase plot reaches -180 degrees. In this case, it is 18.6kHz (I am inferring this from the pictures you sent because the plot covers up this crossing). Then find the gain at that frequency, which is -51.9dB. Invert that number and you have your Gain Margin, which is 51.9dB.

Second find the frequency at which the magnitude plot hits 0dB. This would be 3.54kHz. Then find the phase at that frequency, which is -100.4 degrees.Your Phase Margin will be 180 degrees plus this number, so 79.6 degrees.

The reason we find these two numbers is because we want to find out if the system will increase over time as it oscillates. When your signal hits -180 degrees, it has finished an oscillation (-180 degrees = 180 degrees). If you have positive gain when this happens then you signal increases with every oscillation, to infinity, which is unstable.

So if your Gain and Phase Margin are positive numbers then you know you have a stable system.

2) For the above reason, you closed loop system looks stable, but there might be a part of you bode plot you are not showing/is being covered up.


How Butterworth low-pass filter can be applied on a digital signal (i.e. discrete representation of a signal)? I would like to know how to apply a low pass filter (Butterworth) to a digital signal. So, I have some values of the signal, let's say $S(t)$. Those values are equally spaced in time. I have read the Wikipedia's article, but I do not understand how to apply the transfer function $H$. I assume that I have to use a Z-transformation and apply it on $S$ but I am not certain to know how to do that. By the way, the part Digital implementation of this article is not really detailed.

Can someone help me to understand those things ? I would really appreciate any answers :)

P.S. I know there are functions in Python (e.g. scipy) that handle that issue, but I would like to know the mathematics behind it.

Thanks.

What you need to do is find a digital filter $H(z)$ that behaves similar to the $H(s)$ of your Butterworth filter. That is, you need to read about IIR implementation.

There are several methods described in signal processing books about this subject. I recommend you to check for example "Essentials of Digital Signal Processing" from Lathi and Green.

Two of the most common methods are the "Impulse invariance" and the "Bilinear transform".

The impulse invariance requires you to find the impulse response $h(t)$ of your $H(s)$. After doing this, you sample it replacing $t$ by $nT$, to finally Z-transform your $h[n]$ into $H(z)$. Bare in mind that this method only works for band-limited filters (lowpass and bandpass) because of aliasing issues. You should determine an appropriate small enough T to avoid this.
The bilinear transform, on the other hand, consists on applying the following mapping of your variable $s$ into $z$:

$s=\frac{2}{T}(\frac{z-1}{z+1})$

This method has a drawback as well, which is the frequency warping, you can read about it here https://en.wikipedia.org/wiki/Bilinear_transform#Frequency_warping, but is in general easy to counter by pre-warping your original $H(s)$.


Is my algorithm about blob tracking for vehicle following correct? I'm doing vehicles tracking project. The video is a one-way, no vehicles occlusion.
After doing background subtraction, erosion and dilation, I get bounding box objects for each frame, and I can find $centroid(c_x,c_y)$ from that.
So, the next thing I want to do is blobs tracking. Here is what I'm thinking on how to do the blobs tracking:

At first frame, store all (cx,cy,id)
Next frame, calculate distance each new (cx,cy) with each of old (cx,cy)
If distance < threshold, update the old (cx,cy) which qualify the formula
If distance > threshold, check if the cx,cy is near the entry coordinate (e.g. around x,y=0-640,0-40), if yes store new (cx,cy).

Is it right and possible to apply to tracking?

You might want to consider tracking using MeanShift & Cam Shift:

http://opencv-python-tutroals.readthedocs.org/en/latest/py_tutorials/py_video/py_meanshift/py_meanshift.html#meanshift


Code generation for `sum` function in MATLAB In cross correlation we use the below command to calculate the sum after elementary multiplication of the vectors. But when we try to convert this Matlab code to vhdl using HDL coder, this command is not defined by HDL coder. Is it possible to write code for below line without using 'sum' builtin function please.

xc(i)=sum((y1(j1:k1)).*y2(j2:k2));

where j1, j2, k1 and k2 are the starting and ending limits.
y1 and y2 are two vectors and are being multiplied element wise and the result is being added.

Thanks

Because the size of y1(j1:k1) and y2(j2:k2) must be the same then:
$$
k_1-j_1 = k_2-j_2
$$
and you will also have something like: $j_1+N=j_2$ then you could write a for loop like this:

xc(i) = 0
for n=j1:k1
  xc(i) = xc(i) + y1(n)*y2(n+N)
end


How may samples per symbol in ieee 802.11 g I am doing a baseband simulation for IEEE 802.11 g. The standard uses 64 IFFT length ($3.2 \mu s$) and 20 Msamples per second. This result in one sample per symbol. I am wondering what the structure will be if I use 40 Msample per second, i.e will the IFFT size change, or interpolation will be used?

Thanks

If you oversample the signal by a factor of 2 as you've suggested, then you will have a bunch of empty subcarriers that are outside the specified signal bandwidth. When you modulate the OFDM frame to a time-domain signal, therefore, you'll have to use an IFFT that is twice as long. The result will be a discrete-time signal that is interpolated by a factor of 2 compared to using the nominal sample rate.


Image registration / fusion optimal similarity metrics I have a theoretical question about optimal similarity metrics for comparing data sets.

In reading this linked paper, pp. 488-489 [1], I read the following 2 interesting statements. On page 488:

  In the ... scenario in which $A$ and $B$ differ only by Gaussian noise, then it can be shown that $\rm SSD$ [Sum of Squares Difference] is the optimum measure

And on page 489:

  If the intensities in images $A$ and $B$ are linearly related, then the
  correlation coefficient $\rm CC$ can be shown to be the ideal similarity
  measure

Are there are other conditions under which other similarity metrics [e.g. Mutual Information] can be proven to be optimal, or under which the above 2 metrics ($\rm SSD$ and $\rm CC$) are also optimal? The linked PDF focuses on image registration. I am thinking mostly in terms of image registration/image fusion similarity metrics, but would also be interested in more general similarity between any 2 or more datasets (don't have to be images).

[1] P. A. Viola, Alignment by maximization of mutual information. Ph.D. thesis, Massachusetts Institute of Technology, 1995.

I'm not sure what is meant here by the "optimal" similarity metric for a given situation. The textbook you linked to references the following thesis:

  [87] P. A. Viola, Alignment by maximization of mutual information. Ph.D. thesis, Massachusetts Institute of Technology, 1995.

...but I can't find it anywhere online after a quick search.

From my point of view:

Sum of square differences measures the difference in pixel intensities, so it is useful when pixel intensities are very similar in both images
Normalized cross-correlation measures linear correlation between pixel intensities, so it is useful when image intensities are linearly correlated
Mutual information is a measure of general statistical dependence between pixel intensities, so it is useful when image intensities have a non-linear statistical dependence

If SSD works for a given pair of images, I don't see why NCC or MI wouldn't do the job just as well, since they are more general. They may be slower since the metric takes longer to compute, but I don't understand why one would be more "optimal" than the other.


Question regarding sampling rate of sine wave I am trying to determine my sampling rate, or the delay between each sample, but I am a little unsure if my math is correct.
I have generated a sine wave
f = 1000 Hz
which gives me T = 1 ms
I am sampling the sine wave and imported the data into matlab.
From the graph I can see that I have 51 samples per period.
I then calculated the time between the samples to:
T / 51 = 19.6 us
But I am wondering whether this is the correct way to do so?
Thank you in advance.

The assumption here is the sine wave being samples is a known and assumed accurate reference, and the sampling rate is the unknown.
An algorithm that utilizes every sample with equal weight would provide for the best estimate assuming a white noise process. The approach the OP uses, and any similar "edge-detection" methods where we are making a decision on a cycle boundary and then measuring the distance between two such decisions, are prone to significantly higher errors under lower SNR conditions due to the sensitivity of noise on just the two decisions used; as done the noise effecting the result is the rms of the noise for each of the two decisions so is additive as well.
Such algorithms for determining the frequency of a single tone are well documented (see Question regarding estimation of signal tone parameters (frequency, amplitude, and phase) using Macleods algorithm). The approach would be to use these techniques with a scaled frequency designating the tone frequency to be determined as a frequency ratio relative to the sampleing rate $(f/f_s)$ rather than absolute frequency. Thus once $f/f_s$ is determined with such estimators, and with $f$ known, $f_s$ can be accurately determined.
The FFT is one such algorithm that uses every sample to provide a best estimate (under white noise conditions), but for the case of a single tone, the above references would provide more optimized solutions. See https://www.dsprelated.com/showarticle/1284.php where @CedronDawg worked through the math to derive the precise frequency estimate from two adjacent bins in the FFT result (for all the likely cases when the actual frequency is not centered on a bin), so that could be combined with the Goertzel algorithm which efficiently solves for a smaller subset of all DFT bins for a reasonable solution.


Sampling and Reconstruction of digital signal in Matlab I'm trying to write a program in Matlab that samples (using Nyquist theorem) and recovers signal. However, I cannot write sampling part for sum of 2 signal.

f1 = 30;
T1 = 1/f1;
amplitude1 = 1;
f2 = 60;
T2 = 1/f2;
amplitude2 = 1;
signal1 = amplitude1 * sin(2*pi*t/T1);
signal2 = amplitude2 * sin(2*pi*t/T2);
signal = signal1 + signal2;
plot(t, signal);
grid on;

I write this code.I want to draw Undersampling, sampling at Nyquist
rate and oversampling. When make researching, I find code like that, but i cannot use it for my signal.

fm=input('Enter the Msg frequency(fm):\n');
xa1=amplitude*sin(2*pi*fm*t);
subplot(2,2,1);
plot(t,xa1);

What is fm in my signal ? How can i find it ?

fs1=2*fm + 10
n=0:1/fs1:1;
xa2=amplitude*sin(2*pi*n*fm);
subplot(2,2,2);
stem(n,xa2);
title('Above niquist rate');

I want to write a code like above for all 3 cases. And how can I find the Nyquist sampling rate for signal empirically ? Thank you all.

I know signals are not same. But I just want to do same thing on my signal.

You may be interested in the simple experiment using matlab.
https://poweidsplearningpath.blogspot.com/2020/04/ch4-adcdac-how-to-simulate-adcdac.html

Reconstruction is essentially a kind of interpolation or so called digital to analog conversion (DAC). Detail descriptions are introduced in chapter 4.8.3 of the DSP Bible 1.

However, we all understand the ideal theory. But the problem is how to bridge the theory and practice. Because all the signals we can simulate in Matlab is essentially 'Digital'. So, the so-cal ADC/DAC is just an approaching simulation. Let me briefly justify the difference.

Ideal Reconstruction in DAC: There is typically a zero order hold circuit and reconstruction filter. Please note that the frequency response of zero order hold is not evenly flat so a reconstruction filter comes to an aid.

Matlab Simulation in DAC: There is no need to actually simulate a zero-order hold. All you need are a up sampling and an ideal low pass filter.

Because the t cannot be found in the question. I just post a AD DA practice sample code here. Hope it can help.

This sample code sample a simulated 2 Hz analog signal,x_a, (with sampling rate 1500 Hz) to a discrete time signal,x_d, with 5 Hz sampling rate (nearly Nyquist rate) and then reconstruct this back,x_r.

%% This code simulate the AD/DA processing discussed in Chapter 4.8.3 [1]
close all; clear all;

%% parameters.
% analog
analog_fps = 1500;
analog_window_time = 3; %sec
t = 0:1/analog_fps: analog_window_time-1/analog_fps;

% digital
digital_fps = 5;
n = downsample(t,analog_fps/digital_fps);

% ADC: Quantizer
X_m = 1; % Range
B = 10;% Bit number.

%% Signal generation
freq_hz = 1; % Hz.
x_a_1 = 0.5*cos(2*pi*freq_hz*t+0.1);

% add a small high frequency component as asked.
signal_freq = 2; %Hz
x_a_2 = 0.5*cos(2*pi*signal_freq*t+pi/2);
x_a = x_a_1 + x_a_2;

%% ADC
% Sampling
x_s = downsample(x_a,analog_fps/digital_fps);

% Quantizing (abs of input value should not over 1)
% x_d = Quantizing(x_s,B,X_m);  % A For complete ADC, a quantizing should
% be added here.
x_d = x_s; % For your case, there is no quantizing here.

%% DAC
% up sample / DAC
x_up = upsample(x_d,analog_fps/digital_fps);

% LPF (Reconstruction Filters)
h = intfilt(analog_fps/digital_fps,4,0.9);
%% Important
% please not the parameter 0.9, ideally should be 1 for Nyquist rate.
% 0.9 here is ratio of Nyquist.
% Given known limit band signal, shourter ratio can enhance SNR by oversampling.
%  (i,e, here I filterout the freq larger than 2.5(Nyquist rate) * 0.9 = 2.25Hz)

x_r = filter(h,1,x_up);
x_r(1:floor(mean(grpdelay(h)))) = [];
x_r = [x_r zeros(1,floor(mean(grpdelay(h))))];

%% Display

figure;
plot(t,x_a);
hold on;
plot(n,x_d);
plot(t,x_r);
title('analog signal (1500Hz) v.s. digital signal (5Hz) v.s. Reconstructed signal (1500Hz)');
legend('x_a','x_d','x_r');

You can adjust the digital_fps with sampling rate higher than Nyquist rate (30 for example) and 0.9 in intfilt to 0.5 (shorter group delay).

A. Oppenheim and R. Schafer, Discrete Time Signal Processing 3rd. 2009
As others have said, the book 1 is called the DSP bible.


Poles and zeros in time domain of analogs filters I am currently studying two Butterworth and Chebyshev low-pass filters of order $n =3$ and $n=2$ respectively, whcih are in fact two prototypes to make a bandpass filter. The transfer function that I have is the following one:
$$H(j\omega) = 2\sqrt{\frac{R_G}{R_L}}\frac{1}{R_G(C+\frac{D}{R_L}) + A + \frac{B}{R_L}}$$
where $R_G$ is the impedance in series with the generator and $R_L$ is the load impedance

For the Butterworth filter the data is the following:
$$\begin{align}R_G&=R_L=50\;\Omega\\
A &=1 -0.708\omega^2\\
B &=j59\omega \\
C &=j(0.024\omega -0.008496\omega^3)\\
D &= 1-0.708\omega^2\end{align}$$

and for the Chebyshev filter the data is the following:
$$\begin{align}R_G&=50\Omega\quad R_L=18.8\;\Omega\\
A &=1\\
B &=j34.25\omega \\
C &=j0.036\omega\\
D &=1-1.233\omega^2\end{align}$$
and this is the poles diagram computed with the data above of such circuits:

where:
$$\begin{align}
\omega_{pB_1} &= j1.66\\
\omega_{pB_2} &= -1.45+j0.83\\
\omega_{pB_3} &= 1.45+j0.83\\
\omega_{pC_1} &= -0.9 + j0.55\\
\omega_{pC_2} &= 0.9 + j0.55
\end{align}$$

It is clear that both systems are going to be unstable since they have poles in the right-half plane, but how will these poles look in the time domain? I have read articles like this one or this one and some questions on this website, and I know that if they were complex conjugated poles the response in time domain would be a sinusoid type signal.

I am asking this because the majority of the examples have two conjugated poles, which in this case they are not, and hence the time domain expressions won't be sinusoids, or I am wrong? Could someone answer with the inverse Laplace transform in Matlab of my transfer function?

Note that I am not asking for my job to be done, I have already made it, but some more information about those poles will be appreciated.

Thank you.

I think you calculated the poles incorrectly, namely for $\omega$ instead of the Laplace variable $s=j\,\omega$. You can find the appropriate poles in terms of $s$ by multiplying your poles by $j$. This will give you conjugate poles, with negative real parts.


Estimating Velocity from Fuzzy Position Data with Known Uncertainty I was wondering if someone could maybe clarify or direct me to the best answer for this question.  I want to estimate velocities from position data.  However, the position data is fuzzy but I have an idea of its uncertainty.  My initial sense is that one is supposed to use some iteration of a Kalman filter, but from my limited understanding of it, you need to know at least one parameter with some certainty, i.e., the position data.  So to summarize, my questions are:

Is it possible to estimate velocity values given fuzzy position data?
If so, what are some available techniques for doing so?

Yes it is. That what's Kalman filter is all about. Your measurement includes position data, and there is a hidden velocity parameter. You will need some kind of model, like constant-velocity. You can input your idea of uncertainty into Kalman noise matrixes to improve results.
You can use Kalman filter, which is good to predict the future. If you just want to estimate the velocity in post-processing step (after you have all data), you can use something simpler, like estimating the derivative. This will have the advantage of being smoother and more accurate.


Definition of MIMO (Multi-input Multi-output ),does this output mean antenna or transmitter? I want to ask the definition of MIMO(Multi-input Multi-output )

If there are one transmitter and one receiver,each of them has $N$ antennas,$N >1$.
So i can say it is a MIMO system.

But if now i have one transmitter with $N$ antennas,$N >1$.And $5$ receivers with single antennas for each,can this still be called MIMO?

I think it should be called MISO,because each receiver just has only one antenna.However,in this paper,it called this system "MIMO" system .

Paper:https://arxiv.org/pdf/1903.05758.pdf

I want to ask the definition of MIMO(Multi-input Multi-output ) If there are one transmitter and one receiver,each of them has N antennas,N>1. So i can say it is a MIMO system.

Any communication system comprising of a transmitter with $N >1$ and receiver comprising of $M>1$ antennas is considered to be a MIMO system. So, this is not only limited to $N=M$ in the way you describe it.

  I think it should be called MISO,because each transmitter just has only one antenna.

No, I think it should be called MU-MISO (MultiUser MISO) because it's what it is; you're transmitting to multiple users (each with a single antenna).

  However,in this paper,it called this system "MIMO" system.

In some applications, MU-MISO and MU-SIMO are seen as MIMO.


UWB radar to find heartbeat I have a UWB radar that Im playing with. The problem is that my algorithms at the moment can clearly, reliably and validly pick up 1-2mm movements of a pendulum i've set up. However, when I use it on my chest to find my heart, things doesnt work. The radar should be set up so that the waves penetrate skin and bone, but I am not entirely sure. The antennae are specced specifically to be able to pick up a heartbeat (or penetrate human tissues).

In terms of math, I first do despiking and a bandpass filter around 0.1 - 5 hz, followed by a principle component analysis to reduce the number of samples/sensors. Then I either do a FFT plus some stuff to find the most pronounced frequencies of the top 5 principle components, or I search for the best fit sinuswave (as the heartbeat can be described as such afaik) to the top 5 principle components.

While these methods give similar results when used on a pendulum, they break down when trying to find the heartbeat. The radar can nicely penetrate wood, and plastic and other stuff. Any ideas?

Bonus question: We had some ultrasound gel at the lab, and i tried to see if removing the "air" component would help, but then all signal dropped. Is ultrasound gel absorbing the radar pulses?

Thanks alot!

I am also working on UWB bio radar to detect the heart beats. But the problem is that breath signal is also mixed up with heart signal. I am trying to remove the breath signal (0.25 to 0.7 Hz) and to get the heart frequency signal (1  to 4 Hz). Currently I am also using an ultrasound sensor to detect heart signal.


Effect of noise on high and low frequency components It is given that noise is present in the noisy image in the form of high frequency components. Can any one please explain this concept. Also please tell me that whether it is correct to say that "effect of white Gaussian noise is more on high frequency part of the image as compared to low frequency components" or not.

This is actually due to the characteristics of the actual signal, rather than noise.

You might have heard of pink noise and its spectral characteristics $$S(f)\propto \frac{1}{f^\alpha}$$
Now forget about the term "noise" in here, and imagine a stochastic process whose spectral density vanishes following a $\frac{1}{f^\alpha}$ law, in which $\alpha$ is a constant.

Interestingly, this is the law that nature follows a lot. More specifically, many natural signals follow this law. It has been observed and actually been verified to be valid for instance in human voice, natural images, music, quasar light emissions, heart beat, firings of neurons, and even in financial markets. Note that natural is a key term here and by that, we refer to attributes of a natural environment. For a more detailed discussion see here and here. For 2D (image) signals we should expect $\alpha=2$.

Since according to the $\frac{1}{f}$ law the contribution of the actual [natural] signal will diminish at higher frequencies, we would expect more contribution of noise at these frequencies. Consequently, noise is usually dominant at high frequencies and the signal is more dominant at lower frequencies.


How do I use MATLAB to obtain a continuous state equation from a transfer function and ode45 to solve it? I have the following transfer function that I am trying to model in MATLAB with some additional unity white noise:

  $$G(s)=\frac{1}{s^3 + \beta s^2}$$

I used the following to obtain the continuous state space equation:

beta = 0.2;
b = [0 0 0 1];
a = [1 beta 0 0];
[A,B,C,D] = tf2ss(b,a);

My understanding is that the tf2ss() function provides me with the values A and B which can be used in the following equation:

  dxdt = Ax + Bu, where x is the initial condition and u is the unity white noise.

Upon running the above code snippet, I obtain the following values for A and B:

A = [-0.2,0,0; 1,0,0; 0,1,0]
B = [1; 0; 0]

Here is where I am running into trouble (although I may have run into an issue before with my assumptions regarding tf2ss()).  I do not understand how to write the differential equation function required to use ode45().  I have tried the following:

function dx = mydiffeq(t,x)
    A = [-0.2,0,0;1,0,0;0,1,0];
    B = [1;0;0];

    dx = A*x + B*randn(1);

end

And then I call it as follows:

tspan = [0 10];
xnot = 100;

[t,y] = ode45(@mydiffeq, tspan, xnot);
plot(t,y);

Obviously, I have a problem due to the size of my A and B matrices.  I know that I am doing this wrong, but I do not know how to fix it.

You're doing pretty much everything right. The only problem is that your system is a third order system (s^3) and hence, the state x in

dx = A*x + B*u

is 3x1 but your initial condition xnot is a scalar. Replace

xnot = 100

with

xnot = [0 0 0].';

(or whatever initial condition you have) and it should work.


Why SRS signalling is used in 5g? Since DMRS is also being used as a reference signal in PUCCH and PUSCH channel, so why an extra SRS signalling is required.

Because DMRS is only sent if there is data to transmit and, therefore, narrowband while SRS, as the term sounding suggests, is not associated to any data transmission and also wideband.
For example, SRS can be used by base stations to choose the best portion of the UE channel bandwidth to schedule next PUSCH.


Loop for computing data I have two variables of data (BxHPF and ByHPF), which are equall in number of data. I would like to make a loop, which takes first value of BxHPF and ByHPF compute them, then take next one from BxHPF and ByHPF and so on. For now I have:

So I though about this way:

Computing = (float(BxHPF[0])/float(ByHPF[0]))
Dat1 = math.degrees(Computing)*(-1)

In fact those equations give me correct result. But as I wrote I need to make a loop, which will count every pair from BxHPF and ByHPF with using those two variables Computing and Dat1.

To be exact BxHPF and ByHPF contain 266150 records each.

result = [math.degrees(x/y)*(-1) for x,y in zip(BxHPF,ByHPF)]

It seems that I have found an answer. Maybe it will help someone.


How to know when a sound in a given frequency starts and ends? I need to detect two sounds (913 Hz and 1370 Hz), knowing the time in which each starts and ends.

At this time, I am using a band pass filter (butter), and this is the result:

Here is the first time this sound occurs:

The red line is where the sound is supposed to end (276 ms).

So how can I detect that the first time the tone occurs it starts at 0 ms and ends at 276 ms?

I would recommend reading about onset detection. It's about music signals but nevertheless a few ideas might be extracted from it.

From a more general point of view, what I would do is a frame-based analysis of the signal and define some kind of energy threshold/noise floor with which I could determine if there is activity or not.


Centering zero frequency for Discrete Fourier Transform I am working on a image processing application which uses a discrete fourier transform to implement blurring/sharpening. The application is more or less working, but something about the mechanics is still confusing to me.

In particular, it is how the process of centering the zero frequencies is being done.

The example I have seen preprocesses the input image (of greyscale intensities) by multiplying it with a matrix of size equal to the input image, whose values are $(-1)^{x+y}$, where $x$ is the row, $y$ is the column, so a pattern alternating $1$ and $-1$

According to the notes, this is equivalent of swapping the quadrants of matrix by flipping across the $x$ and $y$ axis.

I understand why this is done, and I would like to stress I understand I have my code/Fourier stuff working, I just don't understand why multiplying the input matrix by 1/-1 ends up centering the zero frequency component around 0.

Thanks

Oh!  What a cool trick!  It works because of the convolution theorem (i.e., multiplication in the spatial/time domain is equivalent to convolution in the frequency domain.)

It's not flipping across the $x$ and $y$ axis, it is rotating the image of the Fourier transform (think of shifting halfway around a cylinder).  The trick here is that alternating -1,1 in the spatial domain is a signal with the highest frequency.  So the Fourier transform of that image is a single point in the frequency domain.  Convolving by a single point is equivalent to shifting (rotating) the image by the offset of the point from the zero frequency.

Here is a test image:
.  It's Fourier transform looks like:

If you take the Fourier transform of the alternating image (), it results in a single point right at the center of the Fourier transform: .  (Recall we haven't done our rotation yet, so the center of the fourier transform is the high frequencies and the low frequencies are at the corners still.)  But this is the "rotation kernel!"  Convolving with this rotation kernel moves everything down and to the right (but things that fall off the bottom right rotate into the top left.)

Convolving the original image with the rotation kernel (in the image domain) gives you: , while convolving the fourier transform image with the rotation kernel (in the frequency domain) gives you: .

And we can check that multiplying the testimage by the checkerboard in the image domain gives , which has a fourier transform of: .


Derivative with respect to complex conjugate I have a real function $C$ of a complex vector $x$. While taking the gradient of the function $C$ for minimising the same, why do we take the derivatives with respect to the complex conjugate of $x$, i.e. $\bar{x}$ and not the actual vector $x$?

I have tried looking up for help but I am directed again and again to  Cauchy–Riemann equations, which I know is not pertinent here.

That's a trick which you will also find in a DSP context, that's why I choose to provide an answer here. It is related to the Wirtinger derivative, and you can find more details about it in this answer over at math.SE.

In practice this trick is often used to compute the extremum (minimum or maximum) of a real-valued function depending on a complex variable (or vector/matrix). In order to find the extremum, you formally take the derivative with respect to the complex conjugate of the variable of interest, set this derivative equal to zero, and from this equation derive the optimum value of the (possibly vector-/matrix-valued) variable.

As a simple example, take the minimization (with respect to the vector $x$) of the following error quantity:

$$\epsilon=||Ax-b||_2^2\tag{1}$$

where the vectors $x$ and $b$, and the matrix $A$ are assumed to be complex-valued. The error can be written as

$$\epsilon=x^HA^HAx-x^HA^Hb-b^HAx+b^Hb\tag{2}$$

where $^H$ denotes the Hermitian (conjugate) transpose. Now you take the derivative of $(2)$ with respect to $x^H$, formally treating this vector as being independent of $x$. The result is

$$\frac{\partial\epsilon}{\partial x^H}=A^HAx-A^Hb\tag{3}$$

Setting this derivative equal to zero and assuming that $A^HA$ is invertible, the optimal value of $x$ is given by

$$x_o=(A^HA)^{-1}A^Hb\tag{4}$$

where the matrix $(A^HA)^{-1}A^H$ on the right hand side of $(4)$ is the Moore-Penrose inverse.


sampling rate Matlab I have a rather strange but i thin interesting question. The idea is to gain a better understanding of sampling rate vs frequency shifting when playing audio. The idea is a little experiment:

%synthetic example

T=1/(10*10^6);%period should be at least 1/(2*10^6Hz) => Nyquist freq if we want to be able to reproduce 10^5Hz max freq
x=0:T:1-T;

f=10^6;%frequency
y=sin(2*pi*f*x);

%i see visually that there is 11 samples constituting 1 period
%plot(y(1:11))
%plot(y(100:111))
%etc
nbPeriods=length(y)/11;%nbtotalsamples/nbsamplesOf1Period
%y contains 10^6 oscillations each of 11 samples
%therefore if i want to reproduce a 1Khz sound, I compute my sampling frequency :
% Fs = nbPeriods/10^3

Fs=909.09;

a=audioplayer(y,Fs)

tic;
play(a)
toc;

my question is: is this reasoning correct? -> that Fs should be 909.09 to get a 1kHz sine wave from  f=10^6;y=sin(2*pifx); when played at Fs=909.09 and considering of course the x range i have taken?

EDIT:
following your answer would this (the comments!) be correct now:

clc;clear all;
%synthetic example

T=1/(10*10^6);%
x=0:T:1-T;%my sampling rate is 1/T= 10*10^6 =10MHz

f=10^6;%frequency
%y: as if i had 1 period T of a sampled signal (sampled at 10MHz)
y=sin(2*pi*f*x);%this signal has a freq of 1 MHz , according to
%Nyquist Shannon sampling thm the MAX freq i could capture with
%my sampling rate of 10MHz would be 5 MHz,
%so i am a bit below that with my 1MHz

%given that i want to create a 1KHz wave from y i must sample at
% (10*10^6 samples)/(10^3 samples)  = 10^4 10 KHz
%idea: cutting my 10M samples into chunks of 10^3 oscillations
Fs=10^4;

a=audioplayer(2*y,Fs/4)

tic;
play(a)
toc;

%indeed it plays now a sound that "sounds" like a reference 1KHz wave,

You reasoning is NOT correct. One period is 10 samples long, not 11. The 11th sample is the first sample of the 2nd period, not the last sample of the 1st period.

So originally you have

sample rate: 10 MHz
sine wave frequency: 1 MHz, which is 1/10th of the sample rate

If you want to get a 1 kHz sine wave you need to set the sample rate to 10 times 1 kHz, i.e. 10 kHz. If you play it back at a sample rate of 909 Hz, your resulting frequency would be 90.9 Hz.


Partial Fraction Expansion for Inverse Fourier Transform In many textbooks, I've seen the application of Partial Fraction Expansion (PFE) to find an inverse Fourier Transform. Let's stick to the discrete time case, and let me give you an example.

Let's say that $$X(e^{j\omega}) = \frac{-\frac{5}{6}e^{-j\omega} + 5}{1 + \frac{1}{6}e^{-j\omega} - \frac{1}{6}e^{-j2\omega}}$$
and we want to find the signal $x[n]$ that has the above DTFT. We can apply PFE and get $$X(e^{j\omega}) = \frac{A}{1 + \frac{1}{2}e^{-j\omega}} + \frac{B}{1 - \frac{1}{3}e^{-j\omega}} \tag{1}$$
where $B$ (and similarly, $A$) is given by
$$B = X(e^{j\omega})\Big(1 - \frac{1}{3}e^{-j\omega}\Big)\Big|_{e^{-j\omega}=3}$$

I know now that we can set $u = e^{-j\omega}$ and turn this into

$$B = X(u)\Big(1 - \frac{1}{3}u\Big)\Big|_{u=3}$$

However, setting $u = e^{-j\omega} = 3$ does not look good to me, because this equation has no real solution for $\omega$.

What do I miss in this application of PFE?

You should indeed replace $e^{j\omega}$ by the more general complex variable $z$ because when talking about poles and zeros of a rational function, you must consider the whole complex plane. The fact that $e^{-j\omega}=3$ doesn't have a real solution just means that the corresponding pole does not lie on the unit circle.

After having performed the partial fraction expansion, there is nothing that keeps you from setting $z=e^{j\omega}$.


Raised-Cosine Filter - phase response graphs for three kinds of filter In the previous post it has been mentioned that the phase response of Raised-Cosine filter is linear and is given with ϕ(ω)=−ωT, where T is the delay. But, I have also read in literature that it is also a minimum phase filter and non-linear. Which one is correct, and if it is minimum phase how I can calculate it?

I have three kinds of Raised-Cosine filters, as per picture below. What is the phase response for each kind of filter? How the phase response graph will look like for each filter?

You need to realize that the term "raised cosine" says absolutely nothing about the related phase response. I think the context of this question of yours is audio equalization, and in that context filters with certain desired magnitude response shapes are sometimes described by piecewise constant values which are smoothly connected by raised cosine shapes.

So those shapes in your figures are idealized magnitude responses. Now you can use any filter design technique to approximate such a response. Still nothing is said about the phase, and it depends on the design technique which phase response you'll end up with. Analog filter designs will usually result in a minimum phase response. But note that this has absolutely nothing to do with the raised cosine shape of the desired magnitude response.

In digital communications raised cosine pulses are used because they satisfy the Nyquist criterion. There also the phase is important and the frequency response of the ideal pulse has a raised cosine magnitude shape and a linear phase.

In sum, it is important to understand that the "raised cosine" property is a description of the magnitude response, and the application and/or the filter design process will determine the phase response, which can be linear, minimum phase, or anything else.


Rayleigh fading in simulink I have a simulink model of zigbee over a Rayleigh fading channel.

The modulation in the transmitter part is O-QPSK modulated baseband.

Most of the examples didn't add the RF frequency after the baseband modulator and I can't understand why.

And what is the effect of cancelling this part in simulink?

However, when I try to simulate the RF band after baseband modulator using up converter, the output of it was real and the Rayleigh multipath fading block didn't accept it.

You're forgetting the most basic truth of baseband signal processing: Any causal, stable, linear channel $h$ can be just as well represented as baseband system as it can be represented as RF passband system. That is a direct result of both the Fourier transform and the Convolution being a linear operation.

That's the reason why we can do any sensible Software Defined Radio, and baseband signal processing, in general.


Recommended Resources / Literature Search Terms for a Solutions to a Specific Kind of Multi Harmonic Signal Structure Hopefully this isn't considered too off-topic. I'm working in industry these days and came up with a solution to a signal processing problem we'd been facing. I'd like to get a sense as to whether said solution has been published before or if there are alternatives I should look at, but I'm having trouble searching for this particular signal scenario. I'll describe the signal structure below, and would appreciate any input on how I should be searching the literature.
So I have two signals: one that is pure noise (noise_only), and one that is very similar noise plus a target signal of interest (noise_plus_target). In each signal, the noise is actually made up of multiple somewhat-frequency-distinct signals, and the same set of said noise signals contribute to each observed signals additively but with different weights between the two. In theory, we should be able to use information from the noise_only signal to help remove the noise from the noise_plus_target signal, and I've come up with a method to achieve this, but I want to know what other solutions might have been already published for this kind of scenario.
In R, code for generating fake signals that match the characteristics of my real signals would be:
library(tidyverse)

# define a function to generate simple sinusoid given time and hz
sine = function(time,hz) sin(time*(2*pi)*hz)

#define a function to scale values to 0:1
scale01 = function(x) (x - min(x)) / diff(range(x))

#specify sample rate
sample_rate = 10 #in Hz
max_time = 30

#construct a tibble
latent_signals = tibble(
    #specify sampling times (in seconds)
      time = seq(0,max_time,1/sample_rate) #30s of data
    #construct some latent noise signals, each at a decently separated Hz
    , noise1 = sine(time,1/11)
    , noise2 = sine(time,1/3)
    , noise3 = sine(time,1)
    #specify a target signal that will be hidden in the noise
    # This could take any shape; here I've chosen a bump midway
    # through the timeseries
    , target = scale01(dnorm(time,mean=max_time/2,sd=3))
)

#show the latent signals
latent_signals %>%
    tidyr::pivot_longer(
        cols = -time
    ) %>%
    ggplot()+
    facet_grid(
        name ~ .
    )+
    geom_line(
        mapping = aes(
            x = time
            , y = value
        )
    )

#combine the latent signals into two observed signals, with different weights
# for each and the latent target only in one
latent_signals %>%
    dplyr::mutate(
        noise_only =
            noise1*runif(1,.5,1.5) +
            noise2*runif(1,.5,1.5) +
            noise3*runif(1,.5,1.5)
        , noise_plus_target =
            noise1*runif(1,.5,1.5) +
            noise2*runif(1,.5,1.5) +
            noise3*runif(1,.5,1.5) +
            target
    ) %>%
    dplyr::select(
        time
        , contains('_')
    ) ->
    observed_signals

#show the observed signals
observed_signals %>%
    tidyr::pivot_longer(
        cols = -time
    ) %>%
    ggplot()+
    facet_grid(
        name ~ .
    )+
    geom_line(
        mapping = aes(
            x = time
            , y = value
        )
    )

```

If I understand this problem correctly you have access to 2 signals:

Noise Signal - $ w \left[ n \right] $. It is composed of a linear combination of harmonic signals. Something like $ w \left[ n \right] = \sum_{i}^{m} {a}_{i} \sin \left[ 2 \pi \frac{ {f}_{i} }{ {f}_{s} } n + {\phi}_{i} \right] $.
Input Signal - $ y \left[ n \right] $ which is composed of the signal of interest $ x \left[ n \right] $ and the same harmonics as in $ w \left[ n \right] $ just with different weights. Something like $ y \left[ n \right] = x \left[ n \right] + \sum_{i}^{m} {b}_{i} \sin \left[ 2 \pi \frac{ {f}_{i} }{ {f}_{s} } n + {\psi}_{i} \right] $.

The objective is to estimate $ x \left[ n \right] $ given those signals.
Since we can model $ y \left[ n \right] = x \left[ n \right] + \left( h \ast w \right) \left[ n \right] $ we can utilize the adaptive filter framework and solve this in highly efficient manner.
Remark
In real world each observed signal would also have some kind of measurement noise we can model as AWGN. But let's leave it a side at the moment.


How might I go about detecting phase changes in a sine wave of a certain frequency? I want to detect occasional instantaneous changes in the phase of a sine wave.

I have a signal that is a sine wave whose frequency is in the kHz range and experiences random changes in its phase. I want to detect these phase changes and count the occurrences. Ideally, I would also be able to extract some approximate phase change in degrees for each occurrence. How might I go about that?

I implemented Jason R's suggestion in Python:

# generating a 30KHz sine wave sampled at 500KHz with a pi/2 phase shift halfway through the signal
F = 30e3
Fs = 500e3
t = np.arange(0, 1e-3, 1/Fs)
S = np.append(np.sin(F*2*np.pi*t[:int(len(t)/2)] + 0), np.sin(F*2*np.pi*t[int(len(t)/2):] - np.pi/2))

S_analytic =  hilbert(S) # calculating the analytic signal

def GetRealImagArray(Array):
    ImagArray = np.array([num.imag for num in Array])
    RealArray = np.array([num.real for num in Array])
    return RealArray, ImagArray

R, I = GetRealImagArray(S_analytic)
phi = np.arctan2(I, R)
f = phi[1:] - phi[:-1]

I then get the following by plotting the signal and FM discriminator:

The large positive spike indicates the change in phase as required in my question. What is the cause of the downward spikes at the minima of the sine wave? Is this expected behaviour or an error in my Python implementation?

The easiest way to do this would be to apply a digital FM discriminator to the signal. At a high level, you would do the following:

Form the analytic signal of your input. This also can be termed as ensuring it is at complex baseband.
Calculate the phase versus time of your analytic signal $x[n]$:
$$
\phi[n] = \text{atan2}\left(\frac{\text{Im}\{x[n]\}}{\text{Re}\{x[n]\}}\right)
$$
The instantaneous frequency of the signal versus time can be estimated by differentiating $\phi[n]$ using any discrete differentiation method (first-order difference is usually used, i.e. $f[n] = \phi[n]-\phi[n-1]$).

Why does this work? If your signal is typically a constant frequency, your FM discriminator's output $f[n]$ should be approximately constant. If there is a phase jump at a particular sample index, this would appear to the discriminator to be a sudden change to a much different frequency. After the jump, the discriminator output $f[n]$ will again track the sinusoid's frequency.

Therefore, you can run a discriminator continuously on your input signal and apply a threshold, looking for large changes in frequency from sample to sample. That would correspond to a phase discontinuity.

Edit: The above is subject to spurious outputs when the phase value output by the arctangent function wraps around $\pm \pi$ (see the comment below). A more robust way of implementing the above would instead be:

$$
\hat{x}[n] = x[n] x^*[n-1] \\
f[n] = \text{atan2}\left(\frac{\text{Im}\{
\hat{x}[n]\}}{\text{Re}\{
\hat{x}[n]\}}\right)
$$

That is, form the product of $x[n]$ with the conjugate of $x[n-1]$. At each sample instant $n$, the result will have a phase that is equal to the difference in phase between $x[n]$ and $x[n-1]$ (which is what you're trying to measure). The subtraction alleviates any wrapping issues when the absolute phase of $x[n]$ is near $\pm \pi$. Just find the phase of the resulting product to get $f[n]$.


Standard deviation of the spectrum of white noise I have a temporal signal which looks like $f(t) = t\eta(t)$, where $\eta(t)$ is a white noise with the mean $\eta_m$ and STD $\sigma$. I want to calculate the corresponding spectrum, $F(\omega)$, take its magnitude, $|F(\omega)|$, and finally compute the STD of $|F(\omega)|$ in terms of the given parameters. I have searched some literatures but I couldn't find anything which helps me directly tackle the problem. So, anyone know how to do this, or know which resources you would refer me to?

The following pictures display $f(t)$ and $|F(\omega)|$.

The one below is the histogram of $|F(\omega)|$ excluding some values around the central peak.

I want to derive a mathematical analysis which can compute the standard deviation of the above histogram.

true white noise has infinite power (because the power spectrum with height $ \frac{\eta}{2}$  has infinite area) which means infinite variance $\sigma^2$ which means infinite standard deviation.  white noise is a conceptual instrument.

to do any analysis of a system with white noise in it, you must somehow determine the (single-sided) bandwidth of the system, $B$, and then the power or variance is $\eta B$.


Is calculus background a must for learn DSP? I'm beginner in DSP and in the literature some Authors suggest Analog (Continuous) Signal Processing a prerequisite for learn DSP. But it demands calculus background, that I do not have. Is calculus and ASP background a must for learn DSP?

With Google's search engine and Google books, I investigated that someone, without calculus background, can learn DSP. I expected learn from you that for avoid laborious searching. I'm very sorry for did not learn that from you. Regards.


Upsampling signal for cross correlation I am recording a 17khz signal at a sample freq of 44.1khz. I want to perform cross-correlation between the received and transmitted signals for calculating TDOA. But when I do xcorr, the results are inconsistent. Is it because my sample freq is too low for my signal freq? If so, how can I overcome this? Will upsampling help?

Edit:
My received signal also consists of noise. I am passing my received signal through a bandpass filter with limits 16khz to 18khz. I am using butterworth filter. Is it better to use FIR filter instead?

A narrow-band filter may remove some of the bandwidth of the signal containing timing information of the signal's envelope modulation.  The narrower the filter, the more any tone burst will be smeared out in time.  Upsampling after the timing information is destroyed won't help much.  Using a wider bandwidth signal, channel, and filter, with perhaps a higher sampling rate to support that bandwidth, will help improve time estimation accuracy.

A chirp, or other large frequency modulation is one common method used to increase a signal's bandwidth to allow better timing estimation.

Also note that any real signal with a frequency near either DC or Fs/2 can face strong interference from it's negative frequency conjugate image in the frequency domain, as well as any of the images sidebands, which can also greatly distort the signal envelope and the envelopes timing information in the received sampled data.  A higher sampling rate can also help solve this problem for high frequency tones.  (And a longer sample window length for near DC signals.)


Cross correlation with FFT and fftshift When I calculate the cross correlation between two vectors with the following code

a = [1 2 3 4 5];
b = [6 9 8 9 10];

% correlation using convolution
conv(a,conj(flipdim(b,2)))
% correlation using FFT
ifft(fft(a,9).*conj(fft(b,9)))
fftshift(ans)

when I use the FFT i need to fftshift the output in order to get the correct cross correlation:

ans = 10    29    56    92   134   116    94    69    30

ans = 134.0000  116.0000   94.0000   69.0000   30.0000   10.0000   29.0000   56.0000   92.0000

ans = 10.0000   29.0000   56.0000   92.0000  134.0000  116.0000   94.0000   69.0000   30.0000

However, if i just calculate the convolution of the two vectors using the FFT I don't need to apply the ffshift

 conv(a,b)
 ifft(fft(a,9).*fft(b,9))

ans = 6    21    44    76   118   124   106    85    50

ans = 6.0000   21.0000   44.0000   76.0000  118.0000  124.0000  106.0000   85.0000   50.0000

Why do I need to do a fftshift when calculating the cross correlation? I guess it must be very obvious but I just don't see it right now and I've been googling for hours.

The (linear or aperiodic) convolution of two vectors $\mathbf x = (x[0], x[1], \ldots, x[N-1])$ and $\mathbf y = (y[0], y[1], \ldots, y[N-1])$ is a vector $$\mathbf z = {\mathbf x}\star {\mathbf y} = (z[0], z[1], \ldots, z[2N-1]).$$
On the other hand, their cross-correlation is a vector
$$\mathbf w = {\mathbf x}\otimes {\mathbf y} = (w[-(N-1)], w[-(N-2)], \ldots,
w[-1], w[0], w[1], \ldots, w[N-2], w[N-1]).$$
Now, the FFT method for computing convolutions computes the cyclic (a.k.a. periodic
or circular) convolution of the two vectors.  To compute the linear or
aperiodic convolution via FFTs, you zero-pad the
vectors to length $2N-1$ so that the cyclic convolution of the zero-padded
vectors is exactly the linear convolution you are looking for.
Thus, after FFTing, term-by-term multiplying, and inverse-FFTing,
what you get is $\mathbf z$. This is exactly what you have done as per
your question. However, when you use the FFT method to compute the
cross-correlation, what you get is the cyclic cross-correlation of
the two zero-padded vectors, and this cyclic cross-correlation is not
$\mathbf w$ but instead
$$\hat{\mathbf w} = (w[0], w[1], \ldots, w[N-1], w[-(N-1)], w[-(N-2)], \ldots,
w[-2], w[-1])$$ so that to get the $\mathbf w$ that you want, you need
to circularly rotate $\hat{\mathbf w}$ so that $w[0]$ appear in the middle
of the vector.


Strange results with all pass compensation I am new to DSP and, as a learning exercise, I am trying to design an all-pass filter in order to compensate the phase of an input impulse response (heqm.mat, which you can download from here).
The MATLAB code is the following; for the phase compensation I am using the time-reversed technique, also called matched-filter:
load heqm.mat; % Load the heqm vector
HEQM      = fft(heqm);
% Design the all-pass filter
HEQMP     = exp(1j.*angle(HEQM));
tmp       = ifft(HEQMP,'symmetric');
% Time-reversed technique
eqpr      = tmp(end:-1:1);
% Figure on the left
figure; plot(unwrap(angle(fft(heqm))));
hold on; plot(unwrap(angle(fft(eqpr))));
% Convolution of the original impulse response with the all-pass filter
o = conv(heqm,eqpr);
% Figure on the right
figure; plot(unwrap(angle(fft(o))));

As you can see from the left picture, the phase of the all-pass filter is the opposite of the original one. But if I convolve the two vectors, why the resulting phase (right figure) isn't flat as it should be? What am I missing here?

UPDATE: As Hilmar pointed out, OP did convolution in the time domain (as is clear in this code). This answer is misleading so will be deleted).
You would not convolve the two responses as you are in the frequency domain. It is the time domain impulse response of the equalization filter that you would convolve with, and convolution in the time domain is the product in the frequency domain. For the product of complex signals, the phases would add, so in your simple example (which may not be realizable due to causality) the phases would add to zero.
For the product of complex signals, the phases would add:
$$K_1e^{j\phi_1}K_2e^{j\phi_2} = K_1K_2e^{j(\phi_1 + \phi_2)}$$


For complex values, why use complex conjugate in convolution? Taken from Adaptive Filter Theory (2014) written by Haykin page 110 :

$$y(n) = \sum_{k=0}^{\infty} w_k^*u(n-k), \quad n=0,1,2,...$$

where $u$ and $w$ are complex values. My question is why use complex conjugate of $w_k$? The answer found in the book says "..., in complex terminology, the term $w_k^*u(n-k)$ represents the scalar version of an inner product of the filter coefficient $w_k$ and the filter input $u(n-k)$ ". I still don't understand, can you elaborate more on this answer?

Turns out that convolution and correlation are closely related.  For real signals (and finite energy signals):

Convolution: $\qquad y[n] \triangleq h[n]*x[n] = \sum\limits_{m=-\infty}^{\infty} h[n-m] \, x[m]$

Correlation: $\qquad R_{yx}[n] \triangleq \sum\limits_{m=-\infty}^{\infty} y[n+m] \, x[m] = y[-n]*x[m]$

Now, in metric spaces, we like to use this notation:

$$ R_{xy}[n] \triangleq \Big\langle x[m], y[n+m] \Big\rangle = \sum\limits_{m=-\infty}^{\infty}x[m] y[n+m] $$

The $\langle \mathbf{x},\mathbf{y} \rangle$ is the Inner product of the vectors $\mathbf{x}$ and $\mathbf{y}$ where $\mathbf{x} =\{x[n]\}$ and $\mathbf{y} =\{y[n]\}$.  Then we also like to define the norm of a vector as

$$\begin{align}
 \| \mathbf{x} \| &\triangleq \sqrt{\big\langle \mathbf{x},\mathbf{x} \big\rangle} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}x[m] x[m]} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}x^2[m]} \\
\end{align}$$

and that looks a lot like the Euclidian length of a vector with an infinite number of dimensions.  All this works very well for the case where the elements $x[n]$ of the vector $\mathbf{x}$ are all real.  The norm $\| \mathbf{x} \|$ is always real and non-negative.

So, if we generalize and allow the elements of $\mathbf{x}$ to be complex-valued, then if the same definition of norm is to be used,

$$\| \mathbf{x} \| \triangleq \sqrt{\big\langle \mathbf{x},\mathbf{x} \big\rangle}$$

then the definition of the inner product needs to be modified a little:

$$\Big\langle \mathbf{x},\mathbf{y} \Big\rangle = \sum\limits_{m=-\infty}^{\infty}x[m] y^*[m]$$

Then if $\mathbf{x}$ has complex-valued elements, the norm comes out as:

$$\begin{align}
 \| \mathbf{x} \| &\triangleq \sqrt{\big\langle \mathbf{x},\mathbf{x} \big\rangle} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}x[m] x^*[m]} \\
 &= \sqrt{\sum\limits_{m=-\infty}^{\infty}\Big|x[m]\Big|^2} \\
\end{align}$$

So, evidently, Haykin is just walking that definition of inner product back to the definition of convolution.


Interpret numpy.fft.fft2 output My goal is to obtain a plot with the spatial frequencies of an image - kind of like doing a fourier transformation on it. I don't care about the position on the image of features with the frequency f (for instance); I'd just like to have a graphic which tells me how much of every frequency I have (the amplitude for a frequency band could be represented by the sum of contrasts with that frequency).

I am trying to do this via the numpy.fft.fft2 function.

Here is a link to a minimal example portraying my use case.

As it turns out I only get distinctly larger values for frequencies[:30,:30], and of these the absolute highest value is frequencies[0,0]. How can I interpret this?

What exactly does the amplitude of each value stand for?
What does it mean that my highest value is in frequency[0,0] What is a 0 Hz frequency?
Can I bin the values somehow so that my frequency spectrum is orientation agnostic?

You actually did a two dimensional Fourier transform there. The 0 Hz frequency is the ,,DC offset'' of your signal, that is, its average value. In computers, it's no wonder that there's a high average value, because signals are often defined using positive numbers.

In a 2D Fourier transform, the coefficients of the transformed correspond to products of 1D Fourier basis elements, that is, signals with a period equal to the length of the signal. You could call these overtones.

2D Products of two such signals are not ,,orientation agnostic''. In fact, they have obvious orientation along the edges of the signal. It is not entirely clear what you are asking for here.


how can I make a simple lowpass filter with this constraints? I'm searching for a simple low pass filter code in C with this constraints

a modulable frequency cutoff (how can I convert all these crazy coefficients into a frequency?) I could change with a slider.
a dB/Octave I can change in my code.
a resonance I can change in my code.

the graphical result will be a simple slider moving from 0-44Khz.
So far I've only managed to cut a bit of the high ends of my sounds by doing

int i = 0;
while (i < x.size - 1)
{
    y[i] = (x[i] + x[i+1]) / 2;
    i++;
}

where x is the input buffer (depends of the settings but it goes from 32 to 2048 samples) and y is the output buffer.

check out the music-dsp source code archive here: http://www.musicdsp.org/showmany.php

I recommend a "Chamberlin" State-Variable Filter (which is 12dB/oct) which gives you other flavors (highpass, bandpass, etc) for free. If you cascade two of them you can get 24dB/oct, etc.


For SNR calculation, should I use FFT|X+iY| or FFT|R|? I use a demodulator that outputs X and Y quadratures, with the R and Phase which can be calculated.

For SNR measurements, when plotting the FFT to find the SNR a signal at a certain frequency, should I be plotting the FFT of |X+iY| or the FFT of R?

My understanding would be to do the FFT of R so that the phase information can be completely canceled out, and R being the amplitude (the signal strength) is what matters to the SNR calculations.

Thank you

To calculate the SNR you need to get the signal's power. One way to do this is in the time domain, but Parseval's theorem says that the power in the time domain is equal to the power in the frequency domain as it should be..it is the same signal after all. The DFT version of the theorem is $\sum_{n=1}^N|x[n]|^2=\frac{1}{N}\sum_{k=1}^N|X[k]|^2$, where $X[k]$ is the DFT of $x[n]$, which means that the power will be determined by its magnitude squared and the phase has no part in the calculation.


Cross Correlation of output with input is equivalent to impulse response of a large bandwidth input (LTI System)? What is the practical way to found the impulse response of a system ? As impulse function has a very large amplitude at zero time (analog signal).

So the cross corelation of output of a large bandwidth input signal (with relatively small amplitude) gives the impulse response of a LTI system .How is these two equivalent?

There are many different ways. Here is one

Excite the system with a periodic broad band noise signal that has a period of larger than the length of the impulse response. Call it $N$ samples.
Capture the output of the system after it's stable and periodic (discard the first $N$ samples). Make sure that excitation and acquisition are using the same clock.
Grab any $N$ samples of your excitation and acquisition
FFT them and divide the acquisition spectrum by the excitation spectrum
Inverse FFT the result: you have an Impulse Response of the systems

Bonus points for:

Using coherent averaging to improve signal to noise ratio
Shape the spectrum of your excitation to match the noise floor in the frequency range of interest so get flat SNR.


Conditions for which the Hilbert transform returns a correct phase I'm quite new to signal analysis, and I'm currently trying to understand under which conditions a Hilbert transform can be used to compute the correct instantaneous phase and enveloppe of a given signal.
Say I start from the example in Python given here (from the scipy website):
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import hilbert, chirp
duration = 1.0
fs = 400.0
samples = int(fs*duration)
t = np.arange(samples) / fs
signal = chirp(t, 20.0, t[-1], 100.0)
signal *= (1.0 + 0.5 * np.sin(2.0*np.pi*3.0*t) )

analytic_signal = hilbert(signal)
amplitude_envelope = np.abs(analytic_signal)
instantaneous_phase = np.unwrap(np.angle(analytic_signal))
instantaneous_frequency = (np.diff(instantaneous_phase) /
                            (2.0*np.pi) * fs)
fig = plt.figure()
ax0 = fig.add_subplot(211)
ax0.plot(t, signal, label='signal')
ax0.plot(t, amplitude_envelope, label='envelope')
ax0.set_xlabel("time in seconds")
ax0.legend()
ax1 = fig.add_subplot(212)
ax1.plot(t[1:], instantaneous_frequency)
ax1.set_xlabel("time in seconds")
ax1.set_ylim(0.0, 120.0)

I get a nice result showing me the enveloppe and desired frequency:

However, if for instance I decide to add a little bit of noise:
signal+=random.rand(len(signal))*0.2

Then the inferred enveloppe and frequencies get not that good (some of the noise goes into the enveloppe, some into the frequencies):

This was quite expectable, and I know I just have to find a way to smooth my signal to solve that. My real problem is, say you add a trend to your signal:
signal+=np.linspace(-1,1,400)

Then the results get completely messed up: *can't add any more links since I'm new to the forum*
And the results are also messed up even for just a vertical shift... So, basically, my question is: what are the conditions needed for the Hilbert transform to return the correct phase and amplitude ?
For now, this is what I suppose:

the signal must not be noisy
the signal must be centered around zero
the signal must not have any trend
amplitude and frequency can vary

Am I right? Thank you.

A single instantaneous phase estimate may or may not make any sense if there is more than one frequency peak in the signal's local spectrum.  So, to get a better single frequency and phase estimate, you may want to first high-pass filter your signal to remove the spectrum of any D.C. bias or slow trend, and also band-pass filter the signal to remove any outside spectral peaks due to noise or interfering signals.

So my off-hand guess for a condition is that the local spectrum has a single clear frequency peak with a monotonicly decreasing spectral magnitude away from that one peak.  All the way to zero at DC for real signals.