Difference in following random proccess Let just say I understand the first process that is white noise,

  Process 1: If $x(n)$ is a Gaussian random variable and a process
  is formed from a sequence of $x(n)$ and all random variables are
  uncorrelated then this process is white noise.

  Process 2: Now if $x(n) = \alpha $ where now $\alpha$ is Gaussian
  random variable then each process in ensamble is equal to constant for
  all $n$

But I am not able to get what difference does it makes when we form Process 2 it seems same to me. Can any one explain the difference? Also how the ensamble became constant in the second case?

In the first case, $x(0)$ and $x(1)$ and $x(2)$ etc are all uncorrelated. (And since they're gaussian, they're independent). This means that, to get a realization of the first random process, you have to take $\infty$ independent samples from a gaussian distribution, and assign each of these $\infty$ samples to a $x(n)$.

The second process is much simpler. You just take one sample from a gaussian distribution (and call it $\alpha$), and then assign all of the $x(n)$ to this $\alpha$. In this way, $\alpha$ may be random, but you will get $x(0)=x(1)=x(2)=\cdots=\alpha$.

I'll post some plots, as you asked. I used the following code in MATLAB:

N = 3;
T = 10;
t = 0:T;
x1 = randn(N,T+1);
x2 = repmat(randn(N,1), 1, T+1);

figure
subplot(3,1,1); stem(t, x1(1,:));
subplot(3,1,2); stem(t, x1(2,:));
subplot(3,1,3); stem(t, x1(3,:));

figure
subplot(3,1,1); stem(t, x2(1,:));
subplot(3,1,2); stem(t, x2(2,:));
subplot(3,1,3); stem(t, x2(3,:));

In this first figure, I plot three independent realizations of Random Process 1:

And in this figure, three independent realizations of Random Process 2. Note that in this case, all the $x$ signals are constant, but the actual level changes from one realization to another:


Is it correct to call a Savitsky-Golay filter of degree 0 a simple moving average? I had this question after seeing that polynomial regressions fit polynomial functions of different degrees to a time-series, since the mean of a time series is a constant and that a constant is also a polynomial of degree 0 i wondered if the SG filter who also fit polynomial functions to a signal over a period window can be called a moving average.

So is it correct to say that a Savitsky-Golay Filter that try to fit a polynomial of order 0 to a signal is considered a moving average filter ?

Yes you can consider a zeroth (or first) order SG filter as a moving average filter. Below MATLAB / Octave code computes the impulse response of a SG filter of order $N$ and length $2M+1$ :
% Savitzky-Golay Filter
%
clc; clear all; close all;

N = 0;                      % a0,a1,...,aN : Nth order polynomial
M = 3;                      % x[-M],...,x[M] : 2M + 1 data

A = zeros(2*M+1,N+1);
for n = -M:M
    A(n+M+1,:) = n.^[0:N];
end

H = (A'*A)^(-1)* A';        % LSE fit matrix

h = H(1,:);                 % S-G filter impulse response (non-causal symmetric FIR)

figure,subplot(2,1,1)
stem([-M:M],h);
title(['Impulse response h[n] of Savitzky-Golay filter of order N = ' num2str(N), ' and window size 2M+1 =  ' , num2str(2*M+1)]);

subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('Frequency response magnitude of h[n]');

figure,plot(linspace(-1,1,1024), 20*log10(abs(fftshift(fft(h,1024)))));
title('Frequency response magnitude of h[n]');

The impulse response for $N=0$ is :

as can be seen, it's a moving average (constant) impulse response.


Bandpass overshoots at the beginning and end of signal: why? The picture below shows the original data in the upper figure and the filtered data in the figure below. The filter applied to the original data is a simple bandpass implementation:

$H(z)  = \frac{0.2 + 0.1z^{-1} - 0.1z^{-3}  - 0.2z^{-4}}{1-0.98z^{-1}}$  (Relative Spectral Processing Filter)

As you can see, at the beginning of the data and at the end the filtered signal overshoots / undershoots sharply.

So far I experimented with a) adding random noise to the original signal and b) removing leading / trailing zeros but those are just my own ideas without any reference for a thesis. So my question is:

How would this problem be handled in digital signal processing?

+++++++++++

If you're interested: Removing leading / trailing zeros reduced the problem (original overshoot is ~40, now it's around 12 which is still the highest value in the filtered signal) and adding random noise reduced the problem even further but the signal is noisy now (which  is a unfavorable tradeoff, isn't it?)

Well what you see is the transient response of your filter.  Your signal looks a bit a like a step signal at the beginning. At the end, you have a negative step, thus you will get another transient that will be the opposite of your first transient.

Plot the step response of your filter and compare the results.

Then ask yourself why does it respond that way?


Response of Linear System to Stochastic Process Somehow I am getting the variance{u(n)} equal to '0' !! This is the case when I take the coefficient 'a' as real. As it is not mentioned in the question I need to find the solution to this question for a complex 'a' also. I would be glad I anybody could help me solve it. I waold also like to clear my concepts about the complex coefficients, as in how does the Z and inverse Z-transform takes place when 'a' is complex.
This is a back problem of the second chapter in the book - Adaptive Filters: Theory and Application by Farhang-Boroujeny

Consider the following LTI system with impulse response $h[n]$

$$ \{v[n]\} \longrightarrow \boxed{H(z)} \longrightarrow \{u[n]\} $$

From the analysis of LTI systems with WSS random inputs, the following results deduced. Given the input $v[n]$ which is zero mean white WSS random process, its ACS (Auto-Correlation Sequence) sequence is:
$$\phi_{vv}[m] = \sigma_v^2 \delta[m]$$ where $\sigma_v^2$ variance of the input.

Then the output ACS is given by:
$$\phi_{uu}[m] = h[m] \star h[-m]^{*} \star \phi_{vv}[m]$$

Similary the input and output PSD (Power Spectral Densities) are related by:
$$S_{uu}(w) = H(w) H(w)^{*} S_{vv}(w) = |H(w)|^2 S_{vv}(w) $$ which follows from the relation that PSD of a WSS RP is the DTFT of its ACS; i.e.,
$$ S_{uu}(w) =\sum_{m=-\infty}^{\infty} \phi_{uu}[m] e^{-j w m} \longleftrightarrow \phi_{uu}[m] = \frac{1}{2\pi} \int_{-\pi}^{\pi} S_{uu}(w) e^{j w m} dw$$

Then the variance of the output is given by the following:
$$ \sigma_{uu}^2 = \phi_{uu}[0] = \frac{1}{2\pi} \int_{-\pi}^{\pi} S_{uu}(w)dw$$

Now using the given information $H(z) = \frac{1}{1-a z^{-1}}$ and $\phi_{vv}[m] = 1 \delta[m]$  you can compute the variance in either way:

$$ \sigma_u^2 = \left( h[m] \star h[-m]^{*} \star \delta[m] \right) |_{m=0} $$ or
$$ \sigma_u^2 = \frac{1}{2\pi} \int_{-\pi}^{\pi} |\frac{1}{1-a e^{-jw}}|^2 dw $$

For example chosing the first approach with $h[n] = a^n u[n]$ yields:
$$
\begin{align}
\sigma_u^2 &= \left( h[m] \star h[-m]^{*} \star \delta[m] \right) |_{m=0} \\
\sigma_u^2 &= \left( h[m] \star h[-m]^{*} \right) |_{m=0} \\
\sigma_u^2 &= \sum_{k=-\infty}^{\infty} h[m-k]h[-k]^{*} |_{m=0} \\
\sigma_u^2 &= \sum_{k=-\infty}^{\infty} h[-k]h[-k]^{*} \\
\sigma_u^2 &= \sum_{k=-\infty}^{\infty} a^{-k} u[-k] (a^{-k})^{*} u[-k]\\
\end{align}
$$

$$ \sigma_u^2 = \sum_{k=-\infty}^{0} (a a^{*})^{-k} = \sum_{k=-\infty}^{0} |a|^{-2k} = \sum_{k=0}^{\infty} (|a|^{2})^k  = \frac{1}{1-|a|^2}$$


Symbol synchronization with M-QAM. M < 16 I working on M-QAM system using a Cortex-M4 processor. Using a alternating symbol as a preamble, I am able to get carrier lock. Symbol synchronization is then achieved by finding the peak. An early late gate is used to adjust for clock differences between modulator and demodulator. An integrator is used to even out the noise from the error detector. However, I can see the peak slipping from the optimum sampling position. Eventually after 500-600 symbols it loses sync. Are there any superior ways of doing symbol sync? Am I missing something?

Are there any superior ways of doing symbol sync?

Yes! As you've noticed, your initial estimation of timing isn't "valid" forever.
This might be caused by clock jitter, or simply remaining sampling rate offset.

It's a very common problem in communication receivers.

What you can of course do is that you simply do the same estimate you did at the beginning periodically (with the already recovered data instead of a preamble, even!), or even simply for every symbol.

A popular way of doing so is basically using the derivative of the pulse shape to track the "flattest" point in the pulse-shaped receive signal even before applying the (typically matched) filter in the receive chain. You might want to watch this video (slides) of GNU Radio's synchronization expert Andy Walls – I like that talk especially since not only is it clear enough to be used as instructional material, it actually also addresses "how it is actually done".


measure angle between two lines, can rotation center constraint be used to refine the initial measured result? I'm trying to measure the angle between two needle positions in an analog meter (gauge).
such as:

After the image processing of several images, I get the processed image below(a image show different needle positions:

Here, the black dot points were the feature points of needle area, so I use the line fitting algorithm(usually by least square method, or Hough transform to find the line) to get several line positions. I can get the L1, L2, so the angle between those two lines was theta1.
But the needle can located in many different places, so I can finally get many positions such as L3, L4...Ln. Those lines should be ideally intersect in one point, which is the rotation center of the needle. I can calculate the equivalent intersection point of all those lines (the red dot), since it is the rotation constraint of the movement.

So, my question is: How can the equivalent point be used to refine the final angle, I mean I want to get a better value of theta1, theta2, theta3....

Thanks.

You could recompute the lines with the constraint that they all go through the red point. If you normally optimize a line

$$y=ax+b$$

by optimizing both parameters $a$ and $b$, you would now only have 1 degree of freedom for each line. If the red point has coordinates $(x_1,y_1)$ then you get
the constrained line equation

$$y=a(x-x_1)+y_1$$

with only one parameter $a$. Now you compute $a$ such that the squared distance to the other dots is minimized, and no matter which value of $a$ you obtain, the line will always pass through $(x_1,y_1)$.


with given Sampling rate, what max frequency human voice can be captured? I'm new to signal processing and sampling rate, I was asked interview question related to this.

This is the exact question asked:

With 8Khz of sampling rate, whats the max frequency that human voice can be captured ?

  Options provided:

a) 4Khz

b) 8Khz

c) 16Khz

d) 32Khz

I have checked few Q/A websites but couldn't find anything helpful. I did find something similar in this SP stack --> link

And I've also found that for human voice range are upto 8Khz for male and 17 Khz for female. I don't know if it's helpful but I thought I should add this in the question.

Any help, would be really appreciated thank you :)

Without more details, when a signal is properly recorded at $F_n$ Hz, one is entitled to expect that all information below the half of it, namely $F_n/2$ Hz, can be recovered, somehow. This is an application of the classical Shannon Sampling Theorem. Twice the maximum frequency is sufficient to sample a signal... in theory. The question says: max frequency that human voice can be captured; don't expect this max frequency (and others) be capture with full accuracy, as suggested in comments by David.

In some cases, this can be extended, as long as the bandwidth can be shrunk. If the actual bandwidth is smaller, like $F_n/4$-$F_n/2$, one might expect to sample it at $2\times (F_n/2-F_n/4) = F_n/2$ without loss, yet I don't think it is the expected answer, cf. Voice Frequency :

  The bandwidth allocated for a single voice-frequency transmission
  channel is usually 4 kHz, including guard bands, allowing a sampling
  rate of 8 kHz to be used as the basis of the pulse code modulation
  system used for the digital PSTN. Per the Nyquist–Shannon sampling
  theorem, the sampling frequency (8 kHz) must be at least twice the
  highest component of the voice frequency via appropriate filtering
  prior to sampling at discrete times (4 kHz) for effective
  reconstruction of the voice signal.


Plotting a density function I am trying to plot the following function in MATLAB

$$f(x) = \lambda e^{-\lambda(x-\Delta)}w(x-\Delta) $$

where
$$w(x) =
         \begin{cases}
           1 &x \geq 0 \\
           0 & {\rm otherwise}\\
         \end{cases}
    $$

I chose a value of $\lambda$ of 1 and $\Delta$ of 2. I am just not sure how to define $w(x)$ and shift the function by delta. Thank you.

You can plot samples of a continuous density function just as you would do with any continuous function. Assuming $w(x)$ represents the continuous-time unit step function,

$$w(x-\Delta) = \begin{cases}{ 1 ~~~,~~~x \geq \Delta \\ 0 ~~~,~~~  \text{otherwise} } \end{cases}$$

then the shifted exponential pdf can be plotted with the following matlab/octave code :

clc; clear all; close all

K1 = -1;        % evaluation interval begins
K2 = 9;         % evaluation interval ends
N = 1000;       % number of points to display
x = linspace(K1,K2, N); % x = domain to use for displaying

lam = 1;        % exponential PDF parameter
delta = 2;      % shift amount
f = lam*exp(-lam*(x-delta));
f = f.*(x>delta);   % to implement w(x-delta)

figure, plot(x,f);
title('shifted exponential pdf for \lambda = 1, \Delta = 2');

The output will be:


Why Runge-Kutta for Quaternion integration in Kalman filter? I'm reading up on Kalman filtering at the moment. In particular, I'm interested in using the "extended" and "unscented" variants for IMU sensor fusion and calibration.

In A comparison of unscented and extended Kalman filtering for estimating quaternion motion, quaternions are used to represent 3d rotation.

I understand unit quaternions can be used to represent a 3d rotation. They suit for representing absolute attitude (a rotation from a universal reference), relative rotation, or angular velocity (a rotation representing rate per second, or some other fixed time period).

However, this papers discuss using Runge-Kutta integration, specifically RK4. It uses RK4 with the quaternions but doesn't seem to furnish details of what this involves or why it is necessary. Here's the part of the paper that mentions it…

  Given the state vector at step k − 1, we first perform the prediction step by finding the a priori state estimate xˆ−k by integrating equation 1 [f = dq/dt = qω/2] through time by ∆t (i.e., 1.0 divided by the current sampling rate) using a 4th Order Runge-Kutta scheme.

I've encountered Runge Kutta before for integrating positions in kinematics. I don't really understand how or why it would be needed here.

My naive approach would be to simply multiply the existing attitude q by the angular velocity ω to get the expected new q – I don't see why numerical integration is necessary here? Perhaps it's to "scale" the unit time ω to the change that occurs in in ∆t, but surely that can be done very simply by directly manipulating of ω (raising it to fractional power ∆t)?

Anyone know?

I think the confusion comes from the authors not parameterizing things clearly.  Furthermore, by switching to geometric algebra rather than quaternions, some additional confusion can be cleared up.

The main difference between normal vector algebra and geometric algebra is that we can multiply vectors. So of $e_x$, $e_y$, and $e_z$ are our (orthonormal) basis vectors, we also have $e_x e_y$ which is not a vector but a bivector which can be thought of as an oriented plane element.  (By "pure vector quaternion" they mean a bivector; you are correct that that means $0$ $w$ part and it's reasonable to say it has a $0$ "real" part.)  A key property is that the basis vectors anti-commute with each other, i.e. $e_x e_y = -e_y e_x$.  This leads to $(e_x e_y)^2 = -1$.  In general, any unit bivector squares to $-1$, which means we can apply Euler's formula: $$R = e^{\theta B} = \cos\theta + B\sin\theta$$ where $B$ is a unit bivector and we call $R$ a rotor.  A unit quaternion is just a 3D rotor.  (Complex numbers are just the even sub-algebra of the 2D geometric algebra.  Here we are looking at quaternions, the even sub-algebra of 3D GA.)

If we want a rotation by $\theta$ in the plane of a unit bivector $B$ we use the rotor $R = e^{\frac{\theta}{2}B}$.  (The half comes in because we rotate a vector $v$ via $Rv\widetilde{R}$; see the link above.)  Of course, we want to allow both the plane and the angle to vary so we define the (non-unit) bivector $\Theta(t) = \theta(t)B(t)$ so if $$R(t) = e^{\frac{1}{2}\Theta(t)}$$ then $$\dot{R}(t) = \frac{1}{2}\Omega(t)R(t)$$ where $\Omega = \dot{\Theta}$.  This is their $f = q\omega/2$; their $q$ is our $R$ and their $\omega$ is our $\Omega$.  Now if we have $R(t_0)$ and we want $R(t_0 + \Delta{}t)$ we need to integrate $\dot{R}$ from $t_0$ to $t_0 + \Delta{}t$.

The scheme you describe is roughly analogous to doing forward Euler integration, which essentially assumes their $q$ and $\omega$ our $R$ and $\Omega$ is constant over $\Delta{}t$.  RK4 is just a better integration method.  I doubt there is any special reason they chose to use RK4 as opposed to other integration methods.  It's just the default choice typically.


Real-valued ringing when zero-padding odd-length FFT So I'm trying to write a frequency-domain interpolator that zero-pads the frequency response of a signal and inverse transforms.  There's two cases I have to deal with:

Even-length response - have to split the $F_s/2$ bin because it's ambiguous.  So I copy the negative part of the spectrum, and add n*(interp-1)-1 zeros in between.
Odd-length response - there is no $F_s/2$ bin so just split positive/negative frequency and insert n*(interp-1) zeros between them.

The code that does the zero-padding can be seen here

// Copy negative frequency components to end of buffer and zero out middle
//  inp    - input buffer of complex floats
//    n    - transform size
//  interp - interpolation amount
void zero_pad_freq(cfloat_t *inp, size_t n, size_t interp) {
    if ((n % 2) == 0) {
        memmove(inp + n*interp - n/2, inp + n/2,     n/2*sizeof(cfloat_t));
        memset (inp + n/2 + 1, 0,       (n*(interp-1)-1)*sizeof(cfloat_t)); // Duplicate Fs/2 so we need one less zero

        inp[n/2]          /= 2.0;
        inp[n*interp-n/2] /= 2.0;
    } else {
        memmove(inp + n*interp - n/2, inp + (n+1)/2, n/2*sizeof(cfloat_t));
        memset (inp + (n+1)/2, 0,         (n*(interp-1))*sizeof(cfloat_t));
    }
}

The first case is working fine, I'm testing it on a chirp signal and it interpolates just fine, there's a little numeric noise, but it's round tripped through an FFT so what can you do (first $50 \mu s$ or so of the signal show):

The problem is with the odd-length transform, I'm getting a pretty heinous transient response on the real samples only ($50 \mu s$ again, real):

The imaginary channel has a small ripple on it, but not nearly as bad:

It's like I've screwed up my $F_s/2$ bin in the odd case, but there is no $F_s/2$ bin, so I'm very puzzled.  Anyone have any thoughts?

By zeroing out high frequency bins, you have effectively multiplied the spectrum of the signal with a rectangular function.  Multiplication in frequency is convolution in time and the Fourier pair of a rect is a sinc.  So what you really did was convolve the time domain signal with a sinc with the width of the main lobe of the sinc inversely proportional to length of the rect.  This is why the numerous filter designing techniques like Parks-McClellan design in what's called a "transition region" or "transition" band so that there is not an instantaneous change in the frequency response of the filter.  These filter design techniques are important because the "ideal" filter like you used has such undesirable effects in the time domain.


How to show Image FFT with origin as DC? The page Here explains how to do the discrete fourier transform on a 2d image.  It mentions this:

  In most implementations the Fourier image is shifted in such a way
  that the DC-value (i.e. the image mean) F(0,0) is displayed in the
  center of the image. The further away from the center an image point
  is, the higher is its corresponding frequency.

I have the implementation working, but I'm not sure how to shift my image to make DC be the center of the image.  Does anyone know how to do that?

The best I can think of is to push everything right 1/2 of width and up 1/2 of height, and wrap the pixels around, but this doesn't seem like the correct thing to do.

Does anyone know how to make this transform?  Right now, DC is in the lower left of my image, at pixel (0,0)

The quadrants of the FFT should be shifted Q1->Q3, Q2->Q4 ...
For a [RxC] matrix X this may be accomplished by shifting as:

tmpX = shift_row(    X, floor(R/2) )
Y    = shift_col( tmpX, floor(C/2) )

where shift_row and shift_col is your shift operator in each dimension and floor is the round towards zero operator.

For the inverse FFT shift you should shift in the same way but with a value -ceil(R/2) and -ceil(C/2), i.e rounded towards infinity


Valid approach to implementing a band pass filter Suppose I want to perform Spectral shaping on a signal, i.e, modify the gains in a band of frequencies. Would there be any difference if I do that using the following two methods? -

Use a band pass filter to apply the gain in the band.
Calculate the bin numbers corresponding to the band, and apply the gain in those bins.

I'm interested to know if the second method is a valid approach, if not, why is it discouraged?

then you'll be doing either "Overlap-add" (OLA) or "Overlap-save" (OLS) which is sometimes more clearly called "Overlap-scrap".  you need to read up about this.

then simply "apply[ing] gain in those bins" is multiplication in the frequency domain, which means convolution in the time domain.  if your frequency domain function is a rectangle function (some number for "those bins" and another number, like 1, for the other bins).  then you need to compute the impulse response that goes with it and make sure that FIR impulse response is not too long for your OLA or OLS.  do you have a good book that describes OLA or OLS?

but a bandpass filter in the time-domain (your option 1.) is simple and can be implemented simply.  maybe more difficult to design your IIR or FIR bandpass filter.


Impulse and step response While trying to determine the unit step and impulse response for a discrete LTI system I had several problems:

$$
\begin{aligned}
(Hx)[n]&=(Hx)[n-1]+\frac{1}{N}(x[n]-x[n-N]) \\x[n]&=(Hx)[n]=0\;\forall\;n\lt0
\end{aligned}
$$

Starting with the Impulse response:

$$h[n]=(H\delta)[n]=(H\delta)[n-1]+\frac{1}{N}(\delta[n]-\delta[n-N])$$

Can I argue that because of the condition $x[n]=(Hx)[n]=0\;\forall\;n\lt0$ the impulse response only exists where $n\gt 1\;\&\;n\gt N?$ But what about the other cases?

Concerning the Step response I am quite lost and am not even sure if this is how it is calculated:

$$a[n]=(H\sigma)[n]\sum_{k=-\infty}^{\infty}h[k]\sigma[n-k]$$

Help as to how to proceed and solutions are greatly appreciated

This is a homework type problem, so I'll give you a few hints to help you solve it yourself (and learn something while doing so).

From what is given, we can assume that the given difference equation describes a causal discrete-time system. Let $h[n]$ denote the impulse response. It must satisfy

$$h[n]=h[n-1]+\frac{1}{N}\big(\delta[n]-\delta[n-N]\big),\qquad n\ge 0\tag{1}$$

with $h[-1]=0$. Note that the second term on the right-hand side of $(1)$ only contributes to $h[n]$ for $n=0$ and $n=N$. Now start at $n=0$ and continue up to $n=N$ to see what happens. It will turn out that the impulse response has a finite length, i.e., the given system is an FIR filter.

The step response is the convolution of the impulse response with the unit step sequence, which simplifies to

$$a[n]=\sum_{k=0}^{\infty}h[k]u[n-k]=u[n]\sum_{k=0}^nh[k]\tag{2}$$

So the step response $a[n]$ is just the cumulative sum of the impulse response $h[n]$.


How to plot magnitude and phase response by hand if I have the Transfer Function? I have the transfer function of the system, which is:

$$H(z) = \frac{1-z^{-1}}{5(1+2z^{-1})}$$

How do I sketch the magnitude and phase response?

I'm sorry for the bad formatting, it's my first time post a question.

Thank you very much in advance for the help!

Use the transformation $z = e^{j\omega}$, you will get

$$H(e^{j\omega}) = \frac{1-e^{-j\omega}}{5\left(1+2e^{-j\omega}\right)}$$

Solving this (using $e^{j\omega} = \cos \omega + j \sin \omega$), you should get something like (please double check):

$$H(\omega) = \frac{\left(\cos\omega-1 \right) + j4\sin\omega}{5\left(5+4\cos\omega\right)}$$

Here,

$$\text{Re}\left(H(\omega)\right) = A = \frac{\left(\cos\omega-1 \right)}{5\left(5+4\cos\omega\right)}$$

and, $$\text{Im}\left(H(\omega)\right) = B = \frac{4\sin\omega}{5\left(5+4\cos\omega\right)}$$

Now, the magnitude response will be

$$\left|H(\omega)\right| = \sqrt{A^2+B^2}$$

and the phase response will be

$$\angle{H(\omega)} = \tan^{-1}\frac{B}{A}$$

You will get both responses as a function of $\omega$, just vary $\omega$, calculate the values and plot the response.


Least Squares Filter Design: Deriving the Objective Function I'm following the derivation in this paper A Comb Filter Design Using Fractional-Sample Delay to obtain the objective function for the least-squares filter design.

N-order FIR filter:

$H(z) = \sum_{n=0}^N h(n)z^{-n}$

Frequency response:

$H(\omega) = \mathbf{h}^T\mathbf{e}(\omega) = \mathbf{e}^T(\omega)\mathbf{h}$

where $\mathbf{h}=[h(0)\quad h(1) \quad\dots \quad h(N)] $ and $\mathbf{e}=[1\quad e^{-j\omega} \quad\dots \quad e^{-jN\omega}] $

Least squares error:

$J(h)=\int_{\omega \in R^+\cup R^-} |H(\omega)-F_d(\omega)|^2 d\omega$

where $R^+=[0,\alpha\omega]$ and $R^-=[-\alpha\omega,0]$

Which is rewritten in the quadratic form: $J(h) =h^TQh - 2h^Tp + c$

$Q$, $p$, and $c$ are given in the paper as follows:

$F_d(\omega)= e^{-jD\omega}$ and $h(n)$ is real.

I've been having trouble with getting these expressions for the Q matrix and p vector. How can I obtain them?

$$
J(h) = \int_{R^+UR-}|e^Th-F_d|^2d\omega\tag{1}\\
= \int(e^Th-F_d)^H(e^Th-F_d)d\omega\\
= \int((e^Th)^H(e^Th) + F_d^HF_d -(e^Th)^HF_d - F_d^He^Th)d\omega\\
= h^H(\int(e^T)^He^Td\omega) h + \int |F_d|^2d\omega -\int (2 Re\{(e^Th)^HF_d\})d\omega
$$
The second term in above integral is the integral of $L_2$ norm of $F_d$ in the region $R^+ U R^-$. Since $F_d(\omega)$ is symmetric about $\omega = 0$, this integral will be $2 \times \int_{R^+}|F_d(\omega)|^2d\omega$ which is $c$ in your questions.

The first term, the argument inside integral is the outer product of column vector $(e^H)^T$ with row vector $e^T$. Again, this will conjugate of each for the regions $R^+$ and $R^-$ (shown later in the appendix). So integral over both regions will result in the real value of integral $e^H e$.  Also, since $H$ is real, $h^H=h^T$. So the first term can be rewritten as (for simplicity dropping $T$ notation as it is assumed to be a column vector)
$$
h^T(\int_{R^+}2Re\{ ee^H\}d\omega )h\\
= h^TQh
$$

The third term is the result of sum of 2 conjugate terms. Again this is symmetric about $\omega=0$. Also, since $H$ is real, $h^H=h^T$. Taking out $2h^T$ from the integral we have
$$
2h^T\int_{R+}2Re\{F_de^H \}d\omega = 2h^Tp
$$
Therefore, summing and rearranging the 3 terms,
$$
J(h) = h^TQh - 2h^Tp + c
$$

Appendix: Showing that $ee^H$ is conjugate of each other in $R^+$ and $R^-$. $E_{mn}$, element of $ee^H$ is $e^{-jm\omega}e^{-jn\omega} = e^{-j(m+n)\omega}$. In $R^-$, this element is $e^{j(m+n)\omega}$ which is $E_{mn}^*$


Centered Fourier transform What is the difference between the non-centered and centered Fourier transforms? In other words, when should you use one instead of the other?

Non-Centered: $\quad \displaystyle X_1(f)=\sum\limits_{n=0}^{N-1} x[n] \, e^{-i 2 \pi f n}$
Centered: $\quad\displaystyle X_2(f)=\sum\limits_{n=0}^{N-1} x[n] \, e^{-i 2 \pi f\left(n-\frac{N-1}{2}\right)}$

I am currently enrolled in a time series-analysis class, and our textbook uses both definitions without explaining the difference or the reason for the two different definitions.

Any help would be appreciated.

the "Non-Centered" is the same as the Discrete Fourier Transform

$$ X[k] \triangleq \sum\limits_{n=0}^{N-1} x[n] \, e^{-i 2 \pi \frac{nk}{N}} $$

when frequency $f = \frac{k}{N}$ is normalized to the sample rate ($f=\tfrac12$ means the Nyquist frequency). And it assumes all of $x[n]$ is for non-negative values of time.  "$n \ge 0$" is "time" in units of the sampling period.  $n=0$ means time at $t=0$.

the "Centered" is the same as the DFT with a linear phase-shift factor

$$ e^{i \pi (N-1) f} X[k] $$

when $f = \frac{k}{N}$ and assumes $x[n]$ is for instances of time that are both positive and negative.  $n=\frac{N-1}{2}$ means time at $t=0$.  if $N$ is odd, then the "time origin" is exactly the time of the sample at $x[\tfrac{N-1}{2}]$.  but if $N$ is even, the time origin (when $t=0$) is midway between adjacent samples $x[\tfrac{N}{2}-1]$ and $x[\tfrac{N}{2}]$.  if you use the DFT to compute your "Centered" Fourier transform, you probably want to remove that linear phase-shift factor by multiplying the DFT result, $X[k]$, with $e^{-i \pi (N-1) f}$ or $e^{-i \pi (N-1) k/N}$.  that will place the time origin (when $t=0$) where $n=\frac{N-1}{2}$.


How to apply a mel-filterbank to a signal? I found this description. According to this, the next step is a "binning" where

  binning means that each FFT magnitude coefficient is multiplied by the corresponding filter gain and the results accumulated.
  What is the exact procedure they are describing?

Binning is an averaging operation on the (squared) magnitudes of the DFT. You would maybe have $256$ DFT bins but only around $20$ outputs of the filter bank. So you need to average groups of DFT bins to reduce the dimension from $256$ to $20$. For a mel-scaled filter bank, the averaging functions (kernels) are usually triangular, i.e. the center DFT bins get more weight than the rest. As an example, let the first kernel that is used for averaging the lowest DFT frequencies be stored in an array $W_1[k]$, $k=0,\ldots,K-1$, then the output of that filter bank channel is computed as

$$y_1=\sum_{k=0}^{K-1}W_1[k]R[k]\tag{1}$$

where $R[k]$ are the (squared) DFT magnitudes. In this way you compute all $20$ filter bank outputs. You just need to know (define) the weights $W_i[k]$ and the range of DFT indices that are averaged by the respective kernel.

This answer to a related question might also be helpful (also have a look at the comments).


Super resolution with aliased signal SKIP to Questions below if you want :

In super resolution I can create a high resolution signal from a low resolution signal.
But it is also possible that the low resolution signal contains aliasing. The problem with aliasing is that it cannot be distinguished between the real low resolution frequencies. But in super resolution you can make multiple low resolution measurements with aliasing and calculate the shift between the signals (afaik) that can help to reconstruct a higher resolution signal. Check some papers like here :
Multi-Frame super resolution,
SR from aliased images

What I am thinking is if I can make multiple shots from one scene with aliasing and construct a higher resolution image out of it, how could I apply this to vibration from a machine.
In the first paper google shows how they use the vibrations from a hand for taking multiple photos from one scene. The vibrations from a human hand help to capture different frequencies with aliasing. Later on, all merged can be used for super-resolution.

Questions:

I have a fixed sensor on a machine which basically means I have a
camera fixed on a stand(compare photo camera wich takes multiple pictures at same time from the same scene with aliasing to construct a higher resolution). How would this affect the super-resolution
in comparison to using different angles in parallel to have multiple measurements with different offsets?
Could I apply this on
vibration with making a long measurement, dividing the measurement
in different blocks and trying to merge them as they where different
captures, instead of using different angles?

I don’t think that you can do multi-frame super resolution in the traditional sense if there is no camera rel scene movement, nor if the camera applies «proper» Nyquistian spatial prefiltering?

I think of good old interlacing as an example of potential super resolution. You get a 2-frame cadence consisting of separate time and spatial samples, usually moderately filtered. If there is no movement, you can get up to double spatial resolution. If there is movement that can not be resolved, you get high temporal resolution.


Radix-2 vs Radix-4 FFT What are the differences between radix-2 FFT and radix-4 FFT,

Are the only differences following two:

In case of radix-2 $N$ is a number that is a power of 2 and in case of radix-4 $N$ is a number that is a power of 4
Incase of radix-2 the butterfly diagram increases or decreases density by a factor of 2, while in case of radix-4 the butterfly diagram increases or decreases density by factor of 4

The reason the Radix-4 FFT is of interest is in the simplicity of multiplying by $\pm j$ in actual implementation. Below shows the Radix-4 4 point DFT core processing element as part of the Radix-4 FFT Butterfly in comparision to the Radix-2 FFT butterfly (with 2 point DFT core processing element) and the resulting decrease in number of operations, applicable when the input signal is of a length that is a power of 4 (or for portions of the signal that are).

To further see how the radix-4 algorithm reduces multiplications consider the number of real multipliers and additions in a complex multiplier:
$$V_1V_2 =  I_{out}+jQ_{out}= (I_1+jQ_1)(I_2+jQ_2) = (I_1 I_2 - Q_1 Q_2) + j(I_1Q_2+I_2Q_1)$$
Where we see from above there are 4 real multipliers and 2 additions (Where the real and imaginary components are maintained on separate paths so not actually added).
Now consider if $V_2$ above was $j$, then $I_2 = 0$ and $Q_2 = 1$:
$$V_1V_2 =  V_1j = (I_1 0 - Q_1 1) + j(I_1 1 +0 Q_1) = -Q_1 + jI_1$$
We see that multiplying a complex number $I_1+jQ_1$ by $j$ is done by simply changing the sign of $Q$ and then swapping $I$ and $Q$. Similarly multiply by $-j$ would be done by changing the sign of $I$ and then swapping $I$ and $Q$.
The radix-4 allows for an increased number of multiplications by $\pm j$ where the radix-2 solution minimizes all complex multiplications including multiplying by $\pm j$.


Why there is error between My DFT and built in Matlab FFT? I'm trying to implement DFT/IDFT and FFT/IFFT in Matlab and C.

I got wrong values while I implement IDFT/IFFT. So I'm trying to find a mistake.

Just for DFT, I used this algorithm to implement it in Matlab (from matlab's fft doc)
But I really don't know why there is very small errors when I run this code.
and even my Radix-2 FFT have same problems.
Is there anything I'm missing or can I ignore this values? (mean of 256 values : -1.6881e-25)

Thank you for reading this.

for i=1:256
    x(i) = sin(2*pi*i*80) + sin(2*pi*i*46);
end

X = zeros(256,1);
for k=1:256
   for n=1:256
       X(k)=X(k) + x(n)*exp((-j*2*pi*(n-1)*(k-1))/256);
   end
end

built_in = fft(x);
error = abs(X) - abs(built_in)';

plot(abs(X)); hold on
plot(abs(built_in));
plot(error);
legend('My code', 'Matlab','error');

with result:

Your error is the yellow line at the bottom.  i see no difference between your DFT and the built-in FFT from MATLAB.


How can I sample a signal at 4 samples per cycle? I am in Chapter 5 of "Understanding Digital signal processing" 3rd edition by Richard G. Lyons.
I am stuck with Section 5.4.

In order to get the bandpass filter centred at $F_s/4 \mathrm{Hz}$, I need to multiply the low pass filter coefficients by a sinusoid of $F_s/4 \mathrm{Hz}$ samples at 4 samples per cycle. I have tried to do this in MATLAB but cannot figure out how to sample at four samples per cycle. Since samplerate was not given, I created a sinusoid:
$$
n=0:30;\\
x=\sin(n\pi/2)%
$$
Is this correct for $F_s/4$ signal?

How do I sample it at four samples per cycle in MATLAB?

I feel like I missing something very fundamental.

You need to understand that $f=f_s/4$ and "$4$ samples per cycle" are two ways of saying the same thing:

$$\text{ # samples per cycle}=\frac{f_s\text{ samples per second }}{f \text{ cycles per second}}$$

The number of cycles per second is equal to the ratio $f_s/f$, where $f_s$ is the sampling frequency, and $f$ is the signal's frequency.

So you just need a sampled sinusoid with a frequency $f=f_s/4$ (or, equivalently, $4$ samples per cycle). The value of $f_s$ is irrelevant.

$$x[n]=\sin(2\pi f\cdot n/f_s)=\sin(n\pi/2)$$

which is exactly what you have already come up with. Multiply this signal with your low pass impulse response, and you'll get a band pass filter centered at $f_s/4$.


Identifying the age of a speaker using speech recognition I am just wondering if it's possible to gather the age group of someone speaking in a sample?

For example, could I tell whether the sample is a child, teenage, adult, etc?

Have a look at the acoustic scale website for some inspiration

http://www.acousticscale.org/wiki/index.php/Main_Page

You can estimate age based on a combination of vocal tract length (vtl) and pitch. Both of these attributes can be extracted from vowel sounds. Children will have short vtls and high pitch, adult males will have low pitch and long vtl, adult women will have fairly long vtls and higher pitch. Dwarfs can have short vtls and low pitch. The castrato will have long vtl and high pitch.


Noise Power Spectrum calculation I want to compute the 1-D noise power spectrum (NPS) for a uniform CT image. The equation of the 2D NPS is as following:
$$
NPS(f_x, f_y) =\frac{1}{N}\sum_{i=1}^N|DFT_{2D}[I_i(x,y)-\bar{I}_i]^2\frac{\Delta_x\Delta_y}{N_xN_y}
$$
whre $DFT$ is the discrete fourier transform, $I_i$ is a ROI, $\bar{I}$ is the mean of the ROI. The 1-D dimensional (or radial) NPS is defined as:
$$
f_r = \sqrt{f_x^2+f_y^2}
$$
Assuming only one ROI and obviating the units factor, I compute the 2D NPS as:
fft = np.fft.fft2(roi-np.mean(roi))
fft = np.fft.fftshift(fft)
NPS_2D = np.abs(fft)**2

How can I compute from here the 1-D radial NPS?

Obtain parallel beam sinogram of your image, then take the FFT of the sinogram (if you consider the Fourier transform in polar coordinate, then Fourier slice theorem tell us value of Fourier transform over a radial line at a given angle is equal to the fourier transform of the projection over that angle), after that average the power over different angles. To obtain sinogram, if there isn't any available function, you could rotate the image by desired angle, then sum over columns or rows.


Changing units in Fourier filtering I'm applying a radial low-pass filter on a $H{\times}W$ pixel 2D image $I$ where $H \neq W$. To implement this filter, I apply 2D DFT to $I$ to produce $I'$, the $H{\times}W$ frequency domain representation of $I$. I apply an element-wise multiplication $I' \odot M$ with a filter matrix $M \in \{0, 1\}^{H{\times}W}$. Let $u \in \{0, H-1\}, v \in \{0, W-1\}$ be indices in $M$. $u$ represents cycles/entire height of image, while $v$ represents cycles/entire width of image. An ideal low-pass filter with radius $R_{img}$ is given by
$$\sqrt{u^2 + v^2} < R_{img},$$
whose units are in cycles/image. I want to change the units of the filter parameter $R_{img}$ from cycles/image to cycles/mm, given that $\gamma$ = mm/pixel.

Here is what I have so far, please let me know if this is correct. Since the filter radius can't exceed $\min(H, W)$, let $f_s = \min(H, W)$ denote the sampling frequency of the image across both dimensions in terms of cycles/image, and $u, v$ denote row and column indices of the $H{\times}W$ representation of the image in frequency domain after applying shifted 2D DFT. We then have

$$\sqrt{u^2 + v^2} < R_{img}\\
\frac{1}{f_s \gamma} \sqrt{u^2 + v^2} < \frac{R_{img}}{f_s \gamma}\\
\sqrt{\left(\frac{u}{f_s \gamma}\right)^2 + \left(\frac{v}{f_s \gamma}\right)^2} < \frac{R_{img}}{f_s \gamma},$$

where the units have changed from cycles/image to cycles/mm. In other words,
$$R_{mm} = \frac{R_{img}}{f_s \gamma}.$$

Short answer: your conclusion seems ok.

Concept of units make sense only in the world of physics. Mathematical objects such as DFT trasforms do not have any kind of units but just their values, quasi exceptions being the planar angle meausure of radian and volumetric angle mesaure of the steradian. However, when a variable represents a physical quantity (as almost always is the case) then it's customary to associate units with those mathematical variables as well.

Let's do a 1D analysis of those units, from which the 2D version follows straightforwardly. Discrete-time variables are indexed sequences such as $x[n]$ or $X[k]$ which may represent samples of a continuous-time signal $x(t)$ and the DFT of $x[n]$ respectively.

We can assume that a sequence $x[n]$ bears a unit which is the unit of the associated physical variable $x(t)$ from which $x[n]$ was sampled. So if $x(t)$ was a voltage in units of Volt,   then the unit of $x[n]$ will be Volt too.

The index $n$ does not have any units which can be seen by the following uniform sampling relation :

$$ t_n = n T_s $$

where $t_n$ is the sampling time with a unit of second (s), and $T_s$ is the sampling period with a unit of second (s) too; hence the index number $n$ is unitless. However, some like to call $n$ with a unit of samples. Yet this is not only unnecessary but also as meaningless as associating a unit with room numbers in a hotel. Nevertheless it's customary to do so and I do also use it.

The 1D-DFT $X[k]$ associated with a given sequence $x[n]$ is defined as follows:

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

The complex exponential function (as well as other trigonometric functions such as sine, cosine, tangent) produce unitless outputs as pure numbers. The trigonometric functions expects angles in units of radians as their input arguments.

Therefore, the unit of the DFT $X[k]$ will be that of the sequence $x[n]$.

The unit of the DFT frequency bin index $k$ is found from the unit of the complex exponential argument :

$$ \frac{2 \pi }{N} k n = \theta $$

In this expression, the sample index $n$ and the total sample count $N$ (or frame length as some would like to call) have no units (or they both have the unit of samples, but it will cancel neverthelesss) and the angles $\theta$ and $2\pi$ have the unit of radians hence the DFT frequency index $k$ is unitless too.

How about calling the frequency variable $k$ as cycles per frame length ? Ok, it is also preferred and some may find it informative to call $k$ as the number of cycles a sine wave fits into the given frame length of analysis.

With the same reasoning, a 2D intensity image sample sequence $f[n,m]$ will have units of intensity and the spatial index variables $m,n$ are either unitless or given a samples unit. And for the 2D-DFT $X[u,v]$, the frequency index variables $u,v$ are unitless. Nevertheless, again some prefer calling them in the units of cycles per image frame length and you can also do it.

Now coming to the units of the expression

$$ \sqrt{ u^2 + v^2 } < R ,$$

the variable $R$ is either unitless as explained above or can be associated with the unit of cycles per frame length as some would prefer. From your expanded questions, it became clear that you wanted to change the scale of this actually unitless parameter from cycles per frame length to cycles per mm.

Now, to find out the continuous-space frequency in cycles per meter associated with a given DFT frequency bin index $u$ or $v$ , you would actually perform the following conversion :

$$ f_u =  \frac{u}{H} F_s $$

where $u$ is the unitless DFT frequency bin index count, $H$ is the number of samples (frame length in samples or image height in your example) and $F_s$ is the sampling frequency in cycles per meter (I assume). In your project you set $f_s = H$ (assuming $H = \min(H,W)$) and you give a pixel-size parameter of $\gamma$ in mm/pixel units. Then this $\gamma$ is the inverse of sampling frequency $F_s$ in samples per mm. And the above formula becomes :

$$ f_u =  \frac{u}{f_s} F_s = \frac{u}{f_s} \frac{1}{\gamma} = \frac{u}{f_s \gamma} $$

Finally since the ideal circular lowpass filter bandwidth radius parameter $R$ is a same kind of variable with DFT frqeuency bin index variable $u$ or $v$, therefore its conversion from cycles per frame length to cycles per mm will be done using the same formula:

$$ r =  \frac{R}{f_s} F_s = \frac{R}{f_s} \frac{1}{\gamma} = \frac{R}{f_s \gamma} $$

where $r$ has the units of cycles per mm.


Does same data inserted to all antenna at beamforming? I wonder do I understand beamforming right.
Question 1
Does same data inserted to all antenna at beamforming? Like figure below? The only different is through different phase shift? If it is right how about Massive MIMO which use 64,128 antennas? The all 64, 128 antennas get same data?

Question 2
Does concept of layer that I understand is right? Like figure below? I mean the two resource blocks transmitted at the same time and same frequency at same antenna?

Question 1
All the antennas might get the same data (besides some phase shift as you noted), but they don't have to. In your picture, you have a blue beam and an orange beam that toggle on and off over time. There is nothing stopping you from using a subset of the antenna elements for each beam and transmit both simultaneously. For example, a 64 element array might use 32 elements to form the blue beam, 32 elements to form the orange beam and then transmit both at the same time while carefully designing the beamforming weights to place nulls (blue beam should have a null in orange beam direction so that UE2 doesn't experience interference from the transmission intended for UE1).
Question 2
The usual flat flading MIMO equation for a single symbol is:
$$ \mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{w} $$
where $\mathbf{y} \in \mathcal{C}^{N_R \times 1}$, $\mathbf{H} \in \mathcal{C}^{N_R \times N_T}$, $\mathbf{x} \in \mathcal{C}^{N_T \times 1}$, and $\mathbf{w} \in \mathcal{C}^{N_R \times 1}$. Notice how there is no mention of number of layers, this equation provides a mapping from transmit to receive antennas. The notion of layer mapping is kind of implicitly built in to this equation and not always shown. The equation including layer mapping is:
$$ \mathbf{y} = \mathbf{H}\mathbf{W}\tilde{\mathbf{x}} + \mathbf{w} $$
where $\mathbf{W} \in \mathcal{C}^{N_T \times N_L}$ and $\tilde{\mathbf{x}} \in \mathcal{C}^{N_L \times 1}$. The matrix $\mathbf{W}$ is called the precoding matrix and does the mapping from layers to transmit antennas. This matrix will decide exactly what data goes on which antenna. As you say, the idea is that this spatially multiplexed data can then be transmitted on the same time/frequency element and it is the job of the receiver to demultiplex the data to get the original symbols intended for that UE back.


IIR issues when order is high Thanks to all the people who helped me here, I finally achieved to implement a notch filter (*), with coefficients given by :
http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html

When I do order=2 (number of poles=4, according to this website), all is working very well with :

y[n] = a0*x[n-4] + a1*x[n-3] + a2*x[n-2] + a3*x[n-1] + a4*x[n]
       + b0*y[n-4] + b1*y[n-3] + b2*y[n-2] + b3*y[n-1]

Great!

When I try to increase order (example : order=4 or even order=3 => number of poles=6), then the filter doesn't work anymore. Very quickly (after n=10 or 100), the signal is saturated to the max value : overflow !

(Nota : I work in double (64 bits) floating point).

Somebody told that the higher the order, the more unstable ! What to do then to solve this problem ?

Thank you very much in advance.

(*) The IIR filter is even ... linear-phase ... thanks to the trick that someone gave me (apply the filter on the signal, then reverse the signal, then apply the filter again, then reverse the signal) !

Break higher order filters down into cascaded second order sections (or biquads). This is covered in every text book on DSP such as https://ccrma.stanford.edu/~jos/filters/


2D Wavelet Transform Explanation I am new to wavelets and for several times I stumbled over these 2D wavelet transformation pictures, but I do not understand how they are created.

Can anybody try to explain in his own words how this wavelet transformation in 2D works and how or what we can see it on the picture?

I took this image from http://graphicon.ru/oldgr/courses/cg_el01/lect_wav/wav_cg/fig1.gif

You can notice that:

the top left looks like the original image, albeit smaller
the top right shows mainly horizontal features
the  bottom left shows vertical trends
the  bottom right is much less clear

What you show "is not wavelets" per se. It is another image, normally with the size of the original one, with four quadrants. However, it illustrates some wavelet features. The top left is a coarse approximation of the image, resulting from filtering and downsampling,  obtained from a scaling function. The three others are details, with at least one high pass component. In the classical version, this image combines low- and high-pass on rows, and the same on columns, each followed by downsampling on both directions, which explains the above items. If you iterate the process on the low-pass image, under some conditions, you get a separable two-dimensional wavelet transformation. The key point is that in your image, all the original information remains.


Analog Square Root Raised Cosine Filter I am wondering if it is possible to realize an analog square root raised cosine filter.
Any idea about how to do that with Matlab?

Here's Octave code for a (digital) RRC filter:
pkg load signal;  % Octave needs this; MatLab doesn't
Fs = 16000;  % sample rate
Rs = 400;    % symbol rate
sps = Fs/Rs; % samples per symbol

%
% Root raised cosine pulse filter
% https://www.michael-joost.de/rrcfilter.pdf
%
r = 0.22; % bandwidth factor
ntaps = 8 * sps + 1; % Pulse duration is 8 symbols

st = [-floor(ntaps/2):floor(ntaps/2)] / sps; % symbol time = t/Ts values
hpulse = 1/sqrt(sps) * (sin ((1-r)*pi*st) + 4*r*st.*cos((1+r)*pi*st)) ./ (pi*st.*(1-(4*r*st).^2));

% fix the removable singularities
hpulse(ceil(ntaps/2)) = 1/sqrt(sps) * (1 - r + 4*r/pi); % t = 0 singulatiry
sing_idx = find(abs(1-(4*r*st).^2) < 0.000001);
for k = [1:length(sing_idx)]
    hpulse(sing_idx) = 1/sqrt(sps) * r/sqrt(2) * ((1+2/pi)*sin(pi/(4*r))+(1-2/pi)*cos(pi/(4*r)));
endfor

% normalize to 0 dB gain
hpulse = hpulse / sum(hpulse);

For practical reasons, the filter impulse response is truncated to 8 symbols in duration, as there isn't much benefit in making it any longer.
I'm not sure what you mean by "analog" RRC filter in Matlab.  The continuous frequency response of the RRC filter is well understood.  It results in an infinitely long time domain impulse response, that is the equation in your question with a scaling factor out front.
In the above code, you can make the sample rate as high as you like, to get more closely spaced samples, to have the RRC filter taps get an ever closer approximation to the continuous function in MatLab.
Matlab itself may have facilities for describing continuous functions in the system toolbox, but I am not familiar with them.
As a technical reality, the infinitely long, continuous RRC filter is not realizable in real-life analog components, as it extends infinitely into the past.  Could you implement a truncated continuous RRC filter with analog components?  Well maybe, but why?  Its use is exclusive to digital modulations and it is easily implemented digitally.


Is there such a thing as a continuous phase BPSK? My knowledge of signals and signal processing is limited, so I apologize in advance for my naïve and possibly confusing question.

In the case of BFSK, the frequencies are altered to represent the data (1's and 0's). The frequency transitions are abrupt. You can use Continuous Phase BFSK to "line up" the transitions in frequency.

Now, for my question: lets say I am using BPSK which uses changes in phase to represent the data. If the phase changes from 0 to 180 (or vice versa), wouldn't that have a similar abruptness as BFSK? Is there something I can do, similar to CP BFSK to smooth the transition out?

You're right that the transitions can be abrupt, but you can always choose to use pulse shaping, i.e., instead of using a rectangular pulse, use some smoother function $p(t)$. The corresponding baseband signal is then

$$s(t)=\sum_kA_kp(t-kT)\tag{1}$$

where $A_k$ are the symbols, and $T$ is the symbol length. A smooth function $p(t)$ will avoid abrupt transitions, and will therefore reduce spectral sidelobe levels. I think this paper uses that idea.


Modulation upsampling method I try to implement a qpsk modulator in Octave.
The chain is the state of the art I guess:
$$\text{Symbol mapper} \rightarrow \text{Upsample 4} \rightarrow \text{Root raised cosine filter} $$
My main question is about the upsampling. I try to do it two ways:

zero hold (I think it is called like that - adding factor-1 last values)
zero padding (adding factor-1 zeros, here the spectrum is duplicated and alias filtered by RRC)

I have not tried yet to do a BER measure with AWGN, but eye-diagram and spectrum looks different depending how I upsample. I have the feeling that the spectrum and eye diagram is more correct with zero hold but get a residual bias at the sampling point.
Am I wrong?

You have two methods. The second one is right, because it's actually upsampling: inserting zeros, and then suppressing images (what you call aliases here is actualyl called images) through the desired filter.
The first one is actually the same, but instead of the RRC, you use a rectangular filter: When you think about it, holding the last value is like inserting zeros, and then convolving with a 1-amplitude rectangle that is as long as the number of zeros you've inserted.


Determining the noise floor of a signal in the frequency domain Is there an accepted way of determining the noise floor of a signal by looking at it in the frequency domain?  Is it a matter of averaging all the bins, or median, or some more complex calculation like those described in the question below?

What's the best criterion for determining a frequency peak?

I want to determine the noise floor to set a threshold for determining whether or not my signal contains a given frequency.

You need to normalize your data based on the type of window you are using to obtain your frequency domain representation of the data.  The normalization differs depending on whether you are measuring a narrow band (your signal peak) or broadband (noise) signal.  Once you have the data properly normalized, the power of the narrow band signal can be read directly from the data.  The noise measurement must be estimated from the "noise floor" of the normalized frequency data.  Your noise power estimate will be 6dB less than the noise floor.  For a detailed discussion,

Go to this link: http://www.fhnw.ch/technik/ime/publikationen

Download the paper ""How to use the FFT for signal and noise simulations and measurements".


How to Remove Gaussian Noise from an Image without Destroying the Edges? What is the best filter for removing Gaussian noise without destroying the edges? I am using the standard Lena images with additive Gaussian noise and I want to denoise before applying anisotropic diffusion. I don't want to median filter because edges become blurred. I tried adaptive filtering but results were not satisfactory.

You might need to consider more advanced techniques. Here are two recent papers on edge-preserving denoising:

Edge-Preserving Image Denoising via Optimal Color Space Projection [in color] This paper preserves edges by decomposing the image into an "optimal" color space and performing wavelet shrinkage. The optimal color space belongs to the luminance/color-difference family (think L*a*b*, or YCrCb).
Edge Structure Preserving Image Denoising From the paper:

  Our method is based on [jump regression analysis], and consists of three major steps, outlined below. First, edge pixels are detected in the whole design-space by an edge detector. Second, in a neighborhood of a given pixel, a piecewise-linear curve is estimated from the detected edge pixels by a simple but efficient algorithm, to approximate the underlying edge segment in that neighborhood. Finally, observed image intensities on the same side of the
  estimated edge segment, as the given pixel, are averaged by the local linear kernel smoothing procedure (cf., [35]), for estimating the true image intensity at the given pixel.

(Jump regression models incorporate discontinuities using step functions. The primary author has a book on this subject.)


Deriving the Matrix Inversion Lemma for RLS Equations vs the Woodbury Derivation Can any one help me in deriving the matrix inversion lemma rule for RLS  algorithm?

I don't know how to start with. Many books have just stated but they haven't derived it.

It is not clear what are you asking but I will try answer both things.

Deriving the Matrix Inversion Lemma

The Matrix Inversion Lemma goes as:

$$ {\left( A + U C V \right)}^{-1} = {A}^{-1} - {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} V {A}^{-1} $$

Deriving it is by utilizing these useful identities:
$$\begin{align}
U + U C V {A}^{-1} U & = U C \left( {C}^{-1} + V {A}^{-1} U \right) = \left(A + U C V \right) A^{-1} U \label{auxEqn001}\tag{1} \\
{\left( A + U C V \right)}^{-1} U C & = {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)} ^{-1} \label{auxEqn002}\tag{2}
\end{align}$$

Where $ \ref{auxEqn002} $ comes from $ \ref{auxEqn001} $ by multiplying the 2 last terms by $ {\left(A + U C V \right)}^{-1} $ on the left and $ {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} $ on the right.

Then it is straight forward:
$$\begin{align}
{A}^{-1} &= {\left( A + U C V \right)}^{-1} \left( A + U C V \right) {A}^{-1} && \text{As $ {\left( A + U C V \right)}^{-1} \left( A + UCV \right) = I $} \\
& = \left( A + U C V \right)^{-1} \left( I + UCVA^{-1} \right) \\
& = \left( A + U C V \right)^{-1} + \left( A + U C V \right)^{-1} U C V {A}^{-1} && \text{Multiplying the terms} \\
& = \left( A + U C V \right)^{-1} + {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} V {A}^{-1} && \text{Utilizing \ref{auxEqn002}} \\
& \Rightarrow {\left( A + U C V \right)}^{-1} = {A}^{-1} - {A}^{-1} U {\left( {C}^{-1} + V {A}^{-1} U \right)}^{-1} VA^{-1}
\end{align}$$

The nice thing is we don't need the Matrix Inversion Lemma (Woodbury Matrix Identity) for the Sequential Form of the Linear Least Squares but we can do with a special case of it called Sherman Morrison Formula:

$$ {\left( A + u {v}^{T} \right)}^{-1} = {A}^{-1} - \frac{ {A}^{-1} u {v}^{T} {A}^{-1} }{ 1 + {v}^{T} {A}^{-1} u } $$

The Sherman Morrison Formula is the MIL where $ C = I $, $ U = u $ and $ V = {v}^{T} $.

Deriving the Sequential form of the Linear Least Squares Estimator

In Sequential Form of the Least Squares Estimator for Linear Least Squares Model I derived the sequential form. To add something new I will extend the derivation to the Weighted Linear Least Squares Estimator.


Impulse response of a continuous system sampled with zero-order hold I've a continuous system $$F(s) = \frac{K}{Ts+1}.$$

I sample it with zero-order hold with sampling period $T_s$. The discrete system transfer function is $$
\begin{aligned}
G(z) &=
%
\frac{z-1}{z} \mathscr{Z} \left\{
  \mathscr{L}^{-1} \left\{ \frac{F(s)}{s} \right\}
\right\}
\\
&= \frac{z-1}{z} \mathscr{Z} \left\{
  \mathscr{L}^{-1} \left\{ \frac{K}{s \left( Ts+1 \right) } \right\}
\right\}
\\
&= \frac{z-1}{z} \mathscr{Z} \left\{
K \left( 1 - e^{-\frac{1}{T} T_s n} \right)
\right\}
\\
&= \frac{z-1}{z}
K \left( \frac{z}{z-1} - \frac{z}{z - e^{-\frac{T_s}{T}}} \right)
\\
&= \frac{z-1}{z} K \frac
  {z \left( 1 - e^{-\frac{T_s}{T}} \right)}
  {(z-1) \left( z - e^{-\frac{T_s}{T}} \right)}
\\
&= K \frac
  {\left( 1 - e^{-\frac{T_s}{T}} \right)}
  {\left( z - e^{-\frac{T_s}{T}} \right)}.
\end{aligned}
$$

From this I can get difference equation of the sampled system: $$
y[n] = e^{-\frac{T_s}{T}} y[n-1] + K \left( 1 - e^{-\frac{T_s}{T}} \right) x[n-1].
$$

I'd like to get  impulse response function of the system in time domain $h[n]$: $$
y[n] = h[n] * x[n].
$$

To get $h[n]$, I'd have to find inverse z-transform of $G(z)$, but I can't figure out what it is. I can't find it in any tables. What is it?

Thank you.

People in ##dsp on Freenode helped me realize that the impuse response is a combination of impulse responses of two simple systems connected in series, that is, convolution of two impulse responses in the time domain.

$$
\begin{aligned}
G(z) &=
K \frac
  {\left( 1 - e^{-\frac{T_s}{T}} \right)}
  {\left( z - e^{-\frac{T_s}{T}} \right)}.
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
z^{-1}
\frac{1}{\left( 1 - e^{-\frac{T_s}{T}} z^{-1}\right)}
\end{aligned}
$$

$$
\begin{aligned}
h[n] &= \mathscr{Z}^{-1} \left\{ G(z) \right\}
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
\mathscr{Z}^{-1} \left\{
z^{-1}
\frac{1}{\left( 1 - e^{-\frac{T_s}{T}} z^{-1}\right)}
\right\}
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
\mathscr{Z}^{-1} \left\{
z^{-1}
\right\}
*
\mathscr{Z}^{-1} \left\{
\frac{1}{\left( 1 - e^{-\frac{T_s}{T}} z^{-1}\right)}
\right\}
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
\delta[n-1] * \left( e^{-\frac{T_s}{T}n} u[n] \right)
\\ &=
K \left( 1 - e^{-\frac{T_s}{T}} \right)
e^{-\frac{T_s}{T}(n-1)} u[n-1]
\end{aligned}
$$

The impulse response is exponential decay shifted to the right by 1.


How to apply the provided wave transforms on an image I want to implement one of the following mappings in MATLAB and verify its performance on an image, but I couldn't. Please help me with the implementation

You want to apply to your original image $I_1[n,m]$ a transform (x,y)=T(n,m) such as a warping effect, to obtain a new image $I_2[n,m]$.

In your example, you want to shift original pixel positions (n,m) to new coordinates based on the calculated position (x,y) from the transform equation (x,y) = T(n,m). As usually the case, this direct transform of indices will possibly result in noninteger positions (x,y) and therefore be quantized to slightly distorted locations, or even worse some pixels of the final image may not be filled, leaving them empty.

One approach to solve these two problems is to apply the algorithm from the reverse, by using the inverse transform of T(,)  provided it exists. As in the pseudo matlab code here: Result is as shown below

As a specific case of your examples:
WAVE1: T(n,m) = (n + 20*sin(2*pi*m/128) , m ) then the inverse transform is
       invT(x,y) =( x - 20*sin(2*pi*y/128) , y) which yields m = y and n =  x - 20*sin(2*pi*y/128)

% NOTE1: This code is not efficient. You shall Vectorize it easily
% and get rid of the for-loops, yielding a substantial boost in efficiency.
% NOTE2: boundary rows/columns may be left untouched which needs a simple fix
% The transform that you want to appply is (x,y)=T(n,m)
% The transform that we will use is its inverse (n,m)=invT(x,y) (you shall invert your transform)
% Original image I1 is of size S1xS2, transformed image I2 is of size N1xN2
for x=1:N1    % let N1 = S1 and N2= S2 for this example WAVE1
  for y=1:N2

     % s1- compute back coordinates in the original image domain
     % [n m] = invT(x,y);
     m = y;
     n = x - 20*sin(2*pi*y/128);

     % s2-s3 Apply (bilinear) interpolation to find the image intensity, if
     % n,m does not fall into integer indices
     d1 = n - floor(n);
     d2 = m - floor(m);

     % s3- define the image pixel
    if ( ((n>=1)&(n<S(1))) & ((m>=1)&(m<S(2))) )
     I2(x,y) = (1-d1)*(1-d2) * I1(floor(n),floor(m)) + ...
               (1-d1)*d2 * I1(floor(n),floor(m)+1) + ...
               d1*(1-d2) * I1(floor(n)+1,floor(m)) + ...
               d1*d2 * I1(floor(n)+1,floor(m)+1) ;
    end%if

  end
end


Confusion in AC and DC signals I have certain doubts in the fundamentals of AC and DC signals.

Most of the time, we deal with DC signals in the form of batteries, whose waveform is a constant signal with peak as 1.5v

However, i came across some resources, that told certain other signals are DC as well. For your reference, i have enclosed the resources and the respective waveforms here.

Figure 1

Figure 2, Please ignore the "bidirectional" word on the image.

So, are these signals AC, or DC? Both the sources have mentioned that they are DC signals (well, the first one is apparently "DC Pulsating" and the second one is a "DC time-varying". This confused me, and i wanted some better clarification and knowledge on these. From what i observe,

1) Any signal that is above or below the time axis (i.e the same polarity all the time) can be called DC?

2) If i have an AC signal, say a sinusoid, and by appropriately DC shifting it above / below the time axis, can i make it a DC signal, if the above point is true?

I would seek your expertise and your insights on these images as well, and tell me why these signals are DC.

Thank you.

Reference 1: https://www.elprocus.com/variable-power-supply-digital-control/

Reference 2: https://www.youtube.com/watch?v=fXBtTI107fQ, time: 2:31

The signal in question is a rectified sine-wave. This signal has a non-zero average, i.e. a DC component. It also has AC components at f, 3f, 5f, etc. where f is the frequency of the non-rectified signal.


What's the Difference between Convolution Kernel and Point Spread Function in the Context of Image Convolution? When you use the deconvolution method to make the blurry image sharper, you will have to estimate the Point Spread Function. Is there a difference between this PSF and an image kernel?

Second question: Is the Point Spread Function in the time domain or frequency domain? Or can it be in both?

The PSF (Point Spread Function) is the system response to Impulse Signal (Point).

If your system model is LSI (Linear Spatially Invariant) then the output image of the system is applying the PSF as a convolution kernel on the input image.

Yet, the PSF is just the response of a system to a certain input.
It doesn't describe the whole system unless it is an LSI system.

If you model for the Deconvolution problem (Often called the Degradation Model) is LSI, then what you estimate is indeed the PSF.


Image Matrix Vector Representation for the Degradation Model I am trying to understand the the degradation model equation but I have doubt that how come y^t.x.h will be equal to  x^t.h^t.y . Aren't they transpose of each other.

Indeed since both expressions are scalars then they are equal to each other since the transpose of a scalar is the same scalar.
See in MATLAB as an example (Calculating $ {x}^{T} H y $ and $ {y}^{T} {H}^{T} x $:
>> vX = randn(10, 1);
>> vY = randn(10, 1);
>> mH = randn(10, 10);
>> vX.' * mH.' * vY

ans = -0.8618

>> vY.' * mH * vX

ans = -0.8618

As you can see, indeed both expressions are scalars.
As expected $ {x}^{T} H y = {\left( {x}^{T} H y \right)}^{T} = {y}^{T} {H}^{T} x $.


Can I search Doppler shift in freqency domain I would like to know if it is possible  to search Doppler shift of a carrier with sufficient oversampling while performing code search in the freqency domain. For instance can I selectively drop every X time domain sample bin from the fft processing? Will that not effectively change the input frequency being compared to the locally generated code?
I would like to perform the Doppler shift search of digital satellite broadcasts without having to change the lo supplying the IF in hardware.
Can this be done this way?

Note my function below showing how I simultaneously searched for code and carrier (Doppler) offset for spread spectrum signals such as GPS in one 2D IFFT computation*. The 2D IFFT is processing intensive, but this approach has a lot of merit if you need very fast ("one-shot") acquisition.    (*Note the total processing involves two 1D FFT's, an element by element matrix multiplication for all Doppler and Code bins, and a 2D IFFT).

The basic premise is using fft's to cross-correlate two sequences (your waveform that has code and carrier offsets, and your reference code):

$\textrm{XCORR} = \textrm{ifft}(\textrm{fft}(a)\textrm{fft}(b)^*)$

Specifically the circular cross correlation between two sequences is equal to the inverse fft of the complex conjugate multiplication of their respective fft's.

If you circularly shift one of the fft's before multiplying, you calculate the result for a different Doppler bin!  It is that easy.  For actual implementation detail I included my Matlab function below. Here is an example plot of the result:

One FFT is computed for the received sequence, and one FFT is computed for the reference code. A matrix of delayed FFT's is created for one of the two that only involves replicating columns with a circular shift of the FFT already computed (which is a Hankel matrix) , and then a matrix of the other is created by simply replicating the FFT computed in each column, matching the number of Doppler columns created. The two matrices are complex conjugate multiplied element by element ($M_{ij}= A_{ij}B_{ij}^*$) Then (this is where all the processing is needed), a 2D inverse FFT is computed of the product matrix resulting in the solution depicted in the surface plot above.

On one axis for the figure above is the normalized frequency offset (Doppler) and the other axis is the relative code delay, and the vertical axis is the correlation where the delay and Doppler were simultaneously determined. So with one capture (whose length depends only on post correlation SNR required of the received sequence) you can complete acquisition processing with the use of two ffts, an element by element matrix multiplication and a 2D IFFT.

function [fdout,freq_axis]= fdcorr(a,b,frange)
%FDCORR Joint frequency and delay correlation
%   Performs an efficient circular cross correlation between
%   two vectors over all possible frequency offsets (Dopplers)
%   between -pi and pi where 2pi is the sampling rate.
%
%[FDOUT, FREQ_AXIS] = FDCORR(A,B)
%
%   A:   vector representing received signal sampled at FS = 2pi
%   B:   vector representing correlation code sampled at FS = 2pi
%
%   FDOUT: 2D Correlation surface vs delay and frequency offset
%   FREQ_AXIS: Frequency values aligned to frequency axis (-pi to pi)
%
%[FDOUT, FREQ_AXIS] = FDCORR(A,B,FRANGE)
%
%   Limits frequency range as defined in FRANGE.    Typically
%   frequency offsets are much smaller than the Nyquist bandwidth, so
%   using a reduced range for FRANGE from -pi to pi will be more
%   efficient.
%
%  FRANGE: Two element vector from lowest to highest frequency to
%          search over within a frequency range from -pi to pi.
%          Example [-pi/10 pi/10]
%
%FDCORR(A,B,<FRANGE>) with no output arguments plots the surface
%       correlation magnitude vs frequency and delay offset in the
%       current figure window.  FRANGE optional.
%
%
% Dan Boschen 11/22/2009

%==========================================================
% VERSION HISTORY
%
% 11/22/2009     Inital Release
%
%
%==========================================================

%==========================================================
% APPLICATION NOTES:
%
% Vector "A" represents a received signal that has been encoded with
% a correlation code (such as GPS and CDMA applications) or a training
% sequence.  It has an unknown frequency offset due to Doppler offsets
% and frequency offsets in the receiver, and an unkown delay shift
% relative to the start of the code sequence represented by vector "B".
% FDCORR will efficiently find the correlation between the received
% sequence and reference code over frequency and delay offsets.
% Vectors A and B are interchangable.
%
% To detect the maximum correlation value and delay and freq location use:
% [temp1, temp2]=max(fdout);
% [maxcorr,freq_index]=max(temp1);
% delay =temp2(freq_index)-1;
% freq= freq_axis(freq_index);
%
% The basic operation uses the Cross-Correlation theorem:
% xcorr= ifft(fft(a).*conj(fft(b)))
%
% Note the detected frequency offset is course and is the
% frequency of the closest Doppler bin where each Doppler bin
% is seperated by 2pi/N
%
%==========================================================

%
%error handling
if nargin==3
    if length(frange)~= 2
        error('FRANGE must be a two element vector')
    elseif abs(frange)> pi
        error('FRANGE values must be between -pi and pi')
    end
end

%force column vectors:
a=a(:);
b=b(:);

%============================================================
%create all Doppler versions of a by shifting the fft:
%============================================================

N=max(length(a),length(b));
fft_a= fft(a,N);

%frequency axis for all doppler bins -pi to pi

%index is the "fft shifted" freq axis. Method below faster than
%using fftshift(0:N-1), but same result.
freq_index=[-floor((N-1)/2):floor(N/2)];
freq_axis= freq_index*2*pi/N;   %convert to -pi to pi range

if nargin==3     %reduced doppler range
    freq_include=find(freq_axis>=frange(1) & freq_axis<=frange(2));
    freq_axis=freq_axis(freq_include);
    freq_index=freq_index(freq_include);
    num_freqs=length(freq_index);
else
    num_freqs=N;
end

%create Hankel matrix for all Dopplers (each column is the fft
%of vector a at a different Doppler offset)
fft_idx=(1:N)';
fft_dopplers = fft_idx(:,ones(num_freqs,1))+freq_index(ones(N,1),:);
fft_dopplers= mod(fft_dopplers-1,N)+1;    %Hankel subscripts
fft_dopplers(:)=  fft_a(fft_dopplers);      %fill data

%alternate method if solving for all Doppler bins (processing intensive):
%fft_all_dopplers= hankel(fft_a,[fft_a(end); fft_a(2:end)]);

%==================================================================
%Circular Correlation using circ_corr= ifft(fft(a).*conj(fft(b))
%==================================================================
fft_code_mtx= fft(b,N)*ones(1,num_freqs);
fdout=(ifft(fft_code_mtx.*conj(fft_dopplers)));

%===============================================================
%plot result
%===============================================================
if nargout ==0
    surf(freq_axis,1:N,(abs(fdout)));
    shading interp;
    xlabel('Normalized Frequency (xpi rad/sample)');
    ylabel('Delay Index');
end


Conversion between DTFT in radians/sample to DTFT in cycles/sample I have found that most commonly the DTFT is defined as:
$X(\omega) = \sum_{n=-\infty}^{\infty} x[n]e^{-j \omega n}$.

However the class I am taking frequently uses the DTFT expressed in "normalized continuous frequency of discrete signals" units i.e. cycles/sample. I am trying to understand how to convert a DTFT pair using one units to other units.

For example, the DTFT pair with $\omega$ [radians/sample]:

$$x[n] = 1  \Leftrightarrow X(\omega) = 2 \pi \sum_{k=-\infty}^{\infty} \delta( \omega-2 \pi k)$$

My understanding that this pair with $\mu$ [cycles/sample] is given by:
$$x[n] = 1  \Leftrightarrow X(\mu)= \sum_{k=-\infty}^{\infty} \delta( \mu -k)$$

I have honestly struggled to find many resources which make any mention of the latter DTFT.

You can use the following :

Consider the discrete time sequence $x[n]$ and the DTFT $X(\omega)$ related by:
$$ X(\omega) = \sum_k x[n] e^{-j \omega n} \longleftrightarrow x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(\omega) e^{j\omega n} d\omega$$

Then define $X'(\mu)$ as the normalized DTFT, which is related to $X(\omega)$ as
$$X'(\mu) = X(2\pi \mu) = \sum_k x[n] e^{-j 2\pi \mu n} $$ and then
$$x[n] = \int_{-0.5}^{0.5} X'(\mu) e^{j 2\pi \mu n} d\mu $$

Applying this to your example case:
$$x[n] = 1 \longleftrightarrow X(\omega) = 2\pi \sum_r \delta(\omega - 2\pi r)$$ then since $X'(\mu) = X(2\pi \mu)$ we have

$$X'(\mu) = 2\pi \sum_r \delta(2\pi \mu - 2\pi r) = 2\pi \sum_r \delta(2\pi [\mu - r]) = 2\pi \sum_r \frac{1}{2\pi} \delta(\mu - r) $$
$$X'(\mu) = \sum_r \delta(\mu - r) $$

as expected, where we've used the property that $\delta(ax) = \frac{1}{|a|} \delta(x)$.


How to convert Third Octave data into single SPL quantity I have a data vector which consists of calibrated pressure at each sampled point in time. I have also taken this signal, applied a filter bank to it using a butterworth filter with passbands according to third octave (i.e. band center frequency at 25, 31.5, 40,...,20,000), and for each band calculated $20\log(p_{rms}/p_{ref})$ where $p_{rms}$ is the root mean squared pressure for the resulting filtered signal in each band, and $p_{ref} =2\times10^{-5}$ Pa. This should give me the third octave "spectrum".

I now wish to calculate the total sound pressure level (SPL). I thought first I could say that $SPL = 20\log(p/p_{ref})$ where $p$ is just the root mean squared value of the entire data vector.

To compare I thought I would then calculate the same quantity but this time starting from my third octave data. I would then hope to get the same answer by adding logarithmically (i.e. adding the rms values from each band and then converting to dB) the dB values from each third octave band. However this did not give matching results.

I therefore suspect I am not converting from third octave to a single SPL quantity correctly, since my first calculation of SPL seems very straightforward. Am I missing something here? For instance does each third octave band need to be scaled by the frequency bin? This would imply I have some sort of density in third octave, but im not sure.

Make sure that you are trying to compare same quantities, i.e. not $SEL$ with $L_{eq,C}$.
If third-octave levels were not measured with any kind of weighting and you are trying to calculate $L_{eq,A}$ for example, then you must apply corrections.
For each band calculate the anti-logarithm using the formula: $p_i=10^{L_i/10}$ (for example for $54.3 \mathrm{dB}$ that would be $218776.2$).
Sum all values you obtained in previous step and calculate the total level as: $10\log_{10}\left(\sum p_i\right)$.

That will produce exact result - tested second ago on some data from SLM.


Reverse engineering a digital filter I have 2 simultaneous signals that both are designed to measure eye movements.
They are sampled at 250 Hz.
We have 12 subject recordings.
For 3 stable periods in each subject, we choose 256 points and did an FFT.
Prior to the FFT, the data were mean-centered and detrended with a 2nd
order polynomial.  They were also windowed with a Hann window.
We are focused on the magnitude spectra plots.
We have a total of 12 X 3 = 36 magnitude spectra, which we average.
These averages are shown in the attached figure.

My hypothesis is that the second signal is a low-pass filtered version of
the first.  The filter has ringing in the passband.

I want to design a digital filter that I would apply to signal 1 that,
after fft, analyses would produce a magnitude spectra like that of
signal 2.

How would I go about this?

Signal 1 was an average of a left and right signal and signal 2 was a binocular signal.

I followed the suggestion of MBaz, and computed the frequency response of the filter for each segment (N = 3) for each subject (N=12).  Then I averaged the frequency responses of the filters. Here is the result:

MBaz and Matt: Thank you so much for solving my problem.
Matt, it looks close to your filter.
Lee


Bandwidth of lowpass filter? This article here on the Window method for FIR filter design speaks of

  the total normalized bandwidth of the lowpass filter in $\textrm{Hz}$ (counting both negative and positive frequencies)

What is it?

If you consider an ideal low-pass filter with cut-off frequency of $f_c$, all frequencies greater than $f_c$ will be removed. Then it's bandwidth is equal to $f_c\;\mathrm{Hz}$ (from $0$ up to $f_c$).

The total bandwidth $B_T$ is simply twice that: $B_T=2f_c$, since we are also considering negative frequencies, from $-f_c$ up to $f_c$.

Here is a picture of that:


FIR Filter Type - MATLAB A question about using fir1/fir2/firpm to design FIR filters in MATLAB.  The first parameter is the 'filter order'. For my filter design I pass in 250 as this argument.  However, when looking at the impulse response of the filter, it is a type 1 FIR filter.  So my question is, is the filter order not the same as the number of taps? Is the number of taps one less, so 249 in this case, which gives the type 1 symmetry?

Cheers

The filter order equals the number of taps minus $1$. The filter order is the order of the polynomial corresponding to the filter's transfer function. E.g.,

$$H(z)=h[0]+h[1]z^{-1}+h[2]z^{-2}\tag{1}$$

has $3$ coefficients (taps) but it is a second order polynomial (in $z^{-1}$), so the filter order equals $2$ (i.e., it has $2$ zeros and $2$ poles).

For type I and type III linear phase FIR filters, the filter order is even, i.e., they have an odd number of taps, whereas types II and IV have an odd filter order (even number of taps).


Calculating the Fourier transform of shifted scaled unit step function I have $x_1(t)$ here. To get $x_2(t)$, I need to differentiate $x_1(t)$. Express $x_2(t)$ as $2u(t+2)-4u(t)+2u(t-2)$.
From Fourier transform definition integral, I got $X_2(j\omega)=\frac{2e^{j\omega 2}}{j\omega}-\frac{4\pi \delta(\omega)}{j\omega}+\frac{2e^{-j\omega2}}{j\omega}$.
Is this correct? It seems weird and would be complicated to calculate the magnitude and phase spectrum. Thanks!

It looks like the idea of that exercise is to compute the Fourier transform of $x_3(t)$, and then from that derive the Fourier transform of the original function $x_1(t)$.

Note that $x_3(t)$ is just a sum of scaled and shifted Dirac impulses, so computing its Fourier transform is trivial. The Fourier transform of $x_2(t)$ is then easily derived by multiplying $X_3(j\omega)$ by $\pi\delta(\omega)+1/j\omega$, and noting that the Dirac delta is cancelled because $X_3(0)=0$. So you just end up dividing by $j\omega$. In the same way, $X_1(j\omega)$ is obtained from $X_2(j\omega)$. Also note that none of the three Fourier transforms $X_1(j\omega)$, $X_2(j\omega)$, and $X_3(j\omega)$ have any Dirac impulses.

It's a good idea to cross-check your result by realizing that $x_1(t)$ can also be represented as the convolution of two rectangular functions, i.e., its Fourier transform must be a squared sinc-function.


Non-zero DFT components where zero is expected? I am implementing DFT in Octave. Here's my code:
function [real_comp, imag_comp] = mydft (samples)
    N = columns(samples);
    for m = 1:N
        real_accum = 0;
        imag_accum = 0;
        for n = 1:N
            real_factor = cos((2 * pi * (n-1) * (m-1)) / N);
            # imag_factor = -I * sin((2 * pi * (n-1) * (m-1)) / N);
            imag_factor = sin((2 * pi * (n-1) * (m-1)) / N);
            real_accum += real_factor * samples(n);
            imag_accum += imag_factor * samples(n);
        endfor
        real_comp(m) = real_acuum;
        imag_comp(m) = imag_accum;
    endfor
endfunction

My goal is to compute the magnitude at each DFT bin and plot it. I am passing in a list of real-only samples and storing the real and imaginary components in different variables. Here's the sample vector that I am using:
samples = [0.3535, 0.3535, 0.6464, 1.0607, 0.3535, -1.0607, -1.3535, -0.3535]

The samples are carefully chosen to prevent spectral leakage at fs = 8KHz. The graph of the DFT magnitude is exactly what I expect. There's nothing wrong with that. But what I don't get is that some components in my 8 point DFT which are supposed be zero, are not exactly zero but rather a very small, close-to-zero value. The built-in fft() function, on the other hand, has zeroed out components where I expect.
For example, I print the elements of output vector produced by my DFT and octave's FFt using the code below:
fft_out = fft(samples);
[real, imag] = mydft(samples);
dft_out = real + (-I * imag);
for i=1:columns(samples)
    disp(fft_out(i));
    disp(dft_out(i));
    printf("\n");
endfor

The point-to-point comparison is given in the table below:

bin_index
fft_out
dft_out

1
-1.0000e-04
-1.0000e-04

2
0.00000 - 3.99988i
4.7184e-16 - 3.9999e+00i

3
1.4141 + 1.4144i
1.4141 + 1.4144i

4
0.0000e+00 - 8.0820e-05i
1.3878e-16 - 8.0820e-05i

5
-1.0000e-04
-1.0000e-04 - 1.4350e-16i

6
0.0000e+00 + 8.0820e-05i
2.3037e-15 + 8.0820e-05i

7
1.4141 - 1.4144i
1.4141 - 1.4144i

8
0.00000 + 3.99988i
-2.0817e-15 + 3.9999e+00i

Now, I wouldn't care too much about the tiny difference but if I compute the phase angles of fft_out and dft_out,
disp(arg(fft_out))
disp(arg(dft_out))

I get the sign flipped in some cases:
fft_out: 3.14159  -1.57080   0.78550  -1.57080   3.14159   1.57080  -0.78550   1.57080
dft_out:-3.14159  -1.57080   0.78550  -1.57080  -3.14159   1.57080  -0.78550   1.57080

So, is there something wrong with my DFT implementation? Is it just a rounding or precision issue?

These are just rounding errors, and you only get a flipped sign for phase values of $\pm\pi$, which is not an issue because the phase value is ambiguous, i.e., adding or subtracting multiples of $2\pi$ doesn't change the value of the complex coefficient.
I assume that your code is just for an educational purpose, because it is neither efficient nor numerically accurate.


Upper bound for realtime application I have developped an application for image inpainting. Now I want to claim that it can be used for realtime purposes.

Does somebody know a paper which gives some boundaries for realtime applications? I want to refer to a paper to prove my assumption.

Real time means that your algorithm can process the video faster than the time it takes to watch the video. In other words, if you were collecting the video live you would be able to finish processing it as soon as you stopped collecting it with no extra work needed to be done after.


ValueError: operands could not be broadcast together with shapes (194,259,3) (600,800,3) I got ValueError: operands could not be broadcast together with shapes (194,259,3) (600,800,3) with the below code:
import numpy
import math
import cv2
original = cv2.imread("original.png")
contrast = cv2.imread("compressed.png",1)
def psnr(img1, img2):
    mse = numpy.mean( (img1 - img2) ** 2 )
    if mse == 0:
        return 100
    PIXEL_MAX = 255.0
    return 20 * math.log10(PIXEL_MAX / math.sqrt(mse))

d=psnr(original,contrast)
print(d)

What was the problem with the implementation?

It sounds like the error is indicating that the two images aren't the same size.  So you can't subtract them point for point.


How to treat noise (in the acceleration) in a Kalman filter when tracking position? I have an equation of motion for which I am wanting to track the position.

The equation looks like this:

$$\ddot{x}(t) = -a\dot{x}(t) - bx(t) + F(t)$$

where $F(t)$ is an external Stochastic force noise (like brownian motion).

I know if $F = 0$ I can build a Kalman filter like so:

$$
\begin{bmatrix}
\dot{x} \\
\ddot{x} \\
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 \\
-b & -a \\
\end{bmatrix}
.
\begin{bmatrix}
x \\
\dot{x} \\
\end{bmatrix}
$$

And then transform $ \begin{bmatrix}
0 & 1 \\
-b & -a \\
\end{bmatrix} $ like so:

\begin{align}
A(t) & = \mathscr{L}^{-1}\left\{\left(s\vec{I} - \begin{bmatrix}
0 & 1 \\
-b & -a \\ \end{bmatrix}\right)^{-1}\right\}\\
& = \mathscr{L}^{-1}\left\{\left(\begin{bmatrix}
s & 0 \\
0 & s \\ \end{bmatrix} - \begin{bmatrix}
0 & 1 \\
-b & -a \\ \end{bmatrix}\right)^{-1}\right\}\\
&= \mathscr{L}^{-1}\left\{\left(\begin{bmatrix}
s & -1 \\
b & s+a \\ \end{bmatrix}\right)^{-1}\right\}
\end{align}

... (at this point I get stuck as I am using Sympy in python to calculate the inverse Laplace transform and it fails for this matrix).

But I want to know how I should treat this Stochastic effect in a Kalman filter.

Uh, this is the whole point of a Kálmán filter.  It has a stochastic measurement model (noisy measurements) and a stochastic process state update model and uses the measurements for estimating the current state of the stochastic process.

If the process were deterministic, you would approach perfect information over time.

So you just plug the variation of your acceleration into the covariance matrix for the state update.

The exact way to do that depends on the conventions your course material employs (lots to choose from).


Why are there comb-like hills in a periodogram? I am playing with the periodogram of MATLAB. I created a simple script to observe how it behaves:

rng(1);  %# initialize the random number generator

Fs = 1000;  %# Sampling frequency
duration = 0.1; %# seconds

A = 1; %# Sinusoid amplitude
f = 150; %# Sinusoid frequency
eps = 0.01;

t = 0:1/Fs:duration;
x = A * sin(2*pi*f*t) + eps * randn(size(t));

periodogram(x,[],1024,Fs);

I have no problem with the code and can write my own periodogram function using the algorithms given in the documentation but I wonder the theoretical reason behind the comb-like hills which are not 150 Hz. What do I get those instead of getting a single spike over 150 Hz? Is there anything special in the distances of the peaks of these hills?

A single spike (as you call it) appears theoretically only for a infinite-length sinusoid. Since your signal is 100 samples length, it is not infinite. You actually multiplied your infinite signal with a window which has a value of 1 over 100 samples, and 0 elsewhere. Since multiplication in time domain is equivalent to convolution in the frequency domain, your spectrum is a convolution of the single spike and the frequency response of the window (btw it is called rectangular window). This is the function you got.

I suggest you read about windowing: http://en.wikipedia.org/wiki/Window_function


2D Fourier transform of Sobel kernel Can someone explain me the highlighted text parts regarding this image ?

Here is a pseudo-code of how it was created:

kSobelX = zeros(100);
kSobelX(1:3,1:3) = [1,0,-1;2,0,-2;1,0,-1];
kSobelXfft = fft2(kSobelX);
clf; scale equal;
image(abs(kSobelXfft)/max(abs(kSobelXfft(:))));
print;

Pixels on the left of the square correspond to low frequencies along the x axis. Because the discrete Fourier transform is discrete in space, it's periodic in frequency: pixels on the right of the square correspond to low negative frequencies along the x axis. And because the spatial data is real, the frequency data is symmetric. The middle of the square correspond to the Nyquist frequency (half the sampling frequency). The regions at 1/4 and 3/4 positions correspond to the highest frequencies the filter operates at. Along the x axis, the light region means the Sobel filter is a highpass filter; along the y axis, the light regions are at low frequencies, i.e. the Sobel filter is a low-pass filter.

Because the discrete Fourier transform is discrete in space, it's periodic in frequency

Let's think of the Fourier Transform of a continuous, periodic function: FT explains the function as a projection of the signal onto a sum of sinusoidal basis functions, i.e. correlates the signal with a bunch of complex sinusoids, each with its own frequency. Thus, any data point contributes to a sinusoid with a particular frequency - data is encoded into the frequency. Thus FT becomes discrete. Now inverse also holds: If the space is discrete, one would require a periodic FT to represent that (non periodic signal on one domain causes continuous signal on the other).

I really like this animation showing the continuous FT.

Mathematically, it is always possible to show that for DFT, $X_K = X_{K+N}$ where $N$ is the period such as in here.

And because the spatial data is real, the frequency data is symmetric. The middle of the square correspond to the Nyquist frequency (half the sampling frequency).

First let's address the symmetry:

I will greatly cite this stackexchange post and for completeness will include some shortened portion here:

Due to its definition, FT will generate symmetric responses for real data. Take the transform equation:
$H(f) = \int h(t)e^{-j2\pi ft}dt$.

Using Euler's formula $e^{jx} = \cos(x) + j*\sin(x)$, we could modify the Fourier transform as follows:
$$
H(f) = \int h(t)(\cos(2\pi ft) - j\sin(2\pi ft))dt
$$

FT for negative frequencies become:
$$
H(-f) = \int h(t)(\cos(2\pi (-f)t) - j\sin(2\pi (-f)t))dt \\
= \int h(t)(\cos(2\pi ft) + j\sin(2\pi ft))dt
$$

Comparing the negative frequency version with the positive frequency version shows that the cosine is the same while the sine has phase inversion of $90^\circ$. Magnitude is the same. Due to this orthogonality, they will both respond to real signals in the same way.

Now for the Nyquist part:

The Nyquist frequency represents the highest frequency that can be represented by the sampling rate (differs from Nyquist rate). The DFT of a real series, i.e. zero imaginary part, results in a symmetric series about the Nyquist frequency - this is where the signal folds. In the image you have given, the mid-part is exactly about where the symmetry is created. This is also where the DC component lies.

I hope these bring in some clarification.


Difference ;Phase and frequency offset can someone explain what a difference is between phase and frequency shift/offset in estimation algorithm?
Honestly I dont understand this difference. If i estimate , i will do it for phase estimation.
EDIT 1 What is a difference in an implementation of FLL and PLL?

The difference in implementation between a FLL and PLL is completely dictated by what device is used as the error detector: if the error detector produces an error that is directly proportional to phase (phase detector) that is then passed into the loop filter, then the implementation will be a PLL and we would treat the VCO (Or NCO) as a phase integrator (since that device producing a frequency in proportion to its input, so produces a phase in proportion to the integral of the input). If the error detector is directly proportional to frequency (such as a frequency discriminator) that is then passed into the loop filter, then that device would be an FLL and we would treat the VCO or NCO as a proportional device producing a frequency in direct proportion to its input. Note that we can have a phase detector that we then take the derivative of its output which would then be a frequency discriminator (and this approach can be used for a combined FLL/PLL) so the actual block diagram and functionality needs to be carefully reviewed for each case.
Frequency is simply the derivative of phase with respect to time. So in a phase estimation algorithm you estimate absolute phase (relative to a reference which is typically the phase of a symbol location on the complex IQ plane), while in a frequency estimation algorithm you estimate the rate of change of phase with respect to frequency which is independent of any offset.
Picture a wheel that is spinning at a constant frequency and we use a control loop to stop the spin; if we measure error using a frequency estimation algorithm to then reverse and ultimately stop the spin, it will eventually stop but at an arbitrary offset like we would typically expect from a roulette wheel. If we use a phase estimation algorithm to measure error which is then used to stop the spin (and assuming we can keep up with and overcome the rate of change which implies sufficient loop order) it will stop at the exact place you desire (and you’ll be a winner every time!)