Interpretation of Clarke's Doppler power spectral density What I understand of Doppler spread is that the relative motion between Transmitter (TX) and Receiver (RX) change the exposing time of signal. In rapport to a constant-distance TX-RX, a moving toward each other TX-RX "compresses" signal in time (signal takes less time to propagate), then signal is "expanded" in frequency domain. Similarly, a moving away RX-TX "expands" signal in time and "compresses" its spectrum. In short, that is scaling Fourier Transform. These two extreme cases set the left and right bounds of spreading an original frequency between $-f_d$  and $+f_d$ where $f_d$ is max Doppler spread.

In looking at the Clarke model, it is just multiple propagation model with rich scattering environment and equal angle of arrival. (link for more details
Clarke model)

If I understand well, there are two assumptions which are reasonale in urban environment:

Rayleigh fading
equal angle of arrival, or equal receiver sensitivity

I have followed the math from the original article, it seems ok. The final Doppler power spectrum is then
$\displaystyle S(f) = \frac{1}{\pi f_d \sqrt{1 - \left(\frac{f}{f_d}\right)^2}}$

What I don't understand is that why energy is concentrated to the two extreme spread frequency $-f_d$ and $f_d$ while angles of arrival are uniform. Is there any physical interpretation ? What am I missing from the famous Clarke model ? Personally, this model seems well-model the typical urban environment.

R. H. Clarke, A Statistical Theory of Mobile-Radio Reception, The Bell System Technical Journal, July/Aug 1968, p. 957ff

  Answers Although the answer of Carlos captures the most fundamental mathematical part, the real answer is in his comment about "mapping
  between angle and frequency". Moreover, the answer of Maximilian is
  interesting too.

A simple, "non-technical" way of thinking of it is the fact that the Doppler frequency is proportional to $\cos\theta$. The amplitudes of cosine, however, are not uniformly distributed, but are heavily weighted towards $\pm 1$.

Example plot to demonstrate, using Python/Pylab code:

theta = linspace(0, 2*pi, 1001)
x = cos(theta)
hist(x)

More rigor can be seen by noting that
\begin{align}
f &= f_d \cos\theta\\
\theta &= \cos^{-1}\left(\frac{f}{f_d}\right)
\end{align}
and the power received at any angle is proportional to a small angle increment $d\theta$:

$$
P(\theta) \propto d\theta = \frac{-1}{f_d\sqrt{1-\left(\frac{f}{f_d}\right)^2}} df
$$

And the total power can be determined by integrating the above quantity, which is identically what defines a power spectral density.


Robust watermarking in change format I watermarked a picture in sound .WAV file, when today i find this sound file after 3 years  try to extract picture embed in sound but i understand format of this chang to .mp3 and thus I can't extract completely picture of watermarked .My question about , find a watermarking algorithm that is resistant to any kind of compression?

Watermarking and steganography are very related to each other. While steganography is interested in capacity and invisibility, watermarking is to provide robustness to protect hidden data against any attacking methods like compressing, cropping, noise addition etc. LSB Coding, Phase Coding, Spread Spectrum, Echo Hiding are basic methods you should search about.

I suggest you to read these articles for the basics:

Audio Watermarking Techniques, Hyoung Joong Kim
Robust and Reliable Audio Watermarking Based on Dynamic Phase Coding and Error Control Coding

First article is very helpful to start. I advice you to try Spread Spectrum method, it was very robust against mp3 compression as I have coded myself.  I had answered my own question here trying to explain steps of Direct Sequence Spread Spectrum Method to embed data inside audio files. It may also help you.


Help with denoising signal and periodogram analysis resources This is a cross posting from the crossvalidated stack exchange as I thought this may be a better forum to ask.
I have a dataset consisting of respiratory time series signals of different lengths obtained from different groups of patients. I want to either classify or cluster the patients using these timeseries by using the commonalities of the time series of each group. However, I have no experience in dsp.

Firstly, I am confused if I am supposed to filter my signals to get rid of any frequencies above the Nyquist frequency. My sampling frequency is 32Hz and my time series is somewhat noisy and has some artifacts. I am also unsure of which filter to select for this.

Secondly, I wanted to look at the periodogram and the average power spectral density at each frequency within a group - but I am not sure if I understand the periodogram very well - if I have different time series lengths then my periodogram length will vary too, so I am not sure how this comparison can be made.

Being from Pure Math, I know Fourier analysis purely from the perspective of functions and using Fourier transforms to obtain the coefficients that describe the projection of these functions onto an orthonormal system. With periodograms however, I noticed that the x-axis represents sample frequencies. I am confused with the distinction between sampling frequencies vs. underlying frequencies of the generating function (say I have $\sin(2\pi x)$ sampled at 10Hz, does the periodogram characterize the 1Hz underlying frequency of the function?)
Any resources on understanding how to analyze and remove noisy components of time signals from a machine learning perspective would be much appreciated! Due to time constraints, I have shied away from long textbooks on digital signal processing. Thanks a lot.

Firstly, I am confused if I am supposed to filter my signals to get rid of any frequencies above the Nyquist frequency. My sampling frequency is 32Hz and my time series is somewhat noisy and has some artifacts. I am also unsure of which filter to select for this.

That ship has sailed.
Let $S=\left\{\left.\alpha e ^{i(\omega t+\varphi)}\right|\alpha > 0, 0\le \varphi <2\pi, \omega \in \mathbb R\right\}$, i.e. the set of all distinct complex sinusoid with an amplitude, frequency and phase.
Then $(S,\cdot)$ is a commutative semigroup (proof trivial).
Introducing the equivalence relation $\sim: a\sim b \iff a\left(\frac{n}{r}\right)= b\left(\frac{n}{r}\right) \forall n\in\mathbb Z$ ("two signals are identical after sampling with rate $r$"), we see that signals $s_l=\alpha e^{i(\omega_l t + \varphi)}, l=1,2,\ldots$ are $s_1\sim s_2$ if $\frac{\omega_1-\omega_2}{2\pi}=nr, n\in\mathbb Z$, i.e. we can't tell signals apart after sampling if their frequency differed by a multiple of the sampling rate.
Let's formalize this: $T=\left\{\alpha e^{i([ft+\phi]\mod 2\pi)}\right\}$ is a quotient semigroup of $S$, i.e. $T\preceq S$, and each of the elements is a leader of a $\sim$-equivalence class – and this homomorpism $S\mapsto T$ is in fact sampling, which, as we can see above, is not bijective.
Hence, all the original frequency components from $S$ got mapped to some component from $T$ with frequency normalized to the sampling frequency. That mapping is called aliasing.
What an anti-alias filter $h$ does is
$$h(s): S\mapsto S, h=\begin{cases}
s, & f< f_\text{nyquist}\\
0 & \text{else,}
\end{cases}
$$
and as you'll figure out when inserting $h(s)$ above is that this yields only the elements that are not aliased to a different frequency.
Thus, everything that "survives" $h$ will also "survive" aliasing without undergoing a change in frequency.
So, if you needed an anti-aliasing filter, it's now too late. Go and do your recording again.

Secondly, I wanted to look at the periodogram and the average power spectral density at each frequency within a group - but I am not sure if I understand the periodogram very well - if I have different time series lengths then my periodogram length will vary too, so I am not sure how this comparison can be made.

In that case, the periodogramm is not a useful mapping on its own – you'll need to add something like a truncation / padding operation to bring all signals to the same duration, for example. At which point the periodogram doesn't seem to be a sensible approach anymore.


What Is the Purpose of Bundle Adjustment? So, my understanding of stereovision is that you essentially have two cameras. Then, offline, you perform camera calibration with checkerboards to get intrinsic camera properties as well as the geometric transformation between cameras. Then, you can take a stereo pair of images, find the disparity map, and using some geometric constraints, can find the 3D positions of points in the scene.

I've read recently about "bundle adjustment" and that it should be performed at the very end? What's the purpose?

If all you want is to reconstruct a scene from a pair of images from a pair of calibrated stereo cameras, and your calibration is sufficiently accurate, then you do not need bundle adjustment.

You do need bundle adjustment if you want to reconstruct a scene from a sequence of images or for a sequence of stereo pairs, where the camera poses in each view are unknown. Search for "structure from motion" and "visual odometry".


When Kalman's assumptions are relaxed From Kalman's seminal paper "A New Approach to Linear Filtering and Prediction Problem", it is clear that Kalman's exposition is based on the following fundamental assumptions:

Measurements that are linearly related to the state.
Measurements are corrupted by white noise:
2a: Serially uncorrelated noise.
2b: Zero mean noise.
Measurements follow the Gaussian distribution.

When these assumptions are met, the Kalman Filter is "the best filter of any conceivable form", according to Maybeck (1979) "Stochastic Models, Estimation, and Control, Vol 1", page 7.

I am looking for an overview of sources that deal with the effect on the estimation error of relaxing the stated assumptions, and I am hoping someone here can help me.

HI: I haven't read it but this discusses the non-linear case so it may be helpful atleast to get an idea of what the issues are. http://sites.utexas.edu/renato/files/2017/04/recursiveUpdate_v3.pdf


Power spectrum RMS or peak I’m a newbie and wondering if I calculate the power spectrum from the peak amplitude or from the amplitude rms ?

This is my matlab script so far:

%-----------------------
% Signal time domain
%-----------------------

Fs = 2000; % Sampling frequency Hz
dT = 1/Fs; % Sampling period
L = 1024; % Length of signal [points]
t = (0:L-1)*dT; % Time vector

% Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and a 120 Hz sinusoid of amplitude 1.

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
% convert row vector to column vector
S = S(:);
%---------------------------
% Windowing
%---------------------------
win = hamming(L); % Window
winS = S.*win; % Windowed signal

% Calc window scaling factor
coherentGain = sum(win)/L;
scaleFactorWin = 1/coherentGain;
%---------------------------
% FFT
%---------------------------
Y = fft(winS);
Y = scaleFactorWin*Y;
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
Ampl = P1;
Ampl(2:end-1) = 2*P1(2:end-1); % Mult by 2 because we discarded the left side of the spectrum (skip first (DC) and last (Ny) values)
%---------------------------
% Power Spectrum : Magnitude^2
%---------------------------
% Calculate with already scaled windowed data
PS = 2 * P1.^2;     % Mult by 2 for discarded left side

% or with rms amplitude spectrum:

PSrms = 2 * (P1.*1/sqrt(2)).^2;

% PLot
f = Fs*(0:(L/2))/L;

plot(f,PS)
title('Single-Sided Power Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('P²')

plot(f,PSrms)
title('Single-Sided Power Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('P²')

This is my matlab script so far.
Im not sure how to calculate the power spectrum, either with PS or Psrms (see below)

Thanks for your help.
Mike

%-----------------------
% Signal time domain
%-----------------------

Fs = 2000; % Sampling frequency Hz
dT = 1/Fs; % Sampling period
L = 1024; % Length of signal [points]
t = (0:L-1)*dT; % Time vector

% Form a signal containing a 50 Hz sinusoid of amplitude 0.7 and a 120 Hz sinusoid of amplitude 1.

S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
% convert row vector to column vector
S = S(:);
%---------------------------
% Windowing
%---------------------------
win = hamming(L); % Window
winS = S.*win; % Windowed signal

% Calc window scaling factor
coherentGain = sum(win)/L;
scaleFactorWin = 1/coherentGain;
%---------------------------
% FFT
%---------------------------
Y = fft(winS);
Y = scaleFactorWin*Y;
% Compute the two-sided spectrum P2. Then compute the single-sided spectrum P1 based on P2 and the even-valued signal length L.
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
Ampl = P1;
Ampl(2:end-1) = 2*P1(2:end-1); % Mult by 2 because we discarded the left side of the spectrum (skip first (DC) and last (Ny) values)
%---------------------------
% Power Spectrum : Magnitude^2
%---------------------------
% Calculate with already scaled windowed data
PS = 2 * P1.^2;     % Mult by 2 for discarded left side

% or with rms amplitude spectrum:

PSrms = 2 * (P1.*1/sqrt(2)).^2;

% PLot
f = Fs*(0:(L/2))/L;

plot(f,PS)
title('Single-Sided Power Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('P²')

plot(f,PSrms)
title('Single-Sided Power Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('P²')

Neither. You calculate the power spectrum from the signal directly. Both RMS value and peak are a single value that quantifies a long period. Hence, for all frequencies higher than the reciproce of that period will be lost.


Finding Interval of Integration If we let :
$$
x(t)=\begin{cases}
1&\text{if $0<t<1$}\\
0&\text{if otherwise}
\end{cases}
$$
and
$$
h(t)=x(t/a)=\begin{cases}
1&\text{if $0<t<a$}\\
0&\text{if otherwise}\end{cases}
$$
where $0<a\leq 1$
I wish to find convolution of $(x*h)(t)$ without graphing. We have that :
$$
\max(0,t-a)<\tau<\min(1,t)
$$

Is the following assertion correct?
$$
\max(0,t-a)=\begin{cases}
0&\text{if $0<t<a$}\\
t-a&\text{if $a<t<1$}
\end{cases}
\quad\text{and}\quad
\min(t,1)=\begin{cases}
1&\text{if $1<t<1+a$}\\
t&\text{if $0<t<1$}
\end{cases}
$$

By definition:
$$x(t)*h(t)=\int_{-\infty}^\infty x(\tau)h(t-\tau) \ \mathrm{d}\tau$$
If any of the functions inside the integral equals 0, then the product does as well. We know that:
$$x(\tau)\neq0  \iff \tau\in(0,1)$$
$$h(t-\tau)\neq0  \iff t-\tau\in(0,a) \iff \tau \in (t-a,t)$$
So the convolution integral can be simplified:
$$x(t)*h(t)=\int_{-\infty}^\infty x(\tau)h(t-\tau) \ \mathrm{d}\tau
= \int_{\max(0, t-a)}^{\min(1,t)} x(\tau)h(t-\tau) \ \mathrm{d}\tau
$$
Just as you had found out.
Notice that this works just because $0<a\leq 1$. If that was not the case, then some other things should be taken into account.


Mysterious notation used in "BLIT" CCRMA paper I'm working on a project involving synthesizing audio, which led me to a promising Stanford CCRMA paper on band-limited synthesis. However, this double-arrow notation keeps popping up throughout some pretty critical portions of the document, and I've been unable to find any mention of it in the paper itself (you can imagine running a Google search for this might be tricky). It seems as though it has something to do with sampling/resampling (as in sampling the continuous variable $t$ at a rate of $nT_s$. But this doesn't seem to explain it's use in $(\leftrightarrow 1)^n$ (for lack of a better symbol). I'm hoping someone here might be able to shed some insight.

Here's a link to the paper for further examples:
https://ccrma.stanford.edu/~stilti/papers/blit.pdf

What a weird typographical error. It's a minus. Get a working copy of the PDF from http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.87.


Does decoding jpeg images uncompresses them? I train a convolutional neural network on compressed jpeg images. I noticed that my code takes pretty much the same amount of space in memory whether I use high quality jpeg images or the same jpeg images but compressed. I know that processing a jpeg image starts by decoding it, so my question is: does 'decoding' a jpeg image uncompresses it?
In other words, does a compressed image regain its original size (in bytes) when decoded?

Image processing is mostly done on frames. The digital image frame is a rectangular raster of pixels, either in an RGB color space or a color space such as YCbCr. So, as you noticed, you can be sure that your compressed images are decoded to rectangular rasters in your application. Each encoded image file stores with its data the source bitmap size (width and height) and the source color space (pixel format and bitness). When decoding, these parameters are restored. So, irrespective of what lossless/lossy algorithms are used, the sizes (in bytes) of recovered raster rectangles are identical, if both lo-res and hi-res compressed image files are received by encoding the same source bitmap data.
EDIT: A "multiresolution" comment fairly corrects my answer, and there is the other explanation why the application may not allocate in CPU memory the entire memory occupied by a source raster: a decoder can use GPU memory. But OP notices that decompressed images take "pretty much the same amount of space in memory" irrespective of encoded image quality. Even with multiresolution, the application can decode the image to source dimensions. Also, code performance dependence on image quality is not part of the question.
Still, the "multiresolution" comment is absolutely necessary for completeness.


Plot the magnitude of coefficients against frequencies of DCT transform When we perform DCT transform using matlab:

 number_of_sample_data = 100;
 f=0.1; %frequency
 fs=1; %sampling rate
 phi_0=0
 t=0:1/fs:((number_of_sample_data-1)*1/fs); %generate 100 samples
 x=sin(2*pi*f*t+phi_0);
dct_transform=dct(x); %perform discrete cosine transform
figure(1), stem(dct_transform);

this piece of code will plot the magnitude of coefficients against frequencies or something else?

Actuually, that code doesn't work without defining axis t variable, and frequency f, regarding to what that code does, is simple:

x=sin(2*pi*f*t);        %x is a sinusoidal vector of frequency f and length t;
dct_transform=dct(x);   %dct_transform is a vector containing the unitary discrete transform of x
figure(1), plot(dct_transform); %plots the the unitary discrete transform of the sinusoid against default axis in Matlab

So you are just plotting the discrete transform vector against sample position, wish is what Matlab does by default.


Documentation about filters in general I already posted this question on StackOverflow. I got a suggestion to go for a better answer here.

To make the question more concise I'm interested in an introduction to digital filters, any resource is good.

But first I will be satisfied if someone could point me to some resources explaining this filter (used in Android):

http://gitorious.org/rowboat/frameworks-base/blobs/671a6ff4be11b3e2d8eb017e0c7a78e6133fb2b8/services/sensorservice/SecondOrderLowPassFilter.cpp

What I'm interested is the way the filter's parameters are chosen. While I can copy it without thinking, I guess I should understand the basic concept/idea behind before using it.

Thanks,

Iulian

You can derive the expression for the coefficients by doing bilinear transformation of the following analog low-pass prototype filter

$$H(s) = \frac{w_0^2}{s^2 + (w_0/Q)s + w_0^2}$$

where $w_0$ is the cut-off frequency.

You can lookup the bilinear transformation on Wikipedia.

The filter used in the Android app is a Butterworth filter because the chosen value of Q is $1/\sqrt{2}$. Note that in the constructor the inverse of Q is computed and assigned to variable iQ which is used in the computation of the coefficients. Note also that the variable K holds the 'frequency-warped' value of the specified cut-off frequency. You can find more information about the frequency-warping phenomenon in the above link.

You can find many examples on digital filter design using bilinear transformation. I found this one, which is quite close to the Android example.


What is "Maximum Quantization Error"? I have an formula for this "Maximum Quantization Error" but i dont know what it is based in. Its just thrown in my study material without further explanation.

It is defined as:

$$Q = \dfrac {\Delta x}{2^{N+1}}$$

where $N$ is the number of bits used for quantization in a analog to digital conversion, and $\Delta x$ is, in portuguese "Faixa de Excursão do Sinal", I don't know what would be the correct translation, but I bet on something like "Signal Excursion Band". I know, its a strange name.

Can someone help me with this? What is this $\Delta x$?
Sorry for my bad english, it isnt my native language.

When you quantize a signal, you introduce and error which can be defined as
$$q[n] = x_q[n]-x[n]$$
where $q[n]$ is the quantization error, $x[n]$ the original signal, and $x_q[n]$ of the quantized signal. The maximum quantization error is simply $max(\left | q \right |)$, the absolute maximum of this error function. Dx in this definition seems to be the range of the input signal so we could rewrite this as
$$Q = \frac{max(x)-min(x)}{2^{N+1}}$$

Let's look at a quick example. Let's assume you have a signal that's uniformly distributed between -1 and +1 and you want to quantize this with 3 bits. You have a total 8 of quantizaton steps which would map to [-1 -.75 -.5 -25 0 .25 .5 .75]. The difference between steps is 0.25. If you round during quantization the maximum error will be half of that (i.e. 0.125).

Now let's try the formula:
$$Q = \frac{max(x)-min(x)}{2^{N+1}} = \frac{1-(-1)}{2^{3+1}}= \frac{2}{16} = 0.125$$


Filtering and peak finding I have an observed signal which is the norm of the sum of two 3-vectors:
$$
S(t) = || \mathbf{A}(t) + \mathbf{B}(t) ||
$$
Now, I have full knowledge of the vector $\mathbf{A}(t)$,
$$
\mathbf{A}(t) = \pmatrix{A_1(t) \\ A_2(t) \\ A_3(t)}
$$
i.e. I know each of the components $A_i(t)$. The signal $\mathbf{B}(t)$ is unknown, and $S(t)$ is what I observe. Now the signal $\mathbf{B}(t)$ has peaks, particularly in the $B_1(t)$ component, and my goal is to determine the position of these peaks. I don't care about the amplitude of the peaks, I only want the position (time).

So what I need to do is somehow filter out the signal $\mathbf{A}(t)$ from the observed signal $S(t)$ and then apply a peak detection algorithm on the residual signal. My question is whether it is possible to do this?

As a simple example, I made a matlab script to construct synthetic data as follows:
$$
\mathbf{A}(t) = \pmatrix{ X(t) \\ Y(t) \\ 0}
$$
$$
\mathbf{B}(t) = \pmatrix{ G(t) \\ 0 \\ 0}
$$
where $X(t) = \sin{(2 \pi (f/4) t)}$, $Y(t) = \sin{(2 \pi f t)}$, and $G(t)$ is a single Gaussian peak. The signals and matlab script are given below. As you can see, it is not easy to tell that there is a peak in my observed signal at $t = 500$ as there is in the unknown signal $G(t)$.

So is it possible to filter $S(t)$ to remove $\mathbf{A}(t)$ and clearly find the peaks in $\mathbf{B}(t)$?

N = 1000;
t = [1:N];

f = 20;
X = sin(2*pi*(f/4)*t/N);
Y = sin(2*pi*f*t/N);

G = gausswin(N, 100.0);
Z = zeros(1,N);

A = [ X ; Y ; Z];
B = [ G' ; Z ; Z ];
C = A + B;

signal = vecnorm(C);

figure

subplot(4,1,1);
plot(t,X);
title('Known X');

subplot(4,1,2);
plot(t,Y);
title('Known Y');

subplot(4,1,3);
plot(t,G);
title('Unknown signal G');

subplot(4,1,4);
plot(t,signal)
title('Observed signal');

Generally, this should not be possible because:

you project a 3D information onto a 1D signal, and the system is underdetermined in general,
you add a non-linearity (the norm), that further hinder restoration.

However, using the triangle inequality, you can bound $B$ from below:

$$ \|A+B\| \le \|A\|+\|B\|$$

hence

$$ B_{\inf} =\max(\|A+B\| -\|A\|,0)\le \|B\|$$

In the example you gave, if the signals $A$ and $B$ are weakly related, outcomes can become interesting:

signal = (sum(C.^2)).^(0.5);
signalInf = max((signal-(sum(A.^2)).^(0.5)),0);
signalG = (sum(B.^2)).^(0.5);

plot(t,[signal',signalInf',signalG'])
title('Observed signal');
legend('A','Binf','B')

From here, you can recover putative peaks from $ B_{\inf}$. However, it is unlikely to work every time.

Could additional filtering help on $S(t)$? Without knowledge about the relative spectra of $A$ or $B$, I don't know (yet).


where we use estimation and detection in communication system? I read from the book that we use estimation when we are dealing with analog signal or continuous random variable and the term detection when the transmitted signal is a discrete value even if a received signal is continuous.
I want to know if in digital communication we use both estimation and detection and in which step we do these process.
Also, when we do equalization in the process.

Don't worry too much about defining these terms too precisely, because they are used in many contexts with slightly different meanings.
In very general terms, "estimation" is the calculation of a signal parameter, for example the phase, the mean, the PSD, etc. In other words, you have a signal, possibly noisy or distorted, and you want to find something about it.
"Detection" is concerned with whether you actually have a signal or not; think about radar, where you send a pulse, and wait to see if you detect the reflection. Since recovering the information from a PAM signal is based on the same principle (i.e. matched filtering), in digital communications it is common to say that transmitted symbols are detected.


Denoising effect in GnuRadio OFDM Serializer block Why does the OFDM serializer has such a strong denoising effect in the flow graph below ? Is that normal ?

The upper constellation is AFTER the OFDM serializer block and the lower is BEFORE the OFDM serializer block.

Here is the OFDM documentation : http://gnuradio.org/doc/doxygen/page_ofdm.html

It does not help to explain why the serializer block has such a strong denoising effect, any idea why this denoising happens ?

Here is the flowgraph (in .grc format): http://pastebin.com/raw/PTY0Q0Ty

Further inspection indicates that the serializer block only does remove non data carriers. It just probably so happens that anything that is non data is super noisy, and data carriers are not noisy at all, but I still wonder how is this possible.

Further inspection shows indicates that the serializer block only does remove non data carriers. It just probably so happens that anything that is non data is super noisy, and data carriers are not noisy at all, but I still wonder how is this possible.

The magic that happens here is in the actual equalizer used in the frame equalizer block. If you'd scroll up in the GRC¹, you'd probably see the payload equalizer object block, holding a "simpledfe" equalizer.

Now, simpledfe stands for simple data feed-back equalizer. It's actually pretty well-documented, but somehow the documentation tool broke and the HTML documentation doesn't actually contain any of the explanation given in the source code.

So, here's a source code excerpt with the docs:

/* \brief Simple decision feedback equalizer for OFDM.
 * \ingroup ofdm_blk
 * \ingroup equalizers_blk
 *
 * \details
 * Equalizes an OFDM signal symbol by symbol using knowledge of the
 * complex modulations symbols.
 * For every symbol, the following steps are performed:
 * - On every sub-carrier, decode the modulation symbol
 * - Use the difference between the decoded symbol and the received symbol
 *   to update the channel state on this carrier
 * - Whenever a pilot symbol is found, it uses the known pilot symbol to
 *   update the channel state.
 *
 * This equalizer makes a lot of assumptions:
 * - The initial channel state is good enough to decode the first
 *   symbol without error (unless the first symbol only consists of pilot
 *   tones)
 * - The channel changes only very slowly, such that the channel state
 *   from one symbol is enough to decode the next
 * - SNR low enough that equalization will always suffice to correctly
 *   decode a symbol
 * If these assumptions are not met, the most common error is that the
 * channel state is estimated incorrectly during equalization; after that,
 * all subsequent symbols will be completely wrong.
 *
 * Note that the equalized symbols are *exact points* on the constellation.
 * This means soft information of the modulation symbols is lost after the
 * equalization, which is suboptimal for channel codes that use soft decision.
 *
 */

Hm, but what about unoccupied carriers?

Now, looking at the implementation in ofdm_equalizer_simpledfe we see:

void
ofdm_equalizer_simpledfe::equalize(gr_complex *frame,
      int n_sym,
      const std::vector<gr_complex> &initial_taps,
      const std::vector<tag_t> &tags)
{
 […]

  gr_complex sym_eq, sym_est;

  for (int i = 0; i < n_sym; i++) {
    for (int k = 0; k < d_fft_len; k++) {
      if (!d_occupied_carriers[k]) {
        continue;
      }
      […]

In other words: those are left untouched from the original FFT.

Now, a typical direct conversion OFDM system will leave the DC carrier unoccupied – that containing the leakage of the receiver LO.
That is typically a relatively powerful FFT bin, shouldn't vary much in magnitude, and its phase should depend on the point in time the DFT was taken, and that's the frame start time as determined by the Schmidl&Cox sync. That would be my explanation of the $|\cdot|\approx 50$ constellation points, and the rest might very well be the unequalized noise in original amplitude, especially since we're looking at a constellation sink that shows 1024 constellation points at once – so we'll see some outliers in noise power.

¹ It's often impossible to fit a whole GRC flow graph on the screen. It's therefore really handy that there's the "save screen capture button" in the toolbar and the menu!


Relationship between signal rate, data rate, minimum bandwidth I need to understand the relationship between bit rate, signal rate and minimum bandwidth. I read in a textbook that

  We define three cases for relationship between bit rate and signal rate for this, they are: the worst, best and average. Worst
  case is when we need maximum signal rate, best case is when we need
  the minimum. In data communication, we usually prefer the average case
  and the relationship between data rate and signal rate is

  $$S = c \times N \times \frac{1}{r} \quad\text{baud}$$

  where $N$ is data rate, $c$ is case factor, $S$ is no number of signal
  elements and $r$ is previously defined ratio.

I don't understand what the above formula signifies.

(i) $S$ is said to be the number of signal elements and not the signal rate, i.e number of signal elements per unit time. How does this then tell the relationship between signal rate and data rate?

(ii) What exactly is $c$ and $r$ in the formula?

(iii) Secondly, there's a formula given for calculating the minimum bandwidth  of a digital signal which is also exactly the same as the formula given above for the relationship between data rate and signal rate, i.e, $B_\text{min} = S = c \times N \times \frac{1}{r}$ . How does this same formula calculates the minimum bandwidth?

Here's an attempt at an answer, by comparing your equations with the more commonly used ones.

We can write the "pulse rate" or the "signaling rate" $R_p$, measured in baud, or equivalently in pulses per second or symbols per second, as: $$R_p = \frac{1}{k} R_b,$$ where $R_b$ is the bit rate, and $k$ is the number of bits carried by a pulse.

To elaborate: a digital communications signal can be written as a sequence of pulses of different amplitudes, each pulse time-shifted an amount $T_p = 1/R_p$, and with the information carried in the pulses' amplitude: $$s(t) = \sum_j a_j p(t-jT_p).$$

Note that:

For each value of $j$, the pulse $a_j p(t-T_j)$ carries the information in its amplitude $a_j$. The amount of information depends on how many distinct values $a_j$ can take. Usually this number is represented by $M = 2^k$ for integer $k$, and we say that each pulse carries $k$ bits.
Since one pulse is transmitted every $T_p$ seconds, the pulse rate is $R_p = 1/T_p$ baud.

Now, (ignoring some more advanced considerations) the bandwidth of $s(t)$ is the same as the bandwidth of the pulse $p(t)$, denoted $B_p$. The pulse that has the minimum possible bandwidth for a given $R_p$ is a sinc pulse, and in that case the bandwidth of $s(t)$ is given by $$B_{min} = \frac{R_p}{2} = \frac{1}{2}\frac{1}{k}R_b.$$ Practical pulses achieve a fraction larger than $1/2$, or in other words, they require a larger bandwidth. If we define this fraction by $c$, we have $$B = c\frac{1}{k}R_b.$$

In conclusion, it seems to me that:

In your equations, $N$ is the bit rate, $r$ represents the number of bits transmitted per pulse, and $c$ represents the bandwidth efficiency of the chosen pulse. The variable name $k$ is more commonly used than $r$.
The first equation in your book is wrong; it should be $S = N \times \frac{1}{r}$.
The equation for the bandwidth is correct.


How to do a spectral fade out? I'd like to do a fade out of a sound in a few seconds, but in a "spectral way", with a special curve :

(here x-axis is time, y-axis is frequency from a few hz to 20khz).

All the part "to be removed" has to be removed !

Do you know how to do this ? (an EQ / simple IIR/FIR filter  won't work for this, etc.)

You want a time-varying low pass filter, automated so that the cutoff frequency follows your curve. For such slow modulations most IIR low pass filters are stable, so you can just find a steep enough low pass filter and update its coefficients while filtering while keeping the filter state.


Sample rates, Samples per Symbol, and Digital Pulse Shaping Having some confusion about digital pulse shaping for complex baseband (passband) signals. The complex baseband linear modulation equation is
$$s(t)=\sum_{m=-\infty}^{\infty}\text{Re}\{a_m\}h(t-mT)+j\sum_{m=-\infty}^{\infty}\text{Im}\{a_m\}h(t-mT)\tag{2}$$ where $a_m$ are the symbols and $h$ is the impulse response of the pulse shaping filter. Consider for now this is a root raised cosine filter with QPSK symbols.
It is often stated in references that the output rate of the pulse shaping filter should be at least twice the input rate - in other words, the pulse shaping filter should be at least a 2x interpolator. For example the following equation holds for the required (positive) bandwidth:
$$BW = (1+a)\frac{R_b}{2\log_2(M)} = (1+a)\frac{R_s}{2}$$ where $R_b$ is the bit rate, $a$ is the excess bw, $R_s$ is the symbol rate, and $M$ is the constellation size. Thus to satisfy the nyquist criterion, we must sample with at least $f_s = (1+a)R_s$. I believe another way of stating this is that the two-sided bandwidth $(1+a)R_s$ is equal to the complex sampling rate.
So considering some values, take $a=0$ which gives the sinc pulse, the digital pulse shaping filter could theoretically operate at 1 (complex) sample per symbol, i.e. the input and output rates are both $R_s$ (disregarding all the other obvious drawbacks of using sinc pulses).
Another common value $a=0.35$, gives us a sampling rate of $1.35R_s$, so we would need at least a fractional interpolating filter (e.g. upsample by 27, downsample by 20) to avoid aliasing. So clearly it seems possible to achieve less than 2 complex samples per symbol. My question is, if this is correct, why do so many references and libraries mandate a minimum of 2 (complex) samples per symbol when the actual limit is 1 complex sample per symbol (e.g. GnuRadio constellation modulator is one). I am assuming this is just due to simplicity and 2 being an easy value to interpolate and cover the range of allowable RRC bandwidths.
A second related question, assume we are doing BPSK, e.g. with $\text{Im}\{a_m\} = 0$. It seems to me in that case if you're in an IQ system still, you'd actually need 2 complex samples per symbol in order to actually give 2 real samples per symbol to the BPSK I channel. Perhaps this is yet another reason for mandating the 2 complex samples per symbol minimum?

Yes the OP is correct in that you can implement pulse shaping in less than 2 samples per symbol for exactly the reasons that was outlined. However importantly we must also keep in mind having excess bandwidth to simplify subsequent filtering required (such as after the DAC on the transmitter side). The Nyquist criteria is the sampling rate must be twice the highest bandwidth of the signal. Therefore slower roll-off signals can be done with less samples per symbol as the OP has outlined. However they will require a much larger time response, meaning more samples overall for each symbol pulse extending well beyond the symbol period, at less than 2 samples per symbol. The reason 2 samples per symbol is commonly stated is for the convenience of implementation given any synchronization between symbol rates and sample rates.
A BPSK system that is complex would be complex for purpose of carrier recovery only; once the received symbol is rotated to the real axis there would be no further need for the imaginary axis. Yet the decision for number of samples per symbol is the same whether the signal is complex or real, exactly as first outlined (must be greater than twice the bandwidth of the signal, and specifically it is the bandwidth of significance- where beyond we are not concerned with the aliasing of residual spectrum that exists beyond the bandwidth given by the sampling rate chosen.
To help see all considerations here are spectral plots for two variants of a root-raised cosine implementation implemented with a 100 Tap FIR. The first is done with 10 samples per symbol while the second is with 2 samples per symbol. Specifically notice the much greater rejection achieved close to the primary bandwidth due to the longer overall time response of the filter (so same complexity but with fewer samples per symbol the time duration of the impulse response is longer, thus closer approximating the ideal infinite long response).  The advantage of the first one with more samples per symbol is in the ability to have more relaxed filtering in subsequent stages, such as required analog filtering after the DAC since the image that would be at the sampling rate is now that much further away from the desired passband (or subsequent digital interpolation stages).

The sampling rate for the lower figure specifically is 2 samples per symbol and here specifically in this example it is 2000 Hz. We see from the plot, as long as we are not concerned about folding noise below -60 dB that we could theoretically decrease this sampling rate to approximately 700 Hz while still maintaining the passband characteristics of the pulse. However consider what would be required for a transmit filter after the DAC to filter out the image of the passband that is centered on the sampling rate (in addition to the natural Sinc filter the DAC provides; which is not nearly sufficient to meet out of band emission requirements!). One solution is to interpolate between the lower rate pulse sampled waveform and the DAC, but even the proper interpolation filter design would be challenged by not having any usable transition band!
This illustrates a fundamental challenge in decreasing the sampling rate to exactly twice the bandwidth (where bandwidth would be the significant bandwidth after which any folding/aliasing is far enough below that dictated by a performance requirement), which is more of a concern on the transmitter side dictated by the subsequent filtering required to meet out of band emissions: given the image frequency response of the passband that would be symmetrically around the sampling rate; if we sampled at exactly twice the significant bandwidth in the transmitter the required filtering would not be realizable given there is no transition band for the filter to reject the images after the DAC that would be positioned at integer multiples of the sampling rate.
On the receiver side, we only need sufficient excess bandwidth to handle frequency offsets due to Doppler and Tx/Rx synchronization, so there is much more of an opportunity/practicality to reduce the sampling rate below 2 samples/symbol. Two samples per symbol is still convenient for common timing recovery implementations such as the Gardner TED (Isn't Gardner's algorithm and Early-Late gate the same thing?), and allows for a simple extraction of the one sample per symbol decisions.
Related to receiver considerations are baud rate equalizers which operate at 1 sample per symbol, with performance implications related to exactly the points the OP makes, but point is they are able to still function. This is detailed further here The benefits of a fractionally spaced equalizer. This is also similar to how the M&M timing recovery algorithm operates at one sample per symbol, but also has inferior performance to the Gardner in the presence of carrier offsets.


How to reconstruct the signal using overlapping frames in MATLAB I understand that we can use the buffer function in matlab to split a long signal into overlapping frames. But is there any functions that we can use to put the frames back to reconstruct the signal?
Here's my current code in Matlab
n = 0.1:0.1:100;
y = sin(n*pi*0.5)';
Nw = 100;  % window length
Nh = 25;   % hop size
Y = buffer(y,Nw,Nw-Nh, 'nodelay');  % Y is the matrix containing the overlapping frames

What I want to do next is to use Y to reconstruct the original signal y. Something like this
yNew = unbuffer(Y, Nw, Nw-Nh)

Of course there is no such function as unbuffer, not that I've known of. So my question is, is there any function available to do the job of unbuffer ?
The situation is similar to this question here: https://www.reddit.com/r/matlab/comments/361y9w/reverse_buffer_function/, which was 5 years ago and still received no answers.

unbuffer in Python; helpful post. Nice visuals in MATLAB docs (bottom).


Oversampled Binary Image Sensor vs Active Pixel Sensor(CMOS) I am a beginner in the field of Image processing.I was designing a camera sensor when I came across these two type of sensors:

CMOS(Active Pixel Sensor)
Oversampled Binary Image Sensor.

Now, I can find many websites for comparing CMOS vs CCD but nowhere can I find a comparison with OBIS. Why is this? Is the OBIS too outdated to be used? I was thinking of designing a camera that captures images using OBIS and later converts them into RGB format (with a variety of processes I have not discussed here). As advantages of Binary Sensor I mentioned Low image size and super fast processing. Are these really going to be an advantage over CMOS sensor? Will my camera be a success or is this a really bad idea trying to use Binary Image sensor? Is using Binary Sensor a good idea or is it a too outdated idea?

Thanks in advance for helping this novice!!

The success of your camera depends on too many factors! and they're out of the scope of this page. If the dynamic range requirements of your application are not too demanding then you could take advantage of the fast sampling rate and try to reconstruct a low quality image (by the same principle by which sigma delta modulation works, eg reconstructing a PCM signal from the output of a 1 bit output digital microphone). There's technique called halftoning, which has been the basis for reproducing images in print media since the invention of the newspaper. It is based on sampling an image with 1 bit and diffusing the error. Determining if using a Binary Sensor instead of a CMOS is a good or bad idea depends on your intended application. You should also pay attention to the size of the objects you want to picture. CMOS technology is widely used and you could probably find a cheap sensor to use in your camera.


Continuous phase for phase delay calculus in FIR filters I would like to estimate the phase delay accurately for any random FIR filter.

The definition of the phase delay is the continuous phase divided by the angular frequency (with a sign change). That means that you can't simply use the wrapped phase, between $-\pi$ and $+\pi$ you can get with a FFT to obtain the phase delay estimation. At least, you need to use an unwrapped version of the phase.

For example, with a linear phase FIR filter, the phase delay is supposed to be a constant, so the phase you have to use to do the calculus has to be monophonic and decreasing.

However, by doing some tests in MATLAB, I have been able to see that the unwrapped phase doesn't have this property everywhere. For a halfband lowpass filter for example, there are a few "jumps" in the stopband with a length not related with $\pi$, so they are untouched in the unwrapped phase, which gives me a wrong estimation of the phase delay, featuring jumps as well instead of a constant value.

Indeed, the "continuous phase" I need isn't the same thing than the unwrapped phase, and is associated with the zero-phase amplitude of the FIR filter. I have not found yet a way to get it, and MATLAB with its zerophase function is giving wrong results as well, because it cheats ! If the filter is linear phase, it uses the formula we all know for the phase delay. If it is a FIR linear phase transformed into a minimum-phase filter for example, then the phase delay displayed features again this low size jumps.

I'm thinking about doing something like detecting the jump length, assuming the phase is supposed to be smooth everywhere, and using the derivative at the sample before the jump to do something which might look good...

But I would like to know if someone knows a rigorous way to get the continuous phase.

EDIT : here is an example. You can find there the phase for a FIR linear phase filter, wrapped and unwrapped. As you can see on the unwrapped, some phase jumps are still here...

EDIT 2 : my problem can be redefined this way. Most of the time, any FIR filter can be studied with its frequency response H($\omega$) = |H($\omega$)| $e^{j \phi( \omega)}$ with |H($\omega$)| the magnitude, always positive or zero, and $\phi(\omega)$ the wrapped phase.

There is an alternative notation : H($\omega$) = A($\omega$) $e^{j \theta( \omega)}$ with A($\omega$) the amplitude or zero-phase amplitude and $\theta( \omega)$ the continuous phase representation, still between -$\pi$ and $\pi$. Apparently, it features only 2$\pi$ jumps, and can be fully continuous with a standard unwrap algorithm. That's because A($\omega$) can be positive or negative. Each time there is a zero in the FIR filter, A($\omega$) changes its sign. That's the zeroes which are causing the lowest jumps. Since A($\omega$) contains the sign information, the low jumps are no more present in the new phase representation $\theta( \omega)$.

So, what I need now, is in fact a way to calculate accurately $\theta( \omega)$, maybe by calculating first A($\omega$). I might detect the low jumps easily and change the sign of |H($\omega$)|, then use it to get the remaining $e^{j \phi( \omega)}$. But there is maybe a better way of doing this, and I'd like to know if someone familiar with this notation has already done that.

The problem is that magnitudes close to zero in the stop band change their sign almost arbitrarily / randomly.  And that's why you're seeing arbitrary jumps of $\pi$ radians --- it's a sign change.

Because of that, I'd be tempted to just ignore the phase of any part of the magnitude response below some low value. What value to use will depend on your application.

Zero padding will not generally help because that just gives you more "random" sign changes around the frequency axis.

I use "random" in quotes as it's not stochastic, but it is hard to predict.

If I interpolate the phase on 100,000 points for an FIR filter of length 126 taps and unwrap the phase I get the picture below.

The reason the unwrapping algorithm fails is because the phase in the stop band is changing very quickly compared with that in the passband. Comparison below of unwrapped phases.  The red lines are $\pm \pi/2$.

Changing the complex_magnitude expression to mimic the fftshift idea from @hotpaw, to:

    complex_magnitude <- polyval(f1[c(63:125,1:62)], exp(-1i*freq))

I get:

which is slightly better, but not ideal.

R Code Below

#29899

library("signal")

f1 <- remez(125, c(0, 0.3, 0.4, 1), c(1, 1, 0, 0))
plot(f1,type="l")
dev.copy(png, 'Q29899-ImpulseResponse.png')
dev.off()

freq <- seq(0,pi,pi/100000)

complex_magnitude <- polyval(f1, exp(-1i*freq))

plot(freq, unwrap(Arg(complex_magnitude)), type="l");
dev.copy(png, 'Q29899-UnwrappedPhase.png')
dev.off()

plot(freq[60000:100000],Arg(complex_magnitude[60000:100000]), type="l")
lines(freq[c(60000,100000)],c(pi/2,pi/2), col="red")
lines(freq[c(60000,100000)],c(-pi/2,-pi/2), col="red")
title('Stop band phase detail')
dev.copy(png, 'Q29899-PhaseDetailWrapped.png')
dev.off()

plot(freq[1:10000],Arg(complex_magnitude[1:10000]), type="l")
lines(freq[c(1,10000)],c(pi/2,pi/2), col="red")
lines(freq[c(1,10000)],c(-pi/2,-pi/2), col="red")
title('Pass band phase detail')
dev.copy(png, 'Q29899-PhaseDetailWrappedPassband.png')
dev.off()


Linear interpolation formula In the following lecture:

http://www.ece.mcmaster.ca/~xwu/interp_1.pdf

the model (formula) for solving the linear interpolation problem (1D) given at p.5 is:

$f(x)= a_1x_1 + a_0x_0$

solve for $a1,a0$

the system:

$f(0)=a_i*0 + a_0*1$

and

$f(1)=a_i*1 + a_0*1$

I don't understand where this formula comes from. Even though deriving in all possible ways I could think of, I cannot get the equation of a line $y=mx+b$  ( also written as: $f(x)= mx+b$, where m is the slope and b is the y-intercept) to look something like $f(x)= a_1x_1 + a_0x_0$. Any idea/proposals/discussions etc.?

To my eye, the equations in the slides look like $f(x) = a_1 x^1 + a_0 x^0 = a_1 x + a_0$. So it's an exponent on the $x$, not an index. As long as you are interpolating linear, the only exponents that you encouter are 1 (linear) and 0 (constant). When you go to quadratic, you'll see $x^2$ terms appearing. Now, $a_1$ is the coefficient of the linear term (your slope $m$) and $a_0$ is the coefficient of the constant term (the y-intercept $b$).

This should explain where the system of equations for $a_1$ and $a_0$ comes from.


Detecting expected patterns in repeated measurements I am new to signal processing and I am trying to learn myself. To understand what i am interested in, let me show this on example (note that this is just a random example to illustrate). Here is some signal detected from a device:

You can notice there are 3 peaks there, and let's say there is a task to detect times of beginning of each peak (where it starts to rise). I then need to evaluate times of these 3 peaks and compare it to expected value. Now I am interested in ways how to detect these "patterns" and I am trying to find more info/learn about (something general, not specifically to this signal).

So far I learnt how to use basic filters to get rid of noise (moving average, exponential smoothing) and some really simple techniques how to find some of these points (checking sudden changes (difference between two successive filtered values) in the filtered value that is above certain threshold, etc...). My question is, where can I get to know more about this area, if you can point me to some resources, direction... Something that is not about finding unknown patterns in signal, but more about being able to accurately detect already known/expected "patterns". I will be thankful for any info/references.

So far I learnt how to use basic filters to get rid of noise (moving average, exponential smoothing)

Those are all low-pass filters!

  checking sudden changes (difference between two successive filtered values)

That is a high pass filter !

To be honest, I don't see "patterns" in your signal, the peaks are all pretty different, but you don't care for "patterns", either:

You really are just interested in the signal that first passes through a low-pass filter, and then a high-pass filter.

Combined, that is a band-pass filter.

So if I'd recommend something: Look into filtering; understand the spectrum of signals, and understand how to build filters that extract specific signals from the original signal.


Converting Gamut3D from RGB colorspace to LAB colorspace I have this code for Gamut3D in RGB, but I'm looking to convert it into LAB colorspace but I'm not really sure what changes I should be making.

rgbImage = imread('pill1.jpg');
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows columns numberOfColorBands] = size(rgbImage);
% Display the original color image.
imshow(rgbImage);
title('Original Color Image');
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Extract the individual red, green, and blue color channels.
redChannel = rgbImage(:, :, 1);
greenChannel = rgbImage(:, :, 2);
blueChannel = rgbImage(:, :, 3);
% Construct the 3D color gamut.
gamut3D = zeros(256,256,256);
for column = 1: columns
        for row = 1 : rows
                rIndex = redChannel(row, column) + 1;
                gIndex = greenChannel(row, column) + 1;
                bIndex = blueChannel(row, column) + 1;
                gamut3D(rIndex, gIndex, bIndex) = gamut3D(rIndex, gIndex, bIndex) + 1;
        end
end
% Get a list of non-zero colors so we can put it into scatter3()
% so that we can visualize the colors that are present.
r = zeros(256, 1);
g = zeros(256, 1);
b = zeros(256, 1);
nonZeroPixel = 1;
for red = 1 : 256
        for green = 1: 256
                for blue = 1: 256
                        if (gamut3D(red, green, blue) > 1)
                                % Record the RGB position of the color.
                                r(nonZeroPixel) = red;
                                g(nonZeroPixel) = green;
                                b(nonZeroPixel) = blue;
                                nonZeroPixel = nonZeroPixel + 1;
                        end
                end
        end
end
figure;
scatter3(r, g, b, 3);
xlabel('R' );
ylabel('G' );
zlabel('B' );

The so-called Gamut in the code, is nothing but a 3d color histogram. Colors from the image are collected into it, and then each unique (r,g,b) is being shown as a scatter point.

In order to show the distribution of the colors in LAB space, you can convert your image from sRGB to LAB, (Assuming that the image is in sRGB space) and then run the same procedure that you have now.

In order to convert it, use:

    cform = makecform('srgb2lab');
    labImage = applycform(cform, rgbImage);

By the way, the code in your question can be optimized into 4 lines:

    im = imread('peppers.png');
    im = reshape(im,[],3);
    im = unique(im,'rows');
    scatter3(im(:,1),im(:,2),im(:,3));


Frequency modulation of audio signals. Does increasing the pulse width of frequency tone result in larger operating distance? I have a simple frequency modulation scheme where in I am sending simple textual data by encoding it into specific frequencies(for example 'a' is converted to ASCII value of 97 and I use 10 frequencies to map decimal values from 0-9. 17 khz for 0,17.2 for 1,17.4 for 2 and so on.) A sound file is created by encoding such data,  and when I playback the tone,  I have a decoder on my phone which decodes that data. Currently each tone has a duration of 500 ms and the phone can get the data even from 10 m from the speaker with 100% accuracy. When I try to reduce the duration to 250 ms in order to increase the speed,  I find that the accuracy decreases drastically for distances >2 m.
Does this imply that if I increase the duration, there will be increase in operating distance?

It is very hard to give a single, correct answer because there are so many variables involved. Some things to keep in mind:

The main variable that determines the number of errors in the receiver is the ratio of energy per digit and the noise (this is the signal to noise ratio, or SNR). When you decrease the pulse duration, you reduce its energy and in consequence reduce the SNR. Increasing the pulse duration will increase the SNR and will allow you to increase the distance, all else equal.
Received power decreases with the square of the distance; for relatively large distances, this will dominate the SNR.
The technique you are using is called frequency shift keying (FSK). It is easy to implement but requires large bandwidth. Your audio setup has limited bandwidth (20 kHz in the most optimistic case). Reducing the pulse duration increases the bandwidth; at some point, your speaker/microphone combination will not have the bandwidth to properly reproduce and detect your signal.
Bandwidth considerations also suggest that you should shift your frequencies to around 10 kHz: at 17 kHz, you will soon hit the system's limit of 20 kHz. Operating at 10 kHz, you have 10 kHz of bandwidth to either side of your carrier.


Translation invariance in max-pooling and cascading with convolutional layer I am reading through this page about convolutional neural networks and I am confused about the part about max-pooling (downsampling a signal/image with the maximum value in a block).

Apparently, one of the advantages of max-pooling is translation invariance:

  It provides a form of translation invariance. Imagine cascading a max-pooling layer with a convolutional layer. There are 8 directions in which one can translate the input image by a single pixel. If max-pooling is done over a 2x2 region, 3 out of these 8 possible configurations will produce exactly the same output at the convolutional layer. For max-pooling over a 3x3 window, this jumps to 5/8.

I kept thinking about this and I can't figure out why this is true.

I am assuming the author means doing a translation to the left/right of the image, passing through a max-pooling layer and convolving. Then comparing the outputs of the convolution with the same process but for a different translation direction

I even did a Python script to illustrate the process of max-pooling -> convolution with 8 different translations by one pixel and the original matrix

from itertools import product

from cv2 import warpAffine
from matplotlib.pyplot import figure, gca, imshow, show, subplot
from numpy import asarray, concatenate, eye, float32, max, min, mod, reshape, zeros
from numpy.random import binomial
from scipy.signal import convolve2d
from skimage.measure import block_reduce

def list_to_2d(x):
    return reshape(asarray(x), (len(x), 1))

def npix_to_m(x):
    return concatenate((eye(2), list_to_2d(item)), axis=1)

def show_images(xl):
    vmin = min(asarray(xl))
    vmax = max(asarray(xl))
    figure(figsize=(4.6,4.6))
    for (x, i) in zip(xl, range(len(xl))):
        subplot(3, 3, i+1)
        gca().xaxis.set_major_locator(plt.NullLocator())
        gca().yaxis.set_major_locator(plt.NullLocator())
        imshow(x, interpolation='None', cmap='gray', vmin=vmin, vmax=vmax)
    show()

x = float32(binomial(3, 0.5, (10,10)))

translations = list(product((-1, 0, 1), repeat=2))

M_list = [npix_to_m(item) for item in translations]

x_trans = [warpAffine(x, M, (10,10)) for M in M_list]
show_images(x_trans)
print 'Original matrices\n'

x_maxp = [block_reduce(x, block_size=(2,2), func=max) for x in x_trans]
show_images(x_maxp)
print 'Max-pooled matrices\n'

filt = float32(binomial(3, 0.5, (2,2)))
x_conv = [convolve2d(x, filt, mode='valid') for x in x_maxp]
show_images(x_conv)
print 'Convolved max-pooled matrices\n'

These are the outputs:

As you can see none of them are exactly the same. So what am I not understanding here?

The idea of complex cell layer is to "pool" a set of simple cells and acquire the same data from those simple cells (e.g invariant translation, rotation invariant, etc). One nice illustration is the computational model of Neocognitron  https://youtu.be/Qil4kmvm2Sw?t=6m9s

As for the example of translating 2x2 and 3x3 matrix by one pixel in deeplearning.net, I think one possible explanation is that the Pooling map returns the same value for the subsequent Convolution map. Considering the following 4x4 image with our significant data in 2x2 region (row = 2, col = 2) to (3,3):

0 0 0 0
0 1 6 0
0 2 3 0
0 0 0 0

The Max-Pooling map always look at the region (2,2) and (3,3) and it returns 6. If our data were translated (moved) one pixel, there would be 8 different ways to do that. Among those there are 3 ways that keep region ((2,2); (3,3)) has 6 as the maximum value:

drow = 1, dcol = 0
0 0 0 0
0 0 0 0
0 1 6 0
0 2 3 0

drow = 0, dcol = -1
0 0 0 0
1 6 0 0
2 3 0 0
0 0 0 0

drow = 1, dcol = -1
0 0 0 0
0 0 0 0
1 6 0 0
2 3 0 0

That said, if our data is represented by the maximum value then the recognition will work for all those 3 cases.

When applying the same translation process (by 1 pixel) for 3x3 matrix

1 2 3
4 5 6
7 8 9

The the maximum value at the corner (1, 3, 7, 9) will be retained in 3 translations, the value in the middle of each edge (2, 4, 6, 8) will be retained in 5 translations. The simple cells are responsible for line/edge detection, so probably the center position (5) is left out.

That's how I understand translation invariance. I'm not completely sure though, since the author stated that "[...]produce exactly the same output at the convolutional layer." It would be great if the author of deeplearning.net could explain a bit more on that.


Sobel filter weird output I am trying to implement a sobel filter to locate horizontal edges in a 160x120 pixel image that is captured from a OV7725 CMOS camera module that outputs RGB565 data.

I've done the RGB565 -> 8-bit grayscale conversion as can be seen in the first picture

However, when I try to do my sobel operation on the grayscale picture it gets a bit funky and the output is multiplied for some reason. I've been looking at my implementation but can't seem to find out what is wrong...

Below is my implementation in C and I don't really care about the pixel at the edges.

void sobelFilter(int8_t *inputImg, int8_t *outputImg)
{
    uint16_t row;
    uint16_t col;
    uint16_t pos;

    int16_t pixel_x;
    int16_t pixel_y;

    const int8_t sobel_h[3][3] = {
            {-1, 0, 1},
            {-2, 0, 2},
            {-1, 0, 1}
    };

    for (row = 1; row < 120 - 1; row++)
    {
        for (col = 1; col < 160 - 1; col++)
        {
            pos = 120 * row + col;

            pixel_x = inputImg[pos - 121] * sobel_h[0][0];
            pixel_x += inputImg[pos - 119] * sobel_h[0][2];

            pixel_x += inputImg[pos - 1] * sobel_h[1][0];
            pixel_x += inputImg[pos + 1] * sobel_h[1][2];

            pixel_x += inputImg[pos + 119] * sobel_h[2][0];
            pixel_x += inputImg[pos + 121] * sobel_h[2][2];

            outputImg[pos] = sqrt( (pixel_x * pixel_x) );) );
        }
    }
}

Is there something funky with my implementation that I've missed or am doing wrong? Like why is it producing 3 copies of the "object" and why is the filtered image not being filtered the whole way down?

I solved it by using a threshold value for the pixel value that I got from sqrt( (pixel_x * pixel_x) );


How to select the sampling frequency when the input signal frequency is not known I am trying to observe the noise produced by UPS present in our lab, under no load (a normal running condition when its power is on). As I don't know its frequency range, I have randomly chosen 4 sampling rates to see which frequency peak remains constant. But surprisingly, peaks obtained in each graph is almost different from one another. So, in this case how to select sampling frequency, so that I can have a rough estimation of noise frequency coming from ups.

I have attached a table which shows my observations. Please have a look.

When the input signal bandwidth in unknown, you can begin with the highest possible sampling frequency, then observe the energy in the high frequency region (close to Nyquist frequency; i.e. half the sampling rate) and try to estimate the bandwidth of the signal and then choose the lowest possible rate still maintainig the Nyquist criterion according to this accepted bandwidth.

If the first trial indicates that there's considerable energy towards the Nyquist frequency, then you should try more complex sampling strategies, such as bandpass sampling at higher frequency regions, unless you can afford an even higher frequency sampling device.


Quantizing a filtered signal In the process of trying to learn dialup, I've managed to learn/figure out everything except how to convert a demodulated/filtered signal back into the original data.

Here I have a random two-bit signal (red) sampled at 16000Hz at 2400 symbols/s (gray lines separate symbols), as well as the filtered signal (blue). The filter is a raised-cosine filter with α=1200Hz, β=0.5, and a length of 27 samples.

Given the filtered signal, it's clear that a simple A/D converter is insufficient; for example, at about 7 symbols in (about a third from the left), the original signal crosses zero four times, whereas the filtered signal only crosses twice. I imagine this is due to higher frequencies being attenuated more than lower frequencies.

I realize that because the signal is being filtered at half the symbol rate, there shouldn't be any loss of information. How do I go about reconstructing the original signal from the filtered version?

In case it helps any: because the dialup training signal consists of known, pseudo-random data, should I use that to construct an adaptive filter/equalizer? Is it as "simple" as "un-attenuating" the higher frequencies?

Your low-pass filter's pass band is smaller than your signal's bandwidth, so it is destroying a significant portion of your signal.  Given that, you cannot reliably reproduce the original signal without using sophisticated techniques like error correction codes.

I would try increasing $\alpha$ to 2400 Hz.  Then the rolloff point will be 2400 * .5 = 1200.  The resulting signal should look much better.  If that filter bandwidth is too large you can find your sweet spot somewhere between 1200 and 2400 Hz.


Derivative of noisy signal My input signal is phase vector. I want to differentiate it to get frequency vector. My input signal is somewhat noisy. Here is the input signal.

This is the derivative of the input signal as calculated using

diff(inputSig)

When i differentiate this signal, I am not getting a smooth curve.
The output looks 'spikey'. I am guessing it is because of the noise in the input signal (is this 'derivative kick'?). How to avoid this and get a smooth derivative curve?

2 point discrete differentiation is bound to produce highly noisy results.
try the 5-points stencil. you can also generate coefficients (i.e. more points) yourself using derivation of Lagrange polynomials.


Amplitude of signal and PSD peak height So, suppose you have a superposition of cosines of the same magnitude:
$\sum_{i}\cos[\omega_i t]$

When you take the Fourier transform, theoretically you'd expect to find $\delta$-like peaks of the same magnitude. But in practice, what I find is this:

How would you explain this variation in the magnitude? Given such variation, would it be possible to recover the amplitudes of each frequency component in a superposition $\sum_{i}A_i\cos[\omega_i t]$

EDIT: I suspect it has to do with the relationship between the sampling frequency and the frequencies of the components. The fact that the theoretical $\delta$ response exists is a consequence of the orthogonality of the FT kernel. So, I suppose that in the discrete case, the frequencies of the DFT lines don't exactly "match" the sampled frequencies of the signal. I don't know how to articulate this properly or if it's even correct though.

For reference, here's the MATLAB code I used to generate this

Fs = 100;
Sp = 1/Fs;
L = 2000;
t = (0:L-1)*Sp;

y2 = cos(5.*t)+cos(25.*t)+cos(10.*t)+cos(1.*t)+cos(100.*t)+cos(200.*t)+cos(300.*t);

Y = fft(y2);
P2 = abs(Y/L).^2;
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;

figure
plot(f,P1,'-k','Linewidth',1);

This question is well-addressed in other posts; I would refer you elsewhere, but since you probably want the answer, here it is.

A couple of notes:

There are three knobs you need to play around with to get the right FFT. The first knob is sampling enough (so the 1/700 parameter); in particular, if you don't do it at least twice the maximum frequency, you run into problems of aliasing (by far the best and easiest to understand resource is here: https://blogs.mathworks.com/steve/category/fourier-transforms/). The second knob is that you need to sample long enough (the 50 number below) - you could do 1/1000 sampling but if you do up to 10 - 1/1000 and do the FFT, you will notice the amplitude of the smallest frequency (1/(2 * pi)) is always a little too high.
The linked blog explains VERY well why this third knob is important, which is the zero-padding, which allows you to get a better DFT with the same data.
Finally, yes I am plotting the two-sided FFT, so usually people will take half the amplitude I am using: $$\text{plot}(\text{fshift},\text{abs}(\text{yshift})/\text{length}(x))$$ instead of $$\text{plot}(\text{fshift},2 * \text{abs}(\text{yshift})/\text{length}(x))$$.

But this recreates the right amplitudes. The idea is you need to normalize by the number of data points you are using (I didn't come up with this code myself; I stole it from here: https://www.mathworks.com/matlabcentral/answers/162846-amplitude-of-signal-after-fft-operation).

Without further ado, here's the code you want:

\begin{equation}
t = 0:1/700:50-1/700; \\
x = cos(5.*t)+cos(25.*t)+cos(10.*t)+cos(1.*t)+cos(100.*t)+cos(200.*t)+cos(300.*t); \\
y = \text{fft}(x, 2^{\text{nextpow2}(\text{length}(t))+2)}); \\
n = \text{length}(y);    \\
\text{fshift} = (-n/2:n/2-1)*(700/n); \\
\text{yshift} = \text{fftshift}(y); \\
\text{figure}; \text{plot}(\text{fshift},2 * \text{abs}(\text{yshift})/\text{length}(x)) \\
\end{equation}


Does the frequency response assume the zero input response to be 0? I'm a bit confused about the frequency response and the state of a system. Is it simply the ratio of the output of the zero state response to the input of the zero state response? Or does it include the zero input response? Any clarifications would be greatly welcomed.

A frequency response describes a linear time-invariant (LTI) system. It is the Fourier transform of the system's impulse response, and only LTI systems are fully characterized by their impulse response.1

A system with non-zero initial conditions is not LTI because its output has a component which only depends on the initial conditions and which is independent of the input signal. Consequently, a system with non-zero initial conditions cannot fully be described by a frequency response. If you assume zero initial conditions, then such a system (if otherwise LTI) can be described by a frequency response.

So to answer your question, a system can indeed only be fully described by a frequency response if its zero-input response (ZIR) is zero, because if its ZIR were not zero it wouldn't be an LTI system.

Also take a look at these two answers to related questions: answer 1, answer 2.

1Also linear time-varying (LTV) systems can be characterized by an impulse response, but it is a two-dimensional function, and LTV systems don't have a frequency response in the conventional sense.


What does it mean in image processing when a filter is called non-linear? In image processing, what does it mean when a filter is called non-linear?

Does it mean the equation of the filter contains derivatives and if it didn't, it would have been called linear?

A filter F is called "linear", iff for any scalars $c_1$, $c_2$ and any images $I_1$ and $I_2$:

$F\left(c_1\cdot I_1+c_2\cdot I_2\right)=c_1\cdot F\left(I_1\right)+c_2\cdot F\left(I_2\right)$

This includes:

Derivatives
Integrals
Fourier transform
Z-Transform
Geometric transformations (rotate, translate, scale, warp)
Convolution and Correlation
the composition of any tuple of linear filters (i.e. applying some linear filter to the output of another linear filter $F(G(I))$)
the sum of the result of any two linear filters (i.e. the output of one filter, added pixel by pixel to the output of another filter $F(I) + G(I)$)

and many others.

Examples of non-linear filters are:

the square, absolute, square root, exp or logarithm of the result of any linear filter
the product of the result of any two linear filters (i.e. the output of one filter, multiplied pixel by pixel with the output of another filter $F(I)\cdot G(I)$)
morphological filters
median filter


removing vocal from audio file I want to create a program in matlab that could act like karaoke software. At this moment i have following lines of code:

clear all;

[y1, fs1] = audioread('miszmasz.wav');
[y2, fs2] = audioread('voice.wav');

%checking if there is an audio file
if y1 == 0
    return
elseif y2 == 0
    return
end

%spliting into left and righ
leftSignal = y1(:,1);
rightSignal = y1(:,2);
voiceSignal = y2(:,1);

leftSignalLength = length(leftSignal);
rightSignalLength = length(rightSignal);
voiceSignalLength = length(voiceSignal);

%furier transform of left, right canal and voice sample
NFFT = 2^nextpow2(voiceSignalLength);
leftSignalFFT = fft(leftSignal, NFFT)/voiceSignalLength;
rightSignalFFT = fft(rightSignal, NFFT)/voiceSignalLength;
voiceSignalFFT = fft(voiceSignal, NFFT)/voiceSignalLength;

for i1 = 1:leftSignalLength
    leftSignalFFT(i1) = leftSignalFFT(i1) - voiceSignalFFT(i1);
end

for i2 = 1:rightSignalLength
    rightSignalFFT(i2) = rightSignalFFT(i2) - voiceSignalFFT(i2);
end

leftOutput = ifft(leftSignalFFT);
rightOutput = ifft(rightSignalFFT);
finalOutput(:,1) = leftOutput;
finalOutput(:,2) = rightOutput;

%wavwrite(y_Out, fs1, 'output.wav');
audiowrite('miszmasz_output.wav', finalOutput, fs1);

unfortunately, it does not work. Any tips?

If my guess is right, your code show you trying remove the vocal from  second audio, the first audio must be the original vocal and the second is some mix (vocal + music)...

So if my guess is true start trying a simple approach to test, here my simple test.

I get just 4096 points of some audio voice and mix with a sine 4096 points in 200hz sampled in 1600hz, all mono.

Plot 4096 voice points:

Now plot a sine in 200 Hz:

Now the Mix (sine + voice):

Now to get the original sine without the voice you need do:

real((ifft(fft(mix) - fft(voice_Signal))))

The final plot after the step above:

course is a very simple test, and in this case we're telling you to have the original clean audio voice before being mixed ...


Why is PSD estimated, and not simply computed? I apologize if this comes as a basic question, but I am struggling to understand why the PSD can only be estimated and not directly computed. For example, this thread discusses several such PSD estimation methods.

If I have a deterministic signal with a fixed number of samples, shouldn't I be able to directly determine its spectral information?  Also, it is well known that PSD can be defined as the Fourier transform of the auto-correlation function. Isn't this calculation deterministic?

Finally, in this paper the authors compute the PSI  directly from the FFT coefficients as shown below. There doesn't seem to be anything stochastic about the computation, and the computation of FFT coefficients is not really an estimation, to the best of my understanding. What am I missing?

Well, because in a lot of real world problems this

  If I have a deterministic signal with a fixed number of samples, shouldn't I be able to directly determine its spectral information?

is just not the case. Very often, measured signals are more of a random process. A simple and common case would be to have the desired signal and some additive noise, very often Gaussian in nature.

While you can capture the signal for a while and then deterministically calculate the PSD of those samples, what you get is the instantaneous PSD of these samples which may or may not be close to the actual PSD of the whole random process that you are sampling. It all boils down to what you want to do but in a lot of cases you want the PSD of the random process that you are sampling and also very often, just capturing some samples and calculating the instantaneous PSD just does not come close to it.

There is a lot of different methods for PSD estimation. Some may work better in some cases than other, maybe some are more computationally expensive than others. Or possibly some reach a better estimate with fewer samples as input data. In the end it boils down to select the most appropriate method for your given scenario.


What is the time unit of the converted Channel Impulse Response? I don't have a lot of background on wireless communication but as far as I understand, applying IFFT on Channel Frequency Response (CFR) would transform the data into time domain Channel Impulse Response (CIR), in which it shows the delay spread of the different multiple paths.

The CSI measuring tool I'm using measures 30 subcarriers resulting in 30 complex number CFRs per a single CSI trace. So when I apply 30-point IFFT on the CFR, the output is a CIR with 30 time-domain values.

My question is,

What is the time unit of the converted CIR?
How much seconds are they (the 30 impulse responses) apart?
Also, are the impulses like the echoes of the transmitted signal?
The number of multi-path is not determined. So I assume that an impulse response does not represent a single path, since the number of impulse response is fixed to 30. Then how should I understand the CIR value with respect to the multi-paths?

Because you don't have background in wireless communications, I will try to answer as simply as possible.

  What is the time unit of the converted CIR?

Just time unit. It can be second, ms, us, etc. They are convertible, aren't they?

  How much seconds are they (the 30 impulse responses) apart?

In the case that CSI-pilots are equally seperated, impulses are $1/(N_p \times \Delta p)$ second apart, where $N_p$ is the number of CSI-pilots, e.g. $N_p = 30$ in your example, and $\Delta p$ is CSI-pilot spacing in Hz.

  Also, are the impulses like the echoes of the transmitted signal?

Impulse reponse is the output if we use a Dirac pulse as input. It represents the effect of scaled-delayed echoes of the TX signal.

  The number of multi-path is not determined. So I assume that an
  impulse response does not represent a single path, since the number of
  impulse response is fixed to 30. Then how should I understand the CIR
  value with respect to the multi-paths?

The generated CIR is the interpolated version of the multipath response. For example, in the figure below, the CIR samples are green and physical paths are in other colors.

Note that if you accept CIR by using iDFT like that, you are implying that your multipath response is invariant during the acquisition time of your 30 CSI-pilots.

Some other resources if you want to learn more or to find the answers by yourself if I have missed something. :)

Multipath fading channel visualization in MATLAB

Statistical Characteristics of Multipath Channels


Cross-correlation or cross-covariance of non-zero mean signals Cross-correlation for uniformly sampled signals is defined as [1]

$$(f \star g)[n]\ \stackrel{\mathrm{def}}{=} \sum_{m=-\infty}^{\infty} f^*[m]\ g[m+n].$$

Cross-covariance for wide-sense stationary (WSS) signals is defined identically [2]. This seems strange to me, because the statistical definition [3] first subtracts the signal means from the signals before correlating it with each other. In other words, cross-covariance is substantially different from cross-correlation if the mean of the signals is substantially different from zero (is the Wikipedia article wrong? Or does it somehow assume zero mean signals?).

So, cross-covariance equals cross-correlation if the signal means are zero, else it is unequal.

Questions:

What are reasons to choose for cross-correlation or cross-covariance when comparing signals with non-zero mean?
Does cross-correlation have any significance in this case, because $f\star g$ has mean $\mu_f\mu_g$?
How would one find out the cross-correlation of signals that are not WSS, e.g. the mean is a linear trend? What would happen if we'd ignore this linear trend and simply subtract the signal means, then calculate cross-correlation? What if we subtracted the linear trend instead of the mean? What are the consequences for such mild violations of WSS-ity?

What are reasons to choose for cross-correlation or cross-covariance when comparing signals with non-zero mean?

Well, part of the issue is that cross-correlation as defined in your equation:

$$(f \star g)[n]\ \stackrel{\mathrm{def}}{=} \sum_{m=-\infty}^{\infty} f^*[m]\ g[m+n].$$

will not exist (or be infinite) if $f$ and $g$ have non-zero mean. So, in order to do the calculation, you really should subtract the mean beforehand.

If you're dealing with finite length signals, it's even worse, because you end up with an underlying triangular offset as per the image below.

This shows the sample cross-correlation between two gaussian processes with non-zero mean.  The black line is the sample cross-correlation calculated without subtracting the mean.  The red line is the sample cross-correlation calculated by first subtracting the mean from both signals.

  Does cross-correlation have any significance in this case, because $f\star g$ has mean $\mu_f\mu_g$?

It doesn't generally, unless you're very interested in the mean value.

And, as I said, as per your definition, any value of $\mu_f$ and $\mu_g$ both different from zero will result in an infinite sample cross-correlation.

  How would one find out the cross-correlation of signals that are not WSS, e.g. the mean is a linear trend? What would happen if we'd ignore this linear trend and simply subtract the signal means, then calculate cross-correlation? What if we subtracted the linear trend instead of the mean? What are the consequences for such mild violations of WSS-ity?

When looking at linear trends, it's usual to detrend before looking at the statistics of the variation around the trend. Sometimes doing this allows the variations to be treated as WSS.

If you keep the trend, you'll end up with a sample cross-correlation that looks like a parabola (the integration of a linear trend).  In this case, simply removing the mean doesn't help. See similar plot to first below, this time with two linear trend noisy signals.

Please see this exposition on the Mathworks site for a Matlab version of the above / below.

R Code Only Below

#34778
#install.packages('SynchWave')
require(SynchWave)

with_mean_ccf <- function(f,g)
{
  length <- length(f) + length(g) - 1

  ff <- c(f, rep(0,length - length(f) + 1))
  gg <- c(f, rep(0,length - length(g) + 1))

  FF <- fft(ff)
  GG <- fft(gg)

  ccf <- fft(FF * Conj(GG), inverse = TRUE)

  return(fftshift(Re(ccf)))
}

T <- 1000
mu_f <- 1
mu_g <- 2

f <- rnorm(T,1) + mu_f
g <- rnorm(T,1) + mu_g

with_mean <- with_mean_ccf(f,g)
no_mean <- with_mean_ccf(f-mean(f), g-mean(g))
plot(with_mean)
points(no_mean, col='red')
title('Sample crosscorrelation with (red) and without (black) mean subtraction')

a <- rnorm(T,1) + seq(1,T)/3
b <- rnorm(T,1) + seq(1,T)/2

with_mean_and_trend <- with_mean_ccf(a,b)
no_mean_and_trend <- with_mean_ccf(a-mean(a),b-mean(b))

plot(with_mean_and_trend, ylim=c(min(with_mean_and_trend,no_mean_and_trend), max(with_mean_and_trend,no_mean_and_trend)) )
points(no_mean_and_trend, col="red")
title('Sample cross-correlation of two linear trend signals')


SFBC Alamouti Scheme LTE A sentence I copied from book.

That a 2$\times$1 Alamouti SFBC is referred to as a rate 1 code, two symbols are sent over two adjacent subcarriers.

This question is regarding this SFBC encoder (Alamouti scheme). I understand how the encoder works for the two transmit antenna case. Below I briefly given an example for convenience.

Assuming a BPSK modulator applied to data input $I$. Then we have

$I$ = 0 0 1 1 0 1 1 1

BPSK symbols=  1 1 -1 -1 1 -1 -1 -1

Then the output of the STBC encoder could be found according to the

\begin{array}{lcr}
\mbox{$$} & Antenna1 & Antenna 2 \\
\mbox{$f_1$} & s_1 & s_2 \\
\mbox{$f_2$} & -s_2^* & -s_1^* \\
 \end{array}

and so applying this to the BPSK symbols I get the following symbols

Antenna1 symbols: 1 -1 -1 1 1 1 -1 1

Antenna2 symbols: 1 1 -1 -1 -1 1 -1 -1

My question is next, what do we do with these symbols for LTE system ie OFDM modulation, do we send the symbols 1 and 2 over antenna 1 and antenna 2 respectively over the same subcarrier ? i.e I don't understand how transmission over subcarriers takes place? More precisely I don't undersatnd the assignement of the symbols to the different subcarriers?

Thank you.

The Alamouti scheme guards against Rayleigh fading channels.  If you were to send two symbols within the same channel within the same coherence time interval, they will experience the same fading.  In order to get diversity gain, you need to send them in different channels (different antennas or different frequencies or distant time slots)


What is meant by "Cubic Difference Frequency"? The cubic difference frequency (or 3rd degree intermodulation product) is given by:

$$D_3 = - 10\log\left(\frac{P_2\left(2 f_{1,2} - f_{2,1}\right)}{P_2\left(f_1,f_2\right)}\right)$$

where $P_2$ is the output power. As far as I understand, the 3rd degree IM products are at frequencies $2f_1 - f_2$ and $2f_2-f_1$. So the formula for the IM product of $2f_1 - f_2$ should be:

$$D_3 = - 10\log\left(\frac{P_2\left(2 f_{1} - f_{2}\right)}{P_2\left(f_1\right)}\right)$$

When the power is given in an logarithmic scale, this would simply be the difference between the two power values.

But what does this actually tell me?
Is this the power of my 3rd degree intermodulation product?
(if it is, $D_3$ should be equal for both frequencies)
And is this formula only valid for an input signal of two different frequencies?

According to a paper D3 can also be expressed in terms of gain compression c and change in AM-PM characteristic kp due to only one carrier:
$$D_3 = -10\log(\frac{c_1}{2}^2 + k_{p1}^2)$$

The gain compression is given by $c = 1 - \frac{\Delta P_2 / P_2}{\Delta P_1 / P1}$
, where $\Delta P_1$ is a constant change of input power and $\Delta P_2$ is a not constant change of output power. I am not sure however, if $P_1$ or $P_2$ refers to the input/output power before the change or after.

After skimming through the paper, I can see more clearly now. The measure $D_3$ quantifies the relative strength of the 3rd order intermodulation product. If two sinusoidal signals with frequencies $f_1$ and $f_2$ are input to a non-linear device (such as a microwave amplifier, as referred to in the paper), there will be intermodulation products at the output. Third-order intermodulation products occur at frequencies $2f_1+f_2$, $2f_1-f_2$, $2f_2-f_1$, and $2f_2+f_1$. Problematic are the ones at frequencies $2f_1-f_2$ and $2f_2-f_1$ because they can end up in the frequency band of interest.

The quantity $D_3$ is defined as the ratio of the power of one of these undesired frequency components and the power of one of the desired frequencies (in dB) at the output of the nonlinear device. There are four possible definitions:

$$D_{3,1}=-10\log\left(\frac{P_2(2f_1-f_2)}{P_2(f_1)}\right)\\
D_{3,2}=-10\log\left(\frac{P_2(2f_1-f_2)}{P_2(f_2)}\right)\\
D_{3,3}=-10\log\left(\frac{P_2(2f_2-f_1)}{P_2(f_1)}\right)\\
D_{3,4}=-10\log\left(\frac{P_2(2f_2-f_1)}{P_2(f_2)}\right)$$

where $P_2(f_x)$ is the output power at frequency $f_x$. These four definitions do not necessarily result in the same value for $D_3$. It's just important to specify which of them is used, and to consistently stick to that definition.

To answer your questions: no, $D_3$ is not the power of the 3rd order intermodulation product, but it is the ratio (in dB) of that power to the power at one of the signal frequencies. This quantity is defined based on the assumption that the input signal is the sum of two sinusoidal signals, but it is used to quantify the strength of 3rd order nonlinearities at the output of a nonlinear system regardless of the input signal.

Concerning the power compression factor, you can view the changes $\Delta P_1$ and $\Delta P_2$ as infinitesimal changes, and so the question if $P_1$ and $P_2$ are the values before or after the change becomes irrelevant. At a given frequency, an ideal linear amplifier would satisfy

$$P_2=g\cdot P_1$$

where $g$ is the gain factor at the given frequency. So you have

$$g=\frac{P_2}{P_1}=\frac{dP_2}{dP_1}\Longrightarrow\frac{dP_2}{dP_1}\frac{P_1}{P_2}=1$$

That's why the power compression factor

$$c=1-\frac{dP_2}{dP_1}\frac{P_1}{P_2}$$

is a measure for the nonlinearity of power amplification.


what are the advanatages of using a scrambler What is the purpose of using a scrambler in communication systems?
The link wiki_for_scrambling says that it is used to randomize data. But the data is already random (like speech,audio,video,etc..), what is the purpose of randomizing again?

Even if we ex-or the sequence for randomizing? does it guarantee a non-repeating pattern?

And regarding the statement "Scrambler removes dependency of a signal's power spectrum on the actual transmitted data." in the article, what is the use of that?.

regards,
phani tej

Scramblers are indeed used to "randomize" the data, i.e. to "whiten" the input sequence. The reason why this is necessary is that a transmission system has no control over the data the user is going to transmit. This causes problems because it violates the assumptions that are usually made when designing transmission systems, such as having independent data symbols. Note that usually the data are not random, as you seem to believe. A frequently occurring problem are long strings of zeros in the data, which can cause difficulties in timing recovery and adaptive equalization. These problematic sequences are removed (or, actually, made much less likely) by the scrambler.


If the mean of a random process is constant, does it imply the process is first order stationary? If a random process is first order stationary, its mean is constant. However, if a random process has a constant mean say $3$ and an autocorrelation equal to $9 + 15e^{|-\tau|}$. The process is clearly wide sense stationary. Is the process first order stationary?

In the usual sense of the term, first-order stationarity means that the first-order distribution of all the random variables is the same: each $X_t$ has the same CDF, and so the same pdf (or pmf) too if the random variables are continuous (or discrete). If the random variables have a mean, then they all have the same mean. But, a sequence of independent
identically distributed Cauchy random variables is a first-order stationary process (as well as a strictly stationary process for that matter) but it cannot be said to have constant mean because none of the variables have a mean at all.

Since your response to my query on the problem statement is that your definition of first-order stationarity is the  standard one, then the answer to the question in the title of your question, viz.

  If the mean of a random process is constant, does it imply that the
  process is first-order stationary?

is NO, and an example of a process that has constant mean but is not
first-order stationary is a sequence of constant-mean random variables whose variances depend on time.

The answer to the question in the text, viz.

  Is a WSS process necessarily first-order stationary?

ia also NO.  Consider the random process
$\{X(t)\colon X(t)= \cos (t + \Theta), -\infty < t < \infty\}$
where $\Theta$ takes on four equally likely values $0, \pi/2, \pi$ and $3\pi/2$ and so $X(t)$ takes on four equally likely values
$\cos(t)$, $\cos(t+\frac{\pi}{2}) = -\sin(t)$,
$\cos(t+\pi) = -\cos(t)$, and
$\cos(t+\frac{3\pi}{2}) = \sin(t)$.
  It is easy to see that in general $X(t)$ and $X(s)$ have different distributions and so the process is not first-order stationary. On the other hand, $E[X(t)] = 0$ for all $t$ while
\begin{align}
E[X(t)X(s)]&= \left.\left.\frac 14\right[2\cos(t)\cos(s) + 2\cos(t+\pi/2)\cos(s+\pi/2)\right]\\
&= \left.\left.\frac 12\right[\cos(t)\cos(s) + \sin(t)\sin(s)\right]\\
&= \frac 12 \cos(t-s)
\end{align}
and so the process is wide sense stationary, but its first-order distribution does depend the choice of $t$ and so the process is not
first-order stationary.

Note: the above proof uses a very commonly accepted definition of wide sense stationarity as found in numerous books, online resources, Wikipedia etc, viz. a random process $\{X(t)\}$ is called a wide sense stationary (WSS) process if

$E[X(t)]$ is a constant, and
$R_{X}(t,s) = E[X(t)X(s)]$ depends only on the difference $t-s$ of the two arguments and not on the individual arguments $t$ and $s$.


2D Image Convolution: Spatial Domain vs. Frequency Domain Convolution in the Computational Complexity Sense I want to convolve 2 images I1 and I2. Let's assume they have the same size.
I know 2 methods:
1) Use direct definition of 2D convolution
2)  i. Computer FFT of I1 (call it I1_FFT) and FFT of I2 (call it I2_FFT)
   ii. Multiply I1_FFT and I2_FFT (call it I1I2_FFT)
  iii. Compute inverse FFT of result I1I2_FFT

I want to know the number of multiplications and additions for each case.  Note I don't want the Big O type of complexity.

It would be great if the answer includes a derivation.

In the end, I want to know for a specific NxM sized image, when I should use method 1 or 2.

Vanilla implementation of each method for image of size m x n and kernel of size k x l will yield:

Spatial Domain Convolution - O(mnkl) as for each pixel in the image we do kl multiplications (Additions are discarded).
Frequency Domain Convolution - O(mn log(mn) + mn) as the complexity of the FFT is mn log(mn) and we add the multiplication (You could add factor of 4 for doing it twice and back and forth).

But in practice there many nuances:

Many times the kernel is Separable or can be approximated by 1-2 separable filters (For instance, Gaussian Filter).
If the kernel is very small compared to the image, better use direct method.
If you apply the same filter on many images you may even consider using Matrix Form (As done in CNN's).
In Frequency Domain you apply Convolution with Circular Boundary Condition. It means that if you're after different boundary conditions you'll need to pad and then complexity is higher and many memory operations are done.

I can go on and on...

You may have a look at Filtering n×n Images by Separable m×m Filters. Computation Time for Filtering Using FFT, 2D Convolution and Two 1D Convolutions.

In order to say which is better, one must be very specific about the case.
In most cases, use efficient spatial domain code. This is my choice in most cases.


Unexpected results for fourier resampling (Python 3) I am trying to resample a signal using the fourier method in python using np.fft.rfft() and np.fft.irfft(). (I'm not using scipy.signal.resample() because in my final application I need to do more to my signal in the frequency domain than what can be accomplished by the window argument of scipy.signal.resample().) I have reduced the problems I see to this simple test case (ipython notebook):

%pylab inline

N = 16
P = 8

S = [0.49, 1.68, 0.78, 0.05, 0.21, 0.29, 1.34, 1.17, 16.73,
     48.78, 16.90, 0.62, 0.40, 1.60, 0.86, 0.57]
assert len(S) == N

figure(figsize=(14, 5))
plot(np.arange(N), S, '.')

dft = np.fft.rfft(S)
plot(np.arange(N), np.fft.irfft(dft), '-')

dft *= P
dft.resize(P*N/2+1)
plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')

xlim(6.5, 11.5)
ylim(-5, 55);

And here is what I get back:

As you can see, the result looks almost correct, but is nevertheless a bit off. (The resampled red curve does not go through the original sample points, even though there are sample points at the same positions in time as in the original signal.)

The greeen line shows that everything is still ok before extending and scaling the dft array. So the problem seems to be in the actual resampling:

dft *= P
dft.resize(P*N/2+1)
plot(np.arange(P*N)/P, np.fft.irfft(dft), '-')

I don't see how this result is mathematically possible because other than scaling of time and amplitude axis, adding zeroes to the high-frequency end of the dft array should only add zero terms to the fourier series expansion of the signal. So the signal should still pass through the original sample values.

But the small error I see does not seem to be related to scaling: sometimes the resampled waveform lags behind the original (t=8) and sometimes it is faster (t=10), sometimes it is larger in magnitude (t=7,9,11) and sometimes it is smaller (t=8,10).

I'm playing with this for hours now and can't find what is wrong with my code.. Maybe I just don't see it because I'm already staring at it for so long now..

(see np.fft documentation for the definitions used by the python FFT implementations)

Edit (in response to comment): When I replace S with an impulse at t=9, or a fast impulse train, I get the following:

If you perform the full DFT of a signal with odd number of samples, the frequency values in the spectrum are a DC coefficient and (N-1)/2 conjugate pairs of sinusoid coefficients.

If you perform the full DFT of a signal with even number of samples, the frequency values in the spectrum are a DC coefficient, (N-1)/2 conjugate pairs of sinusoid coefficients, and a single coefficient corresponding to frequency 0.5.

rfft only keeps the positive frequency part of the spectrum, and the inverse, irfft, presumably recreates the negative part of the spectrum before performing an inverse-DFT.

For example, for an odd size signal with rfft of [0 1 2 3], the full spectrum would be [3 2 1 0 1 2 3].

For an even size signal with rfft of [0 1 2 3 4] it presumably has a flag set somewhere in the code that indicates an even sized signal, and so the full spectrum would be [3 2 1 0 1 2 3 4]. Note the 4 only appears once.

When you pad the rfft result though, e.g. [0 1 2 3 4 0 0 0 0], it probably recreates the spectrum as [0 0 0 4 3 2 1 0 1 2 3 4 0 0 0 0]. That is, the 4 is now repeated twice, and therefore that frequency component is twice as strong.

Therefore, before padding you must halve the last rfft coefficient iff the signal has an even number of samples.


simulink inverse notch filter Is there a way to make an inverse notch filter block in matlab simulink?
I have found the peak-notch block, but I need to amplify the signal instead of attenuate it.

Thank you.

peak-notch block is a two-way component. You can set the response type to peak or notch. Here is the properties page:

int the link: http://www.mathworks.com/help/dsp/ref/peaknotchfilter.html

Using this filter you can obtain desired output by summing input signal and output of the peak filter.