maximum contrast using linear gery-level trasformation this was a question i found online.

given an image with values between [0,1].
how can I get the maximum contrast using linear gray-level transformation?
something like I(x,y)=A*I(x,y)+B , where I is a histogram.
we are also given MAX(I) and MIN(I).

as far as i can understand.
MAX(I) and MIN(I) will give you the values you have now and you wan to stretch the values into 0 and 1.

but how will this give me a linear equation?

If you're subtracting the minimum and stretching the maximum remaining value to be 1, you get the following:

Let $m = \min\{I(x,y)\}$ and $M = \max\{I(x,y)\}$
$$I_{\text{new}}(x,y) = \frac{I(x,y) - m}{M-m} = AI(x,y)+B $$

where $A=\dfrac{1}{M-m}$ and $B = -\dfrac{m}{M-m}$

Hope this clears it up.


Even and odd signal energy property In Signals and Systems by A. V. Oppenheim, A. S. Willsky, S. Hamid Nawab, 2nd Edition, and Signals and Systems, Simon Haykins, Barry Van Veen, 2nd Edition there is a problem related to energy of real-valued even and odd signal.

Energy of an arbitrary real-valued signal $x(t)$ is equal to the sum of the energy of the even component $x_{e}(t)$ and energy of the odd component $x_{o}(t)$ i.e.
\begin{equation}
 \int_{-\infty}^{\infty}x^2(t)dt = \int_{-\infty}^{\infty}x_{e}^2(t)dt +
 \int_{-\infty}^{\infty}x_{o}^2(t)dt\notag
\end{equation}
Which can be proved very easily by expressing $x(t) = x_{e}(t) + x_{o}(t)$ and evaluating its energy as
\begin{align}
 \int_{-\infty}^{\infty}x^2(t)dt &=
 \int_{-\infty}^{\infty}\big(x_{e}(t) + x_{o}(t)\big)^2dt\notag\\
 &= \int_{-\infty}^{\infty}x_{e}^{2}(t)dt + \int_{-\infty}^{\infty}x_{o}^{2}(t)dt
 +2\int_{-\infty}^{\infty}x_{e}(t)x_{o}(t)dt\notag.
\end{align}
The term $x_{e}(t)x_{o}(t)$ in the last integration corresponds to an odd signal (multiplication of even and odd signal results into odd signal) and its integration (joint energy) is equal to zero (area under an odd signal over symmetrical limits is zero). Hence
\begin{align}
 \int_{-\infty}^{\infty}x^2(t)dt =
 \int_{-\infty}^{\infty}x_{e}^{2}(t)dt +
    \int_{-\infty}^{\infty}x_{o}^{2}(t)dt\notag
\end{align}

I was wondering, what if I consider $x(t)$ to be a complex signal, will above property still be valid. I tried as
\begin{equation}
 \int_{-\infty}^{\infty}|x(t)|^2dt = \int_{-\infty}^{\infty}|x_{e}(t)|^2dt + \int_{-\infty}^{\infty}|x_{o}(t)|^2dt\notag
\end{equation}
Expressing $x(t) = x_{e}(t) + x_{o}(t)$ and evaluating its energy as
\begin{align}
 \int_{-\infty}^{\infty}|x(t)|^2dt &= \int_{-\infty}^{\infty}\big|x_{e}(t) + x_{o}(t)\big|^2dt\notag\\
 &= \int_{-\infty}^{\infty}|x_{e}(t)|^{2}dt +
 \int_{-\infty}^{\infty}|x_{o}(t)|^{2}dt
 +\int_{-\infty}^{\infty}[x_{e}(t)x_{o}^{*}(t) + x_{o}(t)x_{e}^{*}(t)]dt\notag
\end{align}

On expanding the term
\begin{align}
 x_{e}(t)x_{o}^{*}(t) &= \big(x_{er}(t) + jx_{ei}(t)\big)\big(x_{or}(t) - jx_{oi}(t)\big)\\
&= x_{er}(t)x_{or}(t) + x_{ei}(t)x_{oi}(t) - jx_{er}(t)x_{oi}(t)  + jx_{ei}(t)x_{or}(t)\\
x_{o}(t)x_{e}^{*}(t) &= \big(x_{or}(t) + jx_{oi}(t)\big)\big(x_{er}(t) - jx_{ei}(t)\big)\\
&= x_{er}(t)x_{or}(t) + x_{ei}(t)x_{oi}(t) - jx_{ei}(t)x_{or}(t)  + jx_{er}(t)x_{oi}(t)
\end{align}
Further on adding this terms results into
\begin{align}
 \int_{-\infty}^{\infty}|x(t)|^2dt &= \int_{-\infty}^{\infty}|x_{e}(t)|^{2}dt + \int_{-\infty}^{\infty}|x_{o}(t)|^{2}dt
 +2\int_{-\infty}^{\infty}\big(x_{er}(t)x_{or}(t) + x_{oi}(t)x_{ei}(t)\big)dt\notag
\end{align}
The terms $x_{er}(t)x_{or}(t)$ and $x_{oi}(t)x_{ei}(t)$ in the last integration corresponds to an odd signal and its integration (joint energy) is equal to zero. Hence
\begin{equation}
 \int_{-\infty}^{\infty}|x(t)|^2dt = \int_{-\infty}^{\infty}|x_{e}(t)|^2dt + \int_{-\infty}^{\infty}|x_{o}(t)|^2dt\notag
\end{equation}
I would appreciate if anyone here can validate the proof. Is it correct to say that the same property too holds for complex signal.

Your result is correct but note that for complex signals, the even and odd parts are defined by

$$x_e(t)=\frac12\left[x(t)+x^*(-t)\right]\tag{1}$$
and
$$x_o(t)=\frac12\left[x(t)-x^*(-t)\right]\tag{2}$$

where $^*$ denotes complex conjugation.

From $(1)$ and $(2)$ it follows that the real part of $x_e(t)$ is even and its imaginary part is odd, whereas $x_o(t)$ has an odd real part and an even imaginary part.

We get

$$|x(t)|^2=|x_e(t)+x_o(t)|^2=|x_e(t)|^2+|x_o(t)|^2+2\,\text{Re}\{x_e(t)x_o^*(t)\}\tag{3}$$

From

$$\begin{align}x_e(t)x_o^*(t)&=\frac14\left[x(t)+x^*(-t)\right]\left[x^*(t)-x(-t)\right]\\&=\frac14\left[|x(t)|^2-x(t)x(-t)+x^*(t)x^*(-t)-|x(-t)|^2\right]\end{align}\tag{4}$$

we see that $x_e(t)x_o^*(t)$ is an odd signal according to definition $(2)$ (it is the odd part of $\frac12 \left[|x(t)|^2-x(t)x(-t)\right]$), and, consequently, its real part is odd. This implies that the integral over the last term in $(3)$ is zero, and the stated property holds.


LTI system response to periodic input I'm trying to find the zero-state response (ZSR) of an LTI system to a one sided periodic input, like a square wave that is equals to zero for $t < 0$.

I know that I can use the Fourier series of said input function to find the steady-state (SS) response, however I'm having trouble understanding how to use the Laplace transform to obtain the ZSR, which contains the SS component plus a transient one.

My guess is that I need to calculate the Laplace transform of the periodic input, and then solve the system for $s$ to obtain the transform of the output.

You can use the Laplace transform, but can also simply use convolution in the time domain. In any case, you will need the system's impulse response $h(t)$. Let the input signal $x(t)$ satisfy $x(t)=0$ for $t<0$, and $x(t+T)=x(t)$ for $t>0$ and $T>0$, as required. Furthermore, let $f(t)$ be the first period of $x(t)$:

$$f(t)=\begin{cases}x(t),&0<t<T\\0,&\text{otherwise}\end{cases}\tag{1}$$

Then

$$x(t)=\sum_{n=0}^{\infty}f(t-nT)=f(t)\star\sum_{n=0}^{\infty}\delta(t-nT)\tag{2}$$

where $\star$ denotes convolution, and $\delta(t)$ is the Dirac delta impulse. If $g(t)$ denotes the convolution $h(t)\star f(t)$, i.e., the system's response to $f(t)$, then the output signal can be written as

$$\begin{align}y(t)&=h(t)\star x(t)\\&=h(t)\star f(t)\star\sum_{n=0}^{\infty}\delta(t-nT)\\&=g(t)\star \sum_{n=0}^{\infty}\delta(t-nT)\\&=\sum_{n=0}^{\infty}g(t-nT)\tag{3}\end{align} $$

According to $(3)$, the output signal can be written as a sum of shifted responses to the finite length signal $f(t)$, which corresponds to the first period of the input signal.

The Laplace transform of the output signal can also be written in terms of the transforms of the input signal and of the function $f(t)$. Note that

$$f(t)=x(t)-x(t-T)\tag{4}$$

The Laplace transform of $(4)$ is

$$F(s)=X(s)(1-e^{-st})\tag{5}$$

If $H(s)$ is the system's transfer function, i.e. the Laplace transform of the impulse response $h(t)$, then the output is given by

$$Y(s)=H(s)X(s)=H(s)\frac{F(s)}{1-e^{-st}}\tag{6}$$


Understanding why (at least) 2 points of the FFT of a real signal are real I am aware that the FFT of a real signal is Hermitian, i.e. $\text{FFT}(f)[i] = \text{FFT}(f)[-i]^*$, where $\text{FFT}(f)[i]$ is the i-th component of the FFT of $f$.
The logical consequence of this is that $\text{FFT}(f)[0]$ is real. Another way of seeing this is that the coefficient with index 0 is the "DC" component, or the "mean", and that this is of course real for a real signal.
Another consequence of course is that one needs to know only "half" the FFT to know it in full (i.e., it should be enough to know the coefficients number $0..N/2$, where N is the number of samples).
I was looking for a "real FFT" algorithm for my microcontroller and found the "official" ARM FFT. Of course they do it effectively, using an output of size N real numbers rather than N complex numbers, and outputing "only half the complex FFT" instead of the full (redundantly hermitian) FFT. This saves 50% space.
However I was surprised to see (in the discussion here https://github.com/ARM-software/CMSIS_5/issues/1091 ) that the FFT of the real signal actually has 2 purely real output values (at least, this is how they pack the output): the coefficient of index 0 (as expected) and also the coefficient of index N/2. I do not understand why the coefficient of index N/2 is real. Any hint why?
Details from the discussion there: "Except the first complex number that contains the two real numbers
X[0] and X[N/2] all the data is complex. In other words, the first complex sample
contains two real values packed."

As pointed out in comments, what I am interested in is actually some DFT properties. So that should be something like ($l$ and $k$ the DFT and function indexes, [] the "index slicing"):
$ DFT(f)[l] = \frac{1}{N} \sum_{k=0}^{N-1} f[k] e^{i 2 \pi k l / N}$

First, "why (at least) 2 points of the FFT is real" is incorrect.  It is why one or two points of the FFT is real.
Each point of a discrete Fourier transform (which the FFT is just an algorithm to compute quickly) is the complex coefficient for its corresponding frequency.  In other words, the DFT starts with the assertion that the signal in question can be fully described by
$$x_k = \sum_N a_k e^{\frac{2\pi}{N}k} \tag 1$$
For $x_k \in \mathbb R$, (1) just carries the constraint that $a_k = (a_{N - k})^*$, where $x^*$ is complex conjugation.
From this, we can show that $a_0$ (and $a_{N/2}$ for $N$ even) are constrained to be real:
Start with the rule I mentioned that in order for $x_k$ to be real, $a_k = (a_{N - k})^*$.  Another property of the DFT is that it's only uniquely defined over an interval $N$ points long, and then it repeats.  So $a_{k + N} = a_k$.  This means that for $k=0$, $a_{N - k} = a_{N} = a_0$.  Applying the conjugate property, this means that $a_0 = a_0^*$ -- so $a_0$ must be real.
For $N$ even, $a_{N/2} = a_{N - N/2}$.  So, again, $a_{N/2} = a_{N/2}^*$, so $a_{N/2}$ must be real.
So -- one (or two, in the case of an even number of points) arguments must be real.


Problem with 1st order Massberg LPF I'm trying to implement this filter in title by following the book "Designing Audio Effect Plug-Ins in C++" By Will Pirkle

Problem I'm facing is the magnitude response error when comparing against analog model of the LP filter. Here is a plot showing the difference.

I get these are values for variables (fc=1000, fs=44100):

>> massberg_test

    w0 =  0.14248
    g1 =  0.045305
    gm =  0.70711
    wm =  6283.2
    Om =  0.071359
    Os =  0.10081
    g0 =  1.1008
    a0 =  0.14612
    a1 =  0.055508
    b1 = -0.89919
Real:
    a = [ 0.132736  0.050424 ]
    b = [ 1.00000  -0.81684 ]

I don't have access to the original paper by M. Massberg so, I can't be sure if the equations in my source book is correct (I've read that the "ac" in g1 formula is just a typo and book errata does not mention other errors for this fileter).

Q:
1. Are the equations found in Pirkle book correct?
2. If 1 is Yes then where to look the issue from?

Problem solved.

Looks like calculation of Os ($\omega_s$) was the culprit for this problem. Proper equation is:

Os=Om*sqrt((gm^2-g1^2)*(1-gm^2))/(1-gm^2);

instead of:

Os=Om*(sqrt((gm^2-g1^2)*(1-gm^2))/(1-gm^2));


How do the English pronunciations of "tiu" and "tou" differ? I'm creating an app (I'm a software developer) that for different sounds (interjections from the user) will do different things.
For example if the user says:
"Tiu", than it will do task A,
if it will say "Tou" than it will do task B.

In digital sound, what's the difference between sound (user interjection) "Tiu" and "Tou"?

I need a fast way to say when the user said "Tiu" and when it said "Tou", I don't want to use word/speech recognition because it is too slow for my needs.

I'm mainly interested in the theoretical part of digital sounds. E.g. how can I put in the "Tiu" and "Tou" the following sounds "Ti Tiau Tiuu Tiiiu Too Toou Taou", what are their common "features".

(1) edit

After reading the answers, I must confess that I'm new to speech recognition and any similar domain. I didn't new in what I'm getting myself with this app. But overall the research is fun because I love algorithms and math.

Where I am now?

Currently I am reading/researching/learning about
CMUSphinx Open Source Toolkit For Speech Recognition http://cmusphinx.sourceforge.net/wiki/tutorial. I don't think that I will use this toolkit but it has a good explanation of the basis in Speech Recognition: Basics of Speech.

Also after reading the responses, I ended to the conclusion that a training phase might be quite helpful, anyway I'm investigating more on it.

New Question:

Do you now any good reading material that might help me in gaining more knowledge in this particular domain?

Can the recording device help me in identifying "tou" and "tiu"? e.g. if the device supports different recording filters - like low frequency filters etc...?

First of all, don't expect an answer like "Compute such and such well-known function on the data and if the result is above such or such value it's tiuuu otherwise it's tou". You won't get anything robust or reliable with approaches like that.

The good representations of sound which can be used for this kind of discrimination task are either the Mel Frequency Cepstral Coefficients (MFCC) or the linear prediction coefficients (Compute the first 20 or so autocorrelation coefficients and use the Levinson-Durbin recursion to get the AR coefficients). Both are somehow independent of pitch (good), moderately robust to noise (good but can be better), but are still somewhat speaker and gender dependent. It's still a vast improvement over trying to guess formants from a spectrum and matching that to a per-phoneme table.

Does your application tolerate a training phase? If so, get your user to record a dozen of utterances of each word ; and use this as a training set for a simple bayesian classifier using a gaussian model. Use this classifier to discriminate new utterances.

Does your application need realtime processing? If not, you can get by with clustering. Record the whole "sentence", split it into "words" ; compute average feature vectors over each word ; and run K-means clustering to get 2 class out of it. One class will be "tou" and the other "tiu". The class whose "words" have the lowest average spectral centroid are likely to be the "tou".

A last word on speech recognition: if you have a vocabulary of two words, the cost of running a decent speech recognition toolchain will be that of:
* Computing the MFCC features (which you have to do anyway no matter what).
* Scoring a dozen of gaussians (which you have to do anyway).
* Finding the lowest cost paths in a FST with less than a dozen states.

This is quite pedestrian and nowhere near the complexity of natural language speech recognition. Since you have no experience in digital signal processing, doing so will probably be faster than whatever custom code you might write. Thus I very strongly suggest you to collect a small training corpus of utterances from people around you and start looking into "off the shelf" speech recognition toolboxes like HTK. Just using the built-in tools you will be able to train a model evaluate it on your data without writing code.


How does constant frequency response relate to equiripple? The teacher says As a consquence, you can think of the frequency response of a filter as a shark. It must always move and it can never be rest. An important case is equiripple.

I am very confused about the abrupt conception..

Why constant frequency response can never be rest?
2.how does constant frequency response relate to equiripple?

Thanks you!!

What your teacher tried to show is that a transfer function which is the ratio of two polynomials can never be constant (excluding the trivial case where both polynomials are of degree zero). This means that if your desired frequency response is a piecewise constant function (as is the case for an ideal frequency selective filter), you can only approximate that piecewise constant function by a realizable transfer function, and one popular approximation is an equi-ripple approximation, where the maximum deviation from the desired constant value is minimized.


Random selection of sub-block in an image with matlab I want to select a random sub-block of 8x8 out of an image 256x256. Is there an instruction that do this thing or i need to calculate it with a different approach. I found how to choose a galois field array out of a  matrix but i think this is not in my case what i need.


Image thresholding while retaining "underneath" contents I have sample images, and I need to extract certain features that depend on color. Images are more like line art, contours, and consist of various line forms, and several colors.

I can extract specific features through thresholding, but if lines intersect I always get discontinuities (holes) that are not acceptable. As example I prepared this:

I have image with two "lines":

I'll threshold to remove black, but underneath line (which of course doesn't exist as information in bitmap representation) gets discontinued:

while I want it like this:

which I'll later equalize, and get wanted context.

I used image processing program, and would appreciate some tip how to do this in Gimp or Photoshop, but any answer is welcomed like Python/ndimage, Matlab/Octave...

This is by no means a complete answer but I think it will get you far enough to achieve what you want. Edit: maybe not.

find main colours using kmeans
assign main colours to different types (brown line, blue line, etc). Some appear in both
clean up with morphology
segment image

What I thought might work was to pick out the most common colours and classify each pixel by the nearest colour. Then classify each pixel by type (brown line, etc), where the colour class can belong to multiple type classes. This would hopefully solve the problem of gaps in the lines - i.e. the colour where blue crosses brown would belong to both type class and fill the gap.

Not sure if it really does a better job. The issue with kmeans is that it uses random starting points so the colour classes were changing every time the script was run. I think I've fixed that now.

I wonder if there is some code out there to connect the broken lines?

% Read file
I = imread('mapdsp.jpg');

% Reshape into 2D matrix of pixel colours
ir = reshape(I, [size(I,1) * size(I,2),3]);

% Choose exact same sample points everytime
rng(0);
ir_indx = randi(size(ir,1),20,1);
ir_startpoints = ir(ir_indx,:);

% Use k-means clustering to find most common colours
[ik,c,sk,sd] = kmeans(double(ir),20,'Start',double(ir_startpoints));

% Show colours
cim = reshape(uint8(c),[1,20,3]);
image(cim); title('Colours'); pause;

% Sort colours
chsv = rgb2hsv(c);
[~,cni] = sort(chsv(:,1),1);
cn = c(cni,:);
cnim = reshape(uint8(cn),[1,20,3]);
image(cnim); title('Sorted colours'); pause;

% Get index of index
[~,cni2] = sort(cni);

% Show pixels classified
in = reshape(cni2(ik,:), [size(I,1),size(I,2)]);
imagesc(in);  title('Pixel classification'); pause;

% Show each class
In = zeros(size(I));
for k = 1:20
    In(:,:,1) = I(:,:,1) .* uint8(in == k);
    In(:,:,2) = I(:,:,2) .* uint8(in == k);
    In(:,:,3) = I(:,:,3) .* uint8(in == k);
    image(uint8(In)); title (['Class ' num2str(k)]); pause;
    im_dist(:,:,k) = reshape(sd(:,k), [size(I,1),size(I,2)]);
    imagesc(-im_dist(:,:,k) .* (in == k));  title('Pixel distance'); pause;
end

% 12,19 and 20 blue lines
im_blue = (in == 13) | (in == 19) | (in == 20);
imagesc(im_blue); title('Blue lines'); colormap gray; pause;
im_blue = imopen(im_blue,strel('diamond',1));
imagesc(im_blue); title('Blue lines opened'); colormap gray; pause;
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_blue); end;
image(uint8(In)); pause;
imagesc(bwmorph(im_blue,'skel',Inf));

% brown lines
im_brown = (in == 1) | (in == 2) | (in == 1) | ...
    (in == 7) | (in == 8) | (in == 11) | (in == 12) | (in == 13);
imagesc(im_brown); title('Brown lines'); pause;
im_brown = imopen(im_brown,strel('diamond',1));
imagesc(im_brown); title('Brown lines opened'); pause;
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_brown); end;
image(uint8(In)); pause;
imagesc(bwmorph(im_brown,'skel',Inf));

% 9 - black line
im_black = (in == 10);
imagesc(im_black); title('Black lines'); pause;
im_black = imopen(im_black,strel('diamond',1));
imagesc(im_black); title('Black lines opened'); pause
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_black); end;
image(uint8(In)); pause;
imagesc(bwmorph(im_black,'skel',Inf));

% 16,17,18 green area
im_green = (in == 16) | (in == 17) | (in == 18);
imagesc(im_green); title('Grean area'); pause;
im_green = imopen(im_green,strel('diamond',1));
imagesc(im_green); title('Green area opened'); pause;
im_green = imclose(im_green,strel('disk',7));
imagesc(im_green); title('Green area opened and closed'); pause;
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_green); end;
image(uint8(In)); pause;


What are advantages of having higher sampling rate of a signal? Being a non signal processing science student I have limited understanding of the concepts.

I have a continuous periodic bearing faulty signal (with time amplitudes) which are sampled at $12\textrm{ kHz}$ and $48\textrm{ kHz}$ frequencies. I have utilized some machine learning techniques (Convolutional Neural Network) to classify faulty signals to the non faulty signals.

When I am using $12\textrm{ kHz}$ I am able to achieve a classification accuracy $97 \pm 1.2 \%$ accuracy. Similarly I am able to achieve accuracy of $95\%$ when I applied the same technique on the same signal but sampled at $48\textrm{ kHz}$ despite the recording made at same RPM, load, and recording angle with the sensor.

What could be the reason for this increased rate of misclassification?
Are there any techniques to spot differences in the signal?
Are higher resolution signals prone to higher noise?

Details of the signal can be seen here, in chapter 3.

Sampling at a higher frequency will give you more effective number of bits (ENOB), up to the limits of the spurious free dynamic range of the Analog to Digital Converter (ADC) you are using (as well as other factors such as the analog input bandwidth of the ADC). However there are some important aspects to understand when doing this that I will detail further.

This is due to the general nature of quantization noise, which under conditions of sampling a signal that is uncorrelated to the sampling clock is well approximated as a white (in frequency) uniform (in magnitude) noise distribution. Further, the Signal to Noise Ratio (SNR) of a full scale real sine-wave will be well approximated as:

$$SNR = 6.02 \text{ dB/bit} + 1.76 \text{dB}$$

For example, a perfect 12 bit ADC samping a full scale sine wave will have an SNR of $6.02\times 12+1.76 = 74$ dB.

By using a full scale sine wave, we establish a consistent reference line from which we can determine the total noise power due to quantization. Within reason, that noise power remains the same even as the sine wave amplitude is reduced, or when we use signals that are composites of multiple sine waves (meaning via the Fourier Series Expansion, any general signal).

This classic formula is derived from the uniform distribution of the quantization noise, as for any uniform distribution the variance is $\frac{A^2}{12}$, where A is the width of the distribution. This relationship and how we arrive at the formula above is detailed in the figure below, comparing the histogram and variance for a full-scale sine wave ($\sigma_s^2$), to the histogram and variance for the quantization noise ($\sigma_N^2$), where $\Delta$ is a quantization level and b is the number of bits. Therefore the sinewave has a peak to peak amplitude of $2^b\Delta$. You will see that taking the square root of the equation shown below for the variance of the sine wave  $\frac{(2^b\Delta)^2}{8}$ is the familiar $\frac{V_p}{\sqrt{2}}$ as the standard deviation of a sine wave at peak amplitude $V_p$.  Thus we have the variance of the signal divided by the variance of the noise as the SNR.

Further as mentioned earlier, this noise level due to quantization is well approximated as a white noise process when the sampling rate is uncorrelated to the input (which occurs with incommensurate sampling with a sufficient number of bits and the input signal is fast enough that it is spanning multiple quantization levels from sample to sample, and incommensurate sampling means sampling with a clock that is not an integer multiple relationship in frequency with the input). As a white noise process in our digital sampled spectrum, the quantization noise power will be spread evenly from a frequency of 0 (DC) to half the sampling rate ($f_s/2$) for a real signal, or $-f_s/2$ to $+f_s/2$ for a complex signal.  In a perfect ADC, the total variance due to quantization remains the same independent of the sampling rate (it is proportional to the magnitude of the quantization level, which is independent of sampling rate). To see this clearly, consider the standard deviation of a sine wave which we reminded ourselves earlier is $\frac{V_p}{\sqrt{2}}$; no matter how fast we sample it as long as we sample it sufficiently to meet Nyquist's criteria, the same standard deviation will result. Notice that it has nothing to do with the sampling rate itself. Similarly the standard deviation and variance of the quantization noise is independent of frequency, but as long as each sample of quantization noise is independent and uncorrelated from each previous sample, then the noise is a white noise process meaning that it is spread evenly across our digital frequency range. If we raise the sampling rate, the noise density goes down. If we subsequently filter since our bandwidth of interest is lower, the total noise will go down. Specifically if you filter away half the spectrum, the noise will go down by 2 (3 dB). Filter 1/4 of the spectrum and the noise goes down by 6 dB which is equivalent to gaining 1 more bit of precision! Thus the formula for SNR that accounts for oversampling is given as:

Actual ADC's in practice will have limitations including non-linearities, analog input bandwidth, aperture uncertainly etc that will limit how much we can oversample, and how many effective bits can be achieved. The analog input bandwidth will limit the maximum input frequency we can effectively sample. The non-linearities will lead to "spurs" which are correlated frequency tones that will not be spread out and therefore will not benefit from the same noise processing gain we saw earlier with the white quantization noise model. These spurs are quantified on ADC datasheets as the spurious-free dynamic range (SFDR).  In practice I refer to the SFDR and usually take advantage of oversampling until the predicted quantization noise is on level with the SFDR, at which point if the strongest spur happens to be in band, there will be no further increase in SNR. To detail further I would need to refer to the specific design in more detail.

All noise contributions are captured nicely in the effective number of bits (ENOB) specification also given on ADC data sheets. Basically the actual total ADC noise expected is quantified by reversing the SNR equation that I first gave to come up with the equivalent number of bits a perfect ADC would provide. It will always be less than the actual number of bits due to these degradation sources. Importantly, it will also go down as the sampling rate goes up so there will be a diminishing point of return from oversampling.

For example, consider an actual ADC which has a specified ENOB of 11.3 bits and SFDR of 83 dB at 100 MSPS sampling rate. 11.3 ENOB is an SNR of 69.8 dB (70 dB) for a full scale sine wave. The actual signal sampled will likely be at a lower input level so as not to clip, but by knowing the absolute power level of a full scale sinewave, we now know the absolute power level of the total ADC noise. If for example the full scale sine wave that results in the maximum SFDR and ENOB is +9 dBm (also note that this level with best performance is typically 1-3 dB lower than the actual full scale where a sine wave would start to clip!), then the total ADC noise power will be +9dBm-70 dB = -61 dBm.  Since the SFDR is 83 dB, then we can easily expect to gain up to that limit by oversampling (but not more if the spur is in our final band of interest). In order to achieve this 22 dB gain, the oversampling ratio N would need to be at least $N= 10^{\frac{83-61}{10}} = 158.5$  Therefore if our actual real signal bandwidth of interest was 50MHz/158.5 = 315.5 KHz, we could sample at 100 MHz and gain 22 dB or 3.7 additional bits from the oversampling, for a total ENOB of 11.3+ 3.7 = 15 bits.

As a final note, know that Sigma Delta ADC architectures use feedback and noise shaping to achieve a much better increase in number of bits from oversampling than what I described here of what can be achieved with traditional ADC's. We saw an increase of 3dB/octave (every time we doubled the frequency we gained 3 dB in SNR). A simple first order Sigma Delta ADC has a gain of 9dB/octave, while a 3rd order Sigma Delta has a gain of 21 dB/octave!  (Fifth order Sigma Delta's are not uncommmon!).

Also see related responses at

How do you simultaneously undersample and oversample?

Oversampling while maintaining noise PSD

How to choose FFT depth for ADC performance analysis (SINAD, ENOB)

How increasing the Signal to Quantization noise increases the resolution of ADC


discrete-time (slew) rate limiting algorithm Consider uniformly-spaced samples of smooth, bandlimited signal observed in noise and subject to some artifacts (small jumps).  Physical restrictions impose a constraint on the maximum magnitude of the derivative of the noiseless, "artifactless" signal.  In an effort to smooth the signal and remove artifacts, I could impose this constraint by simply limiting the size of the sample-to-sample change (e.g., matlab, but this seems like it will introduce other effects to compromise the smoothness of the signal.

The only other option I came across are nonlinear feedback loops (e.g., see this paper.  I suppose that one could also impose derivative inequality constraints on, say, some local polynomial representation of the signal as well.  Any other ideas, references, etc?

Thanks!

Like the old German saying goes:

If you give a man a filter, everything looks like a signal that needs to be low-pass filtered

or so.

Slew rate and signal bandwidth are closely related. Simply low-pass filter your signal; usually, you'd do that with a linear filter, not an non-linear filter, and especially not a non-linear recursive one, because those are pretty hard to get stable.

By the way, if you look at FIR and IIR filters, you'll notice that we model these in z-Domain as polynomials, but by doing so, we can apply them to the signal in a linear way.


Pade Approximation of dead time For the time delay $e^{-sT}$ I shall find the Pade Approximation for $M = 0$ and $N = 1$.

$f(s) = \sum_i^{\infty} a_is^i \approx \dfrac{\sum_{n=0}^{N} b_is^i}{\sum_{m=0}^{M} c_is^i}$

$e^{-sT} = \sum_{i=0}^{\infty}\dfrac{(-sT)^i}{i!}$

Using the Taylor approximation for $e^{-sT}$ I yield

$1 -sT \approx \frac{P_N(s)}{Q_M(s)}$

$1 -sT \approx \frac{b_0 + b_1 s}{c_0}$

with $c_0 = 1$ per definition according to my materials.

This leads me to

$1 -sT \approx b_0 + b_1 s$

$b_0 = 1$

$b_1 = -T$

If I now want to check the step response of $e^{-sT}$ and my approximation MATLAB complains about my function having more zeroes than poles, which is not surprising as I thought I would get my transfer function from $G(s) = \frac{P_N(s)}{Q_M(s)}$.

I'm sure there is a basic misunderstanding of what I am doing here on my side.

There is no misunderstanding at all. The Padé approximant you found is correct.

The "problem" is that you chose $M$ and $N$ such that the transfer function you get to approximate the delay is improper. Namely, the order of the numerator exceeds the order of the denominator.

In Control Theory, improper systems are not too useful because they cannot be realized in real life. A transfer function with more zeros than poles contains pure differentiators, so the transfer function represents a system that is non-causal (assuming stability, as Matt pointed out in the comments below), and non-causal systems cannot be physically realized.


Sine Dipole formation using two hydrophones I am trying to do MATLAB simulations for generation of SINE Dipole using two hydrophones spaced distance 'd' apart for signal coming from direction 'DOA'. The MATLAB code is given below. The confusion is that there is constant phase difference of 90 degree b/w SINE Dipole generated using difference of two hydrophones and simple SINE Dipole generated by multiplication of signal with sin(DOA). There is also a difference b/w amplitude of two Dipoles. The resultant graph of MATLAB code is also attached.
Am i doing right or something is not right in my understanding?
clc
clear
close all

f_sig = 800;%signal frequency in Hz
f_samp = 64000; %Sampling rate
N=81920;%no of data points
c=1500;%sound speed in water in m/s
d=50e-3;%spacing b/w two omni directional hydrophones

DOA=135;%Direction of signal in deg

t=(0:N-1)/f_samp;

lemda=c/f_sig;

Hyd_1=sin(2*pi*f_sig*t+2*pi*(d/2)*sind(DOA)/lemda);
Hyd_2=sin(2*pi*f_sig*t-2*pi*(d/2)*sind(DOA)/lemda);

%Formation of Sinusoidal dipole using two hydrophones Hyd_1 & Hyd_2
Sine_Dipole_1=(Hyd_1-Hyd_2);

%Formation of Sinusoidal by multiplying with sin(DOA)
Sine_Dipole_2=sin(2*pi*f_sig*t)*sind(DOA);

plot(Sine_Dipole_1(1:1000))
hold on
plot(Sine_Dipole_2(1:1000),'r')

As you've defined Hyd_1 and Hyd_2, you are simulating a two element hydrophone array with the array center defined at the middle of the two elements. In that case, Sine_Dipole_1 should be defined as the sum of the two elements, not the difference.
I hope this helps.


Removing mean from signal massively distorts FFT I have a series of signals sampled at 1200 Hz that are 200-300 samples long. I'd like to analyze them in the frequency domain in 1 Hz bins in MATLAB. I'm primarily interested in the frequencies below about 20 Hz, since this is where virtually all of the signal power is contained.

However, I'm finding conflicting opinions on whether it's appropriate to subtract the mean value from the signal before I perform the FFT function.  When I don't subtract the mean, I get a strong 0Hz component, which makes sense given the shape of the original signal. But some sources say you should zero-mean the data so you don't get distortion of the low frequencies, which are the ones I care about.  If I do subtract the mean, the FFT is markedly different, and appears to contain massive sinc-like distortions.

Here is the original signal in blue, and the zero-mean version in red, created using detrend(mySignal, 'constant').

As expected, the detrended signal has a mean of zero. But when I do the FFT and look at the spectral distribution, this is what I get:

Note that X-axis is frequency and I've truncated the axes to the first 50 frequency bins (0-49 Hz).

Is it inappropriate to zero-mean a signal like this? And what's the meaning of this distortion in the zero-meaned signal, if anything? Any help is much appreciated!

I bet you zero-padded the data before FFT. That's not going to work well because the results will be dominated by the discontinuity between your data and the zeros. It is not equivalent to subtract the mean before zero padding vs. after.

Try subtracting the mean, then multiplying the data with a window function, and then zero-padding.

The 0 Hz bin will still typically have a non-zero value. You can get rid of that by calculating the windowed, 0 Hz subtracted signal as (MATLAB/Octave notation):

(signal - mean(signal.*window)/mean(window)).*window

sum can be used in place of mean because the numbers of samples cancel.


OFDM. Cyclic Prefix and ISI Sorry that I created a new question about OFDM, but I donät have enough reputation to comment on other questions.

I have read all the topics here and I have some doubts.

1) I found in the book that If a system has a low symbol rate then --> negligible ISI. Why?
I have thought about it. Maybe because of that the system, in this case, has a longer symbol.

2) A cyclic prefix redices data rate. How can it be improved?

This is because intersymbol interference is caused by dispersion in the channel (assuming properly designed pulse shaping filters such that there is no "native" ISI in the waveform itself). This dispersion can be caused by analog filtering or over the air multi-path as illustrated in my graphic below:

For the case of multi-path, the multiple transmission paths result in a "delay spread".  According the the reference in the next graphic, for the macrocellular systems the typical delay spreads are:

Open areas: 0.2 uS
Suburban: 0.5 uS
Urban: 3 uS

When the symbol duration is significantly longer than the delay spread, there is minimum interference from one symbol to the next, but if the symbol duration approaches the delay spread, the interference will be significant! For a multipath case resulting from multiple replicas of the same waveform arriving at different times, this can be corrected with equalization. For the case of OFDM, we are actually converting the waveform into independent channels of significantly longer duration, thus solving this issue. The cyclic prefix allows a guard interval to eliminate the intersymbol interference from a previous OFDM symbol. To achieve this, the cyclic prefix must be at least as long as the delay spread of the channel.


Turn your guitar to MIDI and Guitar Synthesizer I don't know if this question is related to DSP, but I'll give a chance.

Introduction:
These days there are a Guitar Synthesizers which allow the guitarist
to play another instruments like: brass, drums, strings, pianos and other guitar models
from his guitar.

It possible to turn any guitar with the MIDI Pickup, that sends the all data
that the synth needs.
Now, from what I understand the pickup turns the guitar into MIDI one.

My question is how a signals like those that the guitar makes can be converted into
MIDI instructions ?

There is a sampling before this conversion? I think that it should be, because if not.
How the Synth will know in which level (volume) or duration to produce the sounds?

There is a sensor for each string which individually records the string vibration. This signal is converted to digital, analyzed for estimating its fundamental frequency and its amplitude (envelope detector). Whenever a sharp increase in amplitude is detected, a MIDI note on message is sent, at the detected frequency. Whenever the amplitude falls below a threshold, a MIDI note off message is sent.


Questions Downlink OFDMA in LTE I am looking for a clear and precise explanation of how an OFDMA modulator in downlink LTE operates. I have many books that explain it but each has its own way and do not answer my questions completely. If you believe this questions better fits in another forum please let me know (maybe some PHY layer forum).

My understanding is that an OFDMA modulator consists of the following steps:

a block of bits for each user are first coded using a error correction code.
the stream of bits for each user are then mapped to symbols for example QAM symbols.
There is Serial to Parallel converter, now the QAM symbols for each user are parallel.
Subcarrier mapping ? ( I don't understand this step)
The QAM symbols are input to the L IFFT (L being number of subcarriers)

According to LTE specifications, resources are allocated to users in units of resource blocks over a subframe (12 subcarriers over 14 OFDM symbols) i.e, 168 resource elements.

What does this mean, is there where step 4 comes into place? So is each user's QAM symbols being transmitted over 12 subcarriers and not over all available subcarriers?
What happens in the next subframe? Then why do we perform L point IFFT? how does then one compute the rate for each user over a subframe or frame basis?

In LTE downlink, users are multiplexed in both time and frequency domain. The concept of ressource blocks (RB) describes a certain frequency band in a certain time slot. Your understanding of the steps in the transmitter are correct. Step 4 describes the procedure that assigns a certain frequency band $B_m$ to some user $m$. One OFDM symbol in LTE downlink carries the data of several users. Subcarrier mapping means that the data of user $m$ is allocated to certain subcarriers that represent the frequency band $B_m$ as shown in the following figure (where $M$ is the total number of frequency bands)

This technique corresponds to frequency multiplexing. Time multiplexing means that the the band $B_m$ is reserved for user $m$ only for a certain time duration, i.e. for a certain number of consecutive OFDM symbols (you mention 14 in your question, I haven't checked it). In the next subframe this band might be allocated to a different user.

The bit rate of user $m$ can be calculated as follows. Let $N$ be the number of subcarriers per user, $b$ the number of bits per subcarrier and $T$ the OFDM symbol time including the cyclic prefix. Then the bit rate $R_\mathrm{b}$ of user $m$ is
$$
R_\mathrm{b} = \frac{bN}{T}
$$


Lossless algorithm for hyperspectral imaging I am looking for Lossless compression algorithm for hyperspectral imaging. The only one that I found was LPVQ.
Do you know other algorithms? I am looking for pseudocode or matlab code.

Thanks

I have played with SVD for 3 channel image lossy compression in the past, which prompted me to do a quick search and found the following paper which may help you out.
Randomized SVD Methods in Hyperspectral Imaging

Abstract
We present a randomized singular value decomposition (rSVD) method for the purposes of lossless compression, reconstruction, classification, and target detection with hyperspectral (HSI) data. Recent work in low-rank matrix approximations obtained from random projections suggests that these approximations are well suited for randomized dimensionality reduction. Approximation errors for the rSVD are evaluated on HSI, and comparisons are made to deterministic techniques and as well as to other randomized low-rank matrix approximation methods involving compressive principal component analysis. Numerical tests on real HSI data suggest that the method is promising and is particularly effective for HSI data interrogation.


How to generate realistic echo to sounds? I guess the most straightforward way of adding an artificial echo to a sound signal is to add its time-delayed version to the original signal. However, I am not sure how realistic that would sound like.

I was reading papers on room geometry estimation from impulse responses (see for example this paper). So, perhaps one should have several echoes instead of just one.

I guess something also happens to a sound signal when it is reflected (kind of filtering).

What would you suggest here?

A very simple model is a sum of attenuated and delayed replicas of the signal. Such signal can be be generated using FIR filter. The length of the filter is the "echo time", i.e. how long the echoes last.
The simplest approach would be to generate sparse impulse response. Build a filter with $n$ non-zero coefficient and fading amplitude (for example the amplitude can be generated with Rayleigh distribution).

If you want a more realistic scenario, I suggest using the Room Impulse Response (RIR) generator by Emanuel Habets.


Frequency magnitude distribution of noise By the Wiener-Khinchin theorem, we have a straightforward way to calculate the power spectral density for stationary noise.  If we know the bandwidth of a system, we can further calculate the variance of the noise since it turns out that $v_{noise,\ RMS} = \sigma$ (standard deviation) for zero mean noise.  Therefore, at each point in time we know the probabilistic distribution of noise as shown in this video:

However, can we determine anything about the distribution of noise magnitude at each point in frequency from a noise source with known power spectral density?  I believe this could be related to power spectral density estimation?

This video says that power spectrum realizations are proportional to the square of the underlying PSD, though I am not sure if that is just an artifact of using periodograms.  In any case, 'variance proportional to the square' is still not a precise probabilistic distribution.

Let me try to elaborate a bit on your questions: First, the Wiener-Khintschin-Theorem relates the autocorrelation of a stationary process to its PSD, with the relation of the Fourier Transform, i.e.
$$\mathcal{F}\{r_{nn}(\tau)\}=S_{nn}(f).$$
The autocorrelation of a signal describes the correlation between two noise samples at distance $\tau$. It is just the second moment of the (joint) distribution of both noise samples. Hence, it does not really tell you something about the actual distribution. However, mostly, the noise is assumed to be Gaussian distributed with zero mean. Then, the autocorrelation completely describes the noise process.

Furthermore, in case you have AWGN, then $r_{nn}(t)=N_0\delta(t)$, i.e. the noise samples are uncorrelated and the PSD is a constant with $S_{nn}(f)=N_0$.

Your question is: What is the distribution of the noise magnitude in the frequency domain? So, you want to know $Pr(S(f)|_{f_0})$, where $S(f)$ is a realization of the noise process.

To make this more accessible (at least to me), let me switch to a linear algebra setting (i.e. the finite discrete-time case). The results can be straight-forwardly taken over to the continuos-time case. Let $\vec{n}$ be a realization of the noise process with autocorrelation $R_{nn}=E[\vec{n}\vec{n}^H]$. Let $\mathbf{F}$ be the unitary DFT matrix. So, $\vec{N}=\mathbf{F}\vec{n}$ is one noise realization in the frequency domain. Assuming that $\vec{n}$ consists of Gaussian distributed elements with zero mean, $\vec{N}$ also consists of Gaussian elements (A linear combination of two gaussian random variables remains to be Gaussian). What is the autocorrelation of the frequency-domain noise? It's given by

$$\mathbf{R}_{NN}=\mathbf{F}R_{nn}\mathbf{F}^H.$$

On the diagonal of $R_{NN}$ there is the variance of the Gaussians in the frequency domain.

So, the distribution of noise in the frequency domain remains to be a Gaussian, with variance given by diagonal on the above expression. I believe the diagonal is actually the PSD that follows from the Wiener-Khintschin Theorem. So, the final answer is:

  The noise in the frequency domain is Gaussian distributed, with zero mean and variance given by the PSD of the noise. (Assuming zero-mean Gaussian time-domain noise)


Upsample data using FFTs. How is this exactly done? I am basing my question on this post here, because I would like additional details on it, as I have not had any success in re-creating it.

I would like to simply learn, in no uncertain terms, how to do upsampling via the FFT. There seem to be a lot of details regarding if the signal is even or odd, where exactly to zero-pad, how much exactly, etc, and I am not getting those correct numbers for some reason.

Question: So let us say I have an even or odd length signal, of length $N$. I would like to upsample it to, say, $2N$, using FFTs.

How is this done exactly?

What I tried:

I tried following the post and others on the web but there seems to be very little consistency. It basically amounts to "zero-padding the frequency domain and then IFFTing" but I do not get proper results. In my case, I have a signal of length $100$, and I would like to upsample it to length $200$. So I do the following command:

ifft(fftshift([zeros(1,50) fftshift(fft(signal, 100)) zeros(1,50)]));

This however gives me complex data back, and not a signal that is simply upsampled by a factor of $2$. I am not sure where the problem is.

Thank you.

One trick, for even-length signals, is what to do with the "middle" sample.  Here, I've split it half and half between each side of the FFT.

The other trick is to ensure that you have the right amplitudes in the resampled signal. Here's it's a factor of 2.

Try this in scilab:

x = rand(1,100,'normal');
X = fft(x);
XX = 2*[X(1:50) X(51)/2 zeros(1,99) X(51)/2 X(52:100)];
xx = ifft(XX);

clf;
plot([0:199]/200,xx)
plot([0:199]/200,xx,'go')
plot([0:99]/100,x,'r.')

which appears to generate the right thing. The red dots are the original samples, and the green circles (and blue connecting lines) are the resampled data.

UPDATE

For the case of an odd number of samples, there are fewer fiddly bits (no splitting the last coefficient):

x2 = rand(1,101,'normal');
X2 = fft(x2);
XX2 = 2*[X2(1:51) zeros(1,101)  X2(52:101)];
xx2 = ifft(XX2);

clf;
plot([0:201]/202,xx2)
plot([0:201]/202,xx2,'gx')
plot([0:100]/101,x2,'r.')


phase angle of harmonic I have generated sinsuoids containing fundamental with harmonics using C++. Fundamental frequency is 50 Hz and each of the harmonic has phase offset of 60 degrees. I have applied flattop window before performing and N point FFT. Magnitude repose is correct. Phase reponse is wrong. I should expect 60 degree in all 100,150,200,250 Hz.
I have attached the results.

Sampling rate : 40000
Timeduration : 1s
Maximim amplitude of each of the harmonic : 1

Any help would be highly appreciated.

Regards,
Faisal

Your window is centered around the middle of the time window. This corresponds (roughly speaking) to a delay of half the window size and corresponding additional phase of $e^{-j \cdot \omega \cdot t}$

The "wiggles" that you see is a steep linear phase that wraps around 2*pi occasionally. There a few ways to fix this.

circularly shift the time domain signal so that the max of the window is at 0. This will be better but still isn't great since the windowing will significantly impact the phase.
Make sure that the fundamental fits into your time window and integer amount of times. Then do an FFT without windowing. This should give you an exact result.


Large spike at the center frequency when using Ettus X310 I am using the Ettus X310 SDR and GNU radio. On the GNU radio schema I have a standard UHD USRP source block connected to a GUI FFT sink - basically I'm just displaying the spectrum of the raw signal from the device.

Parameters for the USRP source block:

sample rate = 1.024 MS/s,
center freq = 101 MHz,
gain = 42 dB,
bandwidth = 130 MHz

The problem that I have is that I always have a big spike (10-15 dB) directly on the center frequency (no matter what frequency I set). I am relatively new to all this so I would appreciate any pointers on how to get rid of the spike.

The problem that I have is that I always have a big spike (10-15 dB)
  directly on the center frequency (no matter what frequency I set). I
  am relatively new to all this so I would appreciate any pointers on
  how to get rid of the spike.

That spike is probably nothing surprising – just the LO leakage/DC offset, a very common artifact in direct conversion receivers!

Now, getting rid of it with a USRP is relatively simple.

You have to understand that the USRPs have a two-step tuning process:

There's a physical oscillator, the LO, which gets synthesized by the daughterboard, and is used to mix down the target signal to IQ baseband signals. Typically, synthesizers can only generate a discrete set of LO frequencies $f_\text{RF}$.
The ADC then samples these signals at the master clock rate, and the FPGA digitally shifts that sample stream in the master clock rate (i.e. 200 MHz for the X310 by default) wide baseband by essentially multiplying with $e^{j2\pi f_\text{offset} n}$, to give you your baseband signal from the desired $f_\text{target}$:

By default, $f_\text{RF}$ is set as closely as possible to $f_\text{target}$, minimizing $f_\text{offset}$, and hence, putting the LO into your observed bandwidth $b_\text{sample}$.

You can, however, manually set a "desired" $f_\text{offset}$; as soon as that is more than $\frac 12 b_\text{sample}$, the LO is out of your observer band!

Instead of using 101e6 in the Center Frequency input field, just use

uhd.tune_request(f_target, f_offset)

e.g.

uhd.tune_request(101e6, 4e6)

Let me comment on a few other issues directly, too:

  I am using the Ettus X310 SDR and GNU radio.

The X310 is the platform, indeed, but the actual analog signal processing (generating the LOs for mixing, mixing, amplification, baseband filtering) happens on the daughterboard. In such cases, it's usually important to specify which daughterboard you're using

  On the GNU radio schema I
  have a standard UHD USRP source block connected to a GUI FFT sink -
  basically I'm just displaying the spectrum of the raw signal from the
  device.

Indeed, that should work. Notice that the WX Gui blocks are kind of becoming deprecated; if you're just starting with GNU Radio and don't need to continue an existing application based on WX, I'd recommend changing the "build options" in the option block to "Qt", and replacing the "WX GUI FFT sink" with a "Qt GUI frequency sink". Same thing, more functionality, less CPU usage, usually.

  Parameters for the USRP source block:

  sample rate = 1.024 MS/s,

You probably were using a master clock rate of 184.32 MHz, right?

For future readers, reproducing this (this specific rate really comes up more often then you'd think!):

The X310 doesn't support that sampling rate in the default configuration! The console window should contain a clear warning that UHD decided to use a different frequency, and which frequency that was – in your case, probably 1.02048… MS/s.

UHD Warning:
    The hardware does not support the requested RX sample rate:
    Target sample rate: 1.024000 MSps
    Actual sample rate: 1.020408 MSps

That's because the X310 runs its ADC/DAC at what is called the master clock rate, which can take the values 200 MHz (default), 184.32 MHz (typical for some cellular network standards), 120 MHz (typical for LabView). The sample rate you deliver to the PC is a sample stream sampled at that rate, decimated in the FPGA to the rate you're requesting – and that decimation is done with "normal" decimating FIR filters, which can only give you an integer fraction of the input rate.

So,

$$1.020408… \frac{\text{MS}}{s} = \frac{200\frac{\text{ MS}}{s}}{196}$$

is possible,

$$1.024 \frac{\text{MS}}{s} = \frac{200\frac{\text{ MS}}{s}}{195.3125}$$

is not, when the MCR is set to it's default of 200 MHz.

With an MCR of $184.32\text{ MHz}$, which is $1.024\text{ MHz}*180$, this is indeed possible!

  center freq = 101 MHz,

Which probably means you're using the WBX or UBX daughterboard,

  gain = 42 dB,

which is not a gain setting that UBX or WBX support, their gain ranges are 0 – 31.5 dB and 0 – 31 dB, respectively.

  bandwidth = 130 MHz

That's the analog bandwidth specification; it applies to daughterboards with adjustable baseband low pass filters. However, there's no daughterboard in existence that has an adjustable baseband filter and a maximum bandwidth > 8 MHz.

In fact, this setting does nothing on UBX or WBX.


How to calculate the output image with the following kernel? What is the advantage of using this equation? I guess we may use Taylor series, but  I tried my best and could not get   the equation.

An image is 2D, the kernel seems 1D. So let us double check: from (b), the output image is determined for each fixed $i$, and then only depends in $j$: Equation (1) only uses the second coordinate as the moving index.

So the kernel is effectively a $1\times N $ matrix. Its sum is $1$, and as $N$ is odd, you can interpret it as an equally-weighted moving average filter, applied along lines only. So it can behave like an horizontal motion blur (Question (a)).

The direct convolution requires $N-1$ adds, and one multiply (by $1/N$). Equation (1) tries to make it more efficient, both in terms of operations and calls to samples. It simply uses a recursive formulation for 1D moving average filters (a line of the image, here $x_j = I(.,j)$:

$$ \sum_{j=1}^N x_j = \sum_{j=0}^{N-1} x_j - x_0 + x_N $$

I suspect that in your formula, $K$ should be defined as $2K+1=N$. There is consistent with $K=\mathrm{floor}(N/2)$ for $N$ odd.

In other words: instead of doing sums of pixels each time you go to the right, you can reuse previous computations. The novel mean $I_0(i,j)$ is the last mean $I_0(i,j-1)$  minus $1/N$ the leftmost sample (the one that just left the previous frame) plus $1/N$ the rightmost sample (the one that just arrived).

Almost nothing to do with Taylor. The advantages:

a running buffer of $3$ values only is needed,
its is independent of the kernel size,
some computations can be saved if $N$ is relatively large (say $N\ge 5$).


How to compand an OFDM signal? Im trying to run a simulation of an OFDM system with convolutional coding and u-law compander on MATLAB. The problem is I dont know how to compand the OFDM signal. Do I need to do something to the signal before compressing it? Because value coming out of the ifft block is a complex number.
Xofdm_compand= Xofdm_cp; %OFDM with cp
xofdmcompand=abs(Xofdm_compand);
for j=1:Nframe %Number of symbols
    in=xofdmcompand(j,:);
    v=max(abs(in));
    X_compressed(j,:)=compand(in,mu,v,'mu/compressor');
end

Edit: Here is the snippet of the code where I compand the signal. Am I doing it wrongly?

µ-law is defined for real-valued PCM signals only. You will either need to produce a real-valued signal containing OFDM as passband signal, either

by mixing it to a carrier frequency, or
by modifying your OFDM system to only produce symmetrical spectrum
baseband, or
by u-law coding I and Q separately (but that makes little sense, since a high PAPR happens when both I and Q have a high value at the same time, i.e. coding them separately isn't logical),

but all three options wouldn't be "companding OFDM" anymore, that would be "companding a system which, among other things, involves OFDM".

General remark: Don't forget that you've got channel coding. It's often better to intentionally completely "delete" isolated subcarriers from extreme-peak OFDM symbols than to  "make whole symbols slightly worse", as the former leads to a certain, but isolated, error, whereas the other has a  chance of leading to multiple errors in the same symbol – which would, assuming there's no large interleaver at work, look like a burst error.

I think you'd want to read

X. Wang, T. Tjhung and C. Ng, "Reduction of peak-to-average power ratio of OFDM system using a companding technique", IEEEE Transactions of Broadcasting, vol. 45m, no. 3, Sept. 1999.

on what they did and what their results are.

While you do that, look at Figure 5:

Notice how the second (⛛) curve is the uncompanded one. This looks like for anything but very high SNR, the companded curves are better.

Now, read

A. Mattsson, G. Mendenhall and T. Dittmer, "Comments on reduction of peak-to-average power ratio of OFDM system using a companding technique", IEEE Transaction on Broadcasting, vol. 45, no. 4, December 1999.

and especially this paragraph:

  Which shows that there is a significant increase in power
  for the companded signal. I.e. the SER for the companded
  signal should in our opinion be compared to the SER for
  the non-companded signal amplified by 10.9 dB [for the µ=16 case], 17 dB [for the µ=255 case] etc. It appears that the SNR refered to in Fig. 5 on page 307,
  refers to the signal power of the signal before companding
  and not the actual SNR of the transmitted signal.

In other words: shift the red (uncompanded) curve 11 dB to the left to compare against the µ=16 case, and by 17 dB to the left to compare against the (standard) µ=255 case.

Look how the uncompanded signal (red line) is always better at the same power than the µ=16 companded signal, and how it's better for SNRs > 12 dB for the µ=255 case. I highlighted the amount that uncompanded is better than companded in green for µ=16, and in orange for µ=255.

If your SNR is below 12 dB, you get a SER > 10⁻³, which is a typical design goal for wireless communication; you might want to use a lower modulation than 16-QAM then, which will then typically be an 8-PSK or QPSK. Which is indifferent about  amplitude variation – so instead of companding nonlinearly, you'd just look at the peak in your OFDM symbol and simply either drastically scale a few isolated subcarriers, or the whole symbol slightly. (which you do would depend on how steep the BER curve of your channel decoder is at that point.)

TL;DR: Literature research shows companding only works for low SNR. There, you might not be in a sweet spot for system design, so my expectations for companding to work are low.  Instead, I propose you generally rather design your FEC so that it can live with intentionally broken isolated symbols that you insert on purpose to reduce PAPR.


How to pick coefficients for Fractional Delay Filters? I have a Virtex 6 FPGA running at 200MHz with ADC/DACs on it.  I have been converting a WiFi signal (2462MHz) down to a more reasonable IF of 25MHz, sampling, running the signal through taps, tweaking them, and mixing back up.  This all works fine.

My issue is that since I am mixing down and sampling at a lower rate, I cannot implement RF level multipath by adjusting the delay in the taps for phase changes, instead I had to adjust it at the IF level instead.  This works OK, but it means that timing is off (things like TDOA or stuff with TOF measurements could care about that).  It seems like the best way around this (and probably the right way) is to implement a fractional delay filter (FDF), but I am not sure the best way to do it.

It seems like I need to come up with the coefficients for the taps, but I cannot seem to glean that out of the whitepapers I've read on FDF.  Any ideas?

I usually do this by creating a low pass filter that entirely passes through the signal that I want to delay.  I create the LPF "manually" by creating a windowed sinc function.  Something along the lines of-

filt = sinc(-80:.8:80);
filt = filt .* hamming(length(filt)).';

This gets you a filter that passes about 80% of the nyquist region (the 80% is set by the .8 increment in the sinc indexing).  To get a fractional delay you simply include the fractional delay in the calculation of your filter.  Adding a constant to the sinc indexing delays the filter by $\frac{constant}{increment}$ samples.  For instance, in the case of the filter above, adding 0.1 would create a fractional delay of $\frac{0.1}{0.8}$ samples, or $0.125$ samples.

filt = sinc((-80:.8:80) + .1);
filt = filt .* hamming(length(filt)).';

This works because filtering is a linear operation and delaying the filter is equivalent to convolving your original LPF with a fractionally delayed delta function.  Since linearity implies that the associative property holds, convolving your signal with the delayed LPF is equivalent to convolving your signal with a fractionally delayed delta function and then filtering it with the "undelayed" LPF.

Note: When creating the "undelayed" LPF you must make sure that the indices are symmetric.  The easiest way to do that is to make sure the start and end points are integer multiples of your indexing increment.  That is why I used +/- 80 for an increment of 0.8.


Fourier transform of a noisy signal Let's suppose i have a sensor that measures the output of a process and this cause a not negligible noise that affects my signal. My goal is to analyze the process signal in order to find faults. How can i deal with noise ? Is it possible to apply the fourier transform to the signal with noise or i have to filter it in advance? Maybe is possible to deduce which frequencies does the noise affect...

If i have to filter the noise, which are the techniques to do it ? Thanks

Usually, there are two types of noise that you can eliminate by using the spectrum.

One type is a noise that is in a different frequency band than the signal (it can be a high-frequency noise). In such a case, there is a clear separation between the signal content and the noise in the spectrum and you have to do nothing.

The second case is when the noise is in the same frequency band as the signal and then you can perform Welch's method of averaging the spectrums to reduce the noise component.

In case you have a dynamic model of the signal generator, you can reduce the noise using some kind of adaptive filter such as Kalman-filter prior to the evaluation of the spectrum.


(I)FFT without Bit Reversal It's been a long time since I scratch built an FFT algorithm, and I feel like I should already know this, but I just can't recall and figured I might get a quick answer on here that may be useful to others.  FWIW, I did see one question that was already answered that was very related, but I am hoping for clarification.
Looking at Decimation in Time algorithms, bit reversal is typically shown applied to the input.  For Decimation in Frequency, it is typically shown applied to the output.  My questions are:

What happens if no bit reversal is performed?
If the natural order is important, does it matter if the bit reversal is performed at the beginning or the end?
Are the DIT and DIF algorithms mathematically equivalent?

I am implementing a FFT based convolution.  It would seem that the bit reversal is not needed and I was hoping to figure out why.  I found the figures/analysis here helpful for reference.

If you're using the FFT and iFFT to perform fast convolution, yes, you can do it in-place and you can do without bit-reversing in the sample processing, but you will have to bit-reverse the transfer function $H[k]$.  You would use the DIF for the forward FFT and the DIT for the inverse FFT.
O&S have the clearest treatment, in my opinion.  I am robbing some drawings from O&S.
Note this DIF accepts the data input as normal order in and returns bit-reversed data in the output.  Note that the butterfly has the twiddle factor multiplied in the output of the butterfly.

Here is the DIT.  Note that it accepts the data input in bit-reversed order and returns normal order data in the output.  Also note that the butterfly exactly undoes the DIF butterfly.  It multiplies first at the input to the butterfly.

Now if you use the DIT as the iFFT, then you need to make sure you scale the input by $\frac1N$ (that can be done with the transfer function multiplication) and conjugate the twiddle factors $W_N^r$.  The multiplication in the frequency domain (between the FFT DIF and the iFFT DIT) must be done with the transfer function $H[k]$ pre-bitreversed.
In the 1980s, I wrote a very simple DIF FFT and DIT iFFT in C code.  Do you want it?  I might be able to find it.  Lemme know.
Actually, I posted the C code here before.


Is the following equation established? Is the following equation established?

  $$\int_{-\infty}^{\infty}s(t)r^*(t)dt=\int_W S(f)R^*(f)df $$

  $$s(t)\xrightarrow{\text{Fourier Transform}}S(f) $$

  $$r(t)\xrightarrow{\text{Fourier Transform}}R(f) $$

It's a version of Parseval's Theorem, which can be easily proved by noting that

$$\int_{-\infty}^{\infty}s(t)r^*(t)dt=\mathcal{F}\{s(t)r^*(t)\}\big|_{f=0}\tag{1}$$

where $\mathcal{F}$ denotes the Fourier transform. Note that in general the integration limits should be $-\infty$ and $\infty$ for both integrals. I assume your notation refers to band-limited functions.

From $(1)$, the proof is quite short:

$$\mathcal{F}\{s(t)r^*(t)\}=S(f)\star R^*(-f)=\int_{-\infty}^{\infty}S(\nu)R^*(\nu -f)d\nu\tag{2}$$

where $\star$ denotes convolution. So

$$\mathcal{F}\{s(t)r^*(t)\}\big|_{f=0}=\int_{-\infty}^{\infty}S(\nu)R^*(\nu)d\nu\tag{3}$$

which establishes the equation in your question.


Is there any relationship between the FFT and linear prediction? Is there any relationship between the FFT and linear predictive methods?  Can an FFT result or the input to an FFT be modified to do non-circular prediction/extrapolation from the FFT results.

One thing that comes to mind is that a linear predictor can be implemented using the LMS algorithm, and the LMS algorithm can be implemented in block fashion using FFTs (see Chapter 7 of Haykin's Adaptive Filter Theory). I wouldn't say that implies any kind of relationship, but FFTs can be used to accelerate the implementation of some linear predictors.


What is the relationship between the PSD of a continuous signal and the PSD of its periodically sampled one? In Oppenheim's "discrete time signal processing", it says,

How is the equation (10.65) derived?
PS:
you may try to derive it with Fourier transform of the original and the sampled signals.
But i think it is not appropriate.
Because when you are considering about the PSD of a signal, the assumption will be the fourier transform of the signal doesn't exist...
the definition of the PSD of a signal actually comes with a truncated signal:

the same is for the discrete version:

I guess there may exist something like... "sampling a random process" theory?
Regards.

If you sample a finite-power continuous-time random process $x(t)$ you get a discrete-time random process $y_k$. If $x(t)$ is wide-sense stationary (WSS) you get for the autocorrelation function of $y_k$

$$R_y(k,l)=E\{X(kT)X^*(lT)\}=R_x((k-l)T)=R_x(mT)$$

where $T$ is the sampling period. Obviously, $y_k$ is also WSS (it only depends on the difference $m=k-l$), and its autocorrelation function is a sampled version of the autocorrelation $R_x(\tau)$ of the original continuous-time process. Since the power spectrum is the Fourier transform of the autocorrelation function, the relation between the two power spectra is the same as the relation between the spectrum of a continuous signal and its sampled version:

$$S_y(e^{j\omega T})=\frac{1}{T}\sum_{k=-\infty}^{\infty}S_x(\omega-2\pi k/T)\tag{1}$$

Since the anti-aliasing filter has a cut-off frequency of $\pi/T$ it is sufficient to consider the interval $-\pi/T\le\omega<\pi/T$ (i.e. just consider the term $k=0$ in (1)):

$$S_y(e^{j\omega T})=\frac{1}{T}S_x(\omega)\tag{2}$$

This equation is equivalent to Eq. (10.65) in your question. Just note that I've used $\omega$ in a different way (analog frequency in radians), whereas Oppenheim uses $\omega$ as the normalized angular frequency. So with $\omega'=\omega T$ you finally get

$$S_y(e^{j\omega'})=\frac{1}{T}S_x\left(\frac{\omega'}{T}\right)$$


Why is random noise assumed to be normally distributed? From residual in the linear regression to noise in signal processing are assumed to be normally distributed? By considering them as normally distributed we are kind of telling the pattern in the noise but shouldn't noise be considered random. This seems contradictory to me as on one side it is random then on the other side their distribution is considered normally distributed. Shouldn't the noise distribution be just random?

I believe there is some lacking in my understanding of the concept of statistical distribution which has lead me to this confusion, or I am looking at it all wrong.

One more example- when one augment data by adding Gaussian noise then it is not expected to change the overall distribution of data, why?

Starting at an even more basic level than the other (much smarter) answers, I'd like to pick up on this part of the question:

  This seems contradictory to me as on one side it is random then on the other side their distribution is considered normally distributed.

Perhaps the issue here is what ‘random’ means?

To be clear: ‘random’ and ‘normally-distributed’ do not contradict each other.  ‘Random’ simply means that we can't predict exactly what the next value will be.  But that doesn't mean we can't make probabilistic statements about it.

Consider two experiments:

If you throw a (fair) die, then it could show any number from 1 to 6.  We can't tell which number will come up, but we can say that all numbers are equally likely (i.e. the distribution is uniform).
If you throw two dice and take their sum, that can be any number from 2 to 12.  Again, the sum is still random — we can't predict what it will be — but we can say that those values are not equally likely.  (For example, 7 is six times more likely than 12.)  So in this case it has a non-uniform distribution.  (You can plot all the probabilities; they take on a peaked shape a bit like a normal distribution.)

So there's no contradiction: both cases are random and have a known distribution.

In fact, most things that are random tend to have a non-uniform distribution: electrical noise, weather, the wait for the next bus, voting patterns…  Being able to make general statements about them without being able to predict the exact values is one of the strengths of statistics.

(As for why you often end up with a normal distribution, that's a result of the Central Limit Theorem, which says that when you combine many independent random variables, the result tends towards a Gaussian (normal) distribution.  So you see that crop up a lot.)


Wavelet Filter Coefficients from Scaling Filter Coefficients I am trying to develop a new type of wavelets and I found out a function that following a particular two scale relation. The function at a scale $t$ say $x(t)$ can be related in a finer scale $x(2t)$ by the following relation:

$$x(t)=0.1x\left(2t\right)+ 0.5x\left(2t-1\right)+0.8x\left(2t-2\right)+0.5x\left(2t-3\right)+0.1x\left(2t-4\right),$$

This is found out by trial and error method by simply changing the values. Now I know that the points $\begin{bmatrix}0.1&0.5&0.8&0.5&0.1\end{bmatrix}$ can be considered as  the scaling filter coefficients. Now the problem is finding the wavelet filter coefficients, I am not able to find out a particular relation connecting scaling filter coefficients and wavelet filter coefficients. Is there is a relation connecting the two ??

For perfect reconstruction you need:

$$ g \left[ n \right] * \hat{g} \left[ n \right] + h \left[ n \right] * \hat{h} \left[ n \right] = \delta \left[ n \right] $$

Where $ g \left[ n \right] $ is the Analysis Wavelet Function, $ \hat{g} \left[ n \right] $ is the Synthesis Wavelet Function, $ h \left[ n \right] $ Analysis Scaling Function and $ \hat{h} \left[ n \right] $ is the Synthesis Scaling Function.

Practically, Scaling Function is the LPF and Wavelet Function is the HPF.


Frequency Translation by Band Are there any mathematical transforms which can offset a select band of frequencies across the spectrum? (For sampled digital audio)

(crappy illustration)

What you are describing refers to Pitch Shifting.

A naive approach to pitch shifting is based on the shifting property of the Fourier Transform (and more formally here). According to the shifting property, multiplying a signal with a sinusoid in the time domain results in shifting in the frequency domain.

But, convincing pitch shifting takes a little bit more effort than that because it does not only matter to shift the frequencies of part of the spectrum to a new range, it is also necessary to preserve their relative phase for the sound to be recognisable. If you are trying to pitch shift musical instruments (or vocals), then it gets a little bit more challenging, because of the way pitch is working. In that case, you not only have to shift the components in the frequency domain but also preserve their intervals. In other words, it is not only enough to shift those formants that you show in your illustration, you also need to make sure that the new frequencies they end up into, have the same "distances" between them as the ones they are being pitch shifted from. This is not as straightforward as it...sounds.

Pitch shifting in this way implies a linear pitch scale. But the pitch scale is non-linear and going up in pitch (higher octaves) does not result in equally proportional frequencies.

Because of these reasons, more complex pitch shifting algorithms contain a "processing stage" between analysing and resynthesizing the pitch shifted version of the signal in their input. This processing stage attempts to take into account these relationships between the components of the spectrum according to some approach. For example, directly working out these relationships to the new range as the phase vocoder does, or attempting to fit a model to the spectrum and use it to reproduce it at a different range as spectral modelling does.

Hope this helps.


Output of a stable LTI system (discrete) Consider an LTI (linear and time invariant) system that is BIBO (bounded input
bounded output) stable and is such that $x[n] = 0$ for all $n < 0$ (note: this is sometimes referred to
as a relaxed system). Show that

Suppose that $x[n] = x_0[n] + au[n]$ for some constant a where $x_0[n]$ is a bounded signal with
$$\lim_{n\rightarrow\infty} x_0[n] = 0.$$ Then $x[n]$ is bounded and tends to a constant (namely, $a$). Show that the
output $y[n]$ will also tend to a constant.

Suppose that $x[n]$ has ﬁnite energy. Show that $y[n]$ will also have ﬁnite energy.

For 1, I'm kinda stuck at the step proving $$\sum_{k=-\infty}^{\infty}|x_0[n]h[n-k]| \rightarrow C $$ where $C$ is a constant. I can prove that it's below ($\le$) a constant, but how do I prove it is actually approach to it when $n \rightarrow \infty$?

For 2, can we do it without using the Fourier transform (continuous version here)?

So I'm answering my question here; turned out that I have to solve (2) then solve (1).

To prove $y[n]$ has finite energy, we write

$$\sum_{n=-\infty}^{\infty}y^2[n] = \sum_{n=-\infty}^{\infty}(\sum_{k=-\infty}^{\infty}h[k]x[n-k])^2 = \sum_{n=-\infty}^{\infty}(\sum_{k=-\infty}^{\infty}h[k]x[n-k])(\sum_{l=-\infty}^{\infty}h[l]x[n-l])$$
$$ = \sum_{k=-\infty}^{\infty}\sum_{l=-\infty}^{\infty}h[k]h[l]\sum_{n=-\infty}^{\infty}x[n-k]x[n-l]$$
Now we know that $x[n]$ has finite energy, so
$$ \sum_{n=-\infty}^{\infty}x[n-k]x[n-l] \leq \sum_{n=-\infty}^{\infty}x^2[n] = E$$
(a generalization of $a^2 + b^2 \geq 2ab$)
Then
$$ \sum_{n=-\infty}^{\infty}y^2[n] \leq E\sum_{k=-\infty}^{\infty}|h[k]|\sum_{l=-\infty}^{\infty}|h[l]| $$
For a BIBO system, $$\sum_{l=-\infty}^{\infty}|h[l]| \leq M \leq \infty$$
So we approve $$\sum_{n=-\infty}^{\infty}y^2[n] \leq M^2E \leq \infty$$

To prove $x[n]$ tends to a constant:

$$y[n] = a\sum_{k=-\infty}^{\infty}h[k]u[n-k] + \sum_{k=-\infty}^{\infty}h[k]x_0[n-k] = a\sum_{k=-\infty}^{\infty}h[k]+y_0[n]$$
Using the conclusion of 1, we know that $$\sum_{n=-\infty}^{\infty}x^2[n] < \infty \rightarrow \sum_{n=-\infty}^{\infty}y^2[n] < \infty$$
which implies that $$\lim_{n \rightarrow \infty} |y_0[n]| = 0.$$
So
$$ \lim_{n \rightarrow \infty} y[n] = a \times \lim_{n \rightarrow \infty} \sum_{k=0}^{\infty}h[k] = ac$$
where $c$ is a constant (because of BIBO).


Filtering out low-frequency bias I am trying to extrapolate the bearing of a vehicle using GPS and a gyroscope.

To this end, I first obtain the yaw component using gyroscope readings and the (known) yaw axis of the vehicle, then integrate the result over time and add it to the previous bearing.

I have implemented a real-time visualization of:

the direction of the vehicle yaw axis
the direction of the rotation vector reported by the gyroscope
raw GPS bearing
extrapolated bearing

About half of the time this works very well: the extrapolated bearing smoothly moves as the vehicle turns, and the GPS (which updates a second later) aligns with it almost perfectly. The other half of the time the gyroscope seems to suffer from bias, causing the extrapolated bearing to drift even when the vehicle is stationary.

The visualization shows me that the gyroscope rotation vector changes direction at random as long as I'm not suffering from the bias problem and the vehicle is not turning (indicating I'm getting something close to white noise from the gyroscope, which is what I'd expect). When the bias problem occurs, the gyroscope direction remains constant, corroborating my suspicion that I'm suffering from bias.

Based on the assumption that a vehicle will never turn in the same direction for more than about 1–4 minutes, and on the observation that bias comes and goes at longer intervals (several minutes), I tried using a highpass filter to eliminate bias. I used a 5th-order Butterworth filter.

With an initial cutoff frequency of $\frac{1}{240} \text{Hz}$ (corresponding to a 4-minute period), when I rotate the device by 90° at roughly $\frac{\pi}{2} \frac{\text {rad}}{\text s}$, as soon as the rotation stops, the bearing drifts in the opposite direction. I then started working downward from there: at $\frac{1}{1800} \text{Hz}$ (a half-an-hour period) the effect was still visible but much weaker, and at $\frac{1}{3600} \text{Hz}$ (one hour) it was virtually gone.

A real-world test at $\frac{1}{3600} \text{Hz}$, however, still resulted in random bias, of an even greater magnitude than what I'd experienced with raw values.

I then tried increasing the frequency to $\frac{1}{15} \text{Hz}$ and noticed I am likely suffering from overshoot, as the bearing now bounces back and oscillates around a certain value somewhere between the old and the new bearing.

I then went back to the original $\frac{1}{240} \text{Hz}$ and experimented with the filter order and characteristics: First I tried a 1st-order Butterworth filter, then a 5th-order Bessel filter – to no avail.

Where's the error here, and how can I elimiate bias without introducing any other, even less desirable artifacts?

I have decided to take a different approach: Rather than passing the raw accelerometer data through a filter and choosing between overshoot and ringing versus bias, I am now detecting bias by periodically comparing extrapolated $\Delta \psi _{gyro}$ to $\Delta \psi _{GPS}$:

Initially assume zero bias.
For each extrapolation step, subtract bias from measured $\omega _\psi$
On each new GPS bearing, calculate the change in bias as $\frac{\Delta \psi _{gyro} - \Delta \psi _{GPS}}{t}$
Add this change to the previous bias estimate
Determine the accuracy of the new bias estimate, considering the accuracy of $\psi _{GPS}$, approximate standard deviation of gyroscope output at rest and time elapsed
If we have a previous bias estimate (other than the initial assumption of zero), run the old and new bias estimate through a Kálmán filter; else simply continue with the new bias estimate
Return to step 2.

Results look good so far – this visibly reduces bias without introducing any new artifacts. Responsiveness still needs to be improved, which I can probably tackle by tweaking my accuracy model.


FFT of SIN waves with different phase delays I have come across a peculiarity of FFTs which has got me somewhat baffled.
I've simply summed up 101 sine waves and taken the FFT using this matlab script :

clear all
f=1e9;                                 % Centre Frequency 1GHz
df=2.5e6;                                % Carrier Frequency 2.5MHz
Time=linspace(-100e-9,100e-9,1000);      % Region of time to simulate over
delay=0;
Voltage=Time.*0;                        % Initialise Voltages to zero
for loop=-50:50                         % Sum 101 carrier Frequencies
    Voltage=Voltage+sin(2.*pi().*(f+df.*loop).*(Time-delay));
end
figure(1)                               %Plot Time dependent response
subplot(2,1,1)
plot(Time,Voltage)
subplot(2,1,2)                               %Plot Frequency Content
dt=Time(2)-Time(1);
frequency=linspace(-0.5/dt,0.5/dt,1000);
spectrum=fftshift(fft(Voltage));
plot(frequency,abs(spectrum))

The output is as I had expected, with the correct frequency content :

However, if I simply add a significant time delay (by re-running the script with delay=150e-9;, such that the main constructive interference lobe disappears outside the calculation window) the frequency content of the resulting time trace collapses to two peaks.

However the time trace is still the summation of 101 sin waves albeit now out of phase because of the introduced delay ?? Intuitively I would have expected the absolute frequency content of the trace to be preserved and only the phases modified by the delay. Upon reflection I can perhaps understand that the frequency content must be modified on energy conservation grounds, but can anybody rationalise what is going on here ?

You aren't looking at the frequency content of your sum of sinewave.  You are looking at the frequency content of a rectangular window on your sinewaves, and the window (FFT length) is shorter than the least common multiple of all the sinewave periods.

Signals that are not orthogonal within a window can cancel each other out partially or completely within that window.  What has happened is that you've chosen a window of a length such that each sinewave can get nearly completely cancelled out by the next higher and lower frequency pair of sinusoids of the "right" phase.  The ones on the ends are the only ones that aren't sandwiched, and thus aren't cancelled.


MAD and RMS SNR relation Recently I've found a formula relating mean average deviation (MAD) SNR to root mean square (RMS) SNR:

$$\text{SNR}_\text{MAD}  = \sqrt{\frac{2}{\pi}} \text{SNR}_\text{RMS} $$

which assumes Gaussian statistics. This is for a definition of MAD something like

$$ \text{MovingAverage}[\ |x - \text{MovingAverage}[x]|\ ]$$

for some signal $x$ measured in the time domain. My question is where does the above SNR formula come from? I'm pretty new to signal processing so maybe there's a very basic reference that describes this. Is there some intuition for why MAD SNR is smaller than RMS SNR?

The two main measures of scatter in the Gaussian case are the mean absolute deviation:

$$d_n = 1/n\sum |x_i-\overline{x}|$$

and the mean square deviation:

$$s_n = \sqrt{1/n\sum (x_i-\overline{x})^2}$$

If the random Gaussian noise has deviation $\sigma$, $s_n \to \sigma$ and   $d_n \to \sqrt{\frac{2}{\pi}}\sigma$. I don't have the exact references right now, you can find them in Robust Statistical Procedures, 1996, P. Huber, or texts by Hampel, or Rousseeuw, on robusts statistics.

The derivation is OK, let us use the zero-mean case: the expected value of the absolute value of the random variable:

$$E ( | x | ) = \int_{- \infty}^{ \infty} | x | p ( x ) d x$$

then you derive:

$$E(|x|) =
2\int_0^{\infty} \frac x {\sigma\sqrt{2\pi}}
e^{-\,\frac {x^2} {2\sigma^2} } dx$$

then

$$E(|x|) =
\left[-\,\sigma \sqrt{\frac 2 {\pi}}  e^{-\,\frac {x^2} {2\sigma^2} }\right]_0^{\infty}
=  \sigma \sqrt{\frac 2 {\pi}}$$

But this is for the mean absolute deviation, I am not sure of the equivalence with  the Median Absolute Deviation.


PDF of a Shifted Rectangular Pulse I wanted to determine the PDF of a Stochastic Process. I am familiar with the concept of PDF for a Random Variable which maps the outcomes to its probabilities but I am not able to find the PDF of a Random/Stochastic process. I have the solution but it’s pretty hard to visualize it. The solution which I have contains the PDF of the Stochastic process in each time intervals due to its nature.

Below is the stochastic process:

$x^{(l)}(t)= \operatorname{rect}((t-T_o)/T)$

Where To is uniformly distributed in the range
$|T_o|\leq T$

Now I want to determine the PDF of this random process. Hence I split the Random Process in three intervals:

$|t_o|\geq 3T/2$
$|t_o|\leq T/2$
$T/2\leq |t_o|\leq 3T/2$

The reason why its split into three intervals is that based on the condition of $T_o$, the process cannot occur after $3T/2$. And in between $-T/2$ and $T/2$ most of the processes occur and it has a decreasing slope between $T/2$ and $3T/2$.

Proceeding to the PDF, here is what I have as the solution:

Here is what I guess about the plot:

Region 1:

$|t_o|\geq 3T/2$

Here no processes occur and hence probability of 0 (no processes occur is 1). The probability of occurrence is 0.

Region 2:

$|t_o|\leq T/2$

Here most processes occur. But in this region, given the condition, I can change the process as:

$\operatorname{rect}(t-T_o)\leq T/2$

Here I need to choose the values of To so that the shift lies in this region(2). Hence could choose the value 0 till 1. At 1 it goes to region 1. Since I can choose any values between 0 and 1, I can assign each half the probability.

Region 3:

$T/2\leq |t_o|\leq 3T/2$

In this region, I need to choose the value with the probability of:
$(1/4)+(t_o/2T)$ to get zero shift and $(1-((1/4)+(t_o/2T)))$ to get a shift.

I am not quite sure of the above reasoning behind this graph of PDFs based on which I had made the above conclusions in each regions. Do the graphs of the PDFs need to be interpreted in a different way?

If anyone could provide the reasoning behind any one region (region 2 is highlighted) and how to plot the PDF in that region (region 2), I could formulate the others in a similar fashion.

A random process is a collection of random variables, one for each time instant.  Here, the $t$-th random variable is denoted by $X(t)$ where the assumption is that $t \in (-\infty,\infty)$ and what you call the PDF is actually a whole collection of PDFs, one for each of the (uncountably infinitely many) random variables $X(t)$.  Now, we are told for any given (fixed) value of $t$, say $t=0.313$, the random variable $X(t)$ (specifically $X(0.313)$ in our chosen example) is a function (namely a rect function here) of the random variable $T_0$.
$$X(t) = \operatorname{rect}\left(\frac{t-T_0}{T}\right); \qquad
X(0.313) = \operatorname{rect}\left(\frac{0.313-T_0}{T}\right)$$  Now, the rect function can only take on values $0$ or $1$ and so we immediately can say that regardless of the numerical value of $t$ (doesn't have to be
$0.313$) or whatever value the
random variable $T_0$ might have assumed, $X(t)$ can only take on values $0$ or $1$, that is,

  for each value of $t \in (-\infty,\infty)$, $X(t)$ is a Bernoulli random variable.

As such, $X(t)$ does not have a PDF but has a PMF (probability mass function), unless of course you like to use impulses and define what are called discrete density functions to incorporate both notions into one.

The only question remaining is what is the parameter $p$ (the probability that $X(t)$ equals $1$) of this Bernoulli random variable,
and the value of $p$ depends on what value we assume for $t$. That is, the answer is not something like $p = 0.3$ but rather a multi-line equation like
$$p(t) = \begin{cases}\cdots & \text{if}~ t\in \cdots,\\
\cdots & \text{if}~t\in \cdots,\\
0, &\text{otherwise.}
\end{cases}$$
To figure out the details, note that $\operatorname{rect}(\cdot)$ has value $1$ exactly when its argument has value $\leq \frac 12$. Thus, as a function of $t$, $\operatorname{rect}\left(\frac{t-T_0}{T}\right)$ has value $1$ exactly when $$T_0-\frac 12T \leq t \leq T_0+\frac 12T.$$
"Wait a minute", you splutter, "Didn't you say that $t$ had a fixed value
such as $0.313$? What's all this about as a function of $t$??? $t$ is not a variable; it is a fixed value" Well, we need to turn the thing around a bit. Notice that for a fixed value of $t$ (happy now?) $X(t)$ has value $1$ if and only if $T_0$ is such that
$$t - \frac 12T \leq T_0 \leq t+\frac12T,$$
that is, for the chosen fixed value of $t$, and keeping in mind that
$T_0 \in [-T, T]$
\begin{align}
P\{X(t) = 1\} = p(t) &= P\left(t - \frac 12T \leq T_0 \leq t+\frac12T\right)\\
&= P\left(T_0 \in \left(\left[t- \frac 12T, t+ \frac 12T\right]\cap
\big[-T,T\big]\right)\right)
\end{align}
That intersection on the right depends on what we have chosen $t$ to be,
and so we consider various cases. For brevity, denote the interval
$\left[t- \frac 12T, t+ \frac 12T\right]$ as $\mathcal A$ and the interval
$[-T.T]$ as $\mathcal B$

If $\color{blue}{t < -\frac{3T}{2}}$, then $t+\frac 12T < -T$ and so the right end of the interval $\mathcal A$ is smaller than the left end of $\mathcal B$
and so $\mathcal A \cap \mathcal B = \emptyset$. We get that $p(t) = 0$
for $t \leq -\frac{3T}{2}$. Similarly, for $t > \frac{3T}{2}$, the
left end $t-\frac 12T$ exceeds $T$, the right end of $\mathcal B$. So

  $p(t) = 0$ if $|t| > \frac{3T}{2}$.

If $\color{blue}{-\frac{3T}{2} \leq t < -\frac T2}$, then the right end $t+\frac T2$
is in $[-T,0]$ but the left end $t-\frac T2$ is smaller than $-T$, that
is, $\mathcal A \cap \mathcal B =  \left[-T, t+\frac T2\right]$.  Since
$T_0 \sim U[-T,T]$, we get that
$p(t) = \frac{1}{2T}\times \left(t+\frac T2 - (-T)\right)
= \frac{1}{2T}\times \left(t+\frac {3T}{2}\right)$. Note that
$p(t)$ increases from $0$ at $t=\frac{3T}{2}$ to $\frac 12$ as
$t$ approaches $-\frac T2$.  Similarly,
for $\frac T2 < t \leq \frac {3T}{2}$, $p(t)$ decreases from
$\frac 12$ to $0$ linearly.  We have that

  $p(t) = \frac{1}{2T}\times \left(\frac {3T}{2} - |t|\right)$ if
  $\frac T2 < t \leq \frac{3T}{2}$

If $\color{blue}{|t| \leq \frac{T}{2}}$, then both the left end
$t-\frac T2$ and the right end $t+\frac T2$ of $\mathcal A$ are in
$\mathcal B$, and so

  $p(t) = \frac 12$ if $|t| \leq \frac{T}{2}$.

The process in question is not what is called a stationary process and it does not have a PDF (or PMF) but rather a whole bunch of different PDFs (or PMFs), one for each time instant $t$. Read the first part of
this answer to learn a little more about these ideas.


Fast DCT implementation I'm having trouble figuring out how to follow the fast 8x8 DCT algorithm diagrams found in the following two papers:

(1) A Fast Computational Algorithm for the Discrete Cosine Transform by Chen et al.

and

(2) Practical Fast 1-D DCT Algorithms With 11 Multiplications by Loeffler et al.

In particular, the second diagram that shows the algorithm in (2) looks like the following:

The description of the operations in this algorithm are:

There are a few questions I have about this formulation, and I'm not sure where to find the answers:

(2) suggests that this algorithm generates a DCT that is scaled by some value $C = \sqrt{2}$. It mentions that this $C$ was chosen arbitrarily to avoid any multiplications in computing the DC coefficient. Really the only requirement is that $C_{DCT} * C_{IDCT} = \frac{4}{N^2}$. So my question is this: What is the scaling factor of the output coefficients using this algorithm? It seems like they are different than the original definition of the DCT, but I don't know by how much (mostly because I don't actually see any relationship between this diagram and the original formulation of the DCT):
$$ F(k) = \frac{2c(k)}{N}\sum_{n = 0}^{N - 1}f(n)\cos\left(\frac{\left(2n + 1\right)\pi k}{2N}\right)$$
where $c(k) = \frac{1}{\sqrt{2}}$ for $k = 0$ and $c(k) = 1$ for $k \neq 0$.
The paper states that performing the IDCT can be done using the exact same algorithm but transforming outputs to inputs and vice versa. First, should the DCT coefficients be ordered in bit-reverse order prior to running them through the IDCT? Second, for the rotation blocks (the squares in the diagram), shouldn't the inverse operation be:
$$ \begin{align}
O_0 = I_0 \cdot k \cdot \cos\frac{n\pi}{2N} - I_1 \cdot k \cdot \sin\frac{n\pi}{2N} \\
O_1 = I_1 \cdot k \cdot \sin\frac{n\pi}{2N} + I_1 \cdot k \cdot \cos\frac{n\pi}{2N} \\
\end{align}
$$
My reasoning is this: The inverse of a rotation by $\theta$ is a rotation by $-\theta$. Hence, we just replace the angle by its inverse and use the identities $\cos(-\theta) = \cos(\theta)$ and $\sin(-\theta) = -\sin(\theta)$. Third, what is the scaling factor of the transformed values after the IDCT? (2) says $\frac{2}{N^2}$, but empirically, this hasn't produced correct results.
Suppose after I run the algorithm, I have the result of each lane stored in the values d0 ... d7. Which of the following is correct:
output[0] = d0               or               output[0] = d0
output[4] = d1                                output[1] = d4
output[2] = d2                                output[2] = d2
output[6] = d3                                output[3] = d6
output[7] = d4                                output[4] = d7
output[3] = d5                                output[5] = d3
output[5] = d6                                output[6] = d5
output[1] = d7                                output[7] = d1

If there are any ways to improve this question, or if I should ask elsewhere, please let me know.

Alright, after some days of staring at this problem I hope I can provide a bit of guidance to the next poor soul.

Yes the scaling is different. Compared to scipy.fftpack.dct the DC term is $\frac{1}{2}$ and the other terms $\frac{\sqrt{2}}{2}$. But apparently it all nicely cancels out in the inverse transform.
The inverse input order is exactly as they come out: bit reversed. Literally as if you flip around the figure and connect the lines. And yes, you are correct that the sin gets a minus. I'm seeing a scaling factor of 8, FWIW.
The output order of the inverse is the same as the input order of the forward transform. So $x[n]=\frac{\mathrm{IDCT}(\mathrm{DCT}(x[n]))}{8}$

Also note that there is an error in the graph and it's $\sqrt{2}c6$ for the even-side rotation block.


Finding Fourier Series Coefficients I'm just beginning to learn about Fourier series and I'm trying to figure out how to find the Fourier series coefficients for $$x(t) = e^{j100\pi t}$$

I know that $$x(t) = \sum_{-\infty}^{\infty} a_{k} e^{jk(2\pi/T)t}$$

How do I go about finding the coefficient(s)? I know that I can get T = 1/50, but beyond that I don't even know where to begin. I think I'm supposed to be able to do this just by looking at it without having to solve the integral equation for $a_k$, but I don't know what I'm supposed to be looking for.

The answer is right there in the question. You have expressed the signal as linear sum of scaled complex exponentials.

Now the question is fourier series coefficients of a complex exponential$$e^{j100\pi t}$$

So you can see from the summation of x(t) that you have mentioned in the question that there is only one fourier series coefficient $a_1$ with a fundamental frequency $\omega_0$ because there is only one complex exponential in the input.The value of $\omega_0$ is $$ 100\pi$$. So $f_0$ = $(100\pi)/(2\pi)$ = 50

The value of $a_1$ is 1.

In other words you can express $$x(t)=(1).\ e^{(100 \pi)(1)t}$$ to see that  k=1 and $a_1$=1.


Why do we care about the multiplicity of poles and zeros in rational Z-transforms? In the DSP class I'm taking, a lot of the questions ask me to list the multiplicity of a pole or zero in a rational Z-transform. For example, the multiplicity of the zero at $z=1$ in :$\frac{(z-1)^2}{(z-5)}$ is $2$.

What knowledge do we gain by knowing the multiplicity of some pole or zero?

Of course, if we know the multiplicity of all the poles and zeros, we can describe the z-transform of a rational function, but this is not a very satisfying answer.

Multiple poles have amplitude and phase profile different to that of single poles.