Frequency Analysis (DFT / FFT) of a Signal Without a Constant Sampling Frequency (Non Uniform Sampling in Time Domain) I'm a stack exchange user for some time and now I'm registering to ask a simple question (I think!).

I have a vibration signal with an amplitude and time (sampling frequency not constant) in  a $10000\times 2$ double variable.

The data is available at: https://1drv.ms/x/s!AoCOij4si31tzgY89bhr6XH-_gQq

How can I do some sort of frequency analysis (FFT, DFT) or similar. How can I do it in Matlab?

Sorry if the question is duplicated but I couldn't find any answer for my problem.

The DFT Matrix for Non Uniform Time Samples Series

Problem Statement

We have a signal $ x \left( t \right) $ defined on the interval $ \left[ {T}_{1}, {T}_{2} \right] $.
Assume we have $ N $ samples of it given by $ \left\{ x \left( {t}_{i} \right) \right\}_{i = 0}^{N - 1} $. The samples time $ {t}_{i} $ is arbitrary and not necessarily uniform.

We're after the DFT of the samples $ \left\{ X \left[ k \right] \right\}_{k = 0}^{K - 1} $ as it was samples in a uniform manner (Implicitly means the samples in Frequency domain will be uniform as well).

Deriving the Connection

In the DFT Transform the connection between time and frequency is given by:

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

In $ \eqref{EqnIdft} $ we use $ n $ for modeling the sample index in time. We usually build samples in time as $ x \left[ n \right] = x \left( n {T}_{s} \right) $ where $ {T}_{s} $ is a uniform sampling interval.
Hence we could write:

$$ x \left( n {T}_{s} \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} n {T}_{s}} \tag{2} \label{EqnIdft2} $$

In $ \eqref{EqnIdft2} $ we added explicit scaling of time. This is a known property of Fourier transform family which scales the domain in order to normalize the transform.

Now, there is nothing which blocks us from using arbitrary time:

$$\begin{align*} \tag{3} \label{EqnIdft3}
x \left( t \right) & = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} t} && \text{} \\
& = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k {F}_{s}}{N} t} && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$

As can be seen $ \eqref{EqnIdft3} $ makes sense as it goes through each element according to its frequency and sums to give the output at time $ t $. We can go step farther and generalize it for cases we don't have uniform sampling frequency.
The average sampling frequency is given by $ \bar{F}_{s} = \frac{N}{ {T}_{2} - {T}_{1} } $. Let's define $ T = {T}_{2} - {T}_{1} $ and we'll get:

$$ x \left( t \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$

Which is many ways resembles the DTFT Transform equation which does the same in the other direction, transforming uniform discrete samples in time domain to arbitrary frequency (Within a frequency interval) in Frequency Domain:

$$\begin{align*} \tag{4}
X \left( f \right) & = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi f {T}_{s} n } && \text{} \\
& = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi \frac{f}{ {F}_{s} } n } && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$

We see the same scaling, $ \frac{f}{ {F}_{s} } $ which scales the continuous $ f $ relative to the interval of frequencies $ {F}_{s} $ which is equivalent to $ \frac{t}{ T } $ which scales $ t $ relative to the time interval of the continuous signal.

The Transform Matrix

So, given the set of time indices $ {\left\{ {t}_{i} \right\}}_{i = 0}^{N - 1} $ the transformation matrix, from frequency domain to time domain, is given by:

$$ D \in \mathbb{R}^{N \times K}, \;  {D}_{i, k} = {e}^{ j 2 \pi k \frac{ {t}_{i} }{T} } $$

The Model

In vector form the model is:

$$ x = D y $$

Where $ y \in \mathbb{C}^{K} $ is the vector of the frequency coefficients in uniform grid, $ x $ is the samples in time (Non Uniform, Or at least no assumption of uniformity) and $ D $ as defined above.
Since in our model we're after $ y $ the answer is given by:

$$ y = {D}^{\dagger} x $$

Where $ {D}^{\dagger} $ is the Pseudo Inverse Matrix of $ D $.

Implementation & Results

The code is as following:

subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

samplingFrequency = 101; %<! [Hz]
samplingInterval = 1 / samplingFrequency; %<! [Sec]
startTime = 1; %<! [Sec]
endTime = 4; %<! [Sec]
timeInterval = endTime - startTime; %<! [Sec]

numSamples = round(samplingFrequency * timeInterval);
numSamplesTT = round(1.2 * numSamples);

signalFreq = 2; %!< [Hz]

% The uniform time grid
vT      = linspace(startTime, endTime, numSamples + 1);
vT(end) = [];
vT      = vT(:);

% The non uniform time grid - Reconstruction
vTT = endTime * rand(numSamplesTT, 1);
vTT = sort(vTT, 'ascend');

% The non uniform time grid - DFT
vTD = linspace(startTime, endTime, (10 * numSamples) + 1);
vTD(end) = [];
vTD = vTD(sort(randperm(length(vTD), numSamples)));
vTD = vTD(:);

% The uniform frequency grid
vF      = (samplingFrequency / 2) * linspace(-1, 1, numSamples + 1);
vF(end) = [];
vF      = vF(:);

vK = [-floor(numSamples / 2):floor((numSamples - 1) / 2)];
vK = vK(:);

%% Generate Data

vX  = cos(2 * pi * signalFreq * vT);
vFx = fftshift(fft(vX));

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = subplot(1, 2, 1);
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);

hAxes           = subplot(1, 2, 2);
hStemObj = stem(vF, abs(vFx));
set(hStemObj, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['DFT of the Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - Reconstruction

mD = exp(1j * 2 * pi * (vTT / timeInterval) * vK.') / numSamples;

% Reconstruction according to the model
vY = real(mD * vFx);

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTT, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - DFT of the Non Uniformly Sampled Data

vY  = cos(2 * pi * signalFreq * vTD);

mD = exp(1j * 2 * pi * (vTD / timeInterval) * vK.') / numSamples;
vFy = pinv(mD) * vY;

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTD, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hStemObj    = stem(vF, abs([vFx, vFy]));
set(hStemObj, 'LineWidth', lineWidthNormal);
% hLineSeries     = plot(vTT, vY);
% set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['DFT of the Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

Results are:

Summary

In this post we derived how to estimate the Uniform DFT of a Non Uniform Time Series by solving linear system of equations.

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

Remark: Why Do We Need to Apply fftshift() on the DFT of the Signal?

Indeed in the Reconstruction part we use fftshift(). The shallow answer is easy, we also build the vector vK as symmetric around zero.
But there is a deeper reason for that. In the DFT when we use uniform sampling in Frequency Domain and Time Domain Magic happens without us seeing it explicitly.

When we defined the term $ \frac{k}{ N {T}_{s} } n {T}_{s} $ we replaces $ n {T}_{s} $ with $ t $ hence we prevent the term $ {T}_{s} $ to cancel itself. Now setting $ {F}_{s} = \frac{1}{{T}_{s}} $ means that we multiply by $ k $ and we get frequencies which are out of the Nyquist Frequency.
In most cases when we that happens the Modulo property of the exponent comes in and we get the correct negative value of the frequency in the range $ \left[ -\pi, \pi \right] $. Yet when $ t $ is arbitrary we can think that $ {F}_{s} $ is changing per sample which means when we go farther than $ \pi $ the modulo doesn't bring us to the correct answer.

First, as intuition, always think the DFT is defined on the $ \left[ -\pi, \pi \right] $ interval and it is continuous. So as long as you work on this range things works as intended. This intuition can come from the Fourier Series and Discrete Fourier Series (DFS).

Let's try explaining it using a concrete example. Let's examine the exponent term from the derivation:

$$ 2 \pi \frac{k}{N {T}_{s}} n {T}_{S} = 2 \pi \frac{k}{{F}_{S}} \frac{{F}_{s}}{N} n = 2 \pi \frac{k b}{{F}_{s}} n $$

Where $ b $ is the Bin Resolution in the Frequency domain. Now given the signal is:

$$ x \left( t \right) = \cos \left( 2 \pi f t \right) \Rightarrow x \left( n {T}_{s} \right) = \cos \left( 2 \pi f {T}_{s} n \right) \Rightarrow x \left[ n \right] = \cos \left( 2 \pi \frac{f}{ {F}_{s} } n \right) $$

For $ {F}_{s} = 100 $ [Hz] and $ N = 100 $ (Which means $ b = 1 $) we will have delta at $ k = 2 $ and $ k = 98 $. For $ k = 98 $:

$$ 2 \pi \frac{98}{{F}_{s}} n $$

This is clearly above the Nyquist frequency ($ \frac{{F}_{s}}{2} $) and only for $ {F}_{s} = 100 $ its modulo is $ -2 $ which is correct. But in the model above, since we have arbitrary $ t $ one could think we have changing $ {F}_{s} $ which means we don't get the correct value.

This means the actual equation should be:

$$ x \left( t \right) = \frac{1}{N} \sum_{k = \left \lfloor - \frac{K}{2} \right \rfloor }^{ \left \lfloor \frac{K - 1}{2} \right \rfloor } X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$


How to perform zero-padding to 3 dimensional array Consider an array with $t=10000$ rows, $i=2$ columns and $n=41$ as a 3rd dimension. Now I add the columns of this array to yield a matrix with $t$ rows and $n$ columns. However in MATLAB the size of my new matrix X is:

size(X) = 10000   1   41

I want to add $10000$ zeros to each of the columns of this matrix. This should give a matrix with:

size(X) = 20000   1   41

I would try this code:

X_zeros = [X;zeros(10000,41)];

However, apparently the dimensions of X and zeros(10000,41) are not matching. I assume it is because of my 3rd dimension, although this 3rd dimension is empty. I want to delete this 3rd dimension s.t. :

size(X_zeros) = 20000   41

This is not an option:

X_zeros = [X,zeros(10000,1,41)]

because I assign X_zeros to other matrices with only $2$ dimensions. Is there a way?

A lovely cat ("concatenate arrays along specific dimension") can do the job:

a = randn(10000,1,41);b= cat(1,a,zeros(size(a)));whos b


Compute the two-dimensional DFT Compute the two-dimensional DFT [4x4] for the following 4x4 image $
 \begin{matrix}
  0.5 & 0.5 & 0.5 & 0.5 \\
  0.5 & 0.5 & 0.5 & 0.5\\
  0.5 & 0.5 & 0.5 & 0.5\\
   0.5 & 0.5 & 0.5 & 0.5
 \end{matrix}
 $

I know that DFT is separable by dimensions – one
can calculate 4 vertical transforms first, then 4 horizontal ones

For each row we get [2 0 0 0] in the first and third row
and zero elsewhere.

For each column we get [1 0 1 0] in the first and third column
and zero elsewhere.

How from this two one-dimensional dfts obtain two-dimensional one?

No you are not doing the separation correctly :

The horizontal 1D-DFT of the rows of input will be:

$ H_1 = \begin{matrix}
  2 & 0 & 0 & 0 \\
  2 & 0 & 0 & 0 \\
  2 & 0 & 0 & 0 \\
  2 & 0 & 0 & 0 \\
 \end{matrix}
 $

and the vertical 1D-DFT of the columns of $H_1$ will be:

$
 H_2 = \begin{matrix}
  8 & 0 & 0 & 0 \\
  0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0\\
   0 & 0 & 0 & 0\\
 \end{matrix}
 $

which is equivalent to the 2D-DFT of the original input.


Iterative Kalman filters and system parameters estimation i am working recently on a project in which i want to implement a Kalman filter as being an observer, and i couple this observer with a state feedback controller that produces control actions necessary to control a mobile robot.

The problem is that the governing equations for the mobile robot's motion contain errors, and this will affect the control actions produced, and at the same time the miss modelling error covariance matrix is hard to get as well as the measurement error covariance matrix. I though about using some arbitrary values for these matrices and go on with the filtration process, but i thought that it would be much better if there is anything like a kalman filter formulation that has the ability to estimate closer values for system parameters so better control actions can be obtained.

Thank you in advance.

Regarding the miss modeling, I am not sure how you can handle that because a good mathematical model is important inorder to get kalman filter getting to work. usually the more complex your model is, the better results you can obtain.

Regarding the process and measurement noise, that noise will not effect much of your result but even then if you want to estimate that noise covariance matrix, you can use iterative algorithms like LMS or RLS to fine tune your Noise parameters for some observations and use it further.


Unsupervised Way to Choose K in K-Means Algorithm There are couple of methods listed here about finding the value of k in k means algorithm. How much sound is each of these procedure when used to find out the value of k in unsupervised manner.

The idea underlying the k-means algorithm is to find clusters that minimize the intra-cluster variance and maximize the inter-cluster sum-of-square. Note that the total variance is fixed so minimization of intra-cluster variance is equivalent to maximization of inter-cluster sum of square. Basically various methods aim to achieve such minimization.

However, there is no "correct" K number or "best" method in k means. One clustering may be better than another by one metric, but worse by another metric. And depending on your own project, sometimes one clustering works but sometimes it does not.

You had better try those method to see whether they are suit in your case. If your data is not high-dimensional, you can just try kernel density estimation to observe the number of peaks. Otherwise, the easiest and heuristic one is rule of thumbs, then I would suggest Elbow method (cease increasing the number of clusters when the new generated cluster is closed to some of the existing).

There are many discussions on stackoverflow and cross validation on this topic that you may be interested in as well.


what doest "hit or miss" in morphology do? I searched for hit or miss on the internet and Gonzalez's book but there is just a formula and example. In that example the hit or miss just finds a pixel and no more answers. I want to know by concept what is use and purpose of hit or miss transform.

The Hit-or-Miss transform as the name suggest uses 2 structuring elements (SE) to identify structures which are specific to the foreground (first SE) and background (second SE).

Here, you have a good example for corner detection. The only difference is that the two SE are merged into 1: The 1 are part of the foreground SE in order to determine a specific shape in the pattern, and the 0 are part of the background SE. Consequently, the algorithm hits only if all the conditions on the 1 and 0 are respected, else it's a miss.

This principle is mainly used for skeleton computation, but it can be derived to build all the mathematical morphology basic operations.

Here are more explanations:
 - A complete video of a class.
 - Much more examples in this complete lesson.
 - StackOverflow


How to segment multiple overlapping coins (ellipses) I am doing a program that counts coins (and tells their value by their size) with OpenCV (C++).

Let's say we have an image of some coins.

If there is an overlap of two coins at the same time, it is relatively easy to segment the coins doing the distance transform, thresholding so that coins are sepparated, doing a skeleton of the background, and then drawing that skeleton into the original image.

But what if there is an overlap of 3 coins like in the image below:

What is the best way to detect 3 coins in that case?  Hough circles?

Segmentation is generally a process that is very susceptible to noise. I would better use a detector, especially for geometric shapes like coins. Remember, if you have a good detection, you also ease the segmentation problem dramatically. For the example of coins, a good model would be to use an ellipse: every circle/ellipse appears to be an ellipse under perspective projection, up to a negligible error.

Luckily, there are many good methods that can detect ellipses without segmentation. ELSD is one, Hough transform is another. In a previous post, I have explained how to use ELSD to spot spherical particles. Now I will take the Hough transform approach - this results in a minimal code that works and which I will share below. This might also be more intuitive to understand.

I will now explain a hypothesize-merge algorithm. Assume that we are given a grayscale image $I$, with gaussian smoothed edges $E$. The idea would be to apply a Hough ellipse detector on $E$ to get a large set of ellipse hypotheses $\{\mathbf{\Psi}_{1..n}\}$ and then cluster them to a reduced subset $\{\mathbf{\Psi}_{1..k}\}$ where $k<<n$. If $k$ is known (for instance this image with six coins has $k=6$) we can make use of this knowledge. Though, I will not assume this. Using the referred ellipse detector (see above) and an available mean shift clustering, I have devised the following MATLAB script (all of this can be implemented in OpenCV with reasonable little effort):

%% PARAMETERS
% blur the image to suppress some noise
gaussSigma = 1.0;

% parameters of the canny algorithm
edgeThMin = 0.03;
edgeThMax = 0.3; % usually set to 10*edgeThMin

% override some default parameters - tune them to the size of a coin
minEllipseMajorAxis = 45;
maxEllipseMajorAxis = 110;

% coins are circle like but in case of projectivity set to a lower value
minEllipseAspectRatio = 0.7;

% bandwidth of the meanshift filter. this is a quite critical parameter
% that determines the number of detected ellipses and their grouping.
meanShiftBandwidth = 25;

%%
% some more parameters do exist within the code below, but
% they should be somewhat less important.
I=imread('XqEYO.jpg');

% smooth the image a bit - increases reliability
Ig = imgaussfilt(I, gaussSigma);
G = rgb2gray(Ig);
E = edge(G,'canny', [edgeThMin, edgeThMax]);

params.minMajorAxis = minEllipseMajorAxis;
params.maxMajorAxis = maxEllipseMajorAxis;

params.minAspectRatio = minEllipseAspectRatio;

% there can be multiple correct hypotheses (filter later)
params.numBest = 200;

% pair with k other points N^2 -> k*N
% less randomization means more speed
params.randomize = 5;

% note that the edge (or gradient) image is used
% code from: https://www.mathworks.com/matlabcentral/fileexchange/33970-ellipse-detection-using-1d-hough-transform
bestFitHyp = ellipseDetection(E, params);

% a fit looks like: (ra,rb,ang,x0,y0,C,Nb)
fprintf('Output %d best fits.\n', size(bestFitHyp,1));
%
figure; image(I);
%ellipse drawing implementation: http://www.mathworks.com/matlabcentral/fileexchange/289
ellipse(bestFitHyp(:,3),bestFitHyp(:,4),bestFitHyp(:,5)*pi/180,bestFitHyp(:,1),bestFitHyp(:,2),'g');
title('Hough Ellipse Hypotheses');

% now cluster the ellipses that belong together using x0,y0 only.
% we use mean-shift clustering : https://de.mathworks.com/matlabcentral/fileexchange/10161-mean-shift-clustering
[~,idx,~] = MeanShiftCluster(bestFitHyp(:,1:2)',meanShiftBandwidth, 0);
numClusters = max(idx);
bestFits = zeros(numClusters, 6);
for i=1:numClusters % for each class
    ind = find(idx==i);
    members = bestFitHyp(ind, :);
    for j=1:6 % for each component of an ellipse record the median
        bestFits (i, j) = median(members(:,j));
    end
end

figure; image(I);
%ellipse drawing implementation: http://www.mathworks.com/matlabcentral/fileexchange/289
ellipse(bestFits(:,3),bestFits(:,4),bestFits(:,5)*pi/180,bestFits(:,1),bestFits(:,2),'g');
title('Clustered Hypotheses');

To determine the cluster centers I use a median as this might be more robust than the mean alternative - you can try different things though. Let's look at the outputs of different steps.

The edges of the smoothed image:

Initially we allow for a large number of hypotheses to survive (~120). This initial result, together with the clustered hypotheses (final detections) look like:

Of course this result is a rough detection. In other words, it is an outcome of voting and not an outcome of non-linear parameter refinement. To further improve these results, you can write an optimizer to snap the ellipses on top of the image edges. See here on descriptions of how to do that. You can also use ELSD edges instead of the Canny to ease the job of the Hough transform and to make the algorithm more robust.

It might not be the best but I hope this gives you a direction. Here are more sample outputs of running the method on coin photos fetched from Google Images:


Shannon capacity with distribution different from Gaussian If I understand correctly, in AWGN channel, for a given SNR, there is a code for $M$ input codeword $\mathbf{X}_{i, 1 \leq i \leq M}$ that
$$\lim_{n \to \infty} R=\lim_{n \to \infty}\frac{\log(M)}{n} =\log(1+SNR)$$

(where $n$ is the length of $\mathbf{X}_i$) if $\mathbf{X}_i$ is drawn from multivariate Normal distribution.

However, I don't know if $\mathbf{X}_i$ is not Normal, does the Shannon channel coding theorem hold ?

Wikipedia gives the more general definition. Channel capacity is defined using Information Theory. If you are unfamiliar with that, it will be hard to explain. The channel capacity is always the same for a fixed communications channel (meaning the statistics are fixed).

The more general definition of the capacity $C$ is
$$ C = \sup_{p_X(x)} \mathbf{I}(\mathbf{X};\mathbf{Y}) $$ where $\mathbf{I}$ denotes the Mutual Information, $\mathbf{X}$ is the input to the channel and $\mathbf{Y}$ is the output of the channel, and the supremum is taken over all possible input distributions $p_X(x)$.

The equation that you have provided is specific to the AWGN channel. The plots are demonstrating the performance ceiling that exists because you have chosen not to use a Gaussian input sequence (for convenience usually). Notice that it follows the capacity curve to a certain limit but is then limited by the constellation size.

Codewords are relevant when trying to reach capacity, but not generally applicable in the statement of what capacity is. In other words, the channel capacity is independent of the codewords chosen.


Figure of merit for 'locality of power' I have a raw signal which, if interpreted correctly with various 'tuning' parameters set to their optimum values, can be seen to consist of a relatively small (but a priori unknown) set of discrete frequency components. With the tuning parameters set at sub-optimum values, the power spectrum is spread out - each frequency component is broader or, in the very poorly-tuned case, there can be many 'false' frequency peaks. Some example plots would be:

Good tuning:

Poor tuning:

Awful tuning:

What I'd like is to be to able to auto-tune by identifying a figure-of-merit which is highest for the 'good tuning' case, lowest for the 'awful tuning' case, and somewhere in the middle for the 'poor tuning' case. If I knew a priori that there was a single frequency component, this would be fairly straightforwrad but I haven't been able to think of a good appraoch when the number of components is unknown.

Is there a standard approach to this problem that I'm unaware of?

You could try calculating Shannon entropy of the spectrum. Normalize the Fourier transform $f(x)$ of your signal so that $\int_{-\infty}^\infty |f(x)|^2\, dx = 1$ and calculate Shannon entropy as $- \int_{-\infty}^\infty |f(x)|^2 \log |f(x)|^2\, dx$. You can clamp $|f(x)|^2 \log |f(x)|^2$ to zero for very small values of $f(x)$ if its logarithm blows up and ruins the calculations. Zero-padding before windowing & FFT will give you more frequency resolution and a better approximation, especially if you have narrow peaks in the spectrum.

If you want something more standard, have a look at spectral flatness. But it seems even less stable numerically.


Skel to infinity morphological operations on 3D volume I have a 3D volume image that consists of 2D tomographic images. In MatLab, for 2D images, we can use bwmorph and call for 'skel' operation to infinity to get a pixel thick overall outline of the image. bwmorph cannot deal with 3D images, however, so I am looking around for implementing "skel" operation for such tomographic case.

I see that MatLab has provided an explanation for thinning algorithm (in 2D) but not skeleton algorithm. They are pretty much similar, but I just want to ask the community first: if I want to implement an algorithm for "skeletonize" all voxels of a volume (the end result should be all objects in the volume becoming 1 voxel thin. Ex. A 3D cylinder becomes a 3D line, a 10 voxels thick in diameter spline should become a 1 voxel thick in diameter spline, etc.), what are some of the things should I be looking for?

As a start, you could look at the distance transform. The distance transform will assign each voxel the minimal distance to the volume boundary. Afterwards you could do a non-maximum suppression. The results should look similar to a skeleton.

I bet there are more sophisticated algorithms available for 3d skeleton extraction, but this is not my field of research, so I cannot help you with that.


Inverse $\mathcal Z$-transform when region of convergence goes outwards from the inner pole? I am looking for the inverse $\mathcal Z$-transform of the following:

$$
\frac{1}{1-\frac 12 z^{-1}}+\frac{1}{1+\frac 13 z^{-1}}
$$

When the region of convergence is $z > 1/3$. I have found the $\mathcal Z$-transform for when $z > 1/2$ as a right sided signal:

$$
\left(\frac{1}{2}\right)^n u[n] + \left(-\frac{1}{3}\right)^n u[n]
$$

I can't seem to find a inverse $\mathcal Z$-transform when the ROC goes outwards from the inner pole. Does the inverse $\mathcal Z$-transform not exist for $z>1/3$?

There exist $3$ sequences with the given expression as their $\mathcal{Z}$-transform. These $3$ sequences correspond to $3$ different regions of convergence (ROCs):

$|z|>\frac12$: right-sided
$\frac13<|z|<\frac12$: two-sided
$|z|<\frac13$: left-sided

It's important that you figure out both right-sided and left-sided inverse $\mathcal{Z}$-transforms of the basic transform

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

and their respective ROCs. Then it's easy to find all $3$ inverse transforms for the given example.


Is this HW question even valid? (linear systems) The question is:
Are the following systems linear time invariant?

$y(t)=x(t)$, with $x(t)=u(t)$ (unit step function)
$y[n]=x[n]$, with $x[n]=\delta[n]$ (Dirac delta function)

The reason I am asking is because I don't see the relevance of the values given to the input functions. LTI is a property of the system not its inputs, so I don't understand what is meant!

Is it a restriction?
In that case how would I check for linearity? One probably can't, right?

there is an SE specifically for signal processing.  sorta surprized to see tags for this at the Physics SE.

what you are suspecting is precisely true.  LTI is a property of the system, not of the input.  both 1. and 2. are obviously linear and time-invariant.

one slight misconception.  in DSP we usually use square brackets, as in "$x[n]$", to depict discrete-time signals (as an alternative to "$x_n$").  a discrete-time signal is just a sequence, but we want to leave the subscript for the same purpose we might with continuous-time signals.  i.e. given a vector of continuous-time signals: $\{x_1(t), x_2(t), ... x_N(t) \}$, if we were to sample those signals, the corresponding discrete-time vector would be $\{x_1[n], x_2[n], ... x_N[n] \}$.  the argument (inside of square brackets) to a discrete-time signal is always an integer.

that said, "$\delta[n]$" is the Kronecker delta, not the Dirac delta.  (the latter the math guys like to tell us is not even a "function".)

$$ \delta[n] \ \triangleq
\begin{cases}
1, & \text{if }n = 0 \\
0, & \text{if }n \ne 0
\end{cases} $$


Inverse system in Simulink I have the following inverse system

$$G(s)= s^2 + 2s + 3$$

How do I implement it in Simulink? Note that the transfer function is only accepted if and only if the order of the numerator is less than the order of denominator.

You will have to implement the transfer function in pieces, rather than as a single LTI block.  See http://www.mathworks.com/support/solutions/en/data/1-VPJDU


Does negative frequency actually exist or it is just theoretical? This below questions gives good explanation regarding negative frequency

What is the physical significance of negative frequencies?

But I have one query

Does negative frequency only exist theoretically or does it actually exist in nature and can be measured?

Consider a wheel rotating counter-clockwise at one revolution per second. Its frequency of rotation is 1 Hz. If it rotates clockwise, its frequency of rotation is -1 Hz. It's as simple as that.


Adjusting corner frequency to constrain maximum disturbance in a high-pass filter I have a first-order high-pass filter with transfer function:
$$G(f)=\dfrac{G_0 jf}{jf + f_c}$$

where $G_0$ is the gain at high frequencies.

If I input a sine wave with frequency 1 KHz and I want a maximum disturbance of 0.1% in the amplitude, how can I know the maximum value of the corner frequency ($f_c$) that allows this error?

P.S.: Corner frequency is the frequency when the gain is at -2dB.

You set $|G/G_0| = 0.001$ and $f = 1kHz$ and then solve your equation (in magnitude form) for $f_c$


What would the plot of an two sided digital signal look like? Please, take a look at section, 1.2.3 (Signal Duration), Page-03, of the book Shaum's Outlines: Digital Signal Processing by  by Monson H. Hayes.

According to the definition, I can imagine how a right-sided and a left-sided signal look like.

But, I cannot visualize a two-sided signal.

What would an two sided signal look like?

All I need is to see a double-sided signal. As far as I know, a double sided signal is a complex exponential signal. So, that means it should be a plot of a complex exponential function. Am I correct?

A function of $x$ is right-sided if and only if there is some
$x_0$ such that the value of the function is zero for all $x <
x_0$, and left-sided if and only if there is
some $x_0$ such that the value of the function is zero for all $x >
x_0$.

Likewise for signals with time represented by $x$. A double-sided signal has neither restriction. For example a Gaussian function and $\sin(x)$ (both pictured below) are double-sided.

$f(x) = \frac{e^{-\frac{x^2}{2}}}{\sqrt{2 \pi }}$

$f(x) = \sin(x)$


Why Would Pre Filtering Measurement Data Affect the Least Squares Estimate? In estimating parameters in a discrete time model I've often seen the use of filters applied to the input data, before its applied to least squares processing. I've been told that the filters are useful if one wishes to remove bias or high frequency noise that's not considered a part of the model.

But I'm currently working with a problem where I'm using the Moore Penrose pseudoinverse to calculate the least squares. I'm filtering the input data with a lowpass filter to remove high frequency noise well beyond the model dynamics, but the choice of the filter pole seems to have a significant affect on the outcome of the estimate. And in any case the residues are near zero (they don't seem to be better or worse according to pole selection).

How should one chose to filter input data, and should it affect the estimate?

I'm not sure what's you model is.

Let's say it is something like:

$$ y = H x + n $$

Now, using the Least Squares model is optimal (In the MSE sense) when $ n $ is AWGN (It is the linear optimal estimator if the noise is white).

So unless the noise in your model is colored, no gain by filtering the data before applying the Least Squares method.

Now, what if the noise isn't IID White Noise?

If it is white just with different statistics, then use the Weighted Least Squares.
If it is colored you should apply Whitening filter on the data.

So, if the Whitening filter is given by $ W $ then you have:

$$ z = W y = W H x + W n = A x + v $$

Now you have classic LS model, where your model matrix is $ A = W H $ and your white noise is $ v = W n $.
Now you can just apply the regular Least Squares.


For a 20 MHz channel, with 1200 usable subcarriers(18 MHz occupied), what should be the size of the FFT on the transmitter before IFFT is applied? I am trying to simulate the SC-FDMA transmission and I read H Myung's paper and research on the topic but I am not sure exactly on what should be the size of FFT at the transmitter.

For a 20 MHz channel, with 1200 available subcarriers, what should be the size of the FFT block?

The research talks about IFFT at transmitter sized N = M·Q where Q is the number of users and M is the no. of subcarriers in each block of Q.

So do we perform a M pt. FFT on each Q block or on the entire 1200 symbols at the transmitter before IFFT?

Link to paper: https://ieeexplore.ieee.org/abstract/document/4099344

We perform M point FFT on each Q block, then map the FFT results to intended positions of 1200 available subcarriers (subcarrier mapping), which are finally the input of a 2048-iFFT operation. See the figure below.

Note that this is just one of possible implementations of SC-FDMA which is called DFT spread OFDM. The another is time domain approach. Check this book (Chapter 14) for more details.


Discrete wavelet transform I am unable to understand the discrete wavelet transform on images. I followed Robi Polikar's tutorial and got a brief idea about the theory. But I'm unable to understand w.r.t images.

Using Matlab's ndwt2('chess.jpg', 2, 'haar') function on the chess board , I obtained the other 7 images in the album. (Link to the album given in the end)

ndwt2() returns a structure, whose member dec contains the approximation (A), horizontal (H) , vertical (V) and diagonal (D) details.

This is where I have problem. What does A, H, V and D details of the image mean?

Also, how come is 1st image in the album approximation of the chessboard (assuming approximation means a rough estimate of the image) ? It just has lines where the borders of chess boxes are there. How is that an approximation, or have I understood it wrong?

EDIT

I was doing the mistake of converting the datatype of image returned by ndwt2 (which is the double datatype) to uint8.

these are the images that I now get.

Why does H details image contains some V details also, though the V in H are noisy? Same for V, it also contains some H details.

This was not the case previously, there H had exclusively vertical lines, V had only horizontal details, D was blank.

previous album

One can implement the standard discrete wavelet transform (DWT) on an image (dwt2 in Matlab) with a series of filtering and decimation operations, on the rows and the columns. And the  wavelet by itself results from the iteration at different levels.

Start with the Haar wavelet. In 1D, it can be implemented by a series of sums and differences ($[1,1]$ and $[1,-1]$ filter) on $2$-pixel sets. If you combine these two filters on the rows and the colums, you get $4$ possible combinations. They are illustrated in the four $2\times 2$ matrices.

So take the small $2\times 2$ block from the hat (left image). If you sum the pixels on the rows and the columns, you globally sum all four pixels. That gives you the top-left cyan on the right image. If you repeat the process on all  $2\times 2$ blocks, you get the top-left small image, often called the approximation at level 1 (A1). If you apply $[1,-1]$ on the rows, and sum the columns, you get the top-right small image. The difference along the rows detect some vertical details, as you can see (more or less vetical edges). It is called vertical details at level 1 (V1). Now, if you apply $[1,-1]$ on the columns, and sum the rows, you get the bottom-left small image, which detects
some more of less horizontal  details, called H1. The bottom-left small image computes differences on both rows and columns. It is called D1 for diagonal, though it is not very precise, in general, at detecting $\pm 45°$ edges.

There, you have A, H, V and D. A standard DWT further  decomposes the small top-left image. With the same reasoning, you get approximation, vertical, horizontal and "diagonal" details at a second level (the $2$ in your comand):A2, D2, H2, V2.

Since A1 has been decomposed, you are left with A2, V2, H2, D2,  V1, H1 and D1, hence $7$ images. With the DWT, those  images are smaller (by $4 \times 4$ and $2 \times 2$).

Here, we have taken separate original $2 \times 2$ blocks. Another word for it is decimation, by $2$, after each filter pass.

Wiht ndwt2, the main difference is that you do not decimate, it is like taking blocks with overlap. So, normally, the A, V, H, D images look about the same, except they have almost the same size as the original (up to border effects and extension).

I am surprised by the apprroximation you show. Indeed, I believe there is a read or format issue with it. Open it with an image editor, save it again an a new JPG or PNG seems to fix the problem.


Fourier transform of a sum I have a function :
\begin{equation}
C(t)=\left(1.42*\exp^{-1.192t}- 12.44*\exp^{-1.192t} +11.02 \right) u(t)
\end{equation}
where u(t) is a unit-step function
What is its fourier transform? a step by step breakdown would greatly be appreciated: notably how would you manage the fact that it is a summation, would the linearity of a fourier property apply? How would you compute the power and the phase of the sum?

would the linearity of a fourier property apply?

Of course, why not?

  How would you compute the power and the phase of the sum?

Computing each term in complex numbers, summing and converting to power/phase representation

  a step by step breakdown would greatly be appreciated

\begin{equation}
\begin{split}
C(t)&=\left(1.42 \exp^{-1.192t}- 12.44 \exp^{-1.192t} +11.02 \right) u(t) \\
 & = 1.42 \exp^{-1.192t} u(t) - 12.44 \exp^{-1.192t} u(t) +11.02 u(t)
\end{split}
\end{equation}

\begin{equation}
\mathcal{F}(C(t)) = 1.42 \mathcal{F}\left( \exp^{-1.192t} u(t)\right) - 12.44 \mathcal{F}\left(\exp^{-1.192t} u(t)\right) + 11.02\mathcal{F}\left( u(t)\right)
\end{equation}

Can you go from here? And why aren't you summing your exponential terms?


Audio from multiple microphones mixed to mono and ASR I am trying to do speech recognition on an audio file which has 1 channel (mono). My problem is that the channel is comprised of the signal of 3 microphones which have captured the audio at different distances of a person speaking and the signals have been mixed into mono, which I suspect has added reverberation and noise into the audio, making it impossible for the ASR (automatic speech recognition) to recognize speech.

I am unable to obtain the original signals from each microphone so I need to know in which way I can clean the audio or revert the effect of the downmix.

I have already tried applying the Weighted Prediction Error (Nakatani et al, 2012) filter under the hypothesis that 2 of the 3 microphones would account for reverberation in the audio as the signal has its phase shifted (due to reverberation time) but had no improvements.

I do not have the transcriptions of the audios.

What can I do? I have no more ideas on how to continue.

EDIT: This is an example utterance of the audio file (WAV file, it is in spanish and intelligible for a person)

Yoshioka, Takuya, and Tomohiro Nakatani. "Generalization of multi-channel linear prediction methods for blind MIMO impulse response shortening." IEEE Transactions on Audio, Speech, and Language Processing 20.10 (2012): 2707-2720.

Thank you.

Record the audio from the three microphones separately and simultaneously using three audio cards. Use low sensitivity microphone to avoid audio pilferage. You can use ffmpeg command line to record the audio.


Coefficient weights of Simple Moving Average (SMA) on SMA A simple moving average (SMA: arithmetic mean) is a low-pass FIR-filter. When you cascade 2 SMA's with a window of length $n$, so when you apply the second SMA on the output of the first SMA, and you want to know which coefficient weights a filter would have with the same output result, applied to the original signal, the 2 impulse responses (coefficients/weights) need to undergo a convolution (non-circular): this should return a vector of $(2*n)-1$ coefficients / weights. I refer to the image below: you have a sinusoidal signal (black), an SMA (blue), and an SMA of this SMA (green), and I compared this to a filter with the "convolution weights":

If you do this in R: the command should be convolve(rep(1/n,n),rep(1/n,n),conj = FALSE,type="open"), right? Obviously the red curve does not coincide with the green curve. Is there an analytic solution that returns the weight vector for this new moving average? Not just for 1 SMA on 1 SMA, but for $x$ SMA's on $x$ SMA's?
Here's my pretty basic R code (I assume Matlab-users / dsp-engineers understand):
n<-10
vperiod<-40
vwave<-sin(2*pi*1/vperiod*(1:(600)))
#
SMA<-function(x,n)
{
  out<-c()
  for (i in n:length(x))
  {
    out[i]<-(sum((x[(i-(n-1)):i])*rep((1/n),n)))
  }
  out[1:(n-1)]<-out[n]
  return(out)
}
#
SmaOnSma<-function(x,n)
{
  out<-c()
  cnvweights<-convolve(rep(1/n,n),rep(1/n,n),conj = FALSE,type="open")
  for (i in length(cnvweights):length(x))
  {
    out[i]<-sum(x[(i-(length(cnvweights))+1):i]*cnvweights)
  }
  out[1:(n-1)]<-out[n]
  return(out)
}
#
plot(vwave[(2*vperiod):(3*vperiod)],type="l",lwd=2,main="SMA on SMA vs. convolution weights");abline(h=0,lty=3,col="gray")
legend("bottomleft",inset=.03,c("signal","SMA of signal","SMA of SMA of signal","Filter with convolution weights"),fill=c("black","blue","green","red"),horiz=FALSE,border="white",box.col="white")
#
tempwave<-vwave
lines(SMA(tempwave,n)[(2*vperiod):(3*vperiod)],col="blue",lwd=2)
tempwave<-SMA(tempwave,n)
lines(SMA(tempwave,n)[(2*vperiod):(3*vperiod)],col="green",lwd=2)
lines(SmaOnSma(vwave,n)[(2*vperiod):(3*vperiod)],col="red",lwd=2)

Update:

Thanks to the kind answers here is the code in R for a moving average on moving average:
smavector<-function(n)
{
  return(rep((1/n),n))
}

SmaOnSma<-function(x,n,nit)
{
  if (nit==1)
  {
    cnvweights<-smavector(n)
  }
  if (nit==2)
  {
    cnvweights<-convolve(smavector(n),smavector(n),conj = TRUE,type="open")
  }
  if (nit>2)
  {
    cnvweights<-convolve(smavector(n),smavector(n),conj = TRUE,type="open")
    for (j in 1:(nit-2))
    {
      cnvweights<-convolve(smavector(n),cnvweights,conj = TRUE,type="open")
    }
  }
  #
  out<-c()
  for (i in length(cnvweights):length(x))
  {
    out[i]<-sum(x[(i-(length(cnvweights))+1):i]*cnvweights)
  }
  out[1:(n-1)]<-out[n]
  return(out)
}

I initially thought this could be a method for estimating the instantaneous frequency of a smooth curve, by correcting for the SMA's frequency response, as per my previous question on SE: which was calculated as ($(sin(n*(\pi/p)))/(n*sin(\pi/p))$), (with p = the period of the wave = 1/frequency), but it becomes clear that with every iteration, the length of the weight vector grows by a factor of $(2*n)+1$ the previous weight vector length. Even if you do only half the iterations and multiply with $-1$, the minimum length of required input data (with the same frequency), is still $2$ times the period + $1$.

The green curve in your plot corresponds to the correct result. The red curve is wrong, because the weights in the vector cnvweights are in the wrong order (left and right halves are interchanged). The correct way to compute those weights is
cnvweights <- convolve(rep(1/n,n), rev(rep(1/n,n)), conj = TRUE, type = "open")
According to the R documentation of the function 'convolve', you need to use conj = TRUE for standard convolution.


Complex band pass filter I'm trying to find some mathematical docs about the complex (analytic signal) band pass filtering in discrete time. I've read some text already but most of them describes the problem in real time only.

I would really appreciate some equations.

[EDIT:]

I found this coefficients' computation equation for real signals:

$$
w_{bp}(n) =
\begin{cases}
\displaystyle\frac{\sin\left[2\pi f_{t2}\left(n-\frac M2\right)\right]}{\pi\left(n-\frac M2\right)}-\frac{\sin\left[2\pi f_{t1}\left(n-\frac M2\right)\right]}{\pi\left(n-\frac M2\right)},  & \text{if } n \ne \frac M2 \\[2ex]
2\left(f_{t2} - f_{t1}\right), & \text{if } n = \frac M2
\end{cases}
$$

I know I have to specify the FIR filter order ($M$) which is approximately 4 / Normalised width of transition band. The application of the filer is a simple convolution of the coefficients computed by the equation above and real signal samples.

Now what I need to know is how to achieve this in complex plain. I also found a similar thread here, which describes the low pass filter design.

PS: I already have my complex $f_c = 5 \textrm{ kHz}$ signal sampled at $f_s = 32\textrm{ kHz}$ and I'm familiar with the complex plain problematics.

Take a (real-valued) low pass filter with impulse response $h_{LP}[n]$ and with cut-off frequency $(f_2-f_1)/2$, where $f_1$ and $f_2$ are the lower and upper band edges of the desired band pass filter. Define a center frequency $f_0=(f_1+f_2)/2$. Then a real-valued band pass filter is obtained by

$$h_{BP}[n]=2h_{LP}[n]\cos(2\pi f_0n)\tag{1}$$

Depending on the design of the low pass filter, the filter in $(1)$ can be identical to the one in your question.

If you want a complex-valued band pass filter with a pass band at positive frequencies (but not at negative frequencies), you simply have to modulate the low pass impulse response $h_{LP}[n]$ by a complex exponential instead of by a cosine:

$$h^c_{BP}=h_{LP}[n]e^{j2\pi f_0n}\tag{2}$$


Why are wavelet transforms (Multi Resolution Analysis) used more in practice for compression rather than Fourier series? I know that both Fourier and wavelet can be used for compression of signals.

The Fourier series guarantees that it gets the closest approximation of the original signal in the least squared sense.  It also decomposes the signal into orthogonal components.  So, I wonder why it isn't obvious to always use the Fourier series as a compression algorithm?

The Wavelet transform allows us to observe frequencies at different times, because it is localized in time (at least relative to the fourier transform).  So it is only useful when we care about the interpretability of the components, it seems to me.

Why is the wavelet compression used in practice much more than compression using a fourier series?

In short: wavelets and relatives are pretty good at compacting a sufficiently large class of regular-enough and useful signals and images. What follows are signal properties (and corresponding wavelet features) that make wavelets good (not best) candidates for lossy compression:

piecewise-smooth (vanishing moments, or cancellation of low-order polynomials)
edges or jumps (gradient-like or Laplacian behavior of wavelets)
localized oscillations (zero-average and wiggling wavelet shape)
noise or spurious events (orthogonality and sparsity enhancement)

As said by @hops, the efficiency of wavelets for compression depends of the good matching of a signal class and the chosen wavelets.

Let us  restrict here to non-redundant discrete transformations: discrete Fourier and discrete wavelets. Both are orthogonal, or close enough (bio-orthogonal wavelets) to skip the distinction.

So both, when transform coefficients are discarded, are least-squares approximations. But least-squares are not the best metric for compression: if you double the amplitude of a sample, the energy is multiplied by $4$, but it only adds one bit to the stored data, in a $\log_2$ reasoning.

In a way, a transform will help compression if it reduces the logarithmic cost or bit-budget of a discrete signal; hence, if the coefficients have a power law $1/c^{\beta}$ that decreases fast enough (the highest the $\beta$, the better).

This feature is often called the compressibility of signals. It can be assessed empirically, but also theoretically based on complicated functional analysis (Besov, Sobolev spaces).

However, consider the useful class of $C^\alpha$ piecewise regular signals, $\alpha \ge 1$. They are locally smooth, with jumps (edges). Taking the largest $M$ coefficients for compression (or nonlinear approximation), the mean squared error for the Fourier basis will asymptotically decrease in $1/M$ for Fourier (because of Gibbs ringing). While wavelet approximations can decrease as $1/M^{-2\alpha}$. The smoother the signal, the better the wavelet compressibility. In  other words, Fourier cannot use the regularity of the signal.

In practice, this is not so simple. The complex aspect of Fourier is somewhat difficult to handle.  Quantization should be taken into account. The storage of the highest coefficients location costs too. Perceptual distortion should play a role, at least for lossy compression.

So at low compression rates, local Fourier decompositions like the JPEG DCT can perform as well as wavelets, since asymptotic proofs do not apply. Indeed, local cosine bases, lapped orthogonal transforms that bridge the gap between Fourier and standard wavelets perform well too, and they are used in the MP3 standard (MDCT). For images, the 2D aspect renders separable wavelets less efficient than in 1D.

For strict lossless compression through, (wavelet) transforms are not the state-of-the-art yet.


what is the difference between QPSK and 4QAM? i need to know the main difference between QPSK and 4QAM, i have searched over the internet and still i have not found a good explenation.

There are different types of QPSK: Offset (staggerd) QPSK, $\pi/4$-QPSK, etc.

In its most known form, QPSK is identical to $4$-QAM or $4$-PSK. That is, the same constellation can be referred to by different names. So there is no difference.

This constellation can be seen in the picture bellow (from Wikipedia)

Similarly, $2$-QAM is actually the same as BPSK, as is also pointed out here.


Using morphological operations for color images processing Have you ever seen any articles about using morphological-like technique for denoising color images?

Take a look to this paper from Santiago Velasco, it's about Conditional Toggle Mapping, and it's application to denoising.


Old photographies coloring There is already a succesful question How did they color old black and white movies? on this site. Very general and correspondingly in general answered. Let me please narrow it down to just the photograpy (i.e. no time dimension) and more specified starting point:

What is the algorithm to color the greyscale HD picture given the colors at various points on greyscale original?

Note: Yes, I am aware that the original colors can't be restored any more.

This can be useful to you: Oscar Dalmau Cedeño, Mariano Rivera, Teresa Alarcón: Bayesian Scheme for Interactive Colourization, Recolourization and Image/Video Editing. Comput. Graph. Forum 29(8): 2372-2386, 2010.  ISSN: 0167-7055


Does FFT on two samples from the same audio but with a very small sampling offset return the same results? Using a single audio source, suppose I take two samples of the same length (say, 10s), with the same sample rate (say, 44100Hz), one just 1/3 sample length (7 microseconds) away from the other. Does FFT on the two samples return (roughly) the same shape? Is there a way for me to sync the two samples?

Fourier transforms possess a form of (in magnitude) shift-invariance, also called time-invariance, or shift property (see Properties of the Fourier Transform-shift properties): "a shift in time corresponds to a phase rotation in the frequency domain":
$$F\{x(t−t_0)\}=\exp(−j2πft_0)F\{x(t)\}\,.$$
This property of the continuous Fourier transform has analogues in its discretized versions. So many techniques to compensate for shifts or synchronize signals are based on this property. I'll name a few ones: cross-correlation, time-delay estimation, generalized cross-correlation (GCC), fractional delay filters. They can give you initial pointers.
In Matlab, functions like finddelay or alignsignals can be used.


Inverse fourier transform by using a random phase and an original magnitude (matlab) I want a image. Firstly I want to do fourier transform and get the magnitude and phase information. And then make a new random phase, finally use the changed random phase and original magnitude to make the inverse fourier transform. But I can not achieve it?

The code is shown here.

 img   = imread('flag.jpg');
 imafft = fftn(double(img));
 mag = abs(imafft);
 phi= angle(imafft);
 s= size (img);
 rand_phase= 4*pi*(rand(s(1),s(2))-0.5);
 newfft = mag.*exp(i*rand_phase);
 ima_out = ifft2(newfft);
 figure(2);
 imagesc(abs(ima_out));
 axis image;
 title('Ramdom phase version');

The problem is in the sentence of newfft = mag.*exp(i*rand_phase); I do not know how to modify.

Who can help me? Thank you so much.

You image is not an $N \times M$ matrix. It's actually an $N \times M \times 3$ array of numbers. I assume that's why you're using fftn instead of fft2 to take its Fourier transform.

You have to either

convert your image to grayscale; or
process each component independently.

To work get grayscale, we can naïvely average over the RGB colors like so:

img = mean(double(img),3)

To work with all three components independently, your phase array needs to have an extra dimension:

rand_phase= 4*pi*(rand(s(1),s(2), s(3))-0.5);


How replicas are formed in Frequency domain when a signal is sampled in Time Domain? I know that sampling in one domian (time or frequency) gives raise to replicas in another domain (frequency / time).
How replicas are formed?
What is this Time domain periodicity and frequency domain periodicity here in sampling?

Some people treat sampling in the time domain as the multiplication of a signal $x(t)$ by a periodic impulse train $\sum_n \delta(t-nT)$. Since the Fourier transform of a periodic impulse train is a periodic impulse train $\sum_k \delta\left(f-\frac{k}{T}\right)$ (I might be off by a factor of $T$ here, but the basic idea is correct), and
multiplication in the time domain corresponds to convolution in the frequency domain,
we get that the Fourier transform of the sampled signal is
$$X(f)\star\sum_k \delta\left(f-\frac{k}{T}\right) = \sum_k X(f)\star\delta\left(f-\frac{k}{T}\right) = \sum_k X\left(f-\frac{k}{T}\right)$$
which is a periodic function of $f$ since $X(f)$ is replicated at intervals of
$\frac{1}{T}$ Hz along the frequency axis.
If $X(f)$ is bandlimited to $\left(-\frac{1}{2T},\frac{1}{2T}\right)$, the replicas
do not overlap.

Other people resolutely refuse to consider impulse trains and come up
with other, equally valid, means of arriving at the same result.


Proper method for generating Channel coefficients in MATLAB I am trying to estimate the coefficients of a Moving average model using popular estimation method such as Least Squares (LS). For educational purposes, I am trying to see if LS can work well when the coefficients are complex valued or not. I have not seen any article or text book where LS is applied to complex parameter estimation. In many papers, I have seen that the input data is complex valued, channel is complex valued and the measurement noise is complex valued. When is this setting applicable?

What if only the channel coefficients are complex valued (input and noise are real valued), then where is this setting applicable?
What is the channel coefficients and the input are complex valued, then where is this applied?

Complex channel coefficient is just a way to represent the independent real coefficients. You just need to generate h = [h_0, h_1, h_2] = [hR_0, hR_1, hR_2] + 1i * [hI_0, hI_1, hI_2]. The independence/correlation between coefficients depend on your model. And if I were not wrong, the number of element of h is the order of your MA model.
The idea behind complex representation is that if the orthogonal set of vector space of baseband signal of bandwidth $[0,W]$ is $\left\lbrace \phi_k(t) \right\rbrace$, the orthogonal set of vector space of passband signal $[f_c-W,f_c+W]$ is $\left\lbrace \phi_k(t) \cos(2\pi f_ct), \phi_k(t) \sin(2\pi f_ct) \right\rbrace$. More details can be found at Read chapter 2.
Furthermore, as real signals have symmetric spectrum, it is convenient to represent these signals by its Baseband Equivalent versions. Thus, a complex representation is used.

So, if you work in the baseband version of signal, channel coefficients are generally complex.

The continuos time white noise is a model of infinite-dimensional signal. When you project this signal to the orthogonal set to obtain the discrete time noise samples, you follow the dimension of the orthogonal set, thus the baseband noise samples are generally complex.

The input data can be complex or real depending on how you want to use available dimensions.

After all, yes, you can process your signal as you have written, regardless real or complex assumptions.

for k = 2:N y(k) = x(k) + 0.6*x(k-1) + 0.3*x(k-2);


MMSE- Channel estimated based - DFE Equalization For HF ionospheric channels ( frequency range 3 up to 30 MHz), I Want to make channel equalization in time domain, using MMSE criterion. I have to do it using the channel estimation.

I have some questions:

Can some one elaborate more on the closed-form solution solution for the feedforward and feedback FIR filters ?
And what is the exact relation between the number of taps of the filters and the channel delay spread ?

Any help will be really appreciated

Thank you

It seems that this PhD dissertation covers part of your question:
« Improved receivers for digital High Frequency communications: Iterative channel estimation, equalization, and decoding (adaptive turbo equalization)»
https://ntnuopen.ntnu.no/ntnu-xmlui/bitstream/handle/11250/249838/126139_FULLTEXT01.pdf?sequence=1&isAllowed=y

Edit:
4.2 tapped delay line model, 4.11-4.13 relates the (complex baseband) tapped delay line to physical channel multipath. (p.p 57)


Using FFT to convert to frequency domain, then IFFT back to time domain C++ I've posted the same question on stackoverflow.com with little success, so I thought I would try here!

I'm using C++/C to perform forwards and reverse FFT on some data which is supposed to be the pulsed output of a laser.

The idea is to take the output, use a forward FFT to convert to the frequency domain, apply a linear best fit to the phase ( first unwrapping it) and then subtracting this best fit from the phase information.

The resulting phase and amplitude are then converted back to the time domain, with the ultimate aim being the compression of the pulses through phase compensation.

I've attempted to do this in MATLAB unsuccesfully, and have turned to C++ as a result. The forwards FFT is working fine, I took the basic recipe from Numerical recipes in C++, and used a function to modify it for complex inputs as following:

void fft(Complex* DataIn, Complex* DataOut, int fftSize, int InverseTransform, int fftShift)
{

        double* Data  = new double[2*fftSize+3];
        Data[0] == 0.0;

        for(int i=0; i<fftSize; i++)
        {
                Data[i*2+1]  = real(DataIn[i]);
                Data[i*2+2]  = imag(DataIn[i]);
        }

        fft_basic(Data, fftSize, InverseTransform);

        for(int i=0; i<fftSize; i++)
        {
                DataOut[i] = Complex(Data[2*i+1], Data[2*i+2]);
        }

        //Swap the fft halfes
        if(fftShift==1)
        {
                Complex* temp = new Complex[fftSize];
                for(int i=0; i<fftSize/2; i++)
                {
                        temp[i+fftSize/2] = DataOut[i];
                }
                for(int i=fftSize/2; i<fftSize; i++)
                {
                        temp[i-fftSize/2] = DataOut[i];
                }
                for(int i=0; i<fftSize; i++)
                {
                        DataOut[i] = temp[i];
                }
                delete[] temp;
        }
        delete[] Data;
}

with the function ftt_basic() taken from 'Numerical recipes C++'. The input variable InverseTransform is simply the direction of the FFT: forward or reverse.

My issue is that the form of input seems to affect the output of the Reverse FFT. This could be a precision issue, but I've looked around and it doesn't seem to have affected anyone else before.

Feeding the output of the forwards FFT directly back into the reverse FFT yields pulses identical to the intput.

However taking the power output taken asreal^2+imag^2 of the forwards FFT and copying it to an array such that:

Reverse_fft_input[i]=complex(real(forwardsoutput[i]),imag(forwardsoutput[i]));

and then using this as the input for the reverse FFT yields regular coupled pulses, with the wrong amplitude and periodicity.

And finally, taking the output of the forwards FFT and copying such that:

Reverse_fft_input[i]=complex( Amplitude[i]*cos(phase[i]), Amplitude[i]*sin(phase[i]));

where the Amplitude[i]=(real^2+imag^2)^0.5 and phase[i]=atan(imag/real) yields a really mess output in the time domain, with a peak towards the middle, resembling the power spectrum in the frequency domain.

My question is, is it the precision of the cos and sin functions which cause the output of the reverse fft to become like this? All variables are stored as type 'double'. Why is it that there is such a massive a difference between the different methods of inputting the complex data, and why is it that only when the data is directly fed back into the reverse FFT that the data in the time domain is identical to the original input into the forwads FFT?

Thank you.

( I would have posted pictures but unfortunately I can't do that yet...)

A couple of things:

Your Data[] array seems to start at on offset of "1". That's unusual for C++, where "0" is the first index. Matlab starts at 1. This may be part of the calling convention for fft_basic, but it's a bit strange
Doing the inverse Fourier transform of the Power Spectrum (magnitude squared of the forward transform) gives you the auto correlation. That has nothing to do with phase compensation at all. Wrong algorithm.
You cannot calculate the phase as atan(imag/real) since you'll get half of the quadrants wrong. The correct functions is atan2(imag, real) which takes into account the relative signs of real and imaginary part.
Just adding to the choir: solve it in Matlab first and then port to C++


How can I plot Wiener filter's behaviour for different taps (memory filter)? What's the expected behaviour of Wiener Filter for different taps (memory filter) , when input signal  is a sinusodial tone?

Do you have any idea how I can plot or simulate Wiener Filter's behaviour (estimation basically) for different value of parameter M ?
I hope my imlementation is correct. M is the filter's memory am I right ?

That's my code for Wiener Filter

clc;
close all;
clear all;

T = 1; % total recording time
fs = 500; % sampling frequency
N = fs * T; % samples
tt = (0:N-1) / fs; % time vector
fc = 10; % frequency

d = sin(2 * pi * fc * tt); % desired signal X
figure;
plot(d);

mu_w = 0;
sigma2_w = 5;

w = mu_w +sqrt(sigma2_w) * 0.5 * randn(1, N); % additive gaussian noise W * 0.5 to reduce effect
x = d + w; % signal corrupted with noise Y
l = d - x;

rr = [];
k = 1;
r = xcorr(x);
M = 25; % taps filter memory

for i = 1:1:M
    rr(i) = r(N-i+1);
end

R = toeplitz(rr);
I = inv(R);
p = xcorr(d, x);

for i = 1:1:M
    P(i) = p(N-i+1);
end

W = (inv(R)) * P';
k = 1;
y = zeros(N, 1);

for i = M:N
  j = x(i:-1:i-M+1);
  y(i) = (W)' * (j)';
  ylim([-5 5])
end

e = y' - d; % error signal

Only action needed is to change the value of M and study the plot .


interpolate 1D signal depending on 3D parameter space I have a 1D array of data d(x,p) in which a number of "bumps" or "dips" appear and/or move in the spacial dimension x, depending on 3 parameters p = (p1, p2, p3). I know d(x,p) on a regularly sampled grid in the 3-D parameter space, but now want to interpolate between those samples for arbitrary values of the parameter vector p.

I know how to do a 3D interpolation of 1D data, however this does not really do what I want: suppose a bump moves to a different position upon changing one parameter, but does not change it's amplitude, interpolating the signal would just cause a lower amplitude and wider width of this bump.

I tried fitting the bumps and then interpolating their parameters, however this also does not work since both the number of bumps, their position and amplitude changes.

So I thought that there might be ways of interpolating in Fourier space, however I have very little knowledge about this kind of signal processing. Are there standard methods of achieving this or could someone point me to relevant works or at least keywords to look into?

Here's some made-up data that should illustrate the problem:

What you are trying to do is "guess" the behaviour of $d$ given $p_1,p_2,p_3$ without going through a, potentially costly, full evaluation of the expression that generates it.

Is this interpolation? No.

Inteprolation will be handy, to a certain extent, but the problem sounds more like one of classification. If you treat it as classification then the Discrete Fourier Transform will also be useful.

So, given the values of the set of $p_i$'s, is $d$ spiky, bumpy, smooth, other? "Spiky" now becomes a region of $p_i$. It doesn't have to be continuous, it doesn't have to be smooth. All you need to know is that as $p_i$ tends to certain places in this $\mathbb{R}^3$ space, $d$ is more likely to be "Spiky".

How do you decide if something is Spiky, Smooth or Bumpy? You look at its spectrum using the Discrete Fourier Transform. You use the spectrum as a feature that tells you what $d$ looks like. Specifically, sharp waveforms in the time domain give you wide spectra and vice versa. So, a spiky $d$, has a wide spectrum. A smooth slowly varying $d$ might show a few bumps at the low end of the spectrum but then die down quickly at the top end and of course flat $d$s will only have a DC component.

You will still not be able to do a full blown evaluation of $d$ across all combinations of $p_i$. First of all, some combinations of $p_i$ might be nonsensical, which is a good thing because it reduces the size of your search space. This is where interpolation comes in, by assuming that nearby values of $p_i$ lead to qualitatively similar $d$s. This is a judgement call but of course if the expression that leads to $d$ has $p_i$ involved in subtractions in some denominator, then small changes in $p_i$s can produce wild swings in the values of $d$. So, in the end, you might want to produce a "landscape" of $d$ for a quantised grid of $p_i$s.

Here is such a landscape but from another field:

The figure is reproduced from this paper

The authors take Jansen's model which models the behaviour of a population of neurons and want to see what happens when more than one populations are coupled together.

The model has a number of parameters but there is no function that tells you what parameters to plug to the model so that its output oscillates at a particular frequency. The model is capable of a wild repertoire of outputs from smooth, to fast oscillating to spiking to even just intermittent, a few spikes and then a saturated value.

So, the landscape that you see here was necessary to examine the behaviour of the coupled neuronal population when either of them was supposed to be oscillating in different ways.

(And yes, evaluation of Jansen's model for realistic waveforms (e.g. 20-30 seconds) is costly, the model is defined as a set of coupled differential equations which you have to integrate (carefully) across all the different populations you are trying to simulate. This can be very slow, especially when you try to do it at a very high level language with lots of layers between you and the hardware).

Your application is very similar to this but you are not as lucky as only being interested in oscillatory activity which has a well defined spectrum.

So, the challenge here will be how do you quantify the behaviour of $d$ using the right features to characterise the "bumpiness" (or smoothness) of $d$.

Just out of curiosity, what does the expression that lead to $d$ looks like?

Hope this helps.


In what case do we need to increase the length of the signal compared to zero padding when using FFT? When do we have to increase the total length of the signal considered in time domain compared to just zero padding it when using FFT, this article http://www.bitweenie.com/listings/fft-zero-padding/ has a clear example and explains the difference. I understand the the concept of frequency bin and  FFT resolution but what is the waveform resolution when do we need to increase the length of signal in time domain?

Zero-padding before an FFT results in a form of high quality (Sinc kernel) interpolation of the resulting upsampled spectrum.

The interpolation of frequency peaks works well when the level of the surrounding spectrum, both the noise floor and any adjacent frequency energy, is low enough to not distort the interpolated peak.  (The examples shown in the question have a fairly high S/N).

But an interpolated spectral peak is very rounded after high ratios of zero padding (the result of convolution with the added rectangular window caused by zero padding).  Thus, any surrounding spectral energy can change the shape of this wider spectral hump enough to change the location of its peak in relation to the frequency of any actual underlying sinusoidal or other narrow band energy.  Two peaks close together can either get spread apart or end up appearing to be just one wide spectral peak in the upsampled interpolated FFT spectrum.  This problem becomes even worse when trying to find or measure a spectral notch by interpolation.

So if you are looking for more information about a signal with a relatively high local S/N (signal-to-adjacent-noise ratio) around the signals of interest, zero-padding before an FFT may be helpful (providing a computationally efficient wide band interpolator).

If you are looking for more information about spectrum just above the noise floor, or looking for information about closely spaced sinusoids or spectral peaks or notches, then more samples will be needed to avoid distortions due to interpolation of too little information about the spectrum caused by just zero-padding an FFT.

Also, the Sinc interpolation done by zero-padding a DFT of strictly real data is actually a dual Dirichlet (periodic Sinc) kernel interpolation, which causes even greater distortion in the interpolated spectrum near the DC and Fs/2 bins, or when used with very short amounts of actual data in relation to the size of the FFT.  More actual data is likely to be needed to get a more accurate and detailed spectrum near the DC or Fs/2 result bins of a DFT.


Currently we would employ a Butterworth bandstop of 490-510Hz but this doesn't completely cut-out 500Hz under some conditions

Use a notch filter instead. It will completely remove the 500Hz component and you can adjust the Q to dial in the width vs latency or group delay distortion

  Also if we make the filter too narrow it has a very long latency.

That's unavoidable: the narrower/steeper a filter is in one domain, the wider it is on the other.

  if a similar/same method is exists to subtract/remove the 500Hz signal from the signal?

Yes, many. For example you can use phase locked loop (PLL) to determine both phase and amplitude of the 500 Hz component and then simply subtract it out of the signal. This will cause NO delay or latency to the remaining signal. If your 500Hz has any type of amplitude or phase drift, the PLL will simply track it and you can adjust the speed of PLL controller to match the drift.

If your 500 Hz noise signal is guaranteed to be EXACTLY a 1/4 of the sample rate, that gets really simple. Just heterodyne it down with a cosine and a sine at 500 Hz, that gives you the real part and the imaginary part. Calculate amplitude and phase, create sine wave and subtract.


response of LTI system to a Random Input Signal what is LTI filter? what is the output when x(t) is input?

let x(t) be the input signal to the system and y(t) denote the output signal. The output of the system may be expressed in terms of convolution integral as

Why we are not taking average of h(t) also?

How comes E(X(t-tau)) to mx?

There are some basics that you really need to understand. The expectation operator of a deterministic (non-random) function just equals the function. Since the impulse response $h(t)$ is considered to be deterministic, you have

$$E[h(\tau)X(t-\tau)]=h(\tau)E[X(t-\tau)]$$

Furthermore, since $X(t)$ is known to be stationary, its mean doesn't depend on time, so you have

$$E[X(t-\tau)]=E[X(t)]=m_x$$


Find the minimum sampling rate As the signal $m(t)$ is Baseband signal after modulation it wil be bandpaas signal will look like this

So I need to use the Baseband sampling

Here $f_L=15kHz,f_H=25kHz,B=5kHz,$Now Undersampling formula

$n\leq\frac{f_H}{f_H-f_L}$
,$n\leq2.5,$ so $n=2$ (maximum),

Now $\frac{2f_H}{n}\leq f_s$,$25\leq f_s$,

so minimum sampling rate will be $25kHz$

But the answer given is $10kHz$,What should be the correct answer?

Also i am confused whether will I get back the signal m(t) after bandpass sampling like this?

What is the mistake  I am doing?

PS:Please answer it with some graphical approach if possible i am confused with the theory from this question

EDIT:

Please have a look at this

Also that means modulation will not effect sampling if we want to recover the baseband signal back,so for this minimum sampling can be calcuated by baseband sampling as $2f_m\leq f_s$

I would not use a formula for the understanding of the so called bandpass sampling (or undersampling) operation. Instead try to analyse the situation by yourself considering the signal spectrum, sampling operation and the definition of aliasing which defines the permitted range of sampling frequencies.

First, we state the fundamental principle of sampling: in order to represent a signal $x(t)$ perfectly with a set of samples $x[n]$ taken uniformly at a sampling rate $f_s$ samples per second, there should be no aliasing (spectral overlap) in the sampled signal $x_s(t)$'s spectrum $X_s(\Omega)$.

Then we define the ideally sampled signal as $x_s(t) = x(t) \sum_{k=-\infty}^{\infty} \delta(t - k T_s)$ and its associated CTFT spectrum as:
$$X_s(\Omega) =  \frac{2\pi}{T_s} \sum_{k=-\infty}^{\infty} X(\Omega - k \frac{2\pi}{T_s}) $$

Finally we ask, given the consequences of sampling on the spectrum $X_s(\Omega)$ of the sampled signal $x_s(t)$, which set of frequencies $f_s$ can satisfy the fundamental principle of no-aliasing. Then we try to determine the minimum of this set of valid sampling frequencies.

Lets apply this to your problem:
Given your real bandpass signal with a spectrum as in the figure-1.

First assume that for some $k=m$ the spectrum is shifted to right by $m f_s$ such that the shifted left-piece stands closest from left to the original right piece in its original place. Then consider the next shift at $k=m+1$ so that the shifted left piece jumps over the original right piece and stands closest (from right) to the original right piece. These two conditions are mathematicaly yielding the following constraints on the sampling frequency $f_s$:

For the case one with $k=m$ we have:
$$ -f_L + m f_s \leq f_L $$

For the case $k=m+1$ we have:
$$ -f_H + (m+1) f_s \geq f_H $$

Combining these two yields:
$$ f_s \geq 2(f_H - f_L) $$

In your case $f_H-f_L$ is $10$ kHz, so we conclude that
$$f_{smin} = 20 \text{kHz}$$

Note that this is the minimum sampling rate for real valued samples. If you are allowd to perform complex valued sampling, then you can use only one side of the spectrum (via complex analytic signal generation) and sample it at half the real sampling rate, but then you will be getting twice the number samples per second nevertheless)

Note also that this particular analysis just defines the minimum valid sampling frequency without aliasing. However, unlike the baseband sampling theorem which states that any frequency greater than the minimum Nyquist rate will also be a valid sampling frequency, in the case of bandpass sampling this is not so. For example you can easily verify that for the frequency $f_s = 2 f_L + B/2$ which is can be larger than the minimum we computed above, the resulting spectrum shifts will overlapp; there will be aliasing. So bandpass sampling frequencies are selected with care on the resulting spectral shifts.


Kalman Filter | Difference Between Minimizing the Mean Square Error (MMSE) & Maximizing Likelihood Value in Bayesian Estimation I am going through data assimilation slides on Multi Sensor Data Fusion by Hugh Durrant Whyte and it mentions:

  The Kalman Filter, and indeed any mean-squared-error estimator, computes an estimate which is the conditional mean; an average, rather than a most likely value. (Q: what is the most likely value ?)

I understand what MSQ is, but what does it mean that Kalman Filter estimates mean-square error rather then most likely value? Isnt mean square estimate, the most likely value?

Actually the first section of the notes in the link your provided are about the most likely value in the Bayesian Framework.

So we have a comparison between the Minimum Mean Square Error (MMSE) Estimator and the Maximum a Posterior Estimator.
Both are Bayes Estimator, namely they are a loss function of Posterior Probability:

$$ \hat{\theta} = \arg \min_{a} \int \int l \left( \theta, a \right) p \left( \theta, x \right) d \theta d x $$

Where $ \theta $ is the parameter to be estimated, $ \hat{\theta} $ is the Bayesian estimator, and $ l \left( \cdot, \cdot \right) $ is the loss function. The above integral called the Risk Integral (Bayes Risk).

With the the properties of Bayes Rule it can be shown:

$$\begin{aligned}
\arg \min_{a} \int \int l \left( \theta, a \right) p \left( \theta, x \right) d \theta d x & = \arg \min_{a} \int \int l \left( \theta, a \right) p \left( \theta \mid x \right) p \left( x \right) d \theta d x && \text{By Bayes rule} \\
& = \arg \min_{a} \int \left( \int l \left( \theta, a \right) p \left( \theta \mid x \right) d \theta \right) p \left( x \right) d x && \text{Integral is converging hence order can be arbitrary} \\
& = \arg \min_{a} \int l \left( \theta, a \right) \left( \theta \mid x \right) d \theta && \text{Since $ p \left( x \right) $ is positive}
\end{aligned}$$

So now, the solution depends on the definition of the loss function $ l \left( \cdot, \cdot \right) $:

For $ l \left( \theta, a \right) = {\left\| \theta - a \right\|}_{2}^{2} $ we have the MMSE estimator which is given by the conditional expectation $ E \left[ \theta \mid x \right] $. This is what Kalman Filter estimates.
For $ l \left( \theta, a \right) = {\left\| \theta - a \right\|}_{1} $ we have the Median of the posterior as $ \arg \min_{a} \int \left| \theta - a \right| \left( \theta \mid x \right) d \theta \Rightarrow \int_{- \infty}^{\hat{\theta}} p \left( \theta \mid x \right) d \theta = \int_{\hat{\theta}}^{\infty} p \left( \theta \mid x \right) d \theta $.
For $ l \left( \theta, a \right) = \begin{cases} 0 & \text{ if } \left| x \right| \leq \delta \\ 1 & \text{ if } \left| x \right| > \delta \end{cases} $ (Hit or Miss Loss) we need to maximize $ \int_{\hat{\theta} - \delta}^{\hat{\theta} + \delta} p\left( \theta \mid x \right) d \theta $ which is maximized by the Mode of the posterior - $ \hat{\theta} = \arg \max_{\theta} p \left( \theta \mid x \right) $ which is known as the MAP Estimator.

As you can see above, different estimators are derived from different loss.

In the case the posterior is Gaussian the Mode, Median and Mean collide (There are other distributions which have this property as well). So in the classic model of the Kalman Filter (Where the Posterior is also Gaussian) the Kalman Filter is actually the MMSE, The Median and the MAP Estimator all in one.

Derivation with More Details

To show full derivation we will assume $ \theta \in \mathbb{R} $ just for simplicity.

The $ {L}_{2} $ Loss

We're after $ \hat{\theta} = \arg \min_{a} \int {\left( a - \theta \right)}^{2} p \left( \theta \mid x \right) d \theta $. Since it is smooth with respect to $ \hat{\theta} $ we can find where the derivative vanishes:

$$\begin{aligned}
\frac{d}{d \hat{\theta}} \int {\left( \hat{\theta} - \theta \right)}^{2} p \left( \theta \mid x \right) d \theta & = 0 \\
& = \int \frac{d}{d \hat{\theta}} {\left( \hat{\theta} - \theta \right)}^{2} p \left( \theta \mid x \right) d \theta && \text{Converging integral} \\
& = \int 2 \left( \hat{\theta} - \theta \right) p \left( \theta \mid x \right) d \theta \\
& \Leftrightarrow \hat{\theta} \int p \left( \theta \mid x \right) d \theta \\
& = \int \theta p \left( \theta \mid x \right) d \theta \\
& \Leftrightarrow \hat{\theta} = \int \theta p \left( \theta \mid x \right) d \theta && \text{As $ \int p \left( \theta \mid x \right) d \theta = 1 $} \\
& = E \left[ \theta \mid x \right]
\end{aligned}$$

Which is the conditional expectation as required.


how to sample and modulate a gps gold code I am trying to simulate the data transmission of a GPS satellite. Telecommunications is not my field so I am struggling a bit with some concepts. I managed to generate the Gold code of the GPS transmitter using shift registers. This creates a pseudorandom 1023-vector with [-1 1] values. Now I want to create a sampled version of this signal up to 5.714Hz assuming the clock of the codes work at 1.023 Mchip/s and then modulate this signal with a sinusoidal carrier of 4.309 MHz sampled at the same frequency. I am not sure how to sample the signal.
From what I've understood, I wrote something like this in Matlab but it's not right :

a = [1 1 1 1 1 1 1 1 1 1 ]
b = [1 1 1 1 1 1 1 1 1 1 ]
for k = 1: 1023
    A(k) = a(10) ;
    B(k) = mod(b(2)+b(6),2) ;
    C_A(k) = mod(A(j)+B(j),2) ;
    Af = mod(a(3)+a(10),2) ;
    Bf = mod(b(2)+b(3)+b(6)+b(8)+b(9)+b(10),2) ;
% shifting
    sg_1 = [ Af a(1:n-1) ] ;
    sg_2 = [ Bf b(1:n-1) ] ;
end
code_vect = (C_A-0.5)*2 ;
% UP TO HERE IT'S OK, I GOT THE C/A CODE, how to modulate it?
code_vect = code_vect(index,:) % this is 1x1023 Gold code vector
sa_fq = 5.714e6     %sampling freq
ca_fq = 4.309e6     %carrier freq
cd_fq = 1.023e6   %code frequency (I am right?)
clk = 1.023e6   %clock rate
time = 1/sa_fq : 1/sa_fq :1023/sa_fq   %time vector
for i = 1:1023
%sampled signal, I'm sure this is wrong but don't know how to do it
code_vect_sa((ceil(sa_fq)*i)-(ceil(sa_fq)-1):(ceil(sa_fq)*i))=code_vect(i)*ones(1,ceil(sa_fq))
end
carr = cos(2*pi*(ca_fq+cd_fq)*time) %carrier
sig_mod = code_vect_sa*carr ; modulated signal

Can you guys give a hint on how to solve it. Thanks in advance

Have a look to the definition of the C/A signal here: where it says t you need to put your time, where it says i, you need to put your index i  from the for loop, so that concept is correct.

I think your confusion comes from using the frequency in the equation instead of the sampling time which is your time vector.


Is there a special case in which wave atoms convert to curvelets? I'm new to multi scale transformations. I was wondering if there is a special case where wave atom converts to curvelet transform? Can I use wave atom parameters to have curvelet properties?

In their 2007 paper, Wave atoms and sparsity of oscillatory patterns, Demanet and Ying  draw connections between wave atoms and other types of wavelets:

In plain words:

  We introduce “wave atoms” as a variant of 2D wavelet packets obeying
  the parabolic scaling wavelength

with a better spatial localization, meant for texture analysis. Their choice ($\alpha=\beta=1/2$) was a compromise between two types of sparsification (preservation under warping, sparsity of oscillations). Yet, one could built different  collections of such wave packets, with alternative  choices of  $\alpha$ and $\beta$, bridging the gap between the different transform choices.

Yet, their natural redundancy is 2,  with an orthonormal variant, and a complex one with redundancy 4, different from that of the curvelets. So I do no think they can be used as a special case of curvelets.

A tutorial on 2D geometric transformation can be obtain in A Panorama on Multiscale Geometric Representations, Intertwining Spatial, Directional and Frequency Selectivity, 2011.


Where does the error come from in PID control? Let's say my thermostat is 50 degrees Fahrenheit. I know that 1000 units of kinetic energy are needed to raise the temperature from 80 degrees Fahrenheit. So I calculate how much I need to increase my furnace to produce 1000 units of kinetic energy, and then the internal temperature should be 80 degrees.

What's wrong with what I just described? There's no PID control involved, no error term.

What's wrong with what I just described? There's no PID control involved, no error term.

Ok, even if you knew all this (which you practically never do), then, when do you switch off the heating? The heater's going to be hotter than the target temperature of the room (otherwise, you'd never be able to reach that temperature), so if you keep it on for as long as the room hasn't hit that temperature, then the temperature will overshoot – the heater doesn't suddenly stop being warm when you turn off the heating.

You could alternatively just keep it on for as long as it takes to insert the necessary amount of energy – but then the room is going to take very long to get to its target temperature; you'd be relying on the radiator completely heating up the room / cooling down to room temperature for the heat transfer to complete, that's going to take hours.

so, yes, control theory is necessary, even in this very simplistic example.

So, it's a problem: how to design an oven control that trades off getting it hot enough as quickly as possible (obviously, by maxing out the heating until you've reached target temperature) and not overshooting too much.

Then, of course, none of the things you claim are realistic:

You never know exactly the amount of energy needed to heat up the room, because a room isn't a sealed satellite somewhere in outer space far away from any star – it takes more energy to heat a room with cold walls (winter) than with warm walls (summer, or you just heated before letting in some fresh air). You never know how much energy will escape through walls, windows and doors – that's weather- and usage-dependent.  The list goes on....


Speech Processing applications: stereo data for features I am working on speech databases that have a stereo format. I want to extract spectral information (Mel-filterbanks, MFCC, LPCC) and also some other prosody features like the fundamental Frequency F0.

Is there a standard way to handle the computations of speech features for Stereo data?

This is for speech recognition but also spotting some irregular speech audio events, so it is different from MIR.

The maximum amplitude difference between the two channels is 7% and the mean is 0.3%. Does not tell much about the spectral features, will investigate the difference in the frequency domain now.

Is there a standard way to handle the computations of speech features for Stereo data?

Not really. It really depends on why your recording is stereo in the first place and how exactly the recordings are made.

The preferred method would be to have a single close up microphone recording. Two or more channel can be helpful, if it's far field recording. Google uses two channels, Amazon uses seven. If it's stereo from produced media (movie, song, etc.) then the best extraction method would depend on the way its produced, but typically summing the two channels would be the best.