Optimization of execution time of audio capture call in c/c++ (ALSA) I am working on a Audio Deep Learning Project. In this C++ project I need to capture real-time audio and then do pre-processing and the pre-processed output is fed to the Deep Learning Network.

I am fetching 50 ms audio data in each call at 16000 sampling rate, i.e. 800 samples are fetched in the audio capture call.

Sometimes this call takes more than 50 ms, my processing code is taking ~30 ms.

If the audio-capture + processing takes more than 50 ms, is that mean I lost some audio data?

Can the audio capture call be optimized?

Implementation Info:-
I have the working code on raspberry pi and my processing code is armv7 optimized.

Execution time is measured using std::chrono in c++ using high_resolution_clock.

Capture audio into a FIFO, and process the samples in that FIFO in a separate thread.  That way, timing variations in your processing thread won’t block the audio input.

Very commonly the audio input FIFO is implemented using a lock-free ring buffer so it never blocks.  A separate processing thread can poll the buffer at some suitable fixed rate (video frame rate or network data rate for a given MTU size, etc.)


Adding Salt n Pepper noise into an Image manually When i looked in to the algorithm of Salt n Pepper noise in a Matlab i've found this,

case 'salt & pepper' % Salt & pepper noise
    b = a;
    x = rand(sizeA);
    d = find(x < p3/2);
    b(d) = 0; % Minimum value
    d = find(x >= p3/2 & x < p3);
    b(d) = 1; % Maximum (saturated) value

Which i simplified as this,

% [a, code, classIn, classChanged, p3, p4]
% p3 = density
% p4 = variance
% a = Image
% code = noise Type
% classIn = Image Class i.e unit8 , unit16, double

% -------- ALGORITHM ----------%

%     b = a;   <-- Assign b to the input image
%     x = rand(sizeA); <--- Generate random pixels from the image pixels
%     d = find(x < p3/2); <--- Find the pixels whose values are less than half of the mean value
%     b(d) = 0 <-- Implement minimum noise to them
%     d = find(x >= p3/2 & x < p3) <--- Find the pixels whose values are
%                                        greater than half of the mean
%                                        value & less than mean value
%     b(d) = 1;     <-- Implement maximum noise to them

instead of using imnoise, I am trying to add Salt n Pepper noise manually into the image using the above algorithm,

I = imread('2.jpg');
J = rgb2gray(I);
p3=0.05 %default
b = J;
x = rand(size(b));
d = find(x < p3/2);
b(d) = 0; % Minimum value
d = find(x >= p3/2 & x < p3);
b(d) = 1; % Maximum (saturated) value

I don't know where the final image is stored as when i use imshow for b it show a simple image without salt n pepper noise.
Where the final image is stored ? how can i see the it ?

P.S: What does it mean by image class ? that is unit8 , unit16 ?

It seems that the final image is in the variable "b". Observe that the max (salt) and min (pepper) values are respectively 1 and 0. This indicates that your original image needs to be an intensity image with graylevels normalized to [0,1]. This may be the reason why you don't see a sensible result when you display "b" as an image.

In case your image is grayscale {0,...,255} i.e. 8-bit unsigned integer format (unit8), you can change your code to assign 0 and 255 instead of 0 and 1.

I suggest you read the documentation on various image storage types here: http://www.mathworks.com/help/toolbox/images/f14-13543.html

It can get pretty confusing at times, mainly because there are so many names (intensity image, grayscale image, RGB, HSV, binary image, indexed image, etc.)


Synthetic Aperture Radar (SAR). How to do Range Compression? I am trying to perform range compression to a raw SAR image [fast times, sensor position], Dr in code. I have the transmitted signal (chirp), g in code, and literature says range compression is obtained by convoluting the time reverted conjugate of g along the rows. I use the following MATLAB code.

    h_compressed = conv2(conj(fliplr(g), 1, Dr, 'same');

Where 3003 is the number of columns of the RAW image. Nevertheless my output is the same as my input i.e. no range compression occurs. Any ideas on why?

After range compression I perform 1D matched filtering.

    for k = 1:Nr
    R = sqrt(x_filter.^2 + SAR_data.r_ax(k).^2);
    filter = exp(-1i*4*pi/lambda*R);
    Dfr(:,k) = conv2(Drc(:,k),filter','same');
    end

And obtain the following:

What format is your input data in? It sounds like you have a phase history signal as it's often called. Assuming that the transmitted waveform was a linear FM chirp, has it already been deramped? Some receivers perform deramp-on-receive processing, which removes the chirped nature of the waveform. In that case, range compression is as simple as performing a 1-D FFT on each pulse return (which appears to be one per column in your example).

If you instead have a raw collection from the receiver that hasn't been deramped, then you expect to have an input signal that consists of many reflections of the transmitted waveform, separated in time according to the round-trip delay to the various reflectors in your scene. In this case, in order to perform range compression, you need to generate and apply a matched filter to each pulse return (again, a 1-D convolution of each column of your dataset with the matched filter).


DFT: a function of $n$? I‘m a high school student and I haven’t studied physics or anything.

Why does the DFT depend on an integer, say $k$ or $n$ (it’s usually expressed like $F(n)=...$ or $F(k)$ or $F_k$, etc.) if it is supposed to deliver a frequency information of a sampled signal?

Can the frequency content of the signal be expressed as a multiple of the integer?

Ah okay thank you. I think what irritated me was that I didn’t consider that the sampled function is a periodic or periodically continued function so that the Fourier transform delivers amplitude values for the base frequency and multiples of it only (so the Fourier transform is discrete), which are indexed as n*(base frequency), so that it is obvious why F „depends“ on n. The amplitude value for n can then be converted into an amplitude value of a frequency given by n*(base frequency) (and the base frequency can be obtained by the period since a signal is only sampled between its starting point and endpoint, whose distance is the period).


Compressed / Compressive Sensing - Sensor Placement I have a beam fixed at one end and free at the other. It obviously results in vibrations due to the disturbances in the environment or when induced. This system can be represented using a differential equation. The vibrations are sensed using sensors, analyzed and reduced. Now, the sensing is done through compressed sensing technique.

How is it possible to find the optimal sensor placement (minimal number of sensors along with the details of its spatial position) just with the help of the measurement matrix (that we design in compressed sensing for reconstruction purpose)?

Suppose we want to retrieve $\mathbf x$ (vibrational info) from: $\mathbf{y= Cx}$, where $\mathbf C$ is the fat matrix (measurement matrix * sparse matrix). I want to find a relationship that the measurement matrix can hold with the placement of sensors. It's like an inverse problem, instead of designing a measurement matrix with respect to the sensors placed, I want to place the sensors with the knowledge I have about the measurement matrix.

Suppose your beam has ten vertical degree of freedom. You can only apply CS to compress measurement vector of each DOF. You can not reconstruct all of the responses using limited number of them.
The ten-dof structure would have ten pick- related to ten natural frequency- in frequency domain, so its sparsity has been already guaranteed. Hope this help.


Find the frequency response of the system How can I find the frequency response of the following linear system

OK, in the comments you've shown that you can actually do it by yourself. This answer just summarizes the steps and shows you how to write your (correct) result in the form that you were given as an answer.

First of all, note that you have two identical systems in series. If $G(\omega)$ is the frequency response of the system up to and including the first integrator (from the left), then the total frequency response is

$$H(\omega)=G^2(\omega)\tag{1}$$

The response $G(\omega)$ can be seen as a multiplication of two subsystems, the system up to and including the adder, and the integrator. The first of the two has a frequency response

$$G_1(\omega)=1-e^{-j\omega T}\tag{2}$$

Since $G_1(\omega)$ has a zero at $\omega=0$, we don't need to bother with any problems of the integrator at $\omega=0$ (I'm talking about the delta impulse). So for calculating $G(\omega)=G_1(\omega)G_2(\omega)$, we can simply use

$$G_2=\frac{1}{j\omega}\tag{3}$$

which gives

$$G(\omega)=\frac{1-e^{-j\omega T}}{j\omega}\tag{4}$$

and, from (1),

$$H(\omega)=\frac{(1-e^{-j\omega T})^2}{(j\omega)^2}\tag{5}$$

In order to get from (5) the answer that you were given, use this trick to rewrite $G_1(\omega)$:

$$G_1(\omega)=e^{-j\omega T/2}(e^{j\omega T/2}-e^{-j\omega T/2})=2je^{-j\omega T/2}\sin(\omega T/2)\tag{6}$$

Using (6) you can rewrite $G(\omega)$ as

$$G(\omega)=2e^{-j\omega T/2}\frac{\sin(\omega T/2)}{\omega}=
Te^{-j\omega T/2}\frac{\sin(\omega T/2)}{\omega T/2}\tag{7}$$

Squaring (7), using $\omega=2\pi f$ and $\text{sinc}(x)=\sin(\pi x)/(\pi x)$ will give you the final form of the result.


Bandpass filters with python for low frequencies I have a signal, sampled at 1 kHz, that I would like to analyze in third octave bands. For that purpose I defined Butterworth filters with proportional bandwidth. This works well for higher frequencies, but produces very weird results for low frequencies. What am I missing?
Plotted are the frequency responses of the Butterworth filters.
I find similar behaviour in Hanning filters, and think that I am missing something fundamental.

Code used to define filters:

import scipy
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt
def butter_bandpass(lowcut, highcut, fs, order=5, label=None):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = signal.butter(order, [low, high], btype='band')
    w, h = signal.freqz(b,a,worN=2000)
    plt.plot((fs * 0.5 / np.pi) * w, abs(h), label=label)
    return b, a

center_freqs = np.array([0.5* 2 ** (n * 1 / 3.) for n in range(0, 29)])
center_freqs.sort()
lower_freqs = 2. ** (-1 / 6.) * center_freqs
for lf in lower_freqs:
    butter_bandpass(lf, lf*2**(1/3.), fs=1000, order=5)

plt.show()

I would bet this is just numerical error in the transfer function.  Try using butter_sos = butter(..., output='sos') instead of ba format, and sosfreqz(butter_sos, ...) instead of freqz.  Does that fix it?

https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.signal.sosfreqz.html

Edit: I just tried it and it does.  :D

import scipy
import numpy as np
from scipy import signal
from matplotlib import pyplot as plt
def butter_bandpass(lowcut, highcut, fs, order=5, label=None):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    sos = signal.butter(order, [low, high], btype='band', output='sos')
    w, h = signal.sosfreqz(sos,worN=20000)
    plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label=label)
    return sos

center_freqs = np.array([0.5* 2 ** (n * 1 / 3.) for n in range(0, 29)])
center_freqs.sort()
lower_freqs = 2. ** (-1 / 6.) * center_freqs
for lf in lower_freqs:
    butter_bandpass(lf, lf*2**(1/3.), fs=1000, order=5)

plt.show()


How to undo dynamic Doppler effect by software? I have a signal affected by Doppler effect. In general, the Doppler effect changes the scale of the time and frequency domains. So, to undo this phenomenon, I am using a sampling rate converter. Two possible implementations using polyphase structures are shown in Figures 1 and 2. They both work fine as expected. The recovered signal has no Doppler shift.

.

However, while in presence of dynamic Doppler, i.e. the compression/expansion factor varies over time, I am not sure if this solution is appropriate. I am buffering the signal, and for each buffer I compute the corresponding $L$ and $M$ factors that undo the Doppler effect. Nevertheless, there exists a trade-off between buffer length and the frequency resolution. In other words, in a buffer I still have a varying Doppler shift, so the $L/M$ ratio should not be a constant during a single buffer.

Is there a way of apply a sample-to-sample method to undo dynamic Doppler effect? I am currently interested in software-based solutions, not hardware-based ones. For example, in hardware I would dynamically modify the sampling rate. What about in software? How can I dynamically modify the sampling rate?

References:

Both figures were taken from:

R. E. Crochiere and L. R. Rabiner, Multirate Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1983

If the only reason for the buffer is so your algorithm has something to process, then the simplest way to implement this is to keep your current algorithm but use a shorter buffer. If the larger buffer was present for other reasons (i.e. data arrives in large chunks), then keep that larger buffer, but feed it into a shorter buffer that is given to your processing code. (It seems any algorithm able to recognize Doppler shift in a given snapshot must also be able to process that same snapshot.) In this way you can keep your current algorithm, and tailor the short buffer length as desired. I can't imagine how an algorithm specifically designed to do this per sample would be of benefit, unless it takes less CPU--which seems unlikely.


JPEG compression steps after quantization I have a 3-channel (for colours) a png image that I opened and

I splitted the image into 8x8 blocks
I applied all of the blocks discrete cosine transform
And then applied quantization
I stored the values in an array by zigzag traverse
I do not know what to do to reduce the size in this step

Now I do not really understand what i am supposed to after this steps. How did I compress this image? If I applied inverse of these steps and put the values in a new image, I believe the size will not change at all. I mean I am supposed to save the values in a txt and decode them with my application so that I actually made a JPEG compression?

Ok what's left is the last stage; namely Huffman VLC (variable length code) coding of those quantized DCT coefficients. And indeed it's true that the actual bit savings are obtained in this last stage.

As you may recognize, after the zig-zag scanning stage, there should be long chains of zeros among those 64 quantized DCT coefficients. By taking advantage of those chained zeros, you will encode those zig-zag scanned coefficients according to some CAT and SYM codewords, specifically designed and provided by the JPEG standards committee as a Huffman Table.

The procedure is longer than I wish to lay out here, but not very complex. I highly recommend anyone interested in implementing a standard JPEG the following books

1-) Stanard Codecs_GHANBARI
2-) Introduction to Data Compression_SAYOOD.


Why is the error between the desired signal and estimated signal in the case of LMS filter remaining constant even after n number of iterations I am giving white noise as input to an adaptive filter which is initialized to zero (value of filter coefficients of adaptive filter is 0). I am getting a desired response $d(n)$  by passing white noise through an unknown channel of impulse response $h(n)$.

Now using LMS algorithm, I am trying to update the adaptive filter coefficients by using the equation
$$w(n+1)=w(n)+\mu u(n) e(n)$$
where $\mu$ is the step size or adaptation speed of the algorithm, $u(n)$ is the white noise which is passed to a speaker and $e(n)$ is the difference between desired response and estimated response.

Even after $n$ number of iterations, why is that the error signal, which is the difference between desired response and estimated response, not decreasing but remaining almost constant.

As mentioned in the comment, I modified the code given here and was able to adapt the LMS filter with error tapering to zero. The only assumption I made is that (since I am not an audio expert and do not know how the channel from speaker to the microphone would look like), I assumed a 10 tap channel with only first 3 non-zero values (multi-path reflection from walls). OP is free to use a channel of his/her choice. Using this I generated 'actual_fb_path' by convolving white-noise with assumed channel. If I get access to 'actual_fb_path' it would be great. I would request OP to modify this as per the system model to see if it adapts to the delay.

clc
clear all
close all
training_input = rand(1,100);
h_actual = [0.8 0.3 0.05 zeros(1,7)]; %channel assumption
actual_fb_path = conv(training_input, h_actual);
d = actual_fb_path;
mu = 0.1;
%ha = lms(40,mu);
%[y,e] = filter(ha, training_input, d);
lms = dsp.LMSFilter;
lms.StepSize = mu;
lms.Length = 10;
[y,e]=lms([training_input zeros(1,9)]', d');
%subplot(2,1,1)
N=length(training_input) + length(h_actual)-1;
figure()
plot(1:N,d,'r',1:N,y,'b',1:N,e,'g')
title('System')
legend('Desired','Output','Error')


Intuition behind image derivative using Fourier Transform for edges detection This equation can be shown mathematically:

$\frac{\partial f}{\partial x}=\frac{2\pi i}{N} \mathcal F^{-1}\left(u\cdot \mathcal F(f(x,y)\right)$

I am struggling to understand the intuition behind it when talking about 2D fourier for images and the use of derivatives to detect edges.

How does multiplying the Fourier transform of the image with $u$ help us detect edges in the x direction? Why does it make the "edges" white and the rest of the image darker?

Any help would be appreciated.
Thank you!

In general, the time derivative property of the Fourier Transform is given as

$$\mathscr{F}[\frac{d}{dt}x(t)] = j\omega X(j\omega) $$

Notice that we can simply multiply by the frequency index in the Fourier Transform result.

For the 2D FT result:

$$\mathscr{F}[f(x,y)]= F(u,v)$$

Using the same property results in:

$$\mathscr{F}[\frac{d}{dx}f(x,y)]= uF(u,v)$$

$$\mathscr{F}[\frac{d}{dy}f(x,y)]= vF(u,v)$$

The reason this detects edges is clear if you observe how the differentiation of a 1D step function results in an impulse. Your edges approximate steps. The derivative is the slope, so where the slope is high (edges) the derivative is large.


Period of a quasiperiodic digital signal I have a signal made of a perturbed square wave, sampled so that there are at least six samples per period, for a total of 15 to 50 periods. The sampling frequency and the signal frequency are unrelated.

The perturbations are of several types:

jitter in the period (say 10%),
variable amplitude of the squares (say from 1 to 3),
additional arbitrary signal with an amplitude below 1, often localized.

I am looking for a simple, fast but robust computational procedure to estimate the average period (preferaby without a DFT or autocorrelation).

Typical example:

You can locate each rising edge and use that as the start time of each new cycle. Then compute the time differences between their locations and then take reciprocal of the average (or median, for more robustness against outliers) of these time differences to estimate frequency.

Locating the rising edges can be achieved by checking if the difference in consecutive data points exceeds a predetermined threshold. This threshold can be set adaptively from data. For instance, if you make a scatter plot of all your data samples it should look bimodal: the higher values correspond to the peaks of the square wave and the lower values to the valleys.


Why does the autocorrelation get its peak at zero? I know that zero shifting in the autocorrelation function is equal to its energy, yet, I would like to understand why the peak is at zero.

Are you looking for a formal proof or the intuition behind this? In the later case: "Nothing can be more similar to a function than itself". Autocorrelation at lag $\tau$ measures the similarity between a function $f$ and the same function shifted by $\tau$. Note that if $f$ is periodic, $f$ shifted by any integer multiple of $\tau$ and $f$ coincide, so the autocorrelation has a comb shape - with peaks at the integer multiples of the period with the same height as the central peak.


What Is the Intuition of Convolution in The Signal Processing World It is known that an LTI system is defined by its Impulse Response and the Convolution Operator.

Could some put some intuition behind?

Thank You.

I think this post gives a lot of intuition about Convolution: - What's convolution intuitively?
In the Signal Processing world, an LTI (Linear and Time Invariant) system basically scales and sums delayed versions of the input signal.
This can be mathematically described using the integral known as convolution.


Essential Matrix from line correspondences I'm looking for some library with a function to compute an essential matrix from line correspondences in multiple images. All the computer vision libraries I know use point correspondences to compute the essential matrix.

I've seen some papers about how to calculate E from line correspondences but implementing it myself would be a project on itself, so I prefer to try some existing implementation first. Does anybody know of a library with such an implementation available somewhere?

When I think of lines in correspondence, my mind immediately goes for Pluecker representations, which are very suitable for such estimation tasks. In fact this work exists for camera pose estimation from lines and might be a good starting point:

  Přibyl, Bronislav, Pavel Zemčík, and Martin Čadík. "Camera Pose Estimation from Lines using Pl\" ucker Coordinates." arXiv preprint arXiv:1608.02824 (2016).
  (link to the project page)

As you mentioned, there are various resources, especially when your camera is calibrated (I assume this because you ask for the essential matrix). If you really like to have a quick and dirty solution, you could either use the endpoints of the line segments (if you have them) or the intersection points of the lines.

In fact this just reminded me of one nice work, which generates virtual points on cameras along the lines and uses them for pose computation:

  Line-Sweep: Cross-Ratio for Wide-Baseline Matching and 3D Reconstruction
  S. Ramalingam, M. Antunes, D. Snow, G. H. Lee, S. Pillai
  Computer Vision and Pattern Recognition (CVPR), 2015
  (link to the project page)

Let's keep here posted with updates.


Can't make sense from VHDL butterworth filter implementation I am trying to understand the VHDL implementation of the 3rd order Butterworth filter published on opencores. Supposedly the filter implementation gets away with only division and multiplications by 2 which is cheap when implemented in hardware (only shifting the bits left/right).

I've translated the VHDL to MATLAB to simulate a step response.

x = int16([zeros(1,500) 512*ones(1,1000)]); % step-input

a = zeros(1,3); % accumulator registers
w = zeros(1,5); % analog filter state variables

s = 4; % scaling parameter to adjust cutoff frequency

for t=1:length(x)
  in = x(t);

  a1 = a(1) + w(1) - w(3);
  a2 = a(2) + w(2) - w(4);
  a3 = a(3) + w(3) - w(5);

  w1 = in - w(2);
  w2 = a(1)/power(2,s);
  w3 = a(2)/power(2,s+1);
  w4 = a(3)/power(2,s);
  w5 = w(4);

  % update registers for next cycle
  a = [a1 a2 a3];
  w = [w1 w2 w3 w4 w5];
  y(t) = 2*w(5);
end

The step response looks quite promising; dotted line is the step input, solid red, orange and blue are the filter response for s=3, 4 and 5 respectively.

However, the implementation really does not make much sense to me. Can anyone enlighten me? What are the coefficients of the filter? Which form does the implementation follow?

The filter structure is a digital leapfrog and the structure looks like this picture (note: the picture is a different order than the code):

These filters are discussed at some length on wikipedia and wikibooks.


Determining asymptotic stability using transfer function? In an exam task, I am asked to determine the transfer function of the following direct-time system and decide whether it's stable.

I think this system is canonical and the amplifiers 'on top' correspond to $b_2$, $b_1$ and $b_0$ from left to right in the following formula for the transfer function. Also, the amplifiers 'on the bottom' correspond to $-a_2$ and $-a_1$ similarly from left to right:

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

Therefore the transfer function of this system is

$$ H(z) = \frac{z^{-1} + 0.4 z^{-2}}{1 + (-0.5) z^{-1} + 0.06 z^{-2}} = \frac{z + 0.4}{z^2 + (-0.5)z + 0.06} $$

So it has a single zero at $z = -0.4$ and its poles are $p_1 = 0.3$ and $p_2 = 0.2$

I know I can draw the pole-zero plot for this function that will look like this:

How do I deduce the stability of the system from here?

I have learned things before like given the eigenvalues $\lambda_i$ of the system's $\underline{\underline{A}}$ matrix, a discrete-time system is asymptotically stable if $\forall \lambda_i : |\lambda_i| < 1$. Similarly, in continuous time, the criterion is $\forall \lambda_i : \mathrm{Re}\{\lambda_i\} < 0$.

But now I have poles, not eigenvalues. I have only found that if there is a pole that is to the right of the imaginary axis on the pole-zero plot, the system is unstable. However, here, I think this system should be asymptotaically stable (therefore also BIBO stable), but it has a pole in such a place.

Why is the transfer function given this way for such systems?

Poles in the right half-plane only cause instability if the system is a continuos-time system, because in that case, a pole $s_{\infty}$ contributes a signal

$$k\cdot e^{s_{\infty}t}u(t)\tag{1}$$

to the output signal, with some constant $k$, and with $u(t)$ denoting the unit step function. For $(1)$ to decay we require $\text{Re}\{s_{\infty}\}<0$, i.e., for a (causal) continuous-time system to be stable, all its poles must lie in the left half-plane.

In discrete time, however, a pole $z_{\infty}$ contributes a sequence

$$k\cdot z_{\infty}^nu[n]\tag{2}$$

to the output signal. Consequently, in order for $(2)$ to decay, we require $|z_{\infty}|<1$, i.e., for a (causal and) stable discrete-time system, all poles must be inside the unit circle.

Concerning the equivalence of poles and eigenvalues of the system matrix have a look at this answer.


Right Inverse DFT notation Inverse DFT is given as:
\begin{equation}
y_n=\sum_{k=0}^{N-1}Y_k e^{j\frac{2\pi nk}{N}}\tag{1}
\end{equation}
and knowing that Fourier transform of a real value signal has Hermitian symmetry property, i.e., $Y^*_k=Y_{N-k}$ with $\frac{N}{2}$ positive and $\frac{N}{2}$ negative frequencies, Eq. (1) can be rewritten in which of following ways:

Notation 1:
$$
y_n=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}
$$
Notation 2:
$$
y_n=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=1}^{\frac{N}{2}}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}
$$
Notation 3:
$$
y_n=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}
$$

I can't see any match with the initial $y_n$ in equation $(1)$, (see update). With a slight change to the notation, you have the following results in the first two cases, for Notation 1 and 2 respectively:

Notation 1:
\begin{align}
y_{n_1}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{1a}$}\\
&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_k^* \overbrace{e^{j2\pi n}}^{=1} e^{-j\frac{2\pi nk}{N}}\\
&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\left(\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}\right)^*\\
&=y_{n/2}+y^*_{n/2}\tag{$\scriptstyle{\text{$y_{n/2}=y_n$ for $n=0, \ldots, N/2 -1$}}$}\\
&=2y_{\frac n2}\tag{$\scriptstyle{\text{if $y_n$ is real, then $y_{n/2}$ is real}}$}
\end{align}
Notation 2:
\begin{align}
y_{n_2}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=1}^{\frac{N}{2}}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{1b}$}\\
&=\underbrace{\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+Y_N+\sum_{k=1}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}}_{\textrm{same as notation 1}}-Y_N-Y_{N/2}\\
&=2y_{n/2}-Y_N-Y_{N/2}\\
\end{align}

For Notation 3 however, the equation would be identical to equation $(1)$ if you changed the upper limit of the second summation to $N-1$ instead of $N$, otherwise you have an additional term as a result of the $N^{\rm th}$ term as shown below:
\begin{align}
y_{n_3}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{1c}$}\\
&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\left(\sum_{k=N/2}^{N-1}Y_k e^{j\frac{2\pi nk}{N}}\right)^*+Y_0\\
&=y_{n/2}+\left(y_{n/2}^{'}\right)^*+Y_0\tag{$\scriptstyle{\text{$y_{n/2}^{'}=y_n$ for $n=N/2, \ldots, N-1$}}$}\\
&=y_{n/2}+y_{n/2}^{'}+Y_0\\
&=y_n + Y_0
\end{align}

So, the correct notation should be:

$$
y_n=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\tag{$\mathbf{2}$}
$$

UPDATE:

Regardless of equation $(1)$, to test Notation 1 I will maintain my last notation and go proving if the equality below holds, i.e. if equation $(\mathbf{2})$ equals equation $(\mathbf{1a})$ in Notation 1.

\begin{align}
\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=\frac{N}{2}}^{N-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}&=\sum_{k=0}^{\frac{N}{2}-1}Y_k e^{j\frac{2\pi nk}{N}}+\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\\
\implies\sum_{k=\frac{N}{2}}^{N-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}&=\sum_{k=0}^{\frac{N}{2}-1}Y_{N-k} e^{j\frac{2\pi n(N-k)}{N}}\\
\end{align}

One way could be testing if the following equalities are satisfied:

\begin{align}
Y_{N/2}e^{j\frac{2\pi n(N/2)}{N}}&=Y_Ne^{j\frac{2\pi n(N)}{N}}\\
Y_{N/2-1}e^{j\frac{2\pi n(N/2-1)}{N}}&=Y_{N-1}e^{j\frac{2\pi n(N-1)}{N}}\\
Y_{N/2-2}e^{j\frac{2\pi n(N/2-2)}{N}}&=Y_{N-2}e^{j\frac{2\pi n(N-2)}{N}}\\
\vdots&=\vdots\\
Y_2e^{j\frac{4\pi n}{N}}&=Y_{N/2+2}e^{j\frac{2\pi n(N/2+2)}{N}}\\
Y_1e^{j\frac{2\pi n}{N}}&=Y_{N/2+1}e^{j\frac{2\pi n(N/2+1)}{N}}\\
\end{align}

It's neither the periodicity nor the Hermitian property of the DFT (for a real-valued signal i.e $Y_{N-k}=Y_{-k}=Y_k^*$) that we're seeing in the above equations. But let's check the first two and the last one for instance:

For the first equation:
\begin{align}
&Y_{N/2}e^{j\frac{2\pi n(N/2)}{N}}=e^{j\frac{2\pi n(N/2)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N/2)n}N}\right)=- \sum_{n=0}^{N-1}y[n]e^{-j\pi n}=\sum_{n=0}^{N-1}y[n]=Y_0\\
&Y_{N}e^{j\frac{2\pi n(N)}{N}}=e^{j\frac{2\pi n(N)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N)n}N}\right)=\sum_{n=0}^{N-1}y[n]e^{-j2\pi n}=\sum_{n=0}^{N-1}y[n]=Y_0\\
&\implies Y_{N/2}e^{j\frac{2\pi n(N/2)}{N}}= Y_{N}e^{j\frac{2\pi n(N)}{N}}\\
\end{align}
For the second equation:
\begin{align}
&Y_{N/2-1}e^{j\frac{2\pi n(N/2-1)}{N}}=e^{j\frac{2\pi n(N/2-1)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N/2-1)n}N}\right)= e^{-j\frac{2\pi n}{N}}\sum_{n=0}^{N-1}y[n]e^{j\frac{2\pi n}N}=e^{-j\frac{2\pi n}{N}}Y_1^*\\
&Y_{N-1}e^{j\frac{2\pi n(N-1)}{N}}=e^{j\frac{2\pi n(N-1)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N-1)n}N}\right)=e^{-j\frac{2\pi n}{N}}  \sum_{n=0}^{N-1}y[n]e^{\frac{j2\pi n}N}=e^{-j\frac{2\pi n}{N}}Y_1^*\\
&\implies Y_{N/2-1}e^{j\frac{2\pi n(N/2-1)}{N}}=Y_{N-1}e^{j\frac{2\pi n(N-1)}{N}}\\
\end{align}
For the last equation:

\begin{align}
&Y_1e^{j\frac{2\pi n}{N}}=e^{j\frac{2\pi n}{N}}Y_1\\
&Y_{N/2+1}e^{j\frac{2\pi n(N/2+1)}{N}}=e^{j\frac{2\pi n(N/2+1)}{N}}\left(\sum_{n=0}^{N-1}y[n]e^{-j\frac{2\pi(N/2+1)n}N}\right)=e^{j\frac{2\pi n}{N}}  \sum_{n=0}^{N-1}y[n]e^{\frac{-j2\pi n}N}=e^{j\frac{2\pi n}{N}}Y_1\\
&\implies Y_1e^{j\frac{2\pi n}{N}}=Y_{N/2+1}e^{j\frac{2\pi n(N/2+1)}{N}}
\end{align}

So, it looks like the initial notation in equation $(\mathbf{1a})$ is indeed identical to equation $(\mathbf{2})$. In conclusion, Notation 1 is correct for a real-valued signal. And Notation 3 need a slight change in the summation limit of the second term to be true. That's my two cents.


How is symbol synchronization with OFDM done? I'm trying to understand how symbol synchronization is done in OFDM using pilot tones, cyclic prefixes, or any other technique.

I've read the following answers which provide some explanation, but I still don't totally understand it.

How to demodulate an OFDM signal

How to estimate the number of taps needed for subsequent channel estimation algorithms?

Specific questions:

1) How is a pilot tone found?  What makes it different than the regular data on a sub-carrier?  How can it be used to determine symbol starts and ends?

2) If I understand the answers above correctly, a cyclic prefix can be used to find the symbol start/end because it will auto-correlate with some delay.  However, the cyclic-prefix exists in order to "absorb" ISI.  So if the prefix has been munged with ISI, then how can this auto-correlation be successful?

Regarding you're general question about how symbol sychronization is done in OFDM systems:

One of the most popular and frequently used techniques is the transmission of one or several pilot symbols that are known in the receiver. A pilot symbol is a complete OFDM symbol where the value of each subcarrier is predefined and known in transmitter and receiver. It is repeated with a certain rate that depends on how fast the channel changes. The received signal is correlated with the pilot symbol to detect the OFDM symbol start. It can also be used for channel estimation. Schmidl and Cox have introduced in [1] a pilot symbol based technique where the pilot symbol has a special symmetry so that the pilot symbol need not to be known at the receiver.
As Jason R has noted in his comment, although it is not its initial purpose, the cyclic prefix can also be used for symbol synchronization because its a known repetition of some part of the received signal that can be detected through autocorrelation. It is especially well-suited for fast-changing channels, because the delay time can be updated on a per-symbol basis. Additionaly, it does not add any additional overhead. However, it is more sensitive to noise [2] and presumably also to ISI.
Edit: The maximum delay that can be detected by this method is the lenght of one OFDM symbol. It's therefore only suited for fine synchronization.
There are some more "exotic" techniques. In one of these, for instance, the N-DFT (N = number of subcarriers) of time-shifted versions of the received signal is calculated. If you apply the DFT to the wrong time window, the resulting constellation diagram will be a mess. If you got the correct time window, the constallation digaram shows distinct constellation points. This can be detected by calculating the standard deviation of the DFT output. This method implies a high computational cost.

Regarding you specific questions

  How is a pilot tone found? What makes it different than the regular data on a sub-carrier? How can it be used to determine symbol starts and ends?

Once you have synchronized the received signal the pilot tones are at predefined bins of the DFT. When designing the system the location of pilot tones in the spectrum is fixed. There are more complex schemes, where the location of the pilot tones changes in a predefined pattern to get a good approximation of the channel in both frequency and time domain. Pilot tones cannot be used for synchronization, because the received signal has first to be synchronized before you can even extract the pilot tones in frequency domain. Assume that a wrong time window is used: ortogonality of subcarriers will be lost and the result of DFT is some mixture of two consecutive OFDM symbols. This is a nonlinear effect and the pilot symbols cannot be extracted from this mixture. Pilot tones are used for channel estimation and sometimes phase noise mitigation.
Edit: As Jim Clay has pointed out in his comments, fine synchronization through pilot tones is possible if a coarse value for the delay is known and the residual delay does not exceed the length of the cyclic prefix.

  If I understand the answers above correctly, a cyclic prefix can be used to find the symbol start/end because it will auto-correlate with some delay. However, the cyclic-prefix exists in order to "absorb" ISI. So if the prefix has been munged with ISI, then how can this auto-correlation be successful?

Like all synchronisation techniques this method will suffer from noise and channel dispersion and consequently will only work to some extent of the beforementioned effects. Quantifying to which extent exactly it is still working would require some thorough research that somone has certainly already done.

[1] Schmidl, T.M.; Cox, D.C.; , "Robust frequency and timing synchronization for OFDM," Communications, IEEE Transactions on , vol.45, no.12, pp.1613-1621, Dec 1997

[2] van de Beek, J.J.; Sandell, M.; Borjesson, P.O.; , "ML estimation of time and frequency offset in OFDM systems," Signal Processing, IEEE Transactions on , vol.45, no.7, pp.1800-1805, Jul 1997


Exponential average with time constant of slow, fast and impulse I have been studying about exponential average.  There are enough explanations about this at Internet, but they do not explain about the time constant.

I have one channel with a $T$ seconds time signal with sampling frequency $f_s$.

If I want to do average of this time signal, we need to use either linear or exponential method.  Linear average method is quite simple so there is no trouble for applying.  However, if I try to apply exponential average method, there are some problems.

If time signal varies quickly, we prefer to use fast time constant 125 ms.  Also, the time signal varies slowly, using 1000 ms of slow time constant is better, but in this situation, I do not know how can I apply this time constant with time signal.

Is there any explanation or any example for doing exponential average with time constant?

Exponential averaging is performed using a first order IIR filter according to the difference equation:
\begin{equation}
y[n] = (1-\alpha)\;x[n]+\alpha \; y[n-1]
\end{equation}
where $x[n]$ and $y[n]$ are the input and output signals at time $n$ and the value of the single coefficient $\alpha$ is given by the equation:
\begin{equation}
\alpha = \exp( \frac{-1}{f_s \tau} )
\end{equation}
where $f_s$ is the sampling frequency and $\tau$ is the time constant.

In this light, it can be seen that exponential averaging filters are simply low-pass filters, the roll-off of which is -6dB/oct.  The only free parameter is the cut-off frequency, which is specified by the time constant.  Longer time constants will yield a lower cut-off frequency, as can be see in the Matlab code and plot below.

figure; hold all
fs = 8000;
tau = 10.^(-1:-1:-4);
alpha = exp(-1/fs./tau);
color = 'brgc';
for index = 1 : length( alpha )
    [Hz fVec] = freqz( [1-alpha(index)], [1 -alpha(index)], 8192, fs );
    plot( fVec, 20*log10(abs(Hz)), color(index), 'LineWidth', 2 );
    legText{index} = sprintf( '$\\tau = 10^{%d}$', log10(tau(index)) );
end
hLeg = legend(legText,3);
set( hLeg, 'interpreter', 'latex' );
axis( [8000/8192 4000 -70 10] )
ylabel('Amplitude [dB]')
xlabel('Frequency [Hz]')
set( gca, 'xscale', 'log' )

Notice the parallel -6dB/oct slopes (characteristic of single pole filters) and the falling cut-off frequency with lengthening time constants.

The correct choice of time constant depends on how fast the variation of the desired average is.  In order to capture a quickly varying mean, only signal content above the fastest frequency of variation should be attenuated -- the time constant should be chosen accordingly.  This means that if the signal mean varies very slowly, long time constants are needed, and vice versa.  I find that trial and error works best to determine the optimal value, but formulas to calculate the cut-off frequency do exist, e.g. see this link.

Good luck!


Frequency dependent frequency shift I have the Fourier Transform or FFT of a signal, call it $F(f)$. I would like to apply a "frequency-dependent frequency shift" to it. Specifically, frequency $f$ should be shifted to $3f^2$, for all $f$. Physically, this signal is undergoing dispersion.

Say I have a signal of a one-dimensional wave (sound amplitude, water height, or voltage, etc.) as a function of time at a specific location, and I'm trying to estimate the spatial signal at a single instant in time. Let the dispersion relation be $k=3f^2$, where k is the wavenumber. This procedure will change the temporal signal to the spatial signal

This seems like a linear system to me, so how would I find the frequency response $H(f)$ of the system? Knowing the frequency response, I can simply multiply $F(f)$ by the $H(f)$ to find the resulting signal.

Before I try to answer the question, let's see how the question relates to dispersion.

The solution space of a dispersive (linear) wave equation in $n$ dimensions is spanned by the basis functions

$$f_{\mathbf{k}}(\mathbf{x},t) = \exp( i\, \omega(\mathbf{k}) t )\exp( i\, \mathbf{k}\cdot \mathbf{x} )$$

where $\mathbf{k} \in \mathbb{R}^n$ is the spatial frequency vector or wave vector. The temporal frequency $\omega$ depends on the  wave vector, or just on its magnitude $|\mathbf{k}|$ if the dispersion is isotropic. If $\omega \propto |\mathbf{k}|$ we have the special case of no dispersion. We can also assume in general that $\omega\geq 0$ to simplify calculations.

The basis functions are orthonormal in the Dirac sense with respect to their spatial coordinate:

$$\langle f_{\mathbf{k}},f_{\mathbf{k'}}\rangle := \int_{\mathbb{R}^n} f_{\mathbf{k}}(\mathbf{x},t)\, f_{\mathbf{k'}}(\mathbf{x},t)\, d\mathbf{x} = \delta(\mathbf{k}-\mathbf{k'})$$

And the general solution of the wave equation is:

$$f(\mathbf{x},t) = \int_{\mathbb{R}^n} c(\mathbf{k})\, f_{\mathbf{k}}(\mathbf{x},t)\, d\mathbf{k} $$

If we now observe that general solution at $\mathbf{x}=0$, we get the relationship

$$f(0,t) = \Re\left(\int_{\mathbb{R}^n} c(\mathbf{k})\,\exp( i\, \omega(\mathbf{k}) t )\, d\mathbf{k}^n \right)$$

which looks similar to a Fourier transform of $c(\mathbf{k})$. If we can turn this into a proper Fourier transform, the inverse Fourier transform of the observation will give us the expansion coefficients $c(\mathbf{k})$ required for the reconstruction of the full spatial solution of the wave.

At this point, we need to make a few additional assumptions. The first one is that $n$ equals 1 so that the coefficients are also a function on $\mathbb{R}$. We also have to assume that the dispersion is isotropic and $\omega(|\mathbf{k}|)$ is differentiable and strictly monotonic.

This is not enough however. We also need to assume a direction of travel so that $\omega(|\mathbf{k}|)$ becomes an invertible function on a restriction of $k>0$ or $k<0$ respectively.

We pick a "right" travelling wave, so $k>0$. With this, we can change the integration variable with $dk=\frac{\partial k}{\partial\omega}\,d\omega$ to get

$$f(0,t) = \Re\left(\int_\mathbb{R_{\geq 0}} c(k(\omega)) \frac{\partial k(\omega)}{\partial\omega} \exp(i\,\omega t)\, d\omega \right)$$

where $k(\omega)$ is the inverse function of $\omega(k)$. Identifying

$$C(\omega) := c(k(\omega)) \frac{\partial k(\omega)}{\partial\omega}$$ for $\omega\geq 0$ and the complex conjugate for $\omega<0$, we can see that $f(0,t)$ is the Fourier transform of $C(\omega)$.

With the inverse Fourier transform we get

$$C(\omega) = \int_\mathbb{R}f(0,t) \exp(-i\,\omega t) dt$$

and consequentially

$$c(k(\omega)) = \frac{1}{k'(\omega)} \int_\mathbb{R}f(0,t) \exp(-i\,\omega t) dt$$

or simply

$$c(k) = \omega'(k) \int_\mathbb{R}f(0,t) \exp(-i\,\omega(k) t) dt$$

which we can express in terms of the inverse Fourier transform as

$$c(k) = \omega'(k)\,\mathcal{F_t}^{-1} \left\{ f(0,t) \right\}\left( \omega(k) \right) $$

Using this theoretical background we can now reconstruct the spatial wave function of a free one-dimensional right travelling wave with open boundary conditions.

In the expansion of the general wave solution from above we can substitute the $c(k)$ to get:

$$f(x,t)=\int_\mathbb{R} \omega'(k)\,\mathcal{F_t}^{-1} \left\{ f(0,t) \right\}\left( \omega(k) \right)\, \exp( i\, \omega(k) t )\exp( i\, k x ) dk $$

This can also be written as a Fourier transform in $k$:

$$f(x,t)= \mathcal{F_k}\left\{ \mathcal{F_t}^{-1} \left\{ f(0,t) \right\}\left( \omega(k) \right)\,\omega'(k)\,\exp(i\,\omega(k)\,t)\, \right\}(x)$$

This translates to the following algorithm for the reconstruction of the spatial wave function at time $t$:

Take the inverse Fourier transform of $f(0,t)$
Evaluate the result at $\omega(k)$
Multiply the result with $\omega'(k)$
Evolve the result in time by multiplying with $\exp(i\,\omega(k)\,t)$
Take the forward Fourier transform with respect to k

You can implement this in the discrete domain by taking regular grids for the variables $\omega$ and $k$ and use the discrete Fourier transform.

As you can see, it's a little more complicated than what you described, but your idea was going in the right direction.


Kalman filter after lowpass filter: bad idea? I am new to Kalman filter and am enjoying playing with it. However, I generated some random velocities and acceleration and Kalman filter (with the covariance matrices I have chosen) is comparable with lowpass filtering.

So I thought maybe both lowpass filtering and Kalman filter could be combined, however this completely ruins the assumption of white Gaussian noise. My numerical experiments indicates that Kalman filter does not bring anything interesting after LP filtering for the cases I tried.

What are good preprocessing practices, if any, before applying a Kalman filter?

The Kalman filter is the optimal filter under various assumptions. You need to check whether those assumptions hold in your case:

a) the model perfectly matches the real system,
b) the entering noise is white and Gaussian and
c) the covariances of the noise are exactly known.

Without further detail I can't say whether your statement My experiments indicates that Kalman filter does not bring anything interesting after LP filtering. is true or not.
Use Occam's Razor: the simplest approach that works is the best. That is, it may well be that for your situation a low pass filter is all you need.


I and Q components and the difference between QPSK and 4QAM Both 4QAM and QPSK apparently produce the same waveform, but are they the same mathematically?

In a QPSK constellation, are the mapping points at 45, 135, 225 and 315 degrees while the 4QAM is at 0, 90, 180 and 270?

I also struggle to understand the I/Q components of such a constellation diagram. What does "inphase" and "quadrature-phase" actually mean? Are they just another way to specify the real and imaginary part for this type of usage?

Both QPSK and $4$-QAM constellations have signal points at $45, 135, 225$, and $315$ degrees (note typo in your question).  They arise from amplitude modulation (or, if
you prefer, phase modulation) of two carrier signals
(called the inphase and quadrature carriers) that are orthogonal (meaning that they
differ in phase by 90 degrees.  The canonical representation of a QPSK or $4$-QAM
signal during one symbol interval is
$$s(t) = (-1)^{b_I}\cos(2\pi f_c t) - (-1)^{b_Q}\sin(2\pi f_c t)$$
where $\cos(2\pi f_c t)$ and
$-\sin(2\pi f_c t)$ are the inphase and quadrature carrier signals at
frequency $f_c$ Hz
and $b_I, b_Q \in \{0,1\}$  are the two data bits (called the inphase
and quadrature data bits, naturally, since they are transmitted
on the inphase and quadrature carriers).  Notice that the
inphase carrier $\cos(2\pi f_c t)$ has amplitude $+1$ or
$-1$ according as the inphase data bit has value $0$ or $1$, and
similarly the quadrature carrier $-\sin(2\pi f_c t)$ has amplitude $+1$ or
$-1$ according as the quadrature data bit has value $0$ or $1$.
Some people regard this as an inversion of the normal scheme
of things, didactically asserting
that positive amplitudes must be associated with $1$ data bits and
negative amplitudes with $0$ bits. But if we look at it from the
phase modulation perspective, a $0$ bit means that the carrier
($\cos(2\pi f_c t)$ or $-\sin(2\pi f_c t)$ as the case may be) is
transmitted with no change in phase while a $1$ data bit creates
a change in phase (we will think of it as a phase delay)
of $180$ degrees or $\pi$ radians.  Indeed, another
way of expressing the QPSK/$4$-QAM signal is as
$$s(t) = \cos(2\pi f_c t - b_I\pi) - \sin(2\pi f_c t - b_Q\pi)$$
which makes the phase modulation viewpoint very clear. But,
regardless of which viewpoint we use, during a symbol interval,
the  QPSK/$4$-QAM signal is one of the following four signals:
$$\sqrt{2}\cos\left(2\pi f_c t + \frac{\pi}{4}\right),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{3\pi}{4}\right),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{5\pi}{4}\right),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{7\pi}{4}\right)$$
corresponding to $(b_I,b_Q) = (0,0), (1,0), (1,1), (0,1)$ respectively.

Note that the viewpoint taken here is of QPSK as consisting of
two BPSK signals on phase-orthogonal carriers. The demodulator
thus consists of two BPSK receivers (called the inphase branch
and quadrature branch, what else?).
An alternative view of QPSK as changing the phase
of a single carrier depending on a $4$-valued symbol is developed a little later.

The QPSK/$4$-QAM signal can also be expressed as
$$s(t) = \text{Re}\{B \exp(j2\pi f_c t)\}
= \text{Re}\{[(-1)^{b_I}+j(-1)^{b_Q}] \exp(j2\pi f_c t)\}$$
where $B$ is the complex-valued baseband symbol taking on values
in $\{\pm 1 \pm j\}$ and which, when plotted on the complex plane,
gives constellation points distant $\sqrt{2}$ from the
origin and at  $45, 135, 225$, and $315$ degrees
corresponding to data bits $(b_I,b_Q) = (0,0), (1,0), (1,1), (0,1)$
respectively. Note that complementary bit pairs lie diagonally across the circle from each other so that double bit errors are less likely than single bit errors. Note also that the bits naturally
occur around the circle in Gray code order; there is no need to massage a given data bit pair $(d_I,d_Q)$ (say $(0,1)$) from "natural representation" (where it means the integer $2 = d_I+2d_Q$: $d_I$ is the LSB and $d_Q$ the MSB here) to "Gray code representation" $(b_I,b_Q) = (1,1)$ of the integer $2$ as some implementations seem to insist
on doing. Indeed, such massaging leads to poorer BER performance since the decoded $(\hat{b}_I,\hat{b}_Q)$ must be ummassaged at the receiver into the decoded data bits $(\hat{d}_I,\hat{d}_Q)$ making the single channel  bit error $$(b_I,b_Q) = (1,1) \to (\hat{b}_I,\hat{b}_Q) = (1,0)$$
into the double data bit error
$${(d_I,d_Q) = (0,1) \to (b_I,b_Q) = (1,1)\to (\hat{b}_I,\hat{b}_Q) = (1,0)}\to (\hat{d}_I,\hat{d}_Q) = (1,0).$$

If we delay the four possible signals exhibited above by $45$ degrees or
$\pi/4$ radians (subtract $\pi/4$ radians from the argument of the cosinusoid), we get
$$\sqrt{2}\cos\left(2\pi f_c t + \frac{\pi}{4}\right)\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 0\frac{\pi}{2}\right) = \sqrt{2}\cos(2\pi f_c t),\\ \sqrt{2}\cos\left(2\pi f_c t + \frac{3\pi}{4}\right)\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 1\frac{\pi}{2}\right) = -\sqrt{2}\sin(2\pi f_c t),\\
\sqrt{2}\cos\left(2\pi f_c t + \frac{5\pi}{4}\right)
\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 2\frac{\pi}{2}\right) = -\sqrt{2}\cos(2\pi f_c t) \\
\sqrt{2}\cos\left(2\pi f_c t + \frac{7\pi}{4}\right)\Rightarrow \sqrt{2}\cos\left(2\pi f_c t + 3\frac{\pi}{2}\right) = \sqrt{2}\sin(2\pi f_c t),\\$$
which give the four constellation points at $0,90,180,270$ degrees referred
to by the OP.  This form gives us another way of viewing QPSK signaling:
a single carrier signal whose phase takes on four values depending on
the input symbol which takes on values $\{0,1,2,3\}$. We express this in tabular form.
$$
\begin{array}{|c|c|c|c|c|c|c|}
\hline
(b_I,b_Q) & \text{normal value} ~k & \text{Gray code value} ~\ell & \text{signal as above} &\text{phase-modulated signal}\\
\hline
(0,0) & 0 & 0 & \sqrt{2}\cos(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 0\frac{\pi}{2}\right)\\
(0,1) & 1 & 1 & \sqrt{2}\sin(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 1\frac{\pi}{2}\right)\\
(1,1) & 3 & 2 & -\sqrt{2}\cos(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 2\frac{\pi}{2}\right)\\
(1,0) & 2 & 3 & -\sqrt{2}\sin(2\pi f_c t)
& \sqrt{2}\cos\left(2\pi f_c t - 3\frac{\pi}{2}\right)\\
\hline
\end{array}
$$
That is, we can regard the QPSK modulator as having input
$(b_I,b_Q)$ that it regards as the Gray code representation
of the integer $\ell \in \{0,1,2,3\}$ and produces the
output
$$\sqrt{2}\cos\left(2\pi f_c t - \ell\frac{\pi}{2}\right).$$
In other words, the phase of carrier $\sqrt{2}\cos(2\pi f_c t)$ is
modulated (changed from $0$ to $\ell\frac{\pi}{2}$) in
response to the input $\ell$.

So how does this work in real life or MATLAB, whichever comes first?
If we define a QPSK signal as having value $\sqrt{2}\cos\left(2\pi f_c t - \ell\frac{\pi}{2}\right)$ where the value of $\ell$ is typed in as 0 or 1 or 2
or 3, we will get the QPSK signal described  above, but the
demodulator will produce the bit pair $(b_I, b_Q)$ and we must
remember that the output is $\ell$ in Gray code interpretation,
that is, the demodulator output will be $(1,1)$ if $\ell$ happened
to have value $2$, and interpreting output $(1,1)$ as $3$ is a decoding
error that is not generally discussed in textbooks!


signal processing vs image processing? Originally this question was asked EE SE yesterday but then i saw a comment about DSP SE ,so i am posting my question again here ,in hope to get much better and relevant answers

Link of EE question is
https://electronics.stackexchange.com/questions/497036/signal-processing-vs-image-processing

My question was/is

Although both are domains of EE, but what are differences between them? I am able to draw following conclusions from wikipedia article of dip and also from some other web reources

1) DIP(Digital image processing) is a subset of DSP(digital signal processing)

2)Because of above concept in (1), DIP only deals with 2 dimensional images/signals or in some case 3 dimensional images/signals(forexample RGB images) while in DSP there isn't any restriction on number of dimensions?

Please kindly guide me and correct me if i am wrong in my concepts/assumptions

"General" DSP (and signal processing in general) tends to deal much more with problems that can be solved with linear filters, or at least a preponderance of linear filtering, and there's a lot more content involving 1-D signals -- so you tend to see a lot of emphasis on frequency-domain analysis and FFT (and even some folks who can't seem to distinguish between typing 'fft' in Matlab and understanding DSP).

Image processing (video processing in particular) tends to deal with problems that must be solved with lots of nonlinear processing -- image segmentation, finding edges, etc.  Even super-duper neural networks have significant nonlinearities.  While there are linear operations (i.e., convolution, cross-correlation) these tend to be processor-intensive, and are used sparingly.

One can go through an entire career in DSP and not do image processing, and visa-versa (I'm a case in point, with around 30 years experience in DSP and digital control systems, yet with very little image processing under my belt).

About the only thing that is really common, besides the fact that you're taking in something you could call a "signal" and you're doing something to it you could call "processing", is that "plain old" DSP and image processing both reward the practitioner who's willing to do the hard math (or who's lucky enough to find the hard math easy).

It also helps to be able to make intuitive leaps to potential solutions, and then do the hard math to back-fill to determine that your solution will actually work, to refine it, and to be on a firm mathematical foundation for your next intuitive leap.


Log vs. linear frequency scales of Fourier and wavelet transforms I'm trying to understand the difference between the output of a Fourier transform and a wavelet transform. A Fourier transform is done via the following function:

$$\hat{f}(\xi) = \int^\infty_{-\infty}\ f(t)\ e^{-2\pi i t \xi}\ dt$$

Whereas a wavelet transform is going to use this:

$$F(a,b) = \int^\infty_{-\infty}\  f(x)\psi^*_{a,b}(x)\ dx$$

Now, spectrograms use a windowed Fourier transform. When looking at a signal over time, you can get a graph showing you which frequencies were present at any given time. The frequencies are scaled linearly. However, scalograms use a wavelet transform to obtain the same information. I've frequently seen the output scaled logarithmically (in frequency).

Is there something about wavelets that make their output fundamentally logarithmic? I'm having a hard time seeing it in the above formulas. It seems like you set $\xi$ or $a$ & $b$ give you the frequency you want, and you could calculate results for any frequency you desire.

And FFT can be considered a filter bank. The bandwidth of each filter is inversely proportional to the length of the FFT. The width of each filter sets the spacing such that there isn't either too much overlap or gaps between filters.  In an STFT spectrogram, the window width is fixed for the entire FFT, thus the filter spacing is fixed, which turns out to be a linear spacing.

With wavelets the length of each wavelet, being individually adjustable, can be shorter for the higher frequencies, thus have a wider bandwidth, which allows them to be spaced further apart in frequency without leaving gaps between them in the frequency response, thus allowing logarithmic spacing.


High Order Buttersworth Filter Let's suppose i want to implement a Buttersworth filter of order 12 (N=12).
First of all ,I can decompose the Buttersworth polynomial in elemental polynomials:

Now, if i know how to physically implement the elemental polynomials, how can i construct the original Buttersworth filter of order 12?

The impulse response of two filters in series is the convolution of each filters impulse response. You can multiply two polynomials by convolving their coefficients. Likewise you can factor or combine filters by factoring or multiplying their transfer functions described as polynomials.  For example, a fourth order filter can be factored into two 2nd order filters in series, or two 2nd order filters can be combined into one fourth order filter- it’s simply polynomial multiplication or factoring.

So in your case since you already factored the Butterworth filter into individual filters that you know how to implement, you would then simply cascade those filters.

As @MattL points out in the comments, you must pay careful attention to impedances when cascading analog filter sections. In the analog such structures are quite common with active filter design where it is quite simple to isolate the sections with appropriate design techniques, but not so common with filters constructed with passive components.


How do software equalizers work? Recently, I have been studying IIR and FIR filters, and trying to create a parametric equalizer using a microcontroller. More specifically, I use the ADC on it to sample the audio at about 44 kHz, and in between grabbing samples, I use the time to process the signal and pass it out on the DAC.
What is a bit confusing to me though is that this process seems to be pretty intense, and often takes most of the time in between samples to do all the processing (with 8 2-pole IIR bandpass filters).
Programs like Spotify, or devices like most Galaxy phones have similar equalizers but they are running on non-dedicated hardware, and I am really confused how they are able to pull this off? Does anybody know if they are basically doing the same thing just faster, or is there some secret sauce to this that I am missing?
I have never tested either of the previously mentioned equalizers in those platforms, so maybe they are just really crappy, but at least just by listening to them, it doesn’t sound horrible.

This depends a lot on how you implement it.

A single biquad takes about 10 arithmetic operations. (To be precise
a Transposed Form II takes 4-5 multiplies and 3 adds, depending on how the gain management is done).
Arithmetic operation translates into clock cycles of your processor. That depends a lot on the efficiency of your instruction set and how good yo are at writing code. On a good DSP you can do a single  biquad in 4 clock cycles (ca. 2 ops per cycle) and a general purpose ARM with a decent ALU this may by 15-20 cycles. On a crappy micro controller with no native hardware multiplier this can take up to 50-100 cycles.
Frame based vs sample based processing. Reading data from an ADC, moving it into the processor, and the moving it to the DAC creates a lot of overhead. That's why most implementations will process a frame of data and not one sample at a time. That "amortizes" the overhead across the entire frame. A good processor will also allow to do this with HW support, i.e. it will use DMA channels to move the data from/to DAC/ADC to internal memory without the CPU having to do anything.
The chips inside modern smart phones (even the cheap ones) are quite powerful: they are typically quad core ARMs with some HW accelerators (Neon, SIMD) for media processing. They all clock in the multiple GHz range and have good HW support for peripherals and graphic accelerators. My trusty old Samsung S8 has a 8 core 2 GHz Snapdragon: https://www.zdnet.com/article/qualcomm-snapdragon-835-what-does-the-kryo-280-adreno-540-spectra-180-x16-and-hexagon-682-mean-for/#:~:text=The%20Snapdragon%20835%20is%20comprised,result%20in%20optimized%20power%20savings.

Let's a run a simple example: A well coded biquad on a decent ARM core should take about 15 cycles. For 8 biquads running at 44.1 kHz on a stereo signal (two channels) that would be just about 10M cycles/second. That's about 1% of a single core at 1 GHz or 0.25% of the entire chip. On my Snapdragon 835 it would be less than 0.03%


Efficient implementation of FIR comb/decimating filter I'm working on trying to find an efficient way to implement a FIR comb/decimating filter. For a given "prototype" low-pass filter $
H(z)$, a comb that "repeats" $L$ times in the frequency domain can be constructed by upsampling (stuffing $L-1$ zeros in between each sample), yielding $H(z^L)$ (Proakis and Manolakis 342). We will denote the impulse response of the comb filter by $c[k]$.

I also wish to decimate the comb filter output by a factor of $D$. I do not wish to cascade an integrator to make a CIC filter, and I am aware that aliasing will be introduced. I am okay with this. To implement the decimating FIR filter, I have chosen to use a polyphase implementation. This means that I will have $D$ branches, each of which will contain a downsampled and delayed version of the base FIR filter (Proakis and Manolakis 767), i.e.

\begin{equation}
    p_i[n] = c[nD + i]
\end{equation}

However, since $c[k]$ is $L$-sparse, i.e. only one out of every $L$ samples is nonzero, each polyphase branch, even though they are all downsampled, will still have a lot of zeros in them. In general, $L > D$ for my application, so effectively each branch filter $p_i[k]$ is a resampled version of $h[k]$ by a factor of $L/D$ (with no interpolation/decimation filter in between the upsampling and downsampling steps).

Is there a computationally efficient way to take advantage of the fact that most of the taps are zero? Does multiplying by zero in floating point arithmetic take less computational effort or the same amount as multiplying two nonzero numbers? If so, then perhaps a naive implementation would be sufficient.

Works Cited:

Proakis, John and Dimitris Manolakis. Digital Signal Processing: Principles, Algorithms, and Applications. 4th ed. Upper Saddle River, NJ: Pearson 2006. Print.

"Is there a computationally efficient way to take advantage of the fact that most of the taps are zero?"

is there a pattern to those taps (or to the taps that are potentially non-zero)?  can you skip over those taps that are known in advance to have zero coefficient?


Meaning of AWGN complex and real noise under BPSK Decision Scheme I'm currently studying error correcting codes and decided to replicate some common Bit Error Rate curves under BPSK modulation.

I initially had trouble to convert SNR to a Gaussian $\sigma$ because I was confused why people used diffent $\sigma$ if the channel was real or complex. I initially thought that I could add $\sigma$ for both IQ components. Then I understood from How do I add AWGN to an I and Q representation of a signal? that $\sigma_I = \sigma_Q = \frac{\sigma}{\sqrt(2)}$ for total noise energy consistency.

But after this information regarding how I should add noise to the received signal I don't get how this does not generate a inconsistency under the BPSK decision boundary if I interpret the channel as complex or real.
From what I understand from a given SNR, there is a relationship with the total Gaussian noise $\sigma$ that I should add.

If channel is real (only considering I component), BPSK maps bits to {-1, 1}, noise added in I component is $\sigma_0$ for SNR = 0dB.
If channel is complex (considering IQ components), BPSK maps bits to {(-1, 0), (1, 0)}, noise added for SNR = 0dB has total $\sigma = \sigma_0$, so $\sigma_I = \sigma_Q = \frac{\sigma_0}{\sqrt(2)}$.

Since BPSK should only be concerned with the real values (InPhase values) the Bit Error Rates would be clearly different in both cases.

I guess my lack of understanding must be related to the SNR to $\sigma$ mapping under a real/complex channel. I don't understand how interpreting the problem in 2 dimensions differs from 1 dimension and its relationship with $\sigma$.

It doesn't matter if the channel is complex or real. It's actually always complex, but since you are using BPSK, only the real part matters. How to add it is a mathematical trick. Say you need to generate an SNR of x in linear scale, and your signals have unit energy, then the noise variance must be 1/x. In other words, you don't generate real and complex parts for the noise, because then the SNR will be $1/(2x)$, which is not correct. You generate the real and complex parts only if the signal is complex, at which case the total noise variance gives the desired SNR. Hope this helps.


How to do prediction using frequency domain data? Both linear regression and Kalman filtering can be used to estimate and then predict from a time domain sequence of data (given some assumptions about the model behind the data).

What methods, if any, might be applicable to do prediction using frequency domain data?  (e.g. predict a future step, using the output from suitable FFT(s) of previous data, without just going back to the time domain for the estimation.)

What assumptions about the data, or the model behind the data, might be required for what, if any, quality or optimality of prediction in the frequency domain?  (But assume it is not apriori known whether the data source is strictly periodic in the FFT aperture width.)

An important NOTE: Since you are talking about frequency
  domain, it is implied that entire DFT spectrum is available and hence
  estimation is used for smoothing rather than future prediction.

If the signal is stationary, you can apply wiener filter and the model produced is an FIR filter; in this case, the signal estimation in the time domain will be identical to that of frequency domain.

From wiki:
Wiener's main accomplishment was solving the case where the causality requirement is in effect, and in an appendix of Wiener's book Levinson gave the FIR solution.

Removing noise using wiener filter using deconvolution is called Wiener deconvolution. This works in frequency domain. And is quite well used in Image deconvolution.

I don't know if there is a formulation possible for Kalman filter to be used for given Frequency domain data (assuming DFT) because usual implementations are actually iterative sample by sample.  But kalman smoothing approaches probably can do similar thing.


Fast Recursive 1D Signal Smoothing - IIR / Auto Regressive Implementation of Gaussian Smoothing I have just begun to dive into the field of signal processing, but there is the need to program a digital filter, that has to smooth a realtime signal from a sensor device. As far as I know, in my case a FIR filter is not suitable, because I have a lot of samples (1000) per sec and the convolution that is needed for appling a gaussian smoothing takes to long. Therefore I need some IIR Filter, that is capable of deliver results, that are comparable with gaussian filtering. I have already tested the first order lag filter, but this filter does not smooth out the signal if it was constant over time before.

Do you have any suggestions where to start? What are suitable fast recursive smoothing filters that look like gaussian smoothing? Some hints will be usefull for my to get the keywords for further reading.

Have a look at my Fast Gaussian Blur Project at GitHub.

You will find there implementation of IIR Approximation of Gaussian Blur which implements the following papaers:

Recursive Gabor Filtering.
Recursive Implementation of the Gaussian Filter.
Boundary Conditions for Young - van Vliet Recursive Filtering.

The idea is pretty straight forward.


Cut-off frequencies for fractional sample rate adjustment We have a signal sampled att 22 kHz that we want to interpolate to 40 kHz.

So we can do this by upsampling by a factor 20, then downsampling by a factor 11. My question regards the choice of cut-off frequencies after the upsampling and the downsampling.

The old Nyquist frequency was at $11$ kHz so it would make sense to me to low-pass frequencies above $11$ kHz after upsampling - we're basically replicating our signal 20 times across the new spectra, right? The correct cut-off is apparently 20 kHz, but I can't quite wrap my head around that fact.

We then downsample by a factor 11 to 40 kHz. Now the correct LP cut-off is 11 kHz. I am seriously confused. Can anyone shed some wisdom here?

Don't be confused; you're doing everything correctly here:

$$\begin{align*}
f_{sample,in} &&\overset{\text{interpolate}}{\rightarrow}&& f_{sample,intermediate} &&\overset{\text{decimate}}{\rightarrow}&& f_{sample,out}\\
22 \mathrm{kHz}&&\overset{\uparrow 20}{\rightarrow}&& 440 \mathrm{kHz} &&\overset{\downarrow 11}{\rightarrow}&& f_{sample,out}
\end{align*}$$

To avoid images after interpolating, you will have to suppress all signal above 11kHz, which demands low pass with a cut off of $\frac1{2\cdot20}f_{sample,intermediate}$.
To avoid aliases during decimation, you will have to suppress all signal above 20kHz, which demands a low pass filter with a cut off frequency of $\frac1{2\cdot11}f_{sample,intermediate}$.

Since the filter from 1. implies 2., using one filter is sufficient.

This is a very commonly applied optimization.


What is the role of a LPF in oversampling? Lets assume that a signal $x(n)$ is up-sampled by adding 1 zero between two adjacent samples to form a signal $y(n)$. How does a digital LPF give an oversampled version $z(n)$, with more data points? (As show in the image below) I know that the question might seem stupid but the books that I am using don't really explain that part.

@user29568: Your question is very far from stupid. It's a darned good question. To see why time-domain "zero stuffing" generates the additional spectral images, that must be filtered out as described by Fat32, you're welcome to see my explanation of this effect at: http://www.dsprelated.com/showarticle/761.php


Isolate g(t) starting from $\int_{-\infty}^{+\infty}g^{2}(t)e^{-i2\pi f t}dt$ Is it possible to isolate $g(t)$ starting from
$G^2(f)=\int_{-\infty}^{+\infty}g^{2}(t)e^{-i2\pi f t}dt$ ? Can you show me how?

The integral is the Fourier transform of $g^2(t)$, so an inverse Fourier transform can recover $g^2(t)$. But without any further information you can't recover $g(t)$ from $g^2(t)$. So the problem is not the integral transform but simply the squaring operation.


Periodicity of a constant signal! This can be a very silly question, but I'm quite confused:

If we take the Fourier transform of any constant signal, we get an impulse at zero, which says that its frequency is zero and, hence, it is non-repeating and its period is infinity.
By the definition of a periodic signal, $$F(z+p) = F(z)$$ therefore, for a constant function $$F(z) = c$$ $$F(z+p) = c, \forall p >0$$ Therefore, a constant function is periodic, but its period is undefined or can be defined as anything.

Which of the above arguments are correct and why? Please clarify. Thanks in advance.

As you say, the constant function is periodic. A signal $x(t)$ is said to be periodic with period $p$ or to have a period $p$ if there exists a $p>0$ such that $x(t+p)=x(t)$ for all real numbers $t$. Note that I said "a period" instead of "the period" -- periodic signals have an infinite number of periods, maybe even an uncountable infinite of them!

Often, what is referred to as the period of a signal is actually more formally described as the "fundamental period", which is the smallest among all periods.

So, a constant signal is periodic, it has an uncountably infinite number of periods (since any real number $p>0$ is a period), but it does not have a fundamental period.


Codec Output has "Wavy" Amplitude on Chirp Input I am sending a chirp that sweeps form 20Hz to 500Hz to an audio codec, CS42L52, connected to an MCU, STM32F405, and I am receiving a "wavy" amplitude output.  There is no processing on the MCU, I pass the signal straight from the DMA in to the DMA out.  The DMA services the I2S lines.  If I change the time that the chirp takes to get from 20 to 500 Hz, the amplitude wave changes its frequency.

This picture shows the chirp with a 10 second duration on top and 5 second duration on bottom, and shows that the amplitude frequency is twice as fast in the 5 second chirp vs the 10 second chirp.

(!

Looking at the codec block diagram, it shows that I can bypass the MCU in two different spots. The first being right after the PGA inputs and the next after gain adjust and ALC.

Both of these pass to the codec's DSP and then out to the codec's DAC according to this block diagram.

When I pass the signal through either of these lines, skipping the MCU, the output has a constant amplitude as it should.  The input signal is clean, no "wavy" amplitude. If I put one frequency through the device, no chirp, the amplitude stays constant.  I have played with a lot of the settings on the codec and can't seem to find any that affect the amplitude wave.  Maybe I am overlooking a setting?  This almost makes me think it is an MCU problem, but I am not manipulating the signal at all on the MCU.  Is it a problem that the DMA stores it as an unsigned integer?  What could be the problem?

My sampling frequency is set at 48 kHz, but I measure it at 48.5 kHz, could this affect it?

These could be ripples in the frequency response of an (ill-configured) anti-aliasing filter involved in the acquisition or restitution of the signal. It could be that the code running on your MCU and initializing the codec sends invalid configuration information to the codec; causing the codec to work with filters coefficients defined for a wrong sample rate. One thing worth noting is that the STM32F4 I2S master can only output an MCK at 256x FS. Check that your codec can work at 48kHz with such an MCK frequency.

Another thing I can think of is your software "pass-through" routine on the MCU dropping samples. Depending on the ratio between the block size and the period of the waveform, an amplitude modulation effect could appear.


How to choose a microphone and make sure not to overcharge (clip) it? I'm not sure that's the right forum to ask this question but I'll do it anyway.

In my project, I've to listen to different machines using a MEMS mic by fixing the microphone directly on to the machine. The problem that I get, is  that all of the records that I need, are overcharged (clipped) so useless. My question is how should I choose the microphone for really loud measurements?

Thanks in advance for any hint.

Choose a commercial Sound Pressure Level meter to roughly measure the noise level close to the machine
Look at the data sheet of the MEMS and see if the specification matches the expected level with enough head room (at least 10 dB). IF not, chose a MEMS that meets this spec
Maker sure the MEMS is mechanically de-coupled from the machine, i.e. picks up only sound and NOT vibration.


Understanding SURF Features Calculation Process So, I was reading the paper on SURF (Bay, Ess, Tuytelaars, Van Gool: Speeded-Up Robust Features (SURF)) and I can not comprehend this paragraph below:

  Due to the use of box filters and integral images, we do not have to
  iteratively apply the same filter to the output of a previously filtered
  layer, but instead can apply box filters of any size at exactly the
  same speed directly on the original image and even in parallel
  (although the latter is not exploited here). Therefore, the scale
  space is analysed by up-scaling the filter size rather than iteratively
  reducing the image size, figure 4.

This is figure 4 in question.

PS: The paper has an explanation of integral image, however the whole content of the paper is based on the particular paragraph above. If anybody has read this paper, can you briefly mention what is going on here. The whole mathematical explanation is quite intricate to have a good grasp first up, so I need some assistance. Thanks.

Edit,couple of issues:

1.

  Each octave is subdivided into a constant number of scale levels.
  Due to the discrete nature of integral images, the minimum scale
  difference between 2 subsequent scales depends on the length lo of the
  positive or negative lobes of the partial second order derivative in
  the direction of derivation (x or y), which is set to a third of the
  filter size length. For the 9x9 filter, this length lo is 3. For two
  successive levels, we must increase this size by a minimum of 2 pixels
  (one pixel on every side) in order to keep the size uneven and thus
  ensure the presence of the central pixel. This results in a total
  increase of the mask size by 6 pixels (see figure 5).

Figure 5

I could not make sense of the lines in the given context.

  For two successive levels, we must increase this size by a minimum of
  2 pixels (one pixel on every side) in order to keep the size uneven
  and thus ensure the presence of the central pixel.

I know they are trying to do something with the length of the image, if its even they are trying to make it odd, so that there is a central pixel which will enable them to calculate the maximum or the minimum of the pixel gradient. I am bit iffy about its contextual meaning.

2.

In order to calculate descriptor Haar wavelet is used.

How is the middle region has low $\sum\ dx$ but high $\sum\ |dx|$.

3.

What is the necessity of having an approximate filter?

4.
I have no issue with the way that they found out the size of the filter. They "did" something empirically. However, I have some nagging issue with this piece of line

  The output of the 9x9 filter, introduced in the previous section, is
  considered as the initial scale layer, to which we will refer as scale
  s = 1.2 (approximating Gaussian derivatives with  σ= 1.2).

How did they found out about the value of σ . Moreover how does the calculation of scaling done shown in the image below.The reason I am stating about this image is that the value of s=1.2 keeps on recurring, without clearly stating about its origin.

5.
The Hessian Matrix represented in terms of L which is the convolution of second order gradient of Gausssian filter and the image.

However the "approximated" determinant is said to contain only terms involving second order Gaussian filter.

The value of w is:

My question why the determinant is calculated like that above, and what is the relationship between approximate Hessian and Hessian matrix.

What's SURF?

In order to correctly understand what is going on, you also need to be familiar with SIFT: SURF is basically an approximation of SIFT. Now, the real question becomes: what's SIFT?.

SIFT is both a keypoint detector and a keypoint descriptor.
In the detector part, SIFT is essentially a multi-scale variant of classical corner detectors such as the Harris corner, and that has the ability to auto-tune the scale. Then, given a location and a patch size (derived from the scale), it can compute the descriptor part.

SIFT is very good at matching locally affine pieces of images, but it has one drawback: it is expensive (i.e., long) to compute.
A large amount of time is spent in computing the Gaussian scale-space (in the detector part), then in computing histograms of the gradient direction (for the descriptor part).

Both SIFT and SURF can be seen as difference of Gaussians with automatic scale (i.e., Gaussian sizes) selection. This, you construct first a scale-space where the input image is filtered at different scales. The scale-space can be seen as a pyramid, where two consecutive images are related by a scale change (i.e., the size of the Gaussian low-pass fiéter has changed), and scales are then grouped by octaves (i.e., a big change in the size of the Gaussian filter).

In SIFT, this is done by repeatedly filtering the input with a Gaussian of fixed width until the scale of the next octave is reached.
In SURF, you do not suffer any runtime penalty from the size of the Gaussian filter thanks to the use of the integral image trick. Thus, you compute directly the image filtered at each scale (without using the result at the previous scale).

The approximation part

Since computing the Gaussian scale-space and the histograms of the gradient direction is long, it is a good idea (chosen by the authors of SURF) to replace these computations by fast approximations.

The authors remarked that small Gaussians (like the ones used in SIFT) could be well approximated by square integrals (also known as box blur).
These rectangle averages have the nice property to be very fast to obtain thanks to the integral image trick.

Furthermore, the Gaussian scale-space is actually not used per se, but to approximate a Laplacian of Gaussians (you can find this in the SIFT paper).
Thus, you do not need just Gaussian-blurred images, but derivatives and differences of them. So, you just push a bit further the idea of approximating a Gaussian by a box: first derive a Gaussian as many times as needed, then approximate each lobe by a box of the correct size. You will eventually end up with a set of Haar features.

Increment by 2

This is just an implementation artifact, as you have guessed. The goal is to have a central pixel.
The feature descriptor is computed with respect to the center of the image patch to be described.

Middle region

When going from a black ray to a white ray, you have something like $\sum_{\text{all pix in column}} \partial x = A$. Then, going from white to black, you have the opposite sum: $\sum_{\text{all pix in column}} \partial x = -A$. Thus, you have a small $\sum \partial x$ for the window, but a higher sum of the magnitudes.

Magic number

The first scale is obtained by applying a blur with $\sigma = 1.2$ (or 1.4 in some papers). This is because a natural (real) sharp image can be considered as being the result of the convolution of an ideal (without aliasing) image with a blur kernel of width $\sigma = 1.2$. I can't really remember where it comes from, but it was also explicitly studied in Guoshen Yu's work on A-SIFT, so you might check this page.


MRI Slice image reconstruction I have an MRI K-Space data 320 x 320 x 256 x 8 (4D complex double) from < mridata.org > . The data represents 320 x 320 K-Space of 256 slices from 8 channels. I am trying to reconstruct images for each slice. Here is the matlab code I tried:

% Version - 1
kspacedata= kspacefile.kdata;
imRef = ifftshift(ifftn(kspacedata));
imSOS    = squeeze(sqrt(sum( abs(imRef).^2, 4))); % sum-of-squares to combine all channels

% Version - 2
kspaceSlicedata = kspacedata(:,:,100,:);
imSliceRef = ifftshift(ifftn(kspaceSlicedata));
imSliceSOS    = squeeze(sqrt(sum(abs(imSliceRef).^2, 4))); % sum-of-squares to combine all channels

% Plotting
figure;
subplot(1,2,1); imagesc(imSOS(:,:,100));title('Version - 1');axis image; colormap(gray);
subplot(1,2,2); imagesc(imSliceSOS);title('Version - 2');axis image; colormap(gray);

In Version - 1, I take inverse transform on entire kspace data and plot image data corresponding to slice 100. In Version - 2, I take kspace data corresponding to slice 100 and do inverse transform on this kspace and plot the image data. Here is the output images I obtained.

I thought both versions will reconstruct the same way. But images appeared different. How can I take kspace data corresponding to one slice and perform inverse transform to reconstruct slice image?

Thanks.

The two reconstructions are both valid in general. In MRI images can be acquired in multiple 2D slices and as a 3D volume.

When data is acquired in slices, you would take the 2D Fourier Transform of the k-space of each slice (and each channel, if it is a multi-channel coil). Then you would sum the channel data up in a Sum of Squares fashion as you did. The result would be the reconstructed slices.

When data is acquired in a true 3D mode, you do not have separate k-spaces, but only one big k-space. To reconstruct this data, you would do a 3D Fourier Transformation of your k-space (and do so separately for each channel, if the coil had multiple channels). After a Sum of Squares summation of the channel, you have reconstructed the whole volume.

You can easily check, if you data is acquired in 2D or 3D mode: Take the absolute of your k-space data. If you have only one maximum in your data (one DC component), data was acquired in 3D. If there are multiple maxima in the same order of magnitude (one per slice), then data was acquired in 2D mode.

Note that you cannot replace a 3D FFT by applying 2D FFTs for each slice. That is why your second reconstruction failed. Look at an orthogonal "slice" to what you have shown - you will clearly see the typical k-space modulation with the maximum intensity in the middle.

Hope that helps.


VIbrato (pitch variation) function I'm looking for an algorithm for Vibrato audio effect for a cross-platform C++ program of audio processing. Now I use the realization from Chaotic-DAW sources. The problem is that Chaotic-DAW has a huge class hierarchy, and:
 1. I don't want to use such a huge library
 2. I cannot use it on Mac OS X (cannot compile)
 3. I will not be able to port the program to a microcontroller.

The necessary Vibrato parameters: delay, depth, period (cycle length).

Related topic on Software Recommendations.

Update. I have tried to isolate the Vibrato class, but the hierarchy is so large that it's impossible.

After much searching and asking questions, I have written my own realization for the vibrato effect: https://github.com/Bershov/Vibrato-effect. The code is not hard for understanding and works fine. In fact, the vibrato is a kind of frequency modulation. Its key concepts: variable delay, low-frequency oscillation and intelpolation.


How do audiobook players manage to quickly speedup/slowdown playback speed without using more memory? Many audiobook player apps on Android have a reading speed slider which lets us choose the reading speed e.g. 0.5x 0.75x 1.00x 1.25x 1.50x. This is time-stretching without pitch shifting.

I have some idea about how this could be done using PSOLA algorithm(s), but they would need extra memory equivalent to the memory of original audio file.

ex.

Original audio file -> PSOLA -> Modified audio file

But these Andorid apps don't seem to use double the memory when we try to change reading speed. How do they do an in place time stretching/shortening so quickly?

edit. I found this YouTube engineering blog, which led me to sonic library, which in turn led me to a description of PICOLA algorithm.

I don't understand the algorithm from the images though. Could anyone please explain? What do those triangles mean? How do they decide these A and B points?

Time scale expansion procedure of the PICOLA

WOW long gap to answer this question, maybe if I have a Einstein brain I can understand the image, what is not my case lol.
I haved touched in this algorithm in a long time ago, and I remember that was a surprised for me that almost no one around the Web talks about the Algorithm from Morita Naotaka, I tryed seach his original paper but I cant found online, I think that it is in Japanese, well I dont know ...
So I learned about it looking at among number of patents that make references to the Morita Naotaka, and I freeze when I see a lot of patents around 2000 years that use the same code.
As far as I know this algorithm was derived from TDHS, has advantages about TDHS, this does the job using just Period size Buffer. And believe me it's as easy to code as TDHS.
So lets go, take a look at the Yamaha Patent number US8085953 this is basicly the PICOLA from Morita Naotaka...
Here one image from this Patent that is infinity better to understand how this works to compress :

Looks like one TDHS using a triangular window, but take a look how it are cutting wave A and wave B in the exact Period size(Lp), then they are crossfaded A + B. this image seems to be an compression using rate = 0.6 factor(alpha),  L can be determined using:
$$ \begin{align}
L & = \operatorname{round}\left(\frac{L_p}{\alpha-1}\right)
\end{align} $$
And here how its works to expansion:

Now L can be determined using:
$$ \begin{align}
L & = \operatorname{round}\left(L_p\frac{\alpha}{1-\alpha}\right)
\end{align} $$
Time ago I show a basic code to compress based in TDHS here, now lets try a litle bit more complex and complete example to expand (it will slowdown), so based in this images/equations, here my PICOLA Algorithm:
alpha=0.5;

f=735;
Fs=44100;
signal= 0.9*sin(2*pi*f/Fs*(1:44100)); signal = signal';

period= 60; % is it (Fs/F)

[signal,Fs]=audioread('ederwander.wav');

nsamples = length(signal);
out=zeros(1,floor(nsamples/alpha));
inptr=1;
outptr=1;

minF = 50;
maxF = 900;

minP = floor(Fs/maxF);
maxP = floor(Fs/minF);

winsize = 1024;

while ( inptr+winsize <= nsamples )

    %chunk used to find the period, its need be >= MaxP
    windowedSignal=signal(inptr:inptr+winsize);

    %Basic Autocorrelation to find the period
    Maximum=-Inf(1);
    for P =minP:maxP
        ac = sum(windowedSignal(1:(winsize)-P) .* windowedSignal(P:(winsize -1)));
        if ac > Maximum
            period=P;
            Maximum=ac;
        end
    end

    period = period-1;

    %split all how shown in the patent picture
    waveA = signal(inptr:inptr+period-1);
    waveB = signal(inptr+period:(inptr+period*2-1));
    CrossfadeAB = (waveA .* linspace(0,1,period)') + (waveB .* (1-linspace(0,1,period))');
    L = round(period*alpha/(1-alpha));
    waveB_T0 = signal(inptr+period:inptr+L);
    out(outptr:outptr+period-1)=waveA;
    out(outptr+period:outptr+period*2-1)=CrossfadeAB';
    out(outptr+period*2:outptr+period+L)=waveB_T0';

    %increment input and output pointers
    inptr = inptr + L;
    outptr = outptr + L + period;

end

%not best choice, just append the end to match the out size, this will click, maybe need crossfade here
if outptr < nsamples/alpha

    fim=floor(nsamples/alpha)-outptr;
    out(outptr:floor(nsamples/alpha))=signal(nsamples-fim:nsamples);

end

plot(out);
sound(out,Fs);

I think that standard algorithm of PICOLA just get nice results to 0.5<=alpha<=2
And one more time it's can be considered a proof of concept to show how its can be done!


Implementing the noisy AR(1) Process I'm trying to implement noisy AR(1) process and plot it. The observed noisy sequence x(n) = s(n) + w(n) where variance of w(n) = 0.2.

s(n) is defined as an AR(1) process with s(n) = 0.5 s(n-1) + e(n), for n=1:100, s(0) = 0 and the variance of e(n) = 1. w(n) and e(n) are independent.

I first do the following:

N = 100;                       % number of process samples.
a = [1, 0.5];                  % denominator coefficients, p = 1.
b = 1;                         % numerator coefficient.
s = filter(b,a, randn(1,N));   % generate N sample of AR(1)

nu = sqrt(0.2) * randn(1,N);
x = s + nu;

figure(1)
plot(s,'r');
hold on;
plot(x,'g');
legend('AR(1) Process','Noisy Process');

Can anybody help me whether I do the correct thing. Any help would be appreciated.

Yes, the method you have used generates an AR-(1) random process by inputting Gaussian white noise of unit variance into the LTI filter defined by the coefficeints $a$ and $b$.

Then you are adding some uncorrelated noise to it, which means your random process is not anymore a pure AR-(1) but a noisy one.


Why does my GAIN remain constant after a few cycles? I am assuming that GAIN is the matrix P?

From this example:

function [x,P]=ukf(fstate,x,P,hmeas,z,Q,R)
% UKF   Unscented Kalman Filter for nonlinear dynamic systems
% [x, P] = ukf(f,x,P,h,z,Q,R) returns state estimate, x and state covariance, P
% for nonlinear dynamic system (for simplicity, noises are assumed as additive):
%           x_k+1 = f(x_k) + w_k
%           z_k   = h(x_k) + v_k
% where w ~ N(0,Q) meaning w is gaussian noise with covariance Q
%       v ~ N(0,R) meaning v is gaussian noise with covariance R
% Inputs:   f: function handle for f(x)
%           x: "a priori" state estimate
%           P: "a priori" estimated state covariance
%           h: fanction handle for h(x)
%           z: current measurement
%           Q: process noise covariance
%           R: measurement noise covariance
% Output:   x: "a posteriori" state estimate
%           P: "a posteriori" state covariance

I plot P[0,0] is that correct?

$P$ is the state covariance matrix, not the Kalman gain, though the Kalman gain $K$ is usually calculated using it.

It is normal in simple Kalman filters (ones that do not vary with time) for the state covariance matrix to stabilize after a while so that it is constant.

See, for example, the plots in Figure 3(b) and 3(c) in this PDF which show how both the Kalman gain $K$ and the state covariance $P$ go to a constant value after a relatively short space of time.

Because this happens, the Kalman filter evolution is often ignored in practice and the steady state version is what is implemented.


How can I extract a sinusoid from a very noisy signal for computing the phase shift? I have done some measurements and I got two very noisy signals. I want to extract a certain sinusoid from the signals and I want to compute the phase shift between them.

I made some simulations and I saw that the noise modifies a lot the value of the phase shift.

How can I get a trustful value for the phase shift?

You are correct. Estimating the instantaneous phase of a noisy sinusoid is NOT easy.  I suggest you design a narrow bandpass filter such that your sinusoid-of-interest is in the filter's passband.  (The better the filter the more noise that will be eliminated.) Pass your two signals through the bandpass filter to generate filtered signals $x_1[n]$ and $x_2[n]$.  Next, pass your $x_1[n]$ and $x_2[n]$ signals through a Hilbert transformer to generate $\hat{x}_1[n]$ and $\hat{x}_2[n]$.  Create two analytic (complex) signals as:

$$z_1[n] = x_1[n] + j \, \hat{x}_1[n],$$

and  $$z_2[n] = x_2[n] + j \, \hat{x}_2[n]$$

where

$$ \begin{align}
\hat{x}[n] & = \mathcal{H}\{ x[n] \} \\
           & = \sum\limits_{m=-\infty}^{+\infty} \frac{1 - (-1)^{m}}{\pi \, m}   x[n-m]  \\
\end{align} $$

Next, compute two instantaneous phase sequences:

$$\phi_1[n] = \arg\{z_1[n]\}$$

and $$\phi_2[n] = \arg\{z_2[n]\}.$$

Finally, compare the instantaneous phase difference between the $\phi_1[n]$ and $\phi_2[n]$ sequences.


How to come true the convolution with a 2-D array in frequency domain with python I use the pretty simple example used in many books to understand the convolution in the frequency domain. But I cannot find the real results.
import numpy as np

img = np.array([[0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 1, 0, 0],
                [0, 0, 0, 0, 0],
                [0, 0, 0, 0, 0]])
kernel = np.array([[0, 0, 0, 0, 0],
                   [0, 1, 2, 3, 0],
                   [0, 4, 5, 6, 0],
                   [0, 7, 8, 9, 0],
                   [0, 0, 0, 0, 0]])

fft_img = np.fft.fft2(img)
fft_kernel = np.fft.fft2(kernel)

fshift = (fft_img* fft_kernel)  # Convelution theory at frequency domain

img_back = np.fft.ifft2(fshift)

img_back = np.real(img_back)
img_back = img_back.astype(np.int)

print(img_back)

The output of the img_back is
[[9 0 0 7 8]
 [0 0 0 0 0]
 [0 0 0 0 0]
 [3 0 0 0 1]
 [6 0 0 4 5]]

However, the answer should be
[[0 0 0 0 0]
 [0 9 8 7 0]
 [0 6 5 4 0]
 [0 3 2 1 0]
 [0 0 0 0 0]]

So, is there someone who can help me? I don't know the reason why. And how to come true the ideal result?

This is because the identity for the convolution is
[[1, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0],
 [0, 0, 0, 0, 0]]

You get the correct result by replacing
img_back = np.real(img_back)

by
img_back = np.fft.ifftshift(np.round(np.real(img_back)))

It is important to np.round (to the nearest integer) before casting to np.int, otherwise you will may have wrong results due to rounding errors e.g. np.array([0.999999999999]).astype(np.int) = 0
Notes about fftshift
The the fftshift and ifftshift operations are equivalent to a convolution, we could express fftshift(X) as $S * X$.
Convolutions are commutative and associative $a * (b * c) = (a * b) * c$, and $a * b = b * a$, it is immediate in the frequency from the multiplication operator properties.
Lets say $X$ is your image $H$ is your kernel
you can check that $S^{-1} * X$, or $X = S$, is precisely the term that I mentioned in the beginning of the answer is the identity for the convolution
What you have from my solution above is $S^{-1} * (H * X)$, the solution you metioned in the comments is $H * (S^{-1} * X)$, and there would be one other variation $(S^{-1} * H) * X$ that would represents calling fftshift on the kernel instead of the image. I already saw some doing ffthisft in both inputs and then ifftshift to the output.
Alternative forms...
$(S^{-1} * X) * H = X * (S^{-1} * H) = S * (S^{-1} * X) * (S^{-1} * H)$