Matching LPC frequency response to raw sample data frequency response Background:
The intention is feature/formant extraction from speech input.
I have used an FFT on raw samples in the past and have done simple peak picking from a weighted average of frequencies over time. That gave results that approximated a voiced non-nasal utterance's location in the vowel quadrilateral, but it was not consistent. I heard that LPC would give a smoother response curve and much greater accuracy in feature extraction. The residual is another desired side-effect which will hopefully allow for analysis of other features as well.
Issue:
I have a large disparity in the peak frequencies between an FFT on the raw samples and an FFT on the derived LP coefficients.
I find that the response curves are wildly different in shape and peak frequency locations. Often FFT_LPC (see definitions below) will show large peaks in the 3000/4000Hz range for simple voiced vowel utterances. Perhaps these are pops - maybe I should perform IFFT and play back the output to check. FFT_Raw peaks precisely where I would expect given the source samples.
Process:
Here is my current method:

Obtain samples at 22050Hz.
From the raw sample frame, window with Hamming and perform FFT. Call this FFT_Raw.
In parallel, extract LP coefficients from the raw samples, do NOT window, and perform FFT directly on the LP coefficients padded with zeroes. In effect, every frame would lead with these coefficients. Call this FFT_LPC.
Take RMS of raw samples and apply an (admittedly arbitrary) multiplier to bring FFT_LPC to comparable magnitude compared with FFT_Raw.
Display both frequency response curves with one overlaid on top of the other.

I have tried anywhere from 8 to 32 coefficients with similar results.
Questions:

Is using the coefficients directly as input to the FFT the right approach? The assumption on my end is that leading with the coefficients would approximate a unit impulse. However is there any intermediate processing that I should be doing (such as somehow shifting the starting point of each set of coefficient data to match the fundamental frequency)?
Should the LPC data be windowed prior to FFT? I would think that my current approach of leading with the coefficients every frame means that windowing would reduce the amplitude of the first set of coefficients and FFT would produce different frequencies. However maybe if I window the raw samples and THEN use that as input to the LPC data then I will obtain more accurate results.

I have not used MATLAB or any other third party applications - the desire is to build this myself using Xcode and its accompanying libraries alone. I can display this code but it is lengthy and spans several classes (including adapter classes between data types and math/FFT operations). Ideally I would seek input on a language-agnostic conceptual level.

Is using the coefficients directly as input to the FFT the right approach?

It depends on what are you going to do next. There are many formant extraction algorithms, you can find recent advanced one in COVAREP https://github.com/covarep/covarep/blob/master/formant/formant_CGDZP.m

  and perform FFT directly on the LP coefficients padded with zeroes.

It seems like a crazy step with not much meaning. Why are you doing this? LP coefficients create very rough spectrum approximation, the rest is in residual, but FFT from LPC, it is something strange.

Windowing does not change much, the difference will not be visible. You only can get a few percent more accurate results from it.


Units for a PSD from square of FT? If I am looking to do a PSD plot on a signal $x(t)$ with units $m$ then the units I am expecting in the PSD would be $\frac{m^2}{Hz}$ but I can't seem to get this result from the following reasoning:

If we consider the fourier transform
$$ F(k) = \int^{\infty}_{-\infty} x(t) e^{-2\pi i k t} dt  $$
I think this should have units $ms$ then the PSD should be
$$|F(k)|^2. $$
Which should have units $m^2s^2$ or $m^2Hz^{-2}$ which is not $\frac{m^2}{Hz}$. There is probably an obvious mistake here but I can't see it, can anyone help?

Your reasoning sounds correct (with unit considerations), you can see a discussion on both continuous and discrete aspects of the Fourier transform in What are the units of my data after an FFT?

What is unclear to me is why  you expect the PSD to be in unit $\frac{m^2}{Hz}$? Powers are powers: if a quantity has unit $\textrm{u}$, its square has unit $\textrm{u}^2$. However, when dealing with physics, the quantities you mention are often somehow periodic, and evaluated via averaging over one period, or asymptotically. So basically you had a normalization factor, proportional to the inverse of a time duration, or $\textrm{s}^{-1}$. Hence, possibly, your $\textrm{Hz}$ factor.


Non Periodic, Deterministic Power Signals Any one know of work on non-periodic but deterministic power signals?  Now one member in this class would be the quasi periodic signals.  I wonder if there is a generalized Fourier analysis of non-periodic deterministic power signals.

Consider the signal
$$x(t) = \sum_{n=-\infty}^\infty \operatorname{rect}(t-n)\sin(2\pi nt)$$ which consists of $|n|$ periods of the signal $\sin(2\pi nt)$ in the interval $\left(n-\frac 12,n+\frac 12\right)$ for each integer $n$. Clearly, $x(t)$ is not periodic. Equally clearly, $x(t)$ is a power signal with average power $\frac 12$, and thus it is a nonperiodic but deterministic power signal. I leave it to the OP to determine whether the signal in question has a Fourier transform in the generalized sense (meaning it involves impulses) or not.


Multiband undersampling In Practical Signal Processing (Marc Owen), an exercise explore the topic of undersampling:

  [...] Suppose you know that an audio signal is a sine wave with a frequency that might be in one of two ranges, either between 120Hz and 130Hz or between 150Hz and 170Hz. What is the lowest sample rate you could use, and still be able to reconstruct the signal from (real-valued) samples? What if the second band extends from 150Hz to 172Hz?

The answers are actually provided at the end of the book (60Hz and 86.666Hz) but without further details.

I have apply some of the theory on undersampling and found the different ranges possible for each frequencies. For example, with $n=3$, I can see that the aliases possible for the range 120Hz-130Hz is $86.666 < fs < 120.0$. This theory applies well for one band but I can't find any explanation on the dependency created by another band. I cannot just take the intersection of both sets as it does not cover the case where the ranges are overlapping.

I also did a GNURadio example to see the influence of the sample rate on each range. By moving the sample rate slider around, I can see that 86.66Hz is indeed the solution (for the second case) but is there an analytical method to find this result?

Thanks

I found the answer in Akos, D.M., Stockmaster, M., Tsui, J.B., & Caschera, J. (1999). Direct bandpass sampling of multiple distinct RF signals. IEEE Trans. Communications, 47, 983-988.

The extra condition between the 2 bands is: $|F_{IF1} - F_{IF2}| >= \frac{BW_{I1} + BW_{I2}}{2}$ where $F_{IF}$ is the intermediate frequency and $BW$ the bandwidth.

There are then 4 equations that describe all these constraints (eq $(1)$ to $(4)$) in the paper. Unfortunately, the first constraint contains a non-linear combination which cannot be solved analytically (by something like mixed-integer programming).

This is confirm in the paper: "An iterative procedure, based on these equations, was developed to determine sampling frequencies that would meet the requirement, the results of which are displayed in Fig. 5.". I conclude there is no better solution than iteratively going through the frequencies under consideration.


Spline-Based Adaptive Interpolation Filters? My understanding of interpolation specific to resampling applications is limited to the concept of inserting zeros, then designing a filter to minimize distortion in the passband and reject the images the zero-insert creates (to desired performance levels), such as what is depicted in a simple interpolate by 4 zero-insert shown below (showing the resulting spectrum before filtering).

I am noting that in creating the filter, we are creating the high order polynomial used for interpolating the input samples. I then understand that with a Spline interpolation method (as a user, and perhaps naively as I haven't really dug deeply into the underlying mathematics) we are doing piecewise interpolations using limited polynomials (therefore low order filters) over a smaller span of the sequence, and in many cases due to decreased ripple this will provide better performance.

My question then is if such "adaptive" resampling filters have been designed either using Spline approaches directly, or otherwise piecewise similar in process to a spline, but in a fashion that is best suited for high rate FPGA processing with a FIR based interpolation filter (I say that to avoid the response "just do a spline", but perhaps that would actually result in the solution even in an FPGA FIR based interpolator...). Before digging into that further has anyone seen such an approach, is it common, or are there known pitfalls why this would be a BAD idea?  Thanks for the thoughts and input!

Look up "Farrow filter". It's essential a set of piecewise polynomial interpolators for the coefficients of an even longer (better) FIR interpolation filter (and without the need for a huge polyphase table).  IIRC, it can be implemented in a straightforward arithmetic hardware pipeline (thus is likely quite suitable for an FPGA).


Extending Goertzel algorithm to 24 kHz, 32 kHz and 48 kHz I'm learning to implement Goertzel's algorithm to detect DTMF tones from recorded wave files. I got one implemented in python from here. It supports audio sampled at 8 kHz and 16 kHz. I would like to extend it to support audio files sampled at 24 kHz, 32 kHz and 48 kHz.

From the code I got from the link above, I see that the author has set the following precondition parameters/constants:

self.MAX_BINS = 8
if pfreq == 16000:
    self.GOERTZEL_N = 210
    self.SAMPLING_RATE = 16000
else:
    self.GOERTZEL_N = 92
    self.SAMPLING_RATE = 8000

According to this article, before one can do the actual Goertzel, two of the preliminary calculations are:

Decide on the sampling rate.
Choose the block size, N

So, the author has clearly set block size as 210 for 16k sampled inputs and 92 for 8k sampled inputs. Now, I would like to understand:

how the author has arrived at this block size?
what would be the block size for 24k, 32k and 48k samples?

Choice of $N$

Dual-tone multi-frequency (DTMF) tones should be of certain minimum duration (40 ms) and quality to be detected (reviewed in ITU TABLE A-1/Q.24). The detection algorithm that you cite uses a block duration of $t=N/f_s,$ where $N$ is the block length in samples and $f_s$ is the sampling frequency. The duration is either 210/(16 kHz) = 13.125 ms or 92/(8 kHz) = 11.5 ms. To ensure that the algorithm functions as expected, the block duration should be kept the same. I don't know why the lengths differ for the different sampling frequencies in the cited algo, but it should be okay to calculate $N$ as:

$$N = t\times f_s = 11.5\text{ ms}\times f_s,\tag{1}$$

which should give the same timing and frequency resolution as with the 8 kHz sampling frequency used in the original sources on Wikipedia by an unknown author.

Other changes needed

There are some worrisome large hard-coded constants (4.0e5 and 1.0e9) in the source code. Such hint to a lack of normalization:

#Check for minimum energy
if self.r[row] < 4.0e5:
    msg = "energy not enough"
elif self.r[col] < 4.0e5:
    msg = "energy not enough"
else:

#AT&T states that the noise must be 16dB down from the signal.
# Here we count the number of signals above the threshold and
#there ought to be only two.
if self.r[max_index] > 1.0e9:
    t = self.r[max_index] * 0.158
else:
    t = self.r[max_index] * 0.010

They may require changing if you change $N,$ by the same proportion, I think.

Fine-tuning $N$

Even if the Goertzel algorithm is not restricted to frequencies of a discrete Fourier transform (DFT) of length $N,$ (see this answer) there may be reason for tuning $N$ that is not explained by duration specifications. For sinusoidal input of frequency $f$, the magnitude of the same frequency detected by a Goertzel algorithm fluctuates slightly as function of time depending on the combination of frequency $f$ and $N$. For simplicity of the expressions that follow, frequency is expressed in radians as:

$$\omega=\frac{2\pi f}{f_s}\tag{2}$$

For initial phase $\alpha$ at the start of the window, the unnormalized detected magnitude will be:

$$A(\alpha) = \left|\sum_{k=0}^{N-1} e^{i \omega k}\cos(\omega k + \alpha)\right|\\
= \left|\frac{Ne^{-i\alpha}}{2} + \frac{e^{i(\alpha - \omega+\pi/2)} - e^{i(\alpha + 2N\omega - \omega + \pi/2)}}{4\sin(\omega)}\right|\tag{3}$$

The fluctuation can be seen in this example:

Figure 1. Blue: detected magnitude $A$ as function of initial phase $\alpha$ for sinusoidal input, with $f = 1209\text{ Hz},$ $f_s = 8 kHz,$ $N=92$. Red: extrema of $A(\alpha)$ have been marked.

An interesting observation is that there will no fluctuation of the Goertzel detected magnitude $A$ if the frequency $\omega$ is a multiple of $\pi/N.$ In contrast, DFT bins are at multiples of $2\pi/N.$ (Zaplata & Kasal 2005 discuss the fluctuation phenomenon, but do not recognize this difference between the Goertzel algorithm and DFT.)

For other frequencies, the $\pi\text{-}$periodic function $A(\alpha)$ gets its smallest and largest values (not necessarily in that order) at:

$$\alpha_1 = \frac{\omega(1 - N)}{2}\quad\text{and}\quad\alpha_2 = \frac{\omega(1 - N)}{2} + \frac{\pi}{2}\tag{4}$$

You could tune the integer $N$ to its nearby value which minimizes the maximum of $\left|A(\alpha_1)-A(\alpha_2)\right|$ over the set of DTMF frequencies. Actually, for $f_s = 8\text{ kHz,}$ a better choice than $N=92$ would be $N=93:$

Figure 2. Maximum peak-to-peak fluctuation $\max\left|A(\alpha_1)-A(\alpha_2)\right|$ over DTMF frequencies has a minimum at $N=93$ for $f_s=8000\text{ Hz.}$

A reason why $N=92$ might have been preferred is that in some older version they only allowed DFT bin frequencies. In that case $N=93$ would be a pretty bad choice, as some DTMF frequencies would fall almost just between successive bins, which is perfectly fine for the Goertzel algorithm (or generalized DFT with a frequency shift of $1/2$ binwidth), but almost the worst-case scenario when using DFT:

Figure 3. Maximum absolute frequency error over DTMF frequencies when detecting using DFT of length $N$, for $f_s=8000\text{ Hz}$ and expressed in units of DFT bin width.

The final word would be to test the performance in a test suite and to tune $N$ based on that.

Octave source code for Fig. 3

graphics_toolkit("gnuplot");
fs = 8000; #Sampling frequency
f = [1209, 1336, 1477, 1633, 697, 770, 852, 941]; #DTMF frequencies (same unit as fs)
N = 80:100; #DFT lengths analyzed
omega = 2*pi*f/fs;
y=omega'/(2*pi).*N;
z=round(y)-y;
plot(N, max(abs(z)), "x");
ylim([0, 0.5]);
xlabel("N");
ylabel('max|frequency error| (bin widths)');

References

Zaplata F, Kasal M, "Efficient Spectral Power Estimation on an Arbitrary Frequency Scale", Radioengineering, 2015, vol. 24: p. 178-184. DOI: 10.13164/re.2015.0178


Anisotropic Diffusion Outputs I am providing 2 d grey scale images to the anisotropic diffusion function in MATLAB but not getting desired outputs. How many iterations are ideal and how should i approach this problem?I am keeping kappa values between 20 to 100. I am carrying out the operations on the standard lena image by adding salt and pepper noise.

Input Image:

Expected Output:

The output of 80 iterations of anisotropic diffusion with k=10

As best I can understand, you are interested in achieving edge-preserving denoising of the same "quality" given by the algorithm in the paper:

  Giovani Gómez: Local Smoothness in terms of Variance: the Adaptive Gaussian Filter. Proceedings of the British Machine Vision Confernce, 2000.

The images from that paper that you used above (and in the duplicate of this question that you posted a month later) are available in Section 4 of the link you provided, which seems to be a later version of the paper.

There are two things you should try.

First, as @Jean-Yves suggested to you several times, median filtering is a non-linear filter that preserves edges (although it tends to round off sharp corners).  Whether it does so effectively depends on the amount of noise in the image, and your example images are very noisy.

Second, the currently popular edge-preserving denoising algorithm is the bilateral filter.  This is implemented in photoshop, and in GEGL (which is available through the Gnu Image Manipulation Program in the Tools->GEGL Operation menu).

Here's the pinecone with a radius 1 median filter:

With GEGL's bilateral filter (gaussian radius of 4, and "edge preservation" set at 8%):

With GEGL's bilateral filter (gaussian radius of 4, and "edge preservation" set at 4%):

And with median filter of radius 1 followed by GEGL's bilateral filter with gaussian radius of 4 and "edge preservation" set at 50%:

All of which I prefer to the result of Gómez's adaptive filter:

So, my answer to this question (and to How to remove Gaussian noise without destroying the edges?, which you posted on Feb 7, 2012) is: try median filtering and bilateral filtering.


Amplitude of frequency bin in FFT doesn't match time-domain amplitude I have a sinusoidal current that I am sampling at about 357k SPS. The current signal is about 3A pk-pk @ 750hz (top figure). I am wondering why when I take the fft using numpy/scipy's fft function, the amplitude of the fundamental frequency doesn't match.

I have a sample size of 8000 so I divided the fft by 8000 and multiplied by 2 to get the results below.

Current and fft calculation:

# current sense
df['current'] = (df['Viout'] - 2.5)/20/0.003
df['current_ft'] = np.fft.fft(df['current']) / 8000 * 2

Plot:

plt.subplot(211)
plt.title('current')
plt.xlabel('time [s]')
plt.ylabel('current [A]')
plt.plot(df.t_1, df.current)

plt.subplot(212)
plt.title('current fft')
plt.xlabel('frequency [Hz]')
plt.ylabel('|F(f)| / 8000 * 2')
plt.plot(df.freqz, np.abs(df.current_ft), 'o')
plt.xlim(-ax, ax)

plt.tight_layout()
plt.show()

You seem to know how to handle the DFT/FFT scaling so that its output matches the
sampled signal's power/amplitude levels.

Yet there is a discrepancy between your expectation and the obsevred DFT result? There are a few reasons. Your assumptions on the sampling rate or the signal frequency could be wrong. Indeed as it seems from your figures, the signal is not a pure sine wave. In that case the amplitude looses its crisp definition and spectral leakage and main lobe width smearing becomes an issue. For details please search the site for spectral estimation, spectral display etc...

For your convenience, I have produced an ideal case to your description with MATLAB / OCTAVE code. I'm getting exactly what is expected. So this suggests me that your signal is not a pure sine wave. It's either modulated, or your sampling rate and frequency information is not consistent.

clc; clear all; close all;

T = 0.0225;     % observaton interval in seconds
f0 = 750;      % sine frequency in Hz.
Fs = 357E3;     % sampling frequency in Hz.

t = 0:1/Fs:T;   % analog time

N = length(t);  % number of samples taken...

A = 1.5;        % amplitiude of the sine
x = A*sin(2*pi*f0*t);   % sampled analog signal...

figure,plot(t,x);
title('Analog current waveform');

M = N;          % DFT length (mostly equal to sample size)
X = (2/N)*fft(x,M);

f = linspace(-Fs/2,Fs/2,M); % analog frequenc range in Hz.
figure,stem( f , abs(fftshift(X)));
axis([-1500 1500 -0.1 2]);
title('FFT magnitude of the obtained samples');

The resulting figures are as follows:


Is there something as self-completing data transfer, like Sudoku? I would like to know, since I am not that knowledgeable in this subject, if this idea is feasible or widely implemented:

If a data package is sent, I know parity-bytes are used to detect errors. Could we not use a similar concept to self-complete data packages?
Let's say a client requests an image to be sent from a server. The server first sents a function of each row and column, e.g. the sum of all pixel values in each row and column, so 20 values for a 10x10 image. The server then sends the values of exactly as many pixels as are needed to  let only one possible combination of all pixel values remain.

Would this require less or more data packages to be sent, in the best/average/worst case?

What you are suggesting is indeed rooted on concepts that are important in DSP:

On the channel coding side, you have for instance error correction codes (or error correcting codes, both ECC)
On the source coding side, since you mention Sudoku, techniques like matrix completion, sparse sampling or compressive sensing are possibilities to approximately recover incomplete grids.

However, they work on assumptions on data like positivity, sparsity), and approximate recovery. I have not seen yet such a scheme working for universal lossless compression. An argument is derived from the Dirichlet or pigeonhole principle, another one is easier with recursion. Suppose you have a clever scheme that self-completes data packages, based on some functions that combine inputs. For compression, you expect that it turns all  $B$ bit-stored file on a smaller number of bits, say $B-1$ bits or ledd. Then you could apply it again, and again, until you reach $0$ bits, which does not make sense.

So, this does not appear possible, unless you have additional assumptions that could reduce the space of possible data.


Computer exercises and solutions in signal processing I am very interested in this question and can literally not find anywhere a simple answer, so apologize in advance if it does not fit well this community.

I am looking for either an online resource or a book with computer exercises and solutions in signal processing (can be really anything in this space). I am interested in practicing independently, and while there are countless sources for practice computer exercises, I could not find any source that has solutions as well. Don't mind buying a $100 text book if needed. Appreciate any help in this direction.

"Digital Signal Processing: A Computer-Based Approach"  by Sanjit Mitra is what you need I guess, especially the exercises at the end of each chapter. There is a booklet on the Internet again by Mitra, named Digital Signal Processing Laboratory Using MATLAB.
The other option could be Practical Signals Theory with MATLAB Applications.


How do I find the fundamental period of the given signal? I am having a problem in finding the fundamental period of the signal $x(t)$ given below:
\begin{align}
x(t) &= 2\cos\left(\frac 45 \pi t\right)\sin^2\left(\frac{16}{3} t\right)\\
&= 2\cos\left(\frac 45 \pi t\right)\cdot \frac 12\left[1-\cos\left(\frac{32}{3}t\right)\right]\\
&= \cos\left(\frac 45 \pi t\right) - \cos\left(\frac 45 \pi t\right)\cos\left(\frac{32}{3} t\right)\\
& = \cos\left(\frac 45 \pi t\right)-\frac 12\bigg\{\cos\left[\left(\frac 45 \pi -\frac{32}{3}\right)t\right]+\cos\left[\left(\frac 45 \pi -\frac{32}{3}\right)t\right]\bigg\}
\end{align}

If your top equation is really
$$
x(t) = 2\cos\left(\frac 45 \pi t\right)\sin^2\left(\frac{16}{3} t\right)\tag{1}
$$
You gonna have a hard time getting the fundamental period/frequency as the there isn't an exact integer relating the two periods/frequencies. You can make approximations of the $\pi$ multiplier/divisor but the errors accumulates and this doesn't cut it. However if your equation is:
$$
x(t) = 2\cos\left(\frac 45 \pi t\right)\sin^2\left(\frac{16}{3} \color{magenta}\pi t\right)\tag{2}
$$
You could then easily extract the fundamental frequency $f_0$, and simply get the inverse $\displaystyle\frac{1}{f_0}$  to get the fundamental period $T_0$. For a signal $x(t)$

$$
x(t) = \sum_{k=1}^N a_k\cos\left(2\pi f_k t+ \phi_k\right)
$$

The fundamental frequency is the greatest common divisor (GCD) of all the frequencies; or
$$\mathrm{GCD}\left\{f_k\right\}\quad\big\vert \quad k = 1, \ldots, N$$

With non-integer values of $f_k$, you get the frequencies to the same common denominator or their least common multiple (LCM) and $f_0$ is the GCD of the numerators over this LCM, their common denominator. In your case, you have the form:
$$
x(t) = \prod_{k=1}^N a_k\cos\left(2\pi f_k t+ \phi_k\right)
$$
But the same principle applies. So, looking at the second-last equation you have two distinct frequencies: $\displaystyle f_1 = \frac 25$ and $\displaystyle f_2 =\frac{16}{3}$.

So in the end, you get $\displaystyle f_0 = \frac{2}{15}\implies T_0 = \frac 1f_0 = 7.5$.

Below is a plot illustrating this with the four waveforms below:

\begin{cases}
\color{cyan}{y_1} &= 2\cos\left(\frac 45 \pi t\right)\\
\color{green}{y_2} &= \sin^2\left(\frac{16}{3}\pi t\right)\\
\color{blue}{y_3} &= 2\cos\left(\frac 45\pi t\right)\cdot \sin^2\left(\frac{16}{3}\pi t\right)\\
\color{red}{y_4} &= 2\cos\left(2\pi\frac{2}{15}t\right)\\
\end{cases}

You see that the fundamental period of your signal $y_3$ ($\color{blue}{\rm blue}$) is at exactly $7.5 \ \rm seconds$ in the plot, matching exactly that of a pure sinusoidal wave $y_4$ ($\color{red}{\rm red}$) at $\displaystyle f_0 = \frac{2}{15}$.


What is the difference between upsampling and up conversion (done during signal transmission)? Also is upconversion different from modulation? Up sampling is when you increase the sample rate of a continuous signal. Modulation is when you, by multiplying with a carrier, shift the base frequency of a signal. But what is up conversion? How is it related with the above mentioned terms?

Your first sentence is not correct.  The terminology of 'sample rate conversion' operations is often misused.  I suggest you replace your first sentence with, "Interpolation is when you increase the sample rate of a DISCRETE signal sequence."  I'd also replace your second sentence with, "Modulation is when you, by multiplying with a carrier frequency, shift a signal's spectral components up in frequency by a fixed amount." (The term "modulation" is typically used when talking about continuous signals.)  As far as I understand the terminology, "up-conversion" is when you shift the spectral components of a discrete sequence up in frequency by a fixed amount. [-Rick Lyons-]


How to combine a rotation matrix and a stretch matrix into a single matrix for easy Fourier Transform For full disclosure, this is related to homework.  I have to find the Fourier Transform of a function that I've boiled down to the following.

I have a function $f(x,y)$ that I can think of as another function $g(x,y)$ plus a stretch and rotation.  (For example, think of turning a circle into an ellipse).  So, I have a stretch matrix of:

$$\pmatrix{ \frac{1}{A}& 0\\\ 0& \frac{1}{B}}$$

And I have the standard rotation matrix:

$$\pmatrix{ \cos \theta& -\sin \theta\\\ \sin \theta& \cos \theta}$$

I know that for computing the Fourier transform, $f(Ax)$ is the same as $\frac{1}{det A} F(A^{-T}u)$.  So, if I can combine the stretch and rotation into one matrix A, I'm home free.  How can I combine these two matrices?

Thanks!

You combine those by matrix multiplication.


How can output sequence be equal to the sum of copies of the impulse response, scaled and time-shifted signals? I am sorry, this is a very basic question. But i am having hard time understanding is how it is possible.

I know impulse response is the output of the system when impulse sequence is given as input with initial conditions set to 0.

Scaling is to increase the amplitude of the signal, i.e. if i multiply the input by 2, output will also be multiplied by 2.

Time shifted signal is if I delay the input then the output is also delayed by the same factor.

Now can some one please illustrate this with an example how any sequence can be decomposed into sum of copies of the impulse response, scaled and time-shifted signals?

Thanks a lot in advance.

One interpretation of your question could be as follows:

  Given that a system has the following two properties:

  the scaling or homogeneity property that if the response to input $x(t)$ is output $y(t)$, then for any choice of $\alpha$, the system response
  to scaled input $\alpha\cdot x(t)$ is scaled output $\alpha\cdot y(t)$,
  the time-invariance property that for all choices of $\tau$,
  the response to time-delayed
  input $x(t-\tau)$ is time-delayed output $y(t-\tau)$,

  then why does the system have the additivity or superposition property
  that the response to input $x_1(t)+x_2(t)$ is
  $y_1(t) + y_2(t)$ where
  the system response to $x_i(t)$ is $y_i(t)$, $i = 1,2~$ ???? $~~~~~~~~~~~$
  More generally, why is the system response to input
  $\alpha\cdot x_1(t-\tau_1) +  \beta\cdot x_2(t-\tau_2)$
  given by $\alpha\cdot y_1(t-\tau_1) +  \beta\cdot y_2(t-\tau_2)~$?

The answer is that a system with properties 1 and 2 does not
necessarily have the additivity or superposition property.
If the superposition property also holds, then the system
is called a linear time-invariant system. But this is an
additional assumption that you need to make (or prove).

Commonly, homogeneity and additivity are combined together
into the linearity property that says that the response
to input $\alpha\cdot x_1(t)+\beta\cdot x_2(t)$ (that is,
a linear
combination of inputs $x_1(t)$ and $x_2(t)$) is
$\alpha\cdot y_1(t) + \beta\cdot y_2(t)$
(that is, the same linear combination of outputs $y_1(t)$ and $y_2(t)$).

A couple of points that should be tucked away into the
back of one's mind:

A system can be linear without being time-invariant (e.g.
a modulator $x(t) \to x(t)\cos(\omega t)$, or
time-invariant without being linear (e.g. a square-law
circuit $x(t) \to [x(t)]^2$
A additive system which produces output $y(t) + y(t) = 2y(t)$
in response to input $x(t) + x(t) = 2x(t)$ and so seems
to have the scaling property does not in fact have the
scaling property. Persuade yourself that this is true by
attempting to prove that the response to $0.5x(t)$ is
$0.5y(t)$. In short, scaling and additivity are two
different properties and a system that enjoys one of them
does not necessarily enjoy the other.

A second interpretation of your question could be as follows:

  For a linear time-invariant system, the output is supposed
  to be the sum of scaled and time-delayed versions of the
  impulse response, but I don't see how this is so. For example,
  the standard convolution result (for discrete-time systems)
  says
  $$y[n] = \sum_m x[m]h[n-m]$$
  where $h[\cdot]$ is the impulse (or unit) response of the
  system. But this seems to be completely backwards since the
  impulse response is running backwards in time (as in $-m$
  in the argument of $h$ in the above formula
  compared to $x[m]$ in which time is running forwards.

This is indeed a legitimate concern, but actually the
convolution formula is very successful in concealing
the result that the output is the sum of scaled and time-delayed
versions of the impulse response. What's going on is as follows.

We break down the input signal $x$ into a sum of scaled
unit pulse signals. The system response to the unit pulse signal
$\cdots, ~0, ~0, ~1, ~0, ~0, \cdots$ is
the impulse response or pulse response
$$h[0], ~h[1], \cdots, ~h[n], \cdots$$
and so by the scaling property the single input value $x[0]$,
or, if you prefer
$$x[0](\cdots, ~0, ~0, ~1, ~0,~ 0, \cdots)
=  \cdots ~0, ~0, ~x[0], ~0, ~0, \cdots$$
creates a response
$$x[0]h[0], ~~x[0]h[1], \cdots, ~~x[0]h[n], \cdots$$

Similarly, the single input value $x[1]$ or creates
$$x[1](\cdots, ~0, ~0, ~0, ~1,~ 0, \cdots)
=  \cdots ~0, ~0, ~0, ~x[1], ~0, \cdots$$
creates a response
$$0, x[1]h[0], ~~x[1]h[1], \cdots, ~~x[1]h[n-1], x[1]h[n] \cdots$$
Notice the delay in the response to $x[1]$.  We can continue further
in this vein, but it is best to switch to a more tabular form
and show the various outputs aligned properly in time. We have
$$\begin{array}{l|l|l|l|l|l|l|l}
\text{time} \to & 0 &1 &2 & \cdots & n & n+1 & \cdots \\
\hline
x[0] & x[0]h[0] &x[0]h[1] &x[0]h[2] & \cdots &x[0]h[n] & x[0]h[n+1] & \cdots\\
\hline
x[1] & 0 & x[1]h[0] &x[1]h[1] & \cdots &x[1]h[n-1] & x[1]h[n] & \cdots\\
\hline
x[2] & 0 & 0 &x[2]h[0] & \cdots &x[2]h[n-2] & x[2]h[n-1] & \cdots\\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots & \\
\hline
x[m] & 0 &0 & 0 & \cdots & x[m]h[n-m] & x[m]h[n-m+1] & \cdots \\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array}$$
The rows in the above array are precisely the scaled and
delayed versions of the impulse response that add up to
the response $y$ to input signal $x$.
But if you ask a more specific question such as

  What is the output at time $n$?

then you can get the answer by summing the $n$-th column to get
$$y[n] = x[0]h[n] + x[1]h[n-1] + x[2]h[n-2] + \cdots + x[m]h[n-m] + \cdots
= \sum_{m=0}^{\infty} x[m]h[n-m],$$
the beloved convolution formula that befuddles generations of students
because the impulse response seems to be running backwards in time.


How to determine a very low frequency signal? I have a noisy signal with frequency between 1/3 and 1/5 of the sampling window.

DFT has a very low resolution for low frequency components, what is the best way to find the frequence?

Check out this article by Eric Jacobsen and Peter Kootsookos:
http://www.ingelec.uns.edu.ar/pds2803/Materiales/Articulos/AnalisisFrecuencial/04205098.pdf

I had some success using this method for doing "sub-bin" frequency estimation when only one sinusoid is present.


Calculating a homography based on detected lines I know you can calculate homographies from image to camera plane using correspondence points between a "perfect model" and the image points.

I'm doing it for a football pitch/field, and have used edge detection to find the white lines in the pitch.

But the camera does not (always) cover all of the pitch, so I can't see all the corners... and I only the corners are 100% known points in the model (no other distinguished points).

So the problem is that unless the line intersects with another line and forms a corner, I only know the image points of the line, not it's corresponding "perfect/real-world" coordinates in the model.

Is there some way I can use the detected lines to calculate a homography, or even just a set of candidate homographies, even if the detected lines do not intersect with each other and create a corner?

Example image, showing the pitch, our field of view, and the points of the pitch where I can know the corresponding realworld/model coordinates(green circles), and an example of 2 lines that might be completely useless since in our field of view, I have no clue exactly at which point they start or stop in the corresponding realworld /model of the pitch:

The red lines are examples of lines which I would like to use, but I don't know their realworld coordinates, and it's kind of hard to estimate them because depending on the camera pose, the correspondent points could be "anywhere".

I will explain two approaches for this:

1)
One approach would require a line matching algorithm. After matching the lines, you could simply use the end points of the lines in order to compute homography. To achieve that EDLine or LSD based descriptors are recently proposed in OpenCV. Also, hashing and fast matching of them are also implemented. Check out the videos here:

http://www.youtube.com/watch?v=MqMjvSkM39k

http://www.youtube.com/watch?v=naSWTlbg3To

The recent opencv_contrib repository contains the source code to these methods.

In the case that the line end points are noisy, you could then directly utilize the lines to compute the homographies. Such papers would then read:

  Internal Report: 2005-V04 Computing Homographies from Three Lines or
  Points in an Image Pair
  G. Lopez-Nicolas, J.J. Guerrero, O.A. Pellejero, C. Sagues

  Internal Report: 2003-V01 Robust line matching and estimate of homographies simultaneously G. Lopez-Nicolas

  Probabilistic Matching of Lines for Their Homography Taemin Kim, Jihwan Woo, and In So Kweon

2) There is one method specific to the fields are given here:

  "Using line and ellipse features for rectification of broadcast hockey video.", Gupta, Ankur, James J. Little, and Robert J. Woodham
  Computer and Robot Vision (CRV), 2011 Canadian Conference on. IEEE,
  2011.

and

  "Combining line and point correspondences for homography estimation.", Dubrofsky, Elan, and Robert J. Woodham. International Symposium on Visual Computing. Springer Berlin Heidelberg, 2008.

The idea is as follows: Any line, parametrized by its coefficients $\mathbf{l}_i=(u,v,1)^T$ maps to $\mathbf{l}_i^{'}=(x,y,1)^T$ in the other image using:

$$
\mathbf{l}_i^{'} = \mathbf{H}^T \mathbf{l}_i
$$

In this form the equation can be directly plugged into the DLT method:

$$
A_i =
\begin{bmatrix}
-u & 0 & ux & -v & 0 & vx & -1 & 0 & x  \\
0 & -u & uy & 0 & -v & vy & 0 & -1 & y
\end{bmatrix}
$$

The only difference is the normalization, which you will find in the references above.

Adding Ellipses:
Any point $x$ lies on the conic section $C$ if $x^TCx=0$. This gives rise to the transformation relation:

$$
C' = H^{-T}CH^{-1}
$$

The references above also explain how to insert this constraint to the DLT algorithm.

Using ellipses and lines, it is possible to derive a robust projective relationship.


Identifying whether or not a cyclostationary signal is noisy or not using the cyclic autocorrelation I am trying to determine whether or not a given signal has been corrupted by Gaussian noise, either bandlimited (with a filter) or not. The signal in question is a BPSK or PAM signal that is upsampled and shaped by a Nyquist filter and then frequency modulated. Due to the signal being cyclostationary, I figured the best way was to use the Cyclic Autocorrelation but I am having a hard time understanding how to code it up in Matlab. This paper shows shifting the FT by some frequency $\alpha$ and then doing $F^{-1}(X(f-\alpha)*X'(f + \alpha))$ where $F^{-1}(\cdot)$ denotes the inverse DFT operation. Would I be able to just multiply it in time by a complex exponential at frequency $\alpha$ or is there more to it? Is this a valid approach or am I on a goose chase?

Thanks

You're on the right track. The cyclic autocorrelation can be estimated by inverse Fourier transformation of the cyclic periodogram. If $X_T(t, f)$ is the Fourier transform of your data block (indexed by time $t$), and $\alpha$ is your desired cycle frequency, you can form the cyclic periodogram by shifting the Fourier transform to the left and to the right by $\alpha/2$ (not by $\alpha$ as you state in your question):

$$I_x^\alpha(t, f) = \frac{1}{T} X_T(t, f+\alpha/2) X_T^*(t, f-\alpha/2)$$

Then just inverse Fourier transform that cyclic periodogram. You'll have to be careful about the case where the cycle frequency $\alpha$ is not equal to an integer number of FFT bins.


Multi-taper Coherence, With Tapers=1 I had a question regarding multi-taper estimated coherence. It appears when the number of tapers is set to 1, all coherency values, per frequency, are always 1. It's obvious solving for autospectra $S_{xx}$,$S_{yy}$ and cross spectra $S_{xy}$ between two complex numbers produces this with a single taper, using $C_{xy}^2=\frac{|S_{xy}|^2}{S_{xx}S_{yy}}$. And common matlab functions also reveal this applied to real data. But I struggle to undestand intuitively why having only a single taper completely botches the estimate of coherence, such that all values are one, when expanding the number of tapers gives what could be thought of as a real estimate of the underlying values. Why does the orthogonal basis of tapers have to have $K\ge2$ to obtain a sensible values? Are there any analogies or physical metaphors as to why that is?

For any estimate of the coherence, the cross spectra and spectra need to be averaged before forming the estimate of the coherence.  If there is no averaging, the squared coherence will always be 1.  You can see this by plugging in a single Fourier coefficient for X and Y: $C_{xy}^2=(\alpha_x\alpha_y)^2/(\alpha_x^2\alpha_y^2)=1$.  If you instead use a sum or average of two Fourier coefficients for X and Y, the answer can be less than one.

In multitaper spectral analysis, this averaging is done by averaging spectral estimates made using the same time series but different taper windows.  If you only use one taper window, there is no averaging being done.


FFT of sine wave not coming as expected i.e single point The cyan plot is a spectrum of 50 Hz, and the magenta one is a 50.1 Hz sine wave (having amplitude 0.7). Both are sampled at 1024 samples/s. I performed a 1024 point FFT to get this spectrum.

Why is only the 50Hz spectrum a single value? Why does the  50.1 Hz sine consist of other frequencies apart from 50.1 Hz; where do these new frequencies come from?

I did not do any non-linear processing on the 50.1 Hz signal! Also the 50.1 Hz appears to have smaller maximum amplitude, i.e. it's not 0.7, when in fact the sine wave I generated has an amplitude of 0.7.

Why is this?

Obtained by MATALB command fft();

Acually Matt's answer already gives one view on the problem here: the DFT is implicitly periodic in both time and frequency domain (see this question). From your parameters we can calculate that your observation period is 1 s. That means you observe 50 periods of a 50 Hz tone. Periodically extending that observation interval will always result in a seemless sine wave. If you take the 50.1 Hz tone, you're transforming 50.1 periods of an oscillation. Periodically extending that signal will result in phase jumps that cause additional spectral tributaries.

Another view is the frequency resolution that is defined as $f_\mathrm{s}/N_ \mathrm{DFT} = 1024 \text{Hz}/1024 = 1 \text{Hz}$. As 50 Hz is an integer multiple of 1 Hz, this tone can be reresented by a single DFT bin. As the result of the DFT is discrete values that are space 1 Hz apart, a 50.1 Hz cannot be represented by a single bin and the energy is smeared over the frequency domain. In your example this is the reason for the slightly asymmetric magenta spectrum and for the lower magnitude of the 50 Hz bin compared to the 0.7 magnitude of the 50 Hz spectrum.

Both effects described above contribute to the spectrum you're observing.


Which method should be selected to enhance the compressed image? I have a raw infrared image with 16384(2^14) intensity levels.The raw image must be compressed into one with 256(2^8) intensity levels. In order to retain the original details( small blocks) , the compressed image must be enhanced. Which the best method I should select?Someone advice that "histogram matching" method should be taken.But the both images are different in intensity levels,how can do it?Any advice will be appreciated!

I would advise you to try between two different approaches, then test and see which one performs best.

The first one is intensity clustering. In this approach, you data points are one-dimensional and correspond to the intensity of your raw pixels. Then, run a clustering algorithm such as $k$-means with $k = 256$, and perform the final quantization by replacing the intensity of each pixel by the mean value (or the index) of the corresponding cluster mean. With this approach you can limit the information loss by forming an image of indices and keeping a look-up table of the real (14 bits) cluster mean values.

The second would be to perform some local histogram equalization, i.e., cutting your image into blocks and performing equalization (or any other enhancement procedure) on each block.
You will obtain an image where the features are mstly preserved and that you can remap linearly into $[0,255]$.


What would the target matrix to train Neural Network? I'm new at Artificial Neural Network and I'm using MATLAB developing Facial expression recognition and There are six expressions ; I'm not able to understand about How to create a target matrix? My Question is if I have SampleMatrix containing samples rows*columns. While rows are the features and 29 columns represents a single sample I meant to say 29 columns belongs to the same class and If SampleMatrix dimensions are 21007*174. 174 columns because I have 6 classes 6*29=179 What would the target matrix?

Kindly Help me out!
Thanks.

Example with the following numbers (I use random numbers):

There are 6 possible classes (face expressions in your case).
You have 2000 examples (lets say 2000 photos of faces from which you know the correct face expression).
There are 30 features. From each example you have extracted the 30 features and know the correct class.

With these numbers, the input matrix would be a matrix of $30 \times 2000$, where each column $j$ is an example, and it has the 30 features of it. Hence, the position $(i,j)$ contains the feature $i$ of the example $j$.

The target matrix would be a matrix of $6\times 2000$. Where each column $j$ represents the correct class of the example $j$. The column is formed by zeros, and it has a $1$ at the row that indicates the class. If we have 6 classes, and the first example corresponds to class 3, the first column of the target matrix will be $[0 \ 0 \ 1\  0\  0\  0]$. If the second example corresponds to class 5, the second column of the target matrix will be: $[0 \ 0 \ 0 \ 0 \ 1\  0]$, and so on.

More info here.


Improving Digitization by Adding White Noise Before Oversampling I was reading a physics book and it said that when digitizing a signal white noise was added to improve the process. I don't understand how this works. Therefore my question is:

When digitizing an analogue signal, how can adding noise be beneficial
when combined with oversampling in time?

Adding uncorrelated (i.e. white) noise to an analog signal prior to digitization is called dithering.
To understand why we would do this we need to understand the idea of quantization noise.
Consider an analog system which has signals ranging in amplitude from 0 to 100.
Suppose we digitize this signal with a digitizer whose digital levels are spaced by 1.
In other words, the possible digitized levels are

$$\left\{ -100, -99, -98 \ldots 99, 100 \right\} \, .$$

Now suppose the analog signal $s(t)$ is a DC signal of value 0.8, in other words

$$s(t) = 0.8 \, .$$

If we put this into the digitizer, the digitizer will round it to 1 and our digital samples $s_n$ will be

$$s_n = 1\,.$$

This is not good because now the our digital signal accumulates error as we acquire more signal.
The digitized level is always too high so the longer we average the signal the more we over-estimate the analog level.

Adding white noise helps fix this problem because it pushes the analog level around such that it crosses neighboring digitization levels.
Therefore, as  you average over a set of digitized values you actually get something which is close to the true analog level.
Let's see this via example.

Suppose the noise we add is Gaussian distributed with $\sigma=2$.
Then the distribution of the analog signal is

$$p(x)
\propto \exp \left[ \frac{-(x - s(t))^2}{2\sigma^2} \right]
= \exp \left[ \frac{-(x-0.8)^2}{8} \right] \, .
$$

We now compute the average digital value $\langle s_n \rangle$ by averaging $p(x)$ over the integers.
The normalization constant for the distribution is
$$N = \sum_{m=-100}^{100} \exp \left[ \frac{-(m-0.8)^2}{8} \right]$$
and so we have
$$\langle s_n \rangle =
\frac{1}{N}
\sum_{m=-100}^{100} m \exp\left[ \frac{-(m-0.8)^2}{8}\right] = 0.79999\ldots$$
Thus you can see that adding the white noise caused the average digitized signal to more closely match the true analog value.

Of course, adding the noise makes your signal to noise ratio worse.
That means to to actually have a high probability of measuring the $\langle s_n\rangle$ we just computed, you have to take more samples than you would think in the noiseless case.
This is why you hear about over-sampling and dithering at the same time.
Because the dithering noise is truly uncorrelated, taking more samples always helps improve the signal to noise ratio, even if you're sampling way above the bandwidth of the analog signal coming in.


Frequency modulation audio synthesis: clicks/pops when changing modulator frequency dynamically I'm very new to DSP and working on a FM VST. It's mostly up and running, or at least it generates FM-like sounds. I want to be able to change various parameters dynamically while audio is generated, such as the modulation index, for automations and later for LFOs.

I had a problem with clicks and pops when changing values, and solved it by interpolating the values over a small number of samples. However, this doesn't seem to work when dynamically changing the frequency of the modulator. It just causes strange frequencies during the transition, which is not what I want.

I'm speculating that maybe I could introduce the new frequency immediately just when the carrier has completed a sine cycle, and implemented that, but that doesn't remove the clicks (which might be because of some bug in my code, or because this is the wrong approach completely.)

As you might have noticed, I'm mostly unaware of the mathematics involved.

A very simplified version of my (Rust) code, which should however show the basic method I'm using, is:

let mut time = 0.0; // EDIT: is not reset with each buffer loop
let time_per_sample = 1.0 / 44100.0;

let carrier_frequency = 440.0;
let modulator_frequency = 880.0;
let beta = 2.0;

for sample in buffer.iter_mut(){
    *sample = (
        (carrier_frequency * time * 2 * PI) +
        beta * (modulator_frequency * time * 2 * PI).sin()
    ).sin();

    time += time_per_sample;
}

I'm currently not scaling beta with modulator frequency, if that matters.

In FM synthesis, how do I change the modulator frequency while a sound is playing without clicks or pops?

If you change modulator frequency, do note that it changes a sudden jump in phase, as it changes the result of the equation where the phase is calculated from frequency. This is largely unnecessary, as having a certain frequency means the phase change is constant per sample. Thus, if you change the frequency to be a phase increment per sample, then changing the frequency does not cause a sudden phase jump but just changes the phase increment from that point on. That's how FM synth chips work.


Need help implementing DCT Type 4 algorithm I'm trying to implement Algorithm 2 from this paper:

  Xuancheng Shao and Steven G. Johnson: Type-IV DCT, DST, and MDCT algorithms
  with reduced numbers of arithmetic operations, https://arxiv.org/pdf/0708.4399.pdf

It outlines an algorithm for computing an even-sized DCT type 4 via a DCT type 3 and a DST type 3 of half size.

I've made a few attempts at implementing this, but either the paper has a mistake or I have a serious misunderstanding about what's happening here.

https://pastebin.com/uFiHsMKX Here's my C++ code, which uses FFTW to compute the inner transforms, as well as to compute an "expected output". Aside from FFTW all the code necessary to compile is in this one file. The function compute_dct4 is what does the actual processing -- I did my best to name variables the same as in the paper, so it should be clear what in the code corresponds to what in the paper.

I'd appreciate some extra eyes on what's going wrong. Thanks!

I found the solution. There was a mistake in the paper.

On the first 2 lines of Algorithm 2:

w0 ← x0
vN/2 ← xN−1

x0 and xN-1 should be scaled by 2, giving:

w0 ← x0 * 2
vN/2 ← xN−1 * 2

The following paste has working code (and an open source license, so feel free to copy it)

https://pastebin.com/TvgU4b6d


Filter recognition from transfer function I have a filter with this transfer function in the z-plane $H(z)= z^{-k}$. Which type of filter is this (BP,LP,BS)and why?
I have a delay filter with zero/poles equal to zero and the magnitude constant equal to 1, so can I characterize it as an all-pass filter ?

From one of your comments it appears that you've actually already answered the question. Here are a few remarks and questions that should help you gain some more understanding:

If you have a transfer function $H(z)$ (and if you assume that the corresponding system is stable), then the frequency response is obtained by choosing $z=e^{j\omega}$.
For $H(z)=z^{-k}$ figure out the magnitude and the phase of $H(e^{j\omega})$.
If the magnitude of $H(e^{j\omega})$ does not depend on frequency, what kind of filter have you got?
What is the input-output relation of the given system?


How can I decide proper FFT length(size)? I'm trying to understand about deciding FFT length(size) and relative factors.

what If I got the streaming adc data as the below,

1624
20850
10314
-15042
-18908
3717
21175
8479
-16411
-17907
5796
21308
6493
-17713
-16662
7687
21241
4468
-18724
-15291
9750
20929
2392
-19649
-13833
11508
...

Q1. My first question is from here, I want to process with FFT but I don't know how do we decide proper FFT length (size)?

Q2. If I just decide to use FFT8, then can I use FFT8? or if I want to use FFT8 then what should I have to consider what I need to know?
As I know (I'm slowly studying DSP) FFT8's number of inputs are 8. then FFT8 is working as

First term inputs
1624
20850
10314
-15042      =>       FFT out
-18908
3717
21175
8479

Second term inputs
-16411
-17907
5796
21308           =>       FFT out
6493
-17713
-16662
7687

Third term inputs
21241
4468
-18724
-15291          =>       FFT out
9750
20929
2392
-19649

n terms inputs
-13833
11508
...

is this right?

I've read some articles https://electronics.stackexchange.com/questions/12407/what-is-the-relation-between-fft-length-and-frequency-resolution
https://spectraplus.com/DT_help/fft_size.htm But I didn't get it exactly.

So Would you help me to undersatnd how does FFT size decide and calculate process?

Q3, From http://support.ircam.fr/docs/AudioSculpt/3.0/co/FFT%20Size.html, What does 21,53 Hz. mean?
As their description below.

Reminder : Bins
The FFT size defines the number of bins used for dividing the window into equal strips, or bins. Hence, a bin is a spectrum sample, and defines the frequency resolution of the window.

By default :

N (Bins) = FFT Size/2

FR = Fmax/N(Bins)

For a 44100 sampling rate, we have a 22050 Hz band. With a 1024 FFT size, we divide this band into 512 bins.

FR = 22050/1024 ≃ 21,53 Hz.

But I want to know more about 21,53 Hz. exactly FR.

But I want to know more about 21,53 Hz. exactly FR.

This represents the bandwidth (how wide in Hz a bin is) of a particular bin. The smaller this value is, the better FR you have. In order  to get the bin center frequency, you can use the following formula

bin_center_freq = FR * bin_number;

So, for, say, bin at index 20 , we get

bin_center_freq =   21,53*20 = 430,66 Hz (Note that the bin at index 0 represents DC (0 Hz)

When we have a bin center frequency, we can easily calculate the lower and upper frequencies for that bin.

bin_lower_freq = bin_center_freq-(FR/2) = 419,89 Hz

and

bin_upper_frequency = bin_center_freq-(FR/2) = 441,42 Hz

Note that the real frequency resolution (our ability to separate close-by frequencies) will still depend on the type of the window used for the FFT (Blackman, Gaussian, Hamming, rectangular etc) as well as on the spectral content present (whether we have components that are both close in frequency and amplitude etc) in the signal so you shouldn't assume
that a particular bin contains only frequencies belonging to that bin(due to the spectral leakage, every bin will also have energy contributions from neighbouring bins).

So, what is the right FFT size?

It really depends on what you want to do with your signal. If you have a signal containing 2 sine waves close in frequency and amplitude (e.g. one 220 Hz, and the other 225Hz), you should choose a relatively long FFT length such that
the FR is less then 5 (f1-f2). But since windowing also affects the result, make your FFT length even longer to compensate for the windowing effect. You can read up more on window functions here ( https://en.wikipedia.org/wiki/Window_function )


How to make colored noise white? I am studying linear minimum mean squared error (MMSE) receivers. The technique used in it suggests to make colored noise white by multiplying it with the invertible covariance matrix $K_z^{-\frac{1}{2}}$, with $K$  given by

$$K_{zk} = N_0 I_nr +\sum_{i≠k}^{n_t}P_ih_ih_i^*$$

where $h_i$ and $h_i^*$ are the channel gain and channel gain conjugate transpose, respectively, and $n_t$ is the number of transmit antennas. The colored noise is because of interference from other users (I guess so!).

I am confused as how the multiplication of $y = hx + z$ with $K_z^{-\frac{1}{2}}$ can transform colored noise to white though I can guess that may be multiplying with inverse cancel out the colored part and the noise is white now. Am I correct?

Given the linear model $$y = hx + z$$ where $z \sim \mathcal{N}(0,K_z)$. We say that noise is colored when $K_z \neq \Sigma$, a diagonal matrix. To uncolor the noise, you can pre-multiply $y$ by $K_z^{-\frac{1}{2}}$, so
$$K_z^{-\frac{1}{2}}y = K_z^{-\frac{1}{2}}hx +K_z^{-\frac{1}{2}}z = K_z^{-\frac{1}{2}}hx +\tilde{z}$$
Notice that the noise is now $\tilde{z} = K_z^{-\frac{1}{2}}z$, which is uncolored since $$E(\tilde{z} \tilde{z} ^T) =  K_z^{-\frac{1}{2}}E({z} {z} ^T)K_z^{-\frac{1}{2}} = K_z^{-\frac{1}{2}}K_zK_z^{-\frac{1}{2}} =  K_z^{-\frac{1}{2}}K_z^{+\frac{1}{2}}K_z^{+\frac{1}{2}}K_z^{-\frac{1}{2}} = I$$
So $\tilde{z} \sim \mathcal{N}(0,I)$


Getting blobs from similar colored contours in grayscale image I want to select the staff lines as individual blobs in either of these images:

Is this possible?

Here is a simple steps in c++ to identify the Objects you want to detect! Hope this solves your problem!

void main()
{
    Mat mSource_Bgr,mSource_Gray,mSmoothed,mEdge,mThres;
    mSource_Bgr= imread(FileName_S.c_str(),IMREAD_COLOR);
    cvtColor(mSource_Bgr,mSource_Gray,COLOR_BGR2GRAY);

    medianBlur(mSource_Gray,mSmoothed,5); /*Just to smooth the Edges*/
/*This is for Automatically calculating the Canny Threshold values*/
    double CannyAccThresh = threshold(mSmoothed,mThres,0,255,CV_THRESH_BINARY|CV_THRESH_OTSU);
    double CannyThresh = 0.1 * CannyAccThresh;

    Canny(mSmoothed,mEdge,CannyThresh,CannyAccThresh);

    Mat kernal_E,kernal_D;
    int erosion_size=2,dilation_size=2;
    kernal_E= getStructuringElement( MORPH_RECT,Size( 2*erosion_size + 1, 2*erosion_size+1 ),
                                       Point( erosion_size, erosion_size ) );
    kernal_D= getStructuringElement( MORPH_RECT,Size( 2*dilation_size + 1, 2*dilation_size+1 ),
                                       Point( dilation_size, dilation_size ) );
    dilate(mEdge,mEdge,kernal_D);//Just to make sure all the Edges are closed
    erode(mEdge,mEdge,kernal_E);

    imshow("mCanny_Edge",mEdge);
    imwrite(FileName_S+ "mCanny_Edge.bmp",mEdge);

    vector<vector<Point> > contours;
    vector<Vec4i> hierarchy;
    /// Find contours
    findContours( mEdge.clone(), contours, hierarchy, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_SIMPLE, Point(0, 0) );
    /// Approximate contours to polygons + get bounding rects and circles
    vector<vector<Point> > contours_poly( contours.size() );
    for( size_t i = 0; i < contours.size(); i++ )
    {
        approxPolyDP( Mat(contours[i]), contours_poly[i], 1, true );
    }

    for (int i = 0; i < contours_poly.size(); i++)
    {
        if(contours_poly[i].size()>100)//Just a filter You can add more here
            drawContours(mSource_Bgr,contours_poly,i,Scalar(0,255,0),2);
    }
    imshow("mResult",mSource_Bgr);

    }


Can Kalman Filter be used to track Randomly Moving Target? i want to track random moving object with a camera using kalman filter...i have the following questions...

Randomly moving target means $Corelation(t) = E[ x(T)x(T+t) ]$ is very low...where $x(T)$ is the position of the target at time = $T$ along an axis say $x$ axis..eg...moving hand...moving bird..moving ant..what will be my state transition matrix $F(k)$ in $x(k+1) = F(k)x(k) + w(k)$ because i cannot describe it randomness...
is there any standard model to represent any random moving target?

any suggestion is welcome....i am aware that any object can be determined in a given image frame using its characteristics by image processing techniques and hence tracked....but is it  possible to use Kalman Filter?

It depends on what you mean by a "randomly moving object".  If you are trying to track something that truly moves around in a totally uncorrelated manner from sample to sample (like, say, a laser pointer that flickers on and off and randomly changes position in your camera images) then a linear tracker will not give you insight into the object's state.

The examples you cited (a moving hand, bird, etc.) might exhibit random behavior on "large" time scales, but on "small" time scales their state is highly correlated. For such a system the most important point is to have a high data sampling rate relative to the dynamics of the system.  The higher the sampling and measurement update rate, the more linear the underlying system dynamics will be from one measurement time to the next.  Also if you know the object of interest can be randomly changing direction or accelerating you can incorporate this prior knowledge by increasing the process noise during the prediction phase of your Kalman filter, giving less weight to the predicted position during the "update" stage.


Why this system is linear? Hi guys i'm studying signals and systems, and my professor told us that
$$y(t) = \int\limits_{  t+T }^{t-T/2}  {x(a+T/2)}\mathrm{d} a$$
is a linear system.
But a primitive of $x$ isn't $ x^2$ ? How it's possible that's linear ?

That's not an integral of variable $x$. The notation $x(a+T/2)$ stands for a function $x(\cdot)$ of variable $a$.
So applying the fundamental theorem of calculus, and assuming there exists a function $G(a)$ such that $G'(a) = x(a)$, then you will have :
$$ \int x(a+T/2) da = \int G'(a+T/2)da = G(a+T/2) + C $$
where the constant of integration, $C$, will be omitted in the definite integral :
$$ \int_{t+T}^{t-T/2} x(a+T/2) da = \int_{t+T}^{t-T/2} G'(a+T/2)da = G(a+T/2)|_{t+T}^{t-T/2}  $$
So the system has nothing with a square function.
Coming to its linearity, you can show this in line with the linearity of the integral operator...


Why we need Frequency and Time information at the same time The big disadvantage of a DFT is that it has only frequency resolution and no time resolution. This means that although we might be able to determine all the frequencies present in a signal, we do not know when they are present. So we go for Wavelets.

My question is, What are the applications (other than music & speech processing ) that require both the time and frequency information at same time  ?

How do you relate that we require both spatial (time) and frequency information in an image at a same time ?

Give me some good examples and suggest a book (PDF links) to learn Wavelet Multiresolution Analysis (For Dummies)

For the source, go to end of the answer

Suppose one day you got one note which has some thing written to it, say "Major frequency components are 10 Hz, 25Hz, 50 Hz and 100 Hz". Somehow, you understood that its time-series representation is a very important thing (may be master-piece work of a great musician, or some national security matter, anything). So you decided to find it.

You start by plotting the spectrum which looks like below:

(There can be some noises/disturbances in between these frequencies, but most people neglect it)

Now you want to find its time-series representation. Since we have worked it a lot in our signal processing class in college, most probably, we get below signal:

$$
x(t)=cos(2 \pi 10 t)+cos(2\pi 25 t)+cos(2 \pi 50 t)+cos(2 \pi 100 t)
$$

Great !!!

But is this the only time series representation for those frequencies?

Consider another time series signal which plays a sine wave of frequency 100 Hz for 300 ms, then 50 Hz sine wave for next 300 ms, then 25 Hz sine wave for 200 ms, then 10 Hz sine wave for 200 ms. A total of 1000 ms long signal. If you plot it, it will look like below:

Actually above signal also have similar frequency spectrum (shown below):

Those ripples are because of sudden changes in frequency. Amplitude changes are because of change in how long each frequency was played. They are not of importance to us right now.

So what to do now? How will you pick one? Even if we change the order they were played, still we get the same frequency spectrum. So how many such combinations you have?

  In short words, without knowing when those frequencies were played, it
  is very difficult to find its time series representation, and that is
  where, wavelets come into picture. It provides not only the
  frequencies present in the signal, but also, the time it was present
  in the signal.

UPDATE (17-Feb-14)

Question : How do you relate that we require both spatial (time) and frequency information in an image at a same time ?

Answer : Me too, had the same question before, so after seeing your comment, I thought I should find one example (I haven't studied wavelets in the context of images before wavelets). And below explanation is the summary of an answer in Quora.com (all credits to that author, Ron Reiter) which I felt as a good example.

You might be knowing that we use the DCT for jpeg compression of images. It preserves the low frequency data while throws away some of the high frequency data in the images, based on the image quality we need.Or in other words, this method is good with low frequency data, but not that accurate with high frequency data.

For example, in an image with sky and a tree, the areas of sky will be preserved accurately while some details of tree will be thrown off.

In DCT, we take 8x8 blocks of images and encode it for all frequency of data. For very low frequencies, we could have used higher block size while for very high frequencies, we may have to use lower block size than this.

Now apply our wavelet method as explained in 1D case above. What if we can find both frequency and their location in the image? We can encode the low frequency areas in big chunks while high frequency areas in low chunks. So it produces better results.

Please check the link for more details

SOURCE:

THE WAVELET TUTORIAL by ROBI POLIKAR

All the images, explanation, everything, are extracted from this link. It is a must-read introduction of wavelets to beginners

Amara's Wavelet Page
Wavelet in image compression - Ron Reiter at Quora.com


Multipath fading in wireless channel - parameter needed I’m trying to implement in Matlab the magnitude of the channel taps of the wireless communication system that is described at page 23 (Section 2.1.4) in the classic book of D. Tse, Fundamentals of Wireless Communications.

The $l^{\rm th}$ channel taps at time $m$ is given by eq. (2.51) p. 44:

$$
h_l[m] = \sum_i a_i\left(\frac mW\right)e^{-j2\pi f_c \tau_i\left(\frac mW\right)} \mathrm{sinc}\left[l-\tau_i\left(\frac mW\right)\right]
$$

In this particular case we have $i = 2$ paths and attetuation $a_i(t)$ for these paths are given from eq. (2.15):

$$
a_1(t) = \frac{\lvert\alpha\rvert}{r_0 + vt}, \quad a_2(t) = \frac{\lvert\alpha\rvert}{2d - r_0 - vt}
$$

And propagation delays by eq.(2.16):

$$
\tau_1(t) = \frac{r_0 + vt}{c} - \frac{\phi_1}{2\pi f }, \quad \tau_2(t) = \frac{ 2d -r_0 - vt}{c} - \frac{\phi_2}{2\pi f}
$$

I have calculated everything apart from $\lvert a\rvert $. I don’t really know how to calculate it.

Any suggestions?

You don't need to compute the parameter $|\alpha|$, because it is the antenna gain, which is part of the antenna's specification. This is also mentioned on p.21 of the book you quoted (above Eq. $(6)$).


Alamouti Scheme. MISO. STBC for real constellation I have a real orthogonal design for codeword...

$$
\begin{bmatrix}
    s1 & s2 \\
    -s2 & s1
\end{bmatrix}
$$

$$
s=\begin{bmatrix}
s1\\
s2
\end{bmatrix}
$$

$$
Y^T=\sqrt{\frac{P}{Mt}} h^T X + n^T
$$

...and receive signal:

$$
Y1=\begin{bmatrix}
y1\\
y2
\end{bmatrix}
$$

So for:

$$
X = \begin{bmatrix}
    x1 & x2 \\
    -x2 & x1
\end{bmatrix}
$$

and
$$
Y = \begin{bmatrix}
    y1 & y2
\end{bmatrix}
$$

...i have:

$$
Y=\sqrt{\frac{P}{Mt}} \cdot \begin{bmatrix}h1 & h2\\ h2 & -h1\end{bmatrix}\cdot \begin{bmatrix}x1 & x2\end{bmatrix}^T+N
$$

...right?

Your analysis is correct. I would write it as follows: if you define the transmitted matrix as
\begin{equation*}
X = \begin{bmatrix} x_1 & -x_2 \\ x_2 & x_1
\end{bmatrix}
\end{equation*}
with one column transmitted at each symbol interval (that is, first you transmit $[x_1, x_2]$ and then $[-x_2, x_1]$), then the system can be written as $$Y=hX$$ with $h = [h_1, h_2]$, or equivalently it can be written as $$Y = Hx$$ with
\begin{equation*}
H = \begin{bmatrix} h_1 & h_2 \\ -h_2 & h_1
\end{bmatrix}
\end{equation*}
and $x=[x_1, x_2]^T$.

Of course, you can also transpose all vectors and matrices and arrive at equivalent results.


QPSK demodulation problem due to apparent transmitter clock frequency drift I am using FPGA to generate 10Mbaud QPSK signal in two channels corresponding to two quadratures. Upsampling to 1GSample/s (using raised cos filter for interpolation) and complex upconversion (f_up=250MHz) is done in FPGA. This signal is sent to a DAC and detected back-to-back using Tektronix rsa607a spectrum analyzer. I sent 200k sampes and detected them using the analyzer.

This is a screenshot form the measurement.  On the right side is unzoomed waveform in time where red dots represent sampling points in demodulation. On the left is constellation for the whole stream. As you can see the analyzer synchronizes well with the middle of the captured signal, however it goes off sync with time. If I measure shorter intervals EVM gets smaller. On the left side you can see measured frequency error which is the difference between transmitter frequency and analyzer's internal clock. This error is obviously different with time which tells me the FPGA clock might be unstable. However, I did other measurements when offloading a sine tone where I don't see these frequency drifts.
Does anybody have intuition what can be the cause of this problem?

Cheers,
D

The problem was that the FPGA clock was significantly different than the clock from all other devices. It was not calibrated properly in regard to some universal frequency reference.


Study resource for Digital communication I am in my 5th semester of computer science and we have a digital communication course. All the textbooks of digital comm. available in my library demand a prerequisite of Circuit analysis and Signals and systems.
Is there any book which is self contained or at least less heavy on these background topics?

I highly recommend you the book Schaum's outlines of Analog and Digital Communications.

It will be the fastest refresher for the signal processing and probability concepts which are essential prerequisites for digital communications. It's not using that heavy math also, nevertheless, not a trivial one either ;-)


Pre-envelope of ${\Pi}_{a}(t)\cos(2{\pi}f_{0}t)$ I want to find the pre envelope of

$$x(t) = {\Pi}_{a}(t)\cos(2{\pi}f_{0}t)$$
where

I found the Fourier transform to be
$$  \begin{align} X(f) & \triangleq \mathfrak{F}[{\Pi}_{a}(t)\cos(2{\pi}f_{0}t)] \\ & = a\text{ sinc}(2a(f-f_{0})) + a\text{ sinc}(2a(f+f_{0})) \\ \end{align} $$
So $$X_{+}(f) = 2u(f)X(f)$$
But I don't know how to continue from here.
Is there an easier way to find the pre envelope with another method?

The pre-envelope is also called analytic signal. Its Fourier transform is indeed given by the expression in your question:

$$X_+(f)=2X(f)u(f)\tag{1}$$

where $X(f)$ is the Fourier transform of the original signal, and $u(f)$ is the unit step function. Obviously, $X_+(f)$ has only positive frequency components. The analytic signal $x_+(t)$ with its Fourier transform given by $(1)$ is necessarily a complex-valued signal:

$$x_+(t)=x(t)+j\mathcal{H}\{x(t)\}\tag{2}$$

where $\mathcal{H}\{\cdot\}$ denotes the Hilbert transform. Note that the analytic signal is not the same as the complex envelope. For a band pass signal $x(t)$, the complex envelope is a low pass signal, whereas the analytic signal is a band pass signal.

From $(2)$, in order to compute the analytic signal you need to compute the Hilbert transform of $x(t)$:

$$\hat{x}(t)=\mathcal{H}\{x(t)\}=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau}d\tau\tag{3}$$

For the given signal you get

$$\hat{x}(t)=\frac{1}{\pi}\int_{-a}^a\frac{\cos(2\pi f_0\tau)}{t-\tau}d\tau\tag{4}$$

The result of $(4)$ can be written in terms of the cosine integral $\text{Ci}(x)$ and the sine integral $\text{Si}(x)$:

$$\hat{x}(t)=\frac{1}{\pi}\left\{\cos(2\pi f_0t)\left[\text{Ci}(2\pi f_0(t+a))-\text{Ci}(2\pi f_0(t-a))\right]+\\\sin(2\pi f_0)\left[\text{Si}(2\pi f_0(t+a))-\text{Si}(2\pi f_0(t-a))\right]\right\}\tag{5}$$

I don't think that anybody expected you to come up with that awful expression. Anyway, for large values of $f_0$, the first term in $(5)$ becomes very small, and the second term converges to $\sin(2\pi f_0)\Pi_a(t)$. So for large $f_0$ you get the expected result

$$x_+(t)\approx x(t)+j\sin(2\pi f_0)\Pi_a(t)=\Pi_a(t)e^{j2\pi f_0 t}\tag{6}$$

Obviously, the respresentation

$$x(t)=\text{Re}\left\{\Pi_a(t)e^{j2\pi f_0 t}\right\}\tag{7}$$

is always valid, but the complex-valued signal $\Pi_a(t)e^{j2\pi f_0 t}$ is no analytic signal, it's just a good approximation of the analytic signal for large values of $f_0$.

Note that for $x(t)=m(t)\cos(2\pi f_0t)$ with $m(t)$ a band-limited function, i.e. $M(f)=0$ for $|f|>B$, the complex-valued signal $m(t)e^{j2\pi f_0 t}$ is an analytic signal, as long as $f_0>B$. The problem with the function given in your question is that $\Pi_a(t)$ is not band-limited.


Constraints for Type 3 Linear Phase Filter This might be a strange question - but most of the Type 3 LPF I'm seeing are having this in common for $h[n]$, assuming $h[n]$ is real:

We cannot have outermost element of h[n] lower than middle elements
(or) Have middle element lower than outermost elements.

Eg: if filter coefficients are $[3,2,1,0,-1,-2,-3]$, it there any restriction for a type $3$ filter to not have $[1,2,3,0,-3,-2,-1]$ as a possible LPF $h[n]$?

(Please excuse my ignorance if this is a basic question - I'm quite new FIR Filter design). And I've gone through quite a few design docs before asking this question - none of them mention this, so I'm guessing this is not a req? Yet the diagrams seem to indicate this, hence the question...

There is no such restriction for a type III linear phase filter. Neither is there such a restriction for any other type of linear phase filter. Type III linear phase filters have an odd filter length $N$ and odd symmetry. I.e., as long as the filter's impulse response $h[n]$ satisfies

$$h\left[\frac{N-1}{2}+n\right]=-h\left[\frac{N-1}{2}-n\right]$$

you have a proper type III linear phase filter.

If you look at the impulse responses of type III differentiators or Hilbert transformers then you'll notice that in both cases the filter taps are decreasing in magnitude as you move away from the center tap. So these common applications of type III filters are counterexamples to your assumption.


Basis signals for the discrete-time Fourier Series I am using Alan Oppenheim's Signals and Systems and I am a bit confused by the notion of discrete-time periodic exponentials as basis signals for the discrete-time Fourier Series (and later the DTFT). For the discrete-time Fourier Series we use the basis signals $$\phi_k [n] = e^{jk\frac{2\pi}{N}n} = e^{jk\Omega_0n} ,$$ where $\Omega_0 = \frac{2\pi}{N}$, to build general discrete-time periodic signals. Unlike the continuous-time case where we use complex exponentials as basis functions for the Fourier  Series representation, in this case the fundamental frequency of the $k$th basis signal is not $k\Omega_0$. This is because it can be shown that for a discrete-time periodic complex exponential signal $e^{jk\frac{2\pi}{N}n}$, the fundamental period $N_0$ is not $N$, but rather $$N_0 = \frac{N}{\gcd(k,N)},$$ where $\gcd(k,N)$ is the greatest common divisor of the integers $k$ and $N$. From that we conclude that the fundamental frequency of the $k$th basis signal $\phi_k[n]$ is not $k\frac{2\pi}{N}$, but rather $\gcd(k,N) \frac{2\pi}{N}$. So it is still an integral multiple of the frequency $\Omega_0$, except $\Omega_0$ is multiplied by the greatest common divisor of $k$ and $N$, rather than $k$. Of course in the case where $k$ divides $N$ the frequency would be $k \Omega_0$, but we cannot know a priori whether that's the case. Am I correct in my understanding?

Assuming that what I said above is true, when we look at the discrete-time Fourier Series representation of the sequence $x[n]$, we have $$x[n] = \sum_{k = <N>} a_k \phi_k[n],$$

which is a linear combination of the basis signals I defined above (and what the author means by the notation above is that the sum is over any $N$ values of $k$). Now given the fact the fundamental frequency of the $k$th component in this sum is $\gcd(k,N) \Omega_0$, it might turn out that $k$ and $N$ have no common factors at all, in which case the fundamental frequency of the $k$th basis signal is simply $\Omega_0$, which is the same as the fundamental frequency of the first basis signal. Therefore, those two basis signals seemingly have the same fundamental frequency and period (within the range of values of $k$ defined earlier). Does that make those signals identical, or are they different? Are the signals guaranteed to be different? And if they are not, doesn't that make the set of basis signals $\{\phi_k\}$ linearly dependent? And does the Fourier Series representation break down in this case?

Thank you in advance

The complex exponentials

$$\phi_k[n]=e^{j\frac{2\pi k}{N}n},\qquad k=0,1,\ldots, N-1\tag{1}$$

form an orthogonal basis, i.e.

$$\sum_{n=0}^{N-1}\phi_k[n]\phi_l^*[n]=0,\qquad k\neq l\tag{2}$$

So clearly they are all different, and they are even orthogonal (and hence not linearly dependent).

It's straightforward to show that Eq. $(2)$ is true:

$$\sum_{n=0}^{N-1}\phi_k[n]\phi_l^*[n]=\sum_{n=0}^{N-1}e^{j\frac{2\pi (k-l)}{N}n}=
\begin{cases}
\sum_{n=0}^{N-1}e^{j\frac{2\pi (0)}{N}n}=\sum_{n=0}^{N-1}1=N,\qquad & k = l \\
\\
\frac{1-e^{j2\pi (k-l)}}{1-e^{j\frac{2\pi (k-l)}{N}}}=\frac{1-1}{1-e^{j\frac{2\pi (k-l)}{N}}}=0,\qquad & k\neq l \\
\end{cases}\tag{3}$$


Filter pixels with respect to specific wavelength frequencies I'm not sure if this is possible, but let's say you had access to the original raw sampled digital pixels from a sensor, captured with nothing between the sensor and the lenses, i.e.  Anti aliasing filter, colour filter array etc etc removed. Would it be possible to filter the raw pixels in relation to the original wavelength frequencies?

Thanks

Aliasing irreversibly transforms frequencies above the Nyquist limit to appear as frequencies below the limit. Almost no mathematical algorithms are effective compared with a physical anti aliasing filter in front of the sensor.

However, if the sampling frequency of the sensor itself is higher than the numerical aperture of the lens, it may not need the low pass anti aliasing filter any more since lens itself would serve as a low pass filter.

Regarding the color filter, I think it depends on your application. If you are try to image fluorescence, it is possible without a color filter if the emission spectrum is distinguishable with excitation spectrum. Background subtraction operations will recover the pixel information on the fluorescence. Even if several emission spectra are multiplexed, there is still possibility to differentiate them by solving a linear matrix equation under the assumption that the sensor system is a linear system with different wavelengths.

But color correction filters in digital photography may still needed to get more uniform exposure in the desired channels. Otherwise, the denoising process, though maybe possible through some image processing algorithms, is still challenging in complete noise removal without the prior knowledge of the noise distribution.


Under what conditions is there a one-to-one mapping between continuous-time and discrete-time signals? As the sampling theorem dictates that the uniform sampling frequency must be at least twice the maximum frequency present in the bandlimited signal (Nyquist rate), a question arises about the uniqueness of this representation of the analog signal to digital one.
Apart from practical errors related to the truncation of the infinitely-long energy signal when discretized to digital, can we find a UNIQUE mapping between the two signal forms? If yes, then what are the conditions that limit or control this unique representation?
Note: I am not assuming by default that the sampling topology should be uniform. Random sampling also could be considered as well.

Apart from practical errors related to the truncation of the infinitely-long energy signal when discretized to digital, can we find a UNIQUE mapping between the two signal forms?

Yes.
Say you're starting with $x(t)$, which is perfectly bandlimited to less than $f_B = \frac{1}{2 T_s}$, and you take samples at every $t = T_s k$: $x_k = x(T_s k)$.
All you need to do is to treat the  as a train of Dirac impulses (i.e., let $x_p(t) = \sum x_k \delta(t - T_s k)$), and then filter $x_p(t)$ with a brick-wall low-pass filter.  I.e., let $x_f(t) = h(x_p(t))$, where $h$ is that brick-wall filter.  $h$ will have an infinitely long impulse response, so you cannot realize this system in the real world - but not only will the mapping from $x_k$ to $x_f$ be exact, but if $h$ has zero phase delay then you'll find that $x_f(t) = x(t)$.


Why would I leave a signal oversampled? I can't think of a better way for asking this question so I will start with an example. Suppose that I have an input signal with a max frequency of 50Hz (sampled at 100Hz). Now the signals of interest lie in the range 0-5Hz, so I can add a low-pass filter with a cut-off of 5Hz, and use the resulting signal for further processing. My understanding is that now I can downsample the filtered signal by a factor of 10 and hence reduce processing load. Am I right?
If yes, why is downsampling not ALWAYS performed after filtering because it seems to me as the obvious way to go?
And if I am wrong in my assumption, where am I mistaken?

You are correct that if your signal is bandlimited to <5 Hz, then you can perfectly represent it with a 10Hz sampling rate.  This is the well-known sampling theorem

But ... there may be practical considerations for why one would not be able and/or inclined to use critically sampled data.

One reason is the difficulty of making a signal critically sampled. Any operation you perform to change the rate of the signal is going to have some filter with a non-zero transition bandwidth.  In your example, this limits the unaliased frequency content to 5-ftrans  This transition bandwidth can be made very narrow with long impulse response filters but this has costs both in terms of processing and in transients (ringing) at signal start and end.

Another reason is the efficacy of algorithms that work on the resulting signal. If you need to work with a blackbox component that can only choose the nearest sample, then you'll be better off feeding it oversampled data.

Most (all?) non-linear operations will behave differently with critically sampled vs oversampled data.  One example is squaring a signal, a well known method of BPSK carrier recovery. Without a 2x oversampled condition, the multiplication of the time domain signal with itself causes wraparound garbage aliasing when the frequency domain convolves with itself.


Lifting scheme versus filter banks I am having some trouble understanding the intuitive connection between the two. I worked through all the theory and learned that Sweldens showed how to derive lifting steps when having admissible filter taps (Daubechies, etc.).
The opposite is simple mathematically, i.e. deriving the filter taps when having discrete lifting steps. However, I dont get how one would come up with 'admissible' lifting steps in the first place, if not by getting them from valid filter taps like Sweldens showed.

Traditional (analysis) filter banks are banks of linear filters followed by downsampling. This  is not so efficient, because you convolve and then subsample: possibly a waste of operations.  People have found ways to invert the operations, downsample and convolve, with other filters involve the polyphase components (even and odd ones in the two-channel case) of the signal (polyphase matrix). Such matrix can be factored into a product simpler matrices, like diagonal, upper and lower diagonal ones, acting more separately on the polyphase components.
A reinterpretation consists in using the upper matrix on some polyphase component (say, odds) to predict and combine outcome to other components (say, even) to get what is called an update, possibly with scaling factors.
Properties of polyphase matrices tell that this is always possible, yet not unique. And one can save some operations, compute inverses more easily (upper and lower matrices have simpler inverses). So all invertible linear filter banks can be converted in lifting.
But more importantly, one can insert non-linear blocks and still be invertible and efficient. So the lifting is more generic that "linear" filter banks (as long as they have finite support). An example of the above is known as the S-transform: given a sequence of integers  $c[n]$, you can decompose it into ($\lfloor \cdot\rfloor$ denotes rounding down):

$l[n] =\Large\lfloor \frac{c[2n]+c[2n+1]}{2}\Large\rfloor$
$h[n] = c[2n]-c[2n+1]$

How to design a lifting scheme from scratch requires to:

understand how to split the data into several subsets
design proper predict and update operators
choose how to iterate them on different levels.

I do like a couple of papers from people at KU Leuven:

Integer wavelet transfroms using the lifting scheme, Geert Uytterhoeven, Dirk Roose, Adhemar Bultheel
Wavelet Transforms Using the Lifting Scheme, Uytterhoeven, Geert  and Roose, Dirk  and Bultheel, Adhemar
WAILI: Wavelets with integer lifting, G. Uytterhoeven, F. Van Wulpen, M. Jansen, D. Roose and A. Bultheel


Are there any real world applications for complex-valued signals or impulse responses? I was just curious...
$$x[n] {\longrightarrow}  \boxed{h[n]} {\longrightarrow} y[n]$$
I've never seen a real world filter where the coefficient of $h[n]$ were complex, or where $x[n]$ was a complex sequence.
But, the DSP book always makes a big deal about conjugate symmetric sequences, which implies that either $x[n]$ or $h[n]$ is composed of complex numbers.
What would the real world use case be for the following three scenarios:

complex $x[n]$, real $h[n]$
real    $x[n]$, complex $h[n]$
complex $x[n]$, complex $h[n]$

Do complex numbers exist in real world implementations?

Absolutely! Conjugates are mentioned in textbooks because conjugation has no effect on real signals, but it does on complex ones. This way, formulations are more general and apply to both real and complex valued signals. Complex numbers don't exist themselves, they are a mathematical construct.
Having said that, their mathematical properties can be replicated using real systems. You can separate the real and imaginary parts and treat them individually as real signals, but you must use additional hardware in order to do so. This manifests as needing more wires to handle both components, as well as additional memory to store complex values.
This is especially straight forward in the digital domain. However, I'm going to use continuous time signals to avoid introducing sample rates and is cleaner to present.
1. Complex $x(t)$, Real $h(t)$  - Moving Average Filter
Let's say we have a complex signal which is noisy and you want to smooth it out. One way to do this would be to employ a moving average filter on the signal. This requires that the filter be applied separately to both real and imaginary parts. In this example the input signal $x(t) $ is an arbitrary complex triangular signal with noise added. The moving average filter is given by
$$h(t) = \frac{1}{L}$$
So the output is
$$y(t) = x(t)*h(t) =  \frac{1}{L}\int_{-\infty}^{\infty}x(t - \tau)d\tau$$
The plot below shows the filter smoothing out the complex input signal.

We have processed a complex signal with a real-valued system. This is a very specific example as there are many types of real-valued systems that operate on complex (quadrature) inputs.
2. Real $x(t)$, Complex $h(t)$  - Bandpass Filter Design Using Lowpass
Using the frequency shift property of the fourier transform, you can yield a bandpass filter design given a lowpass filter $h(t)_{LP}$. Using this property, we can move the lowpass filter to be centered around a desired frequency $f_0$ and is given by
$$h(t)_{BP} = h(t)_{LP} \space e^{j2{\pi}f_0t}$$
Doing this makes $h(t)_{BP}$ complex and can be used to filter a signal.
Let's say we have an input signal that contains frequency components at $f_0 = 200 kHz$ and $2f_0 = 400 kHz$ but we only want $f_0$. We can start with an appropriate lowpass filter $h(t)_{LP}$ and apply the frequency shift to yield the new filter and process the signal.
$$x(t) = cos(2{\pi}f_0t) + cos(2{\pi}(2f_0)t) $$
$$h(t) = h(t)_{LP} \space e^{j2{\pi}f_0t} $$
Below we can see the initial lowpass filter and the bandpass filter we design using frequency shifting. The new filter performs as intended and we are left with the sinusoid at $f_0$.

We have processed a real signal with a complex-valued system.
3. Complex $x(t)$, Complex $h(t)$  - Radar LFM Pulse Compression
In pulse-Doppler radar systems, a technique called pulse compression is employed to achieve both good pulse widths (better energy on target) while maintaining fine range resolution. This is usually done via matched filters to achieve the highest SNR possible for a given target return.
A popular modulation scheme is linear-frequency modulation (LFM). The complex LFM signal transmitted with a chirp bandwidth $\beta$ and pulse width $\tau$ is
$$s(t) = e^{j{\pi}\frac{\beta}{\tau}t^2}$$
For the following example, we'll be using a bandwidth of 10 MHz and a pulse width of 10 $\mu$s. Below shows the real and imaginary parts of the LFM pulse.

The matched filter for this waveform is given by
$$h(t) = s(-t)^* = e^{-j{\pi}\frac{\beta}{\tau}t^2}$$
The return signal from a target arrives at a delay of $t_d$, so the signal we will process with the matched filter is
$$x(t) = s(t-t_d)$$
Convolution with a matched filter yields the cross-correlation output as
$$y(\tau) = x(t) * h(t) = \int_{-\infty}^{\infty}x(t)h(t+\tau)dt$$
Using the delay $\tau$, we can determine the range of the target since we know our pulse travels at the speed of light. Below we see the output of the matched filter for the nominal zero-delay case and for a target at 300 m.

We have processed a complex signal with a complex-valued system.


Contour Integral and Residue Theory for Inverse $z$-Transform I'm aware that the inverse $z$-transform can be evaluated using contour integration which leads to the use of Residue Theory as a corollary and I do know of the two definitions. My question is how does one lead to the other being used? It's been many years since I studied the two.

Thanks.

It's a consequence of Cauchy's Residue Theorem, which can be used to evaluate contour integrals of analytic functions. Take a look at this answer for an example.


Circular and Linear Convolution What is the difference between circular and linear convolution? When would I choose one over the other? In image processing where a filter is applied to an image with a mask which type of convolution should I choose?

If you have a vector of data, $d$, that is composed of elements $d_1, d_2,... d_N$, then linear convolution operates on them in order, starting with $d_1$ and ending with $d_N$.

Imagine that the data vector $d$ is represented by a slip of paper with the $N$ elements written in order.  Now, imagine forming the slip of paper into a circle by touching the end (where $d_N$ is written) to the beginning (where $d_1$ is written).  Convolving that is circular convolution.  In practice linear convolution and circular convolution are nearly the same, the difference happening at the beginning and the end of linear convolution.  In linear convolution you assume that there are zero's before and after your data (i.e. we assume that "$d_0$" and "$d_{N+1}$" are 0), while with circular convolution we wrap the data to make it periodic (i.e. "$d_0$" is equal to $d_N$ and "$d_{N+1}$" is equal to $d_1$).

The same principles hold for multi-dimensional arrays.  For linear convolution there is a definite start and end for each axis, with zeros assumed before and after.  For circular convolution the data wraps around in each axis.

When would I choose one over the other?

With a few very rare exceptions we don't "choose" circular convolution.  We almost always want linear convolution.  The reason that circular convolutions pop up as much as they do is because convolutions via FFT's (FFT, multiply, inverse FFT) are circular convolutions, not linear.


Apply AWGN noise to QPSK-OFDM symbol I'm working on an OFDM-QPSK based system and I'm having some trouble to match the BER of my system with the theoretical curve. My OFDM implementation fills all the subcarriers and have t_symbol s duration and t_cp cyclic prefix duration. I'm trying this to validate my system implementation before go on with my system development. I'm using the following python code for apply noise:

def awgn(self,data,noise_power):
    Es = sum(abs(data)**2.)
    Eb = Es/(self.OFDM_SIZE*2.)
    EbN0_dB = noise_power + 10*log10(self.t_symbol/(self.t_symbol - self.t_cp))
    EbN0 = (10.)**(EbN0_dB/10.)
    N0 = Eb/EbN0
    noise = sqrt(N0/2.)*(random.standard_normal((data.size,)) + random.standard_normal((data.size,))*1j)
    return (data + noise)

I just realize that the normalization done here is not correct, since this signal is post OFDM modulation. But the problem remains.

The equation for the theoretical curve used is:

BER_t   = 0.5*erfc(n_sqrt((10.**(EbN0_vt/10.))/2.))

I had used similar code before on Matlab and it worked fine, but now It's not.

I guess it's not working now because my system don't have the constellation with (1+1j),(1-1j)... I have to use values that are less than one, since my system is implemented in fixed point. Now I'm using 0.75 as an amplitude but it's arbitrary chosen. I tried to correct the energy symbol to get a mean energy equal to one but using this take my signal out of range. I'm thinking about correct the EbN0 before apply on the theoretical BER equation,  but I'm just wondering how.

Can anyone give me the correct way to go?

My central question

Trying to improve my question:
What I need to know is how to apply awgn on ofdm signal when it's subcarrier are modulated with QPSK and how to match it with the theoretical BER plot. My QPSK symbol have an amplitude of 0.75. The above code is to show how I'm doing it now.

update

I run the code changing the FFT subsystem to the scipy.fftpack and the amplitude of the QPSK symbol to 1 instead of 0.75. And the curve match. I also had corrected the theoretical curve to:

BER_t   = 0.5*erfc(n_sqrt((10.**(EbN0_vec/10.))))

So I still think it's a matter of correct the EbN0 in the above equation.

update2

After several testing of the above functions I found thet they are correct. The source of the problem is that I'm dealing with fixed point implementation of IFFT/FFT and I have to handle the channel model with this situation in mind. I'm trying to figure out how to handle this but without a solution yet.

Your measurement of the symbol and noise energy is not correct.  Please note that the following is pseudo-code, not valid Python code.  EbN0_dB should be an input to this code, not an output.

Esymbol = sum(real(data*conj(data))) / numSymbolsInData
Eb = Es/numPayloadBitsPerSymbol
Eb_dB = 10*log10(Eb)
N0_dB = Eb_dB - EbN0_dB
N0 = (10.)**(N0_dB/10.)


OFDM Short Training Field clarification I want to clarify my understanding of OFDM training fields. Please correct me if I am wrong as I am getting 2 different information from online resources.

For 802.11a standard the preamble contains 10 short training Field (STF) and 2 Long Traning Field (LTF). Each of the STF (A1,...A10) contains 52 subcarriers that contain symbols as shown in the picture.

And In some other document, I am seeing that there is 2 short training field instead of 10, the 2 LTF is there as same. Which one is true, or both are true in a different perspective? I have also seen in many implementation 2 STF are used.

In addition:

a) Is the same packet detection method used for both 802.11a and 802.11n? (i.e with STF)?

b) The standard says each of the STF symbols has 16 samples. Can anyone explain this? How this 16 sample is achieved or comes from?

Typically both are used, coming from an 802.11p perspective. The short preamble (STF) is used for packet detection and coarse frequency offset. The long training preamble (LTF) is used for timing synchronization and fine frequency offset. The LTF can also be used for channel equalization. Both are used from my understanding to accomplish total synchronization.


OFDM Phase compensation I need to compensate a OFDM-16_qam constelation rotation but I could not find out how to implement it with MATLAB
I have a Matrix [48 x 25] where each colum is a OFDMsymbol without CP , pad and pilots, this matrix was passed through a Multipath channel using the matlab command Ch=rayleigh and filter(CH,OFDMsymbol) and so the pilots modulated in bpsk .
With this pilots I've constructed a Pilot_matriz with each colum (an OFDM symboll) has an angle to be compensate.
The OFDMsymboll  is a complex number as  +/-a+/-ib and I have the angles (theta) in rad from the pilot matrix but I'm not find a way to turn bach the theta angle from the complex number +/-a+/-ib.
any suggestions
Regards
Luiz

It's not completely clear from your question, but I interpret your situation as:

You have a set of OFDM symbols that have been passed through a multipath channel.
You have generated an equalization matrix of phase shifts (in radians) that you want to apply to each subcarrier of each symbol.
You're not sure how to apply the phase shifts to the OFDM symbols.

Being able to implement an equalizer via simple multiplication is one of the benefits of OFDM, since one can assume that fading is flat across each subcarrier's channel. To apply a phase shift of $\theta$ radians to a complex number $x = a + jb$, do the following multiplication:

$$
\tilde x = e^{j\theta}x
$$

$$
\begin{align*}
\tilde x &= \left(\cos(\theta) + j \sin(\theta)\right) \left (a + jb\right) \\
&= \left(a \cos(\theta) - b \sin(\theta)\right) + j \left(a \sin(\theta) + b \cos(\theta) \right)
\end{align*}
$$