How to avoid leakage in the following FFT I have written a code for FFT for a signal which time domain plot and FFT is shown bellow.
I am not sure about what FFT length to do in order to avoid leakage of energy?
I know the number periods has to be integer, but i cant see how to implement this principle in here practically.
Also I get an amplitude which is very far from 5 5 10 as written in the time domain signal.
Where did i go wrong?

    clc
    clear all
    Fs=200e3;
    Ts=1/Fs;
    dt=0:Ts:5e-3-Ts;

    f1=1e3;
    f2=20e3;
    f3=30e3;

    y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+5*sin(2*pi*f1*dt)+10*sin(2*pi*f3*dt);
    %plot(dt,y)
    nfft=length(y);
    nfft2=2.^nextpow2(nfft);
    fy=fft(y,nfft2);
    f_half=fy(1:(nfft2/2));
    xfft=Fs.*(0:nfft2/2-1)/nfft2;
    plot(xfft,abs(f_half));

fixed final:
clc
clear all
Fs=200e3;
Ts=1/Fs;
dt=0:Ts:5e-3-Ts;

f1=1e3;
f2=20e3;
f3=30e3;

y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+10*sin(2*pi*f3*dt);
%plot(dt,y)
nfft=length(y);
nfft2=1000;
fy=fft(y,nfft2);
f_half=2*fy(1:(nfft2/2));
xfft=Fs.*(0:(nfft2/2)-1)/nfft2;
plot(xfft,abs(f_half)/(nfft));

Updated FFT conv FIR
clc
clear all
Fs=200e3;
Ts=1/Fs;
dt=0:Ts:5e-3-Ts;

f1=1e3;
f2=20e3;
f3=30e3;

y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+10*sin(2*pi*f3*dt);
%plot(dt,y)
nfft=length(y);
nfft2=1000;
fy=fft(y,nfft2);
f_half=2*fy(1:(nfft2/2));
xfft=Fs.*(0:(nfft2/2)-1)/nfft2;
%plot(xfft,abs(f_half)/(nfft));

cut_off=25e3/(Fs/2);
order=64;
h=fir1(order,cut_off);
con=conv(y,h)
%plot(con)
fz=fft(con,nfft2);
f_half_z=2*fz(1:(nfft2/2));
plot(xfft,abs(f_half_z)/(nfft));

Updated with filter response plot:
clc
clear all
Fs=200e3;
Ts=1/Fs;
dt=0:Ts:5e-3-Ts;

f1=1e3;
f2=20e3;
f3=30e3;

y=5*sin(2*pi*f1*dt)+5*sin(2*pi*f2*dt)+10*sin(2*pi*f3*dt);
%plot(dt,y)
nfft=length(y);
nfft2=1000;
fy=fft(y,nfft2);
f_half=2*fy(1:(nfft2/2));
xfft=Fs.*(0:(nfft2/2)-1)/nfft2;
%plot(xfft,abs(f_half)/(nfft));

cut_off=25e3/(Fs/2);
order=32;
h=fir1(order,cut_off);
fh=fft(h,nfft2);
plot(abs(fh(1:nfft/2)));
con=conv(y,h)
%plot(con)
fz=fft(con,nfft2);
f_half_z=2*fz(1:(nfft2/2));
%plot(xfft,abs(f_half_z)/(nfft));
```

Since your frequencies f2 and f3 are a multiple of f1=1k, and since your time vector is until 5 ms, that means your signal fits for a whole 5 periods, harmonics and all, so it's just ripe for a FFT. The problem you're having is that you chose the FFT to be of a number of points (1024) different than the number of samples (1000), so now one bin will be Fs/nfft = 2e5/1024 = 195.31 Hz. If you check the first values of the frequency vector xfft(1:10), you'll see that the first frequency falls in between two bins:
            0
   1.9531e+02
   3.9062e+02
   5.8594e+02
   7.8125e+02
   9.7656e+02 % f1 is between this
   1.1719e+03 % ...and this
   1.3672e+03
   1.5625e+03
   1.7578e+03

The other two have the same symptom. This means the amplitude (which is not normalized, that's why it's ~5k) can't be the real one. When the equivalent DFT will calculate the butterflies, the indices will not match f1 at the same spot. If you would have performed the FFT for 1000 points, the bins would have been 2e5/1e3 = 200 Hz in width. That would have made the frequency vector display these frst 10 values:
      0
    200
    400
    600
    800
   1000 % spot on
   1200
   1400
   1600
   1800

Using a 1024 length FFT means now there are more than one butterfly with different values, and they all add up. Which means you have a mild case of spectral leakage, i.e. the frequency that should have been present in one bin, only, is spreading across its neighbours. Also see the Wikipedia article.


Nonlinear wavelets transform? Is wavelet a Nonlinear transform, or Not?
specifically, continuous wavelet transform with morlet function.
I am studying behavior of a dynamic system, and it has nonlinear behaviour. can I employ wavelet transform?

A transform being linear has very little to do with its ability to analyze linear or nonlinear systems.

The wavelet transform $W[s(t)]$ of a signal $s(t)$ is linear because $$W[a s_1(t) + b s_2(t)]=a W[s_1(t)]+b W[s_2(t)]$$ for real or complex $a$ and $b$.

The signal you're analyzing is just a signal, it has no concept of linearity. However, if you try to come to conclusions about system properties of a nonlinear system, then you cannot break the analysis down to just a set of base signals to understand the system. In the worst case you would have to look at every possible intput/output pair. Often this can be simplified using known system properties like symmetries (i.e. time invariance).


Active noise cancellation using kalman filter I am doing signal processing on audio data sampled at 8Ksps in matlab but it is corrupted with random noise.  Therefore,  I decided to use LMS and RLS ANC algorithms to remove overlapped frequency noises and I have found RLS performance was better than LMS and NLMS. But now I want to use kalman filter so that I can achieve better result. I know kalman filter and have used in many basic applications but don't know how to use it for ANC. I would be grateful if anyone give me some idea .

I guess you use the LMS and RLS filters because they are adaptive, and you can do the same thing with the Kalman filter, if you augment the state-vector with the parameters of the model. You then get a non-linear system, and you can solve that using any number of non-linear Kalman-filter implementations, the most common being the extended Kalman filter (EKF), then second most common the unscented Kalman filter (UKF), and then there are a fair few more exotic variants, such as the second order Kalman filter. The convergence properties of any of these filters have not been proved, but they seem to work really well in practice. The UKF is typically regarded as better than the older EKF, since it does not require the inversion of the Jacobian, and it behaves as a higher order approximation -- 3rd order for Gaussian noise I seem to remember (compared to first order for the EKF).

The EKF is perhaps a bit more straight forward to get started with. You use the EKF for parameter adaptation by modeling unknown parameters as Wiener processes. Consider the linear system
\begin{align*}
 \dot{x} &= A(\theta) x + B(\theta) u + w \\
 y &= C(\theta) x + v
\end{align*}
with unknown parameters $\theta \in \mathcal{R}^{p}$. By augmenting the state vector to include the unknown parameters, $\chi^{\text{T}} = [x^{\text{T}}, \theta^{\text{T}}]$, we obtain the non-linear system $\dot{\chi} = f(\chi,u) + w$, $y = h(\chi) + v$, where
\begin{align*}
 f(\chi,u) &=
 \left[ \begin{array}{c}
  A(\theta) x + B(\theta) u \\
  0
 \end{array} \right] \\
 h(\chi) &= C(\theta) x .
\end{align*}
The Jacobians $F$ and $H$ for this system are found as
\begin{align*}
 F &=
 \left[ \begin{array}{cc}
  A(\theta) & \frac{\partial}{\partial \theta}[A(\theta) x + B(\theta) u] \\
  0 & 0
 \end{array} \right]_{\hat{x}, \hat{\theta}} , \\
 H &=
 \left[ \begin{array}{cc}
  C(\theta) & \frac{\partial}{\partial \theta}[C(\theta) x]
 \end{array} \right]_{\hat{x}, \hat{\theta}} .
\end{align*}

The general implementation of the EKF and UKF, as well as for parameter adaption, can be found in several publication. I personally like "Optimal State Estimation: Kalman, H$_{\infty}$, and Nonlinear Approaches" by Dan Simon, 2006 (ISBN: 978-0-471-70858-2).


Conditions for analytical signal reconstruction with the Hilbert Transform (HT) I'm new to the Hilbert Transform, I was reading the paper: "SISAR Imaging - Radio Holography Signal Reconstruction Based on Receiver-Transmitter Motion", and the next sentence was written in the paper:

  The HT states that an analytical complex signal can be
  reconstructed from the signal’s real part under the condition
  that the real amplitude and frequency modulated signal in the
  form:
  $$
S(t) = I(t)\cos(\theta(t))
$$
  has the spectrum of $\cos(\theta(t))$ outside the spectrum of its envelope $I(t) \in [−f_0, f_0]$.

I have checked some documentation on the Hilbert Transform and I have not been able to find that condition, so, why is this condition needed? What is the reason behind it?

Let me use the notation from Bedrosian's research report "The analytic signal representation of modulated waveforms" in a slightly simplified version:

$$s(t)=c(t)m(t)\tag{1}$$

where $s(t)$ is the modulated signal, $c(t)$ is a (not necessarily sinusoidal) carrier function, and $m(t)$ is the message signal. The carrier is assumed to be a relatively narrow-band band pass signal, and the message signal has a low pass character.

Let $c_a(t)$ be the analytic signal corresponding to $c(t)$. By definition, $c_a(t)$ has only positive frequency components. The analytic signal $s_a(t)$ corresponding to $s(t)$ can then be written as

$$s_a(t)=c_a(t)m(t)\tag{2}$$

if and only if the spectra of $c(t)$ of $m(t)$ do not overlap. Because only then does $s_a(t)$ as defined in $(2)$ contain positive frequency components only. If the spectra overlap, then $s_a(t)$ is not analytic because it will have negative frequency components. This is illustrated by the figure below, taken from Bedrosian's report:

The condition that the spectra of $m(t)$ and $c(t)$ do not overlap is equivalent to the requirement that $s(t)$ be a true band pass signal, i.e., $s(t)$ has no DC component.

I'm sure that this is what is alluded to by the sentence you quoted. I do agree that this is not obvious from the chosen formulation.


Difference between the DTFT and DFT I know this question has been asked before. It is, however, so confusing, that I'd like to give this another try:
I have come across the following 2 definitions of the DTFT:

So, the first line is the definition given in my lecture. The second line was given in a tutorial on the difference between DTFT and DFT. The first, obvious difference is the limit of the sum. While the first goes from -infinity to +infinity, the second treats a finite signal. I think DTFT are applied to infinite signals. Now, ,the second difference appears in the exponential function.  In the second line, we are using small omega (2*pi*f) multiplied with n, which should be the sample index. In the first line we are using 2*pi*f divided by fs. I think the signal's highest frequency component divided by fs gives the number of samples. So... no, I still don't see the link to n. Can somebody please help me understand this? Although this question appears a lot in this forum?

A DFT produces a finite number of results (frequency bins), so the complex result vector can be indexed by 0..(N-1).

A DTFT produces a result over a continuous range of frequencies (more than a finite number of results) so must be specified by a real number as opposed to an integer.

As MBaz says, you can sample the DTFT at a fixed integer number of points.


MFCCs and mean normalization I'm reading a blog about extracting MFCCs features for Machine Learning applications, but I didn't understand the following points about the mean normalization:

  To balance the spectrum and improve the Signal-to-Noise (SNR), we can simply subtract the mean of each coefficient from all frames.

mfcc -= (numpy.mean(mfcc, axis=0) + 1e-8)

What does he mean by "balance the spectrum and improve the Signal-to-Noise (SNR)"?
Before feeding the MFCCs features to my neural network, I'm doing min-max normalization for the whole features in the dataset instead of normalizing each file alone (like what the author does above), and I'm not sure which method gives more accurate results.
Why did he add the epsilon value 1e-8, and didn't divide by the standard deviation?

I'm reading a blog about extracting MFCCs features for Machine Learning applications, but I didn't understand the following points about the mean normalization:

Do not read random blog posts, read books instead. At least they are reviewed by the publisher.

  What does he mean by "balance the spectrum and improve the Signal-to-Noise (SNR)"?

Not spectrum is balanced but features. Signal to noise ratio remains the same after this operation.

  Before feeding the MFCCs features to my neural network, I'm doing min-max normalization for the whole features in the dataset instead of normalizing each file alone (like what the author does above), and I'm not sure which method gives more accurate results.

Per-file normalization is often better if per-band amplitude of files in the database is different. At the same time if files are small per-file estimate might be inaccurate, then per-dataset normalization might be better.

  Why did he add the epsilon value 1e-8, and didn't divide by the standard deviation?

He just does mean normalization, not variance normalization, so he doesn't divide. As per 1e-8, it is totally meaningless.


How to send / receive a short impulse in an extremely noisy signal with high timing accuracy? I am pretty new to DSP. I need a computer to transmit data to a receiver through audio signals in a loud outdoor environment. (I am writing the software for both sender and receiver so I can design the signal however I want.) I am using cheap speakers and microphones operating at 44.1Khz.

I would like for the sender to emit a signal that can be recognized by the receiver with extremely high accuracy on the time domain. This signal does not need to contain any information, the receiver just needs to be able to figure out exactly when the signal started, ideally down to the exact sample.

My first idea was for the sender to emit a 2kHz wave for 250ms. The receiver would continuously perform FFT to identify that wave, and when present, it would look at the FFT history to figure out the sample at which the wave had started. However this approach feels clunky because of the number of samples needed to satisfy the FFT. I could never achieve 1-sample accuracy without doing an FFT for each sample, and I don't think it would even give me good results.

I have been learning about PSK for unrelated reasons, and I am getting the impression that there are lots of ways to detect properties of waves – without the FFT. So my question is – what methods could I use to identify an instantaneous change in my input signal, in an extremely noisy environment?

Using a sinusoid is a bad idea. In order to precisely detect the time of arrival, you want a signal with wide bandwidth. Here's a brief description of a good approach:

Transmit a known direct-sequence spread spectrum (DSSS) waveform (typically a phase-shift-keyed digital signal, either BPSK or QPSK). You should pick a spreading code that has good autocorrelation properties: its autocorrelation should look close to an impulse. A sequence of random symbols will often work, or you can use something more well-defined like a maximal-length sequence.
At the receiver, apply a matched filter to the known transmitted pattern.
You will see a spike in magnitude at the output of the filter that corresponds to the time of arrival of the signal of interest. Apply a threshold to the filter output to identify the locations of these spikes.
The matched filter will have a delay of $N-1$ samples, where $N$ is the length of the sequence. Subtract this duration from the detected pulse locations to yield an estimate of the pulse arrival time.

Optionally, take the samples surrounding the peak and interpolate to get a finer (sub-sample) estimate of the actual peak location. Polynomial, sinc-based, or other methods may be used here.

This is a very simplified description of what you might want to do, but it should get you started. One consideration to keep in mind when designing the length of the sequence: there is a tradeoff between the total integration period (i.e. the sequence length) and the frequency selectivity of the matched filter. If you have appreciable frequency error between your transmitter and receiver, the peaks at the matched filter output will become attenuated. The amount of attenuation depends on the frequency offset's relationship to the matched filter length. This is a common attribute of DSSS systems that must be handled appropriately.


How to compute dBFS from rms-values if the samples exceed 32bits For a first approach on audio-sampling I'm currently working with the raspberry pi's pcm pins and MEMS-microphones, using the ALSA library.
A problem arises when I want to calculate dBFS from a buffer's rms-value because i have to assume that a sample can be as big as (2^32)/2. In order to obtain rms, I would have to square every single sample, some of which might not be computable due to their size.
What would be a good way of obataining the dBFS of a finite array of discrete values without causing overflows?

There are a few different ways:

Use a fractional multiply or data types. Unfortunately that's not a standard "C" or "C++" instruction but most embedded processors do support it and have some sort of non-standard way to call it.
Use a "long int" 64-bit integer data type, do a double precision multiply.
Convert everything to floating point (float or double) before multiplying.
Scale the numbers by $2^{-16}$ (arithmetic right shift by 16-bits) before squaring.

Option 1 is typically the "best" if your processor supports it as a native instruction. Option 4 is often "good enough" if you have half-way decent signal to noise ratio but that may depend on the specifics or you application and data.
Note on "fractional" data types: This is a fixed point type that is scaled so that the data range is assumed to be  from -1 to +1, hence a multiply can never overflow. That's actually the most "natural" type for fixed point signal processing. Multiplying two 32-bit numbers yields a 64 bit result. The difference between an integer and fractional multiply is that the integer multiply keeps the lower 32 bit and the fractional keeps the higher 32 bit of the result.
Both integer and fractional format are subsets of the "Q" number format. See https://en.wikipedia.org/wiki/Q_(number_format).  Loosely speaking integers are Q0 and fractionals are Q31(signed).


What is the highest relevant stopband attenuation for a 1-D digital filter? I have the impression there is little added value for more than 100 dB of attenuation in a  stopband of significant width (i.e. not just a specific frequency). What is the point at which there is little practical application for more attenuation (80 dB, 100 dB, 120 dB, ...) over a stopband of significant width? Also what factor drives this point of no added value? My guess is that the highest useful attenuation is determined by the limited precision that we can measure voltage or any parameter that we might want to filter.

That really depends on the application. For example this processor: http://www.analog.com/en/processors-dsp/sharc/adsp-21364/products/product.html
has 140 dB SNR in it's sample rate converter. That means the stop band attenuation of the interpolation filter must be at least as high but it's more likely to be 150dB - 160dB. You need extra head room to manage the fractional indices into the poly-phase filter bank in irrational sample rate converters


Why does the FFT of a complex variable create negative frequencies? So I have a demodulator that streams out X and Y values. I use a spectrum analyzer within this demodulator instrument which plots the |FFT(X+iY)| against the frequency domain, which shows up with negative frequencies. There is an equal number of frequency points in both positive and negative axis.

Whereas if I take the data of X and Y directly into MATLAB, perform the FFT of this, I only get the positive half of the frequencies as seen from the spectrum analyzer.

Excuse me if my question is vague, but I am still getting to grips with DSPs and would like to understand this more.

x = dev2305.demods.sample_x{1, 1}.value; %reading X values
y = dev2305.demods.sample_y{1, 1}.value; %reading Y values
data = complex(x, y); %Create complex variable X+jY
Y = fft(data);
P2 = abs(Y/length(data));
P1 = P2(1:length(data)/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = 109900*(0:(length(data)/2))/length(data); %109,900 is the sample rate

plot(f, mag2db(P1)); %dbV
%plot(f, P1);
title('Manual MATLAB FFT')
xlabel('f /(Hz)')
ylabel('FFT|xiy| / dBV')

It's about the type of signal being analysed.

It's known from signal analysis that a real valued signal's Fourier transform is conjugate-symmetric $H(w) = H^*(-w)$ ; i.e., its positive frqeuencies are conjugate's of its negative frequencies. Hence it's not necesary to plot the whole spectrum, but just the positive half ot it is sufficient.

However, a complex signal's Fourier transform does not posses this symmetry property and the whole spectrum would contain relevant information about the signal, therefore you should consider not only the positive frequencies but also the negative frequencies when you analyse a complex valued signal.


On the Use of OMP Algorithm to Estimate Sparse Vector As known, Orthogonal Matching Pursuit (OMP) Algorithm is to recover the sparse channel after convolution with another vector. But when I implement that in MATLAB, I don't get the sparse vector correctly. I provided below a simple example of multiplication of sparse vector with a matrix and then I tried to get back the sparse vector, but I don't get the correct results. Could you please guide me how to correct that.  Below is my example code:

   clear all; clc; clear;

    h = [1 0 0 1 0 1 1 0];    %Sparse channel
    V = randn(8,8);           %Matrix
    Yy = V*h.';

   h_esti=OMP(Yy,V,4);             %Perform OMP algorithm

Regarding the OMP function, it's related to the algorithm of OMP as here OMP

NP:  The above code gives sometimes the exact results of the sparse vector but it's not always. !!   I don't know why that, maybe matrix V must have special conditions or something wrong in my code.

Well, in your example, the channel isn't exactly sparse.

It has been shown that $\ell_0$ minimization can recover any $K$-sparse vector $x$ from observations $\Phi x$ as long as $2K < {\rm spark}(\Phi) \leq M+1$, when $\Phi$ is $M \times N$ (so that $x$ is $N\times 1$), i.e., $K<M/2$ more or less. This is a necessary condition which means that if $K$ exceeds $M/2$ it can be shown that there is no unique solution anymore (sidenote: for an MMV problem, i.e., if you have multiple snapshots, this condition can be relaxed to $K<M-1$).

The example you picked has $M=N=8$ and $K=4$. This means that you are operating right at the limit of $K=M/2$. To have an $\ell_0$ algorithm perform well, you would want to have $K$ a bit below this limit, even more so if you choose a heuristic one like OMP. Try $K=2$ or $3$ or increase $M$ and $N$ a bit, it should work better.

Further notes:

Results  like these are typically asymptotic in $M$ and $N$ and tend to become tight only as $M$ and $N$ are somewhat larger. Performance for small-scale systems is much harder to analyze.
OMP is a sub-optimal technique. It is quite sensitive to correlations in the measurement matrix $\Phi$. For $M=N$, your could use orth(randn(n,n)) to make $\Phi$ orthonormal. OMP would appreciate that.
As you asked whether there are conditions on the matrix $\Phi$. Yes, there are. It should satisfy the restricted isometry condition (RIC). If you draw it randomly, chances are high that it does, even more so the bigger it is (see first point).

*edit: Regarding your comment: yes, OMP is just not such a great algorithm. And your randn design is also not the greatest. Success rates for M=8, N=8 with your design (empirical): K=1: 100%, K=2: ~88%, K=3: ~65%, K=4: ~45%. If you use orth(randn(n,n)) instead, you should see 100% (though in this case you wouldn't need OMP). Going from 8 to 16 also improves things. The truly interesting case for sparse signal recovery is $M<N$. Here you can make $\Phi$ at least row-orthogonal to help OMP a little. Normalizing columns would be another good idea.

*edit2: To emphasize more clearly what I said above, here is a numerical example varying $K$ and showing the empirical success rate. Comparing pure random Gaussian design (red line) to one where the rows are made orthogonal to each other (blue line). It helps. Note that here $M=8 < N = 16$.


Why do digital filters work? So I just started reading on FIR and IIR filters and am amazed by how "simple" the theory seems, so far.

But what confuses me is, why does filtering work by creating a weighted sum of previous samples?
What intuition makes one think this can produce desired filtering effects?
It seems a bit unintuitive to me eventhough anyone can verify that summing delayed signals together does produce comb filtering. But desirable filtering? Why?

Consider a discrete-time input signal of the form:
$$ x[n] = \cos(\omega_0 n) ~~~,~~~-\infty < n < \infty, ~~~~~ n\in \mathcal{Z}$$
where the radian frequency $\omega_0$ is set between 0 and $\pi$ radians per sample.
Now, consider two basic discrete-time (digital) filters which are defined through addition and subtraction of their input $x[n]$ for producing their outputs $y_1[n]$ and $y_2[n]$ as:
$$y_1[n] = (x[n]+x[n-1])/2 ~~~,~~~ \text{sum filter}$$ and
$$y_2[n] = (x[n]-x[n-1])/2 ~~~,~~~ \text{difference filter}$$
Lets make a qualitative analysis of these two filters by setting their input frequency $\omega_0$ to low (close to $0$), and high (close to $\pi$) values, and then observing the corresponding outputs respectively;
First, assume that $\omega_0$ is set to low frequencies: Then the consecutive input samples $x[n]$ and $x[n-1]$ will have highly similar values, as a low frequency sine wave will not change much from one sample to the other. When this is the case, their sum will add up, whereas their difference will cancel.
Therefore, $y_1[n]$ will be approximately equal to the input's value $x[n]$, while the output $y_2[n]$ will be close to zero, due to cancelling of $x[n]$ by $x[n-1]$ when subtracted. Then we conclude that the sum filter passes the low frequencies while the difference filter attenuates (blocks) them.
For the second part of the analysis, set $\omega_0$ to high frequencies: then the values of the input samples $x[n]$ and $x[n-1]$ will be of opposite polarities, as the cosine will be rapidly changing from sample to sample. Then, their sum  will cancel, whereas their difference will add up.
Therefore, $y_1[n]$ will be approximately zero, while $y_2[n]$ will be following the input. Then we see that the sum filter blocks high frequencies while the difference filter passes them.
Combining these two analyses we conclude that the sum filter is a lowpass filter, which passes low frequencies and attenuates high frequencies, while the difference filter is a highpass filter, which attenuates low frequencies and passes high frequencies.
Following this basic basic setting, more complex filters are realized through using more samples at farther delays, and are weighted properly. Then the passband and stopband cutoff frequencies, the transition bandwidth and ripples are all determined by those delays and weights applied into the summing (or differencing) samples and the number of samples (length of filter) used in the sums (or differences).
Those weights are then called as the filter coefficients (or its impulse response $h[n]$ for the FIR filter) that characterize the filter.


Decreasing Sample Rate in DFT (FFT) for Audio Analysis I'm working with audio data analysis through FFT algorithm.
My example audio is a sine wave at 440 Hz and 44100 Hz sampling rate.

FFT methods in programming, like scipy.fft(y), require a vector of samples and a number of points/samples to analyze.
Sample rate is used to determine the width of each frequency range/bin in the resulting data vector (FS/N).

I've read that to increase frequency resolution of FFT results one should dicrease sampling rate and increase window size (number of samples).

I understand that I can change the N value (samples count) that I pass to the algorithm.

But if the audio is at 44100 Hz how can I dicrease it?
Would the resulting frequency bands reflect incorrect information about magnitudes of the FFT data?

Or should I change the input audio data somehow (manipulating directly the samples in the vector) before FFT so the data is passed with lower rate (say at 22050 Hz)?

Thanks.

I've read that to increase frequency resolution of FFT results one should dicrease sampling rate and increase window size (number of samples).

To increase frequency resolution, you increase the window size (increase the number of samples).  The sample rate does not change.

For simplicity, I'll pretend your signal is sampled at 800 Hz.

If you were to take 4 samples of the audio, and run an FFT on it, you would get 4 samples in the frequency domain.  They would represent 4 frequency bins:

[0 Hz (DC), 200 Hz, ±400 Hz (Nyquist), -200 Hz]

Each of these bins is 200 Hz wide.  If instead, you had taken 8 samples of the audio, and run an FFT on it, you would get 8 samples in the frequency domain, representing 8 frequency bins:

[0, +100, +200, +300, ±400, -300, -200, -100]

Each of these bins is 100 Hz wide.  Does that make sense?  So to get narrower bins, you need to take more samples of your audio.  But the maximum frequency in the FFT (the Nyquist frequency of ±400 Hz) does not change; it's the same regardless of the number of samples.

Note that you can get these values like this:

scipy.fftpack.fftfreq(4, 1/800)
Out[8]: array([   0.,  200., -400., -200.])

scipy.fftpack.fftfreq(8, 1/800)
Out[9]: array([   0.,  100.,  200.,  300., -400., -300., -200., -100.])


Sampling Audio Signal I generate a 19kHz audio signal with sampling frequency of 44100 using Matlab. Due to the sampling, some other frequencies are generated in addition to the 19 kHz signal. When I play this audio, I hear additional voices. How I can remove these additional sides? I think it needs aliasing filtering, but I cannot filter signal after DAC, because I generate the audio signal and play it with phone or computer, I cannot put any filter after A/D.

If your code is working correctly, than there should be no aliasing.

Most likely this is a problem with your playback system. Most sound cards in computers are terrible and will create a lot of noise and other artifacts. Since the 19 kHz tone is inaudible for most people, you just hear the artifacts.

It's also possible that you turn up the volume too high. The 19 kHz may easily clip your playback hardware. You can't tell if it's "too loud" since your auditory sensitivity  at 19 kHz is basically zero. You may have also clipped the sine wave when sending it to the sound card.

There is also the software stack: On Windows, most sound will go through the kernel mixer, which often applies sample rate conversion and other processing. It's default behavior is quite bad, so this is also likely to create a lot of artifacts and noise.

Playing sine waves cleanly requires decent equipment. Typically a good external sound card and a dedicated SW driver that bypasses the operating system. You can try saving it as a wave file and play it over some "known good" hardware.


Difference between 'conventional' and 'adaptive' beamformers? This might be a terminology question but I am not sure.

Basically, what is the difference between conventional beamformers, and adaptive beamformers? I thought that all beamformers were inherently adaptive to some criteria, like minimization of distortion or variance, or some other spatial filtering criteria etc. So what is the difference between them?

A beamformer is basically a spatial filter.  It can be passive, just like a temporal filter.

Instead of samples separated by time, they are separated by space.  A passive temporal filter can be a bandpass that is "aimed" or "steered" at a particular frequency.  For passive spatial filters (i.e. beamformers), the filter can be steered towards a particular angle of arrival, instead of temporal frequency.

Adaptive filters/beamformers can be incredible because they can "steer" a null towards a frequency or angle where there is an interfering signal.

Just like a temporal adaptive filter, a spatial adaptive filter (i.e. beamformer) is constantly adjusting the filter weights/coefficients to optimize some criteria that usually involves "nulling out" or "rejecting" an interferer.

Here's a diagram of a CBF that is called a "K-Omega" Beamformer.

Here's a diagram to reinforce the idea that a passive beamformer is possible.

I realize that these are a bit random, but hopefully you can follow the logic of the above.
I'll see if I can find a better diagram that more clearly shows what's happening.  To be clear, the windows referred to above are Hamming/Hanning type windows, and this process is basically returning a 2-d matrix where temporal frequency is on the x-axis and a special spatial variable is on the y-axis.  This special variables makes the math easier, and it's a one step process to convert the special spatial variable and temporal frequency to an angle of arrival.

The green lines below are lines of constant angles.

Keep in mind that all of this is for CBF (conventional beamformers) and the above covers all spatial and temporal frequency (within the limits of Nyquist.
Some common ABF techniques are:

MPDR - Minimum Power Distortionless Rejection BF (Described by Van Trees)
DMR - Dominant Mode Rejection BF (Abraham and Owsley)
R-DMR - Robust Dominant Mode Rejection BF (Cox and Pitre)
EBAE - DMR BF with Eigenvector Beam Association and Excision (Kogon)

Instead of using a spatial FFT as the spatial filtering operation, these techniques typically involve constructing an updated spatial autocorrelation matrix of the incoming signal, and then using that matrix or eigenvectors of that matrix to adaptively influence the spatial filter.

Update for @Mohammad:
Below is the list that of beamforming texts I received from my professor:

Van Veen, B.D.; Buckley, K.M.; , "Beamforming: a versatile approach to spatial filtering," ASSP Magazine, IEEE , vol.5, no.2, pp.4-24, April 1988
Efficient digital beamforming in the frequency domain
Brian Maranda, J. Acoust. Soc. Am. 86, 1813 (1989)
Array Signal Processing by Don Johnson and Dan Dudgeon, Prentice Hall, 1993
Optimum Array Processing by Harry L. Van Trees, Wiley, 2002.

Also, after searching around, this looks very interesting.  It seems to be more interested in the practical/implementation side than the theory.  I don't have it but I'll probably buy a copy:

Practical Array Processing by Mark Sullivan, McGraw Hill, 2008 ( http://www.practicalarrayprocessing.com/index.html )


Estimating confidence intervals for the mean value Say I conducted a numerical experiment -a time-dependent fluid simulation- and extracted a time-series of some variable. It looks like this:

As a result of the simulation, I am supposed to report the mean value of this signal along with some estimate how accurate this mean value is from a purely statistical point of view. The Idea is: If I conduct the same experiment several times, I get different signals with different mean values. These mean values would be normally distributed. The mean value and standard deviation of this distribution are unknown.

Now if I could perform the experiment several times, estimating the variance of the mean would be straightforward. Then I could estimate some confidence intervals using Student's t-distribution. Yet since compute resources are scarce, this is not an option. All I have is this one sample.

tl;dr:
Is there a way to estimate a confidence interval for the mean when all you have is one sample?

I was thinking about making use of the fact that the signals autocorrelation function vanishes rather quickly. But I haven't been able to turn this into a reliable method.

To the extent you can assume your process is a white stationary ergodic process, the variance of the mean is the variance of the process divided by N where N is the number of samples.

Assuming that you can estimate the variance from your sample using

$$\sigma_x = \frac{1}{N-1}\sum_{i=1}^N(x_i-\mu_x)^2$$

NOTE: This is the unbiased estimator since the variance is being estimated from a sample of N elements using a mean that is estimated from the sample itself. (see Zwillinger, D. (Ed.). CRC Standard Mathematical Tables and Formulae. Boca Raton, FL: CRC Press, 1995.). Otherwise if the mean was known (which it isn't) there would be a $\frac{1}{N}$ in front of the summation above.

You should then be able to derive confidence intervals on the mean itself with the variance of the mean as $\frac{\sigma_x}{N}$:.

This is seen in the expression for the mean and variance of $Y$ where $Y$ is the average of $X$ over $N$ samples, as follows where $\mu_x$ is the mean of x and $\sigma_x$ is the variance of x:

\begin{align}
E(Y) &= E\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N}E\left[(X_1+X_2+ \ldots +X_N)\right]\\
\textrm{Var}(Y) &= \textrm{Var}\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N^2}\textrm{Var}\left(X_1+X_2+ \ldots+X_N\right)
\end{align}

\begin{align}
E\left[(X_1+X_2+ \ldots +X_N)\right] &= E[X_1]+E[X_2]+ \ldots E[X_N] = N \mu_x\\
\textrm{Var}(X_1+X_2+ \ldots +X_N) &= \textrm{Var}(X_1)+\textrm{Var}(X_2) + \ldots \textrm{Var}(X_N) = N\sigma_x^2
\end{align}

Therefore
\begin{align}
E(Y) &= \frac{1}{N}E[(X_1+X_2+ \ldots +X_N)] =\frac{1}{N}N\mu_x =  \mu_x\quad\text{and}\\
\textrm{Var}(Y) &= \frac{1}{N^2}\textrm{Var}(X_1+X_2+ \ldots +X_N) = \frac{1}{N^2}N\sigma^2 = \frac{\sigma^2}{N}
\end{align}

How to make use of the Autocorrelation Function

Note that this reduction by N is valid as long as the sequence is white. Once samples are correlated, there will be no further reduction in the variance of the estimate. You can therefore make use of your autocorrelation to determine the number of samples that will reduce the variance. I do not have an exact calculation for this, but to provide a rough order magnitude I would use the number of samples that results in the normalized autocorrelation dropping to below 0.5. For example, in your data the total number of samples is 30,000. If the autocorrelion immediately dropped below 0.5 after just one sample, then all samples are independent and your variance estimate would be the $\sigma_x/30,000$. However if the autocorrelation does not drop below 0.5 until 100 samples, then the variance estimate would be $\sigma_x/300$.


Brightness equalization between different images I am dealing with a proprietary whole slide image format from 3DHistech, namely MRXS. This image is produced by scanning a glass slide biopsy. The scanner takes pictures at certain positions which eventually make up one big image (in the range of hundreds of thousands pixels per side).

Subsequently, the software produces downscaled versions of the original image, effectively generating a pyramid, in order to make viewing the image faster.

The camera pictures, as well as the downscaled versions, are saved as small rectangular JPEG images, usually around 300x300 pixels.

I have written a parser for the image and I am able to render these digital slides appropriately. However, I came across some slides that have different brightness (?) among tiles. It seems that the vendor's software does compensate for this, while my software obviously doesn't.

I have opened some of the JPEG streams in Photoshop (among other tools), and I was able to manually fix the brightness difference and make the image look smooth. My question therefore, is what kind of algorithm could I incorporate to do this automatically as I render these images?

I am attaching a sample original tile produced by the scanner software and an edited version I did in Photoshop.

ORIGINAL

EDITED

You could do Histogram equalization on the histogram of brightness values across all tiles.

Alternatively, you can loop over all tiles, accumulate the minimum and maximum brightness value of each tile, and then scale all brightness values in all tiles against the overall minimum/maximum.


DSP project suggestion I apologize in advance if this question is not suited to this forum.

I have a DSP background from school(graduate EE courses). I am good on theory, but have never worked on a real project. After graduating and spending 2 years working as firmware engineer, I come to realize that I want to pursue DSP and algorithms for my career rather than pure firmware coding. I want to develop a portfolio of projects that can show the future employer that I can do this stuff. But I am a little short on ideas. Do you have any project suggestions?

Thank you very much.

There are soo many projects, field of applications, platforms, architectures and tools to mention about but let me put here a few.

1- Begin with Audio projects, as probably every DSPeng. did so.
There's an open source program AudaCity which involves very rich set of audio functions implemented. So write your own, add your own, develop your own.

Again, audio dsp effects is an area rich for new and better algorithms. Try developing your own.

Digital audio synthesizers is another rich field with immense set of applications.

2- Do image video processing applications. If you are into data compression as well I highly suggest you FFMPEG project. Get into such a project and develop your own applications. There are many filters, effects, algorithms etc that are implemented and you can definietely add your own.

All these require nothing but a PC to develop and implement.

3- In addition, you can search hardware based projects too but they will require a more elaborate set-up. Try implementing some embedded systems using development boards. ADI, ST, TI, NI, NXP provide very convenient, versatile application development boards using their own microcontrollers, microprocessors and DSP procesors, multimedia processors and sensors.

4- Look for some SDR (software defined radio) for most affordable digital communications applications. RTL-SDR is its home. Get a usb-sdr dongle (a hack of DVB-T receiver) and write your own pure C-code decoding algorithm there.

5- A raspberry pi 3 with a pi-camera can provide you a real opportunity to deal with raw camera processing. You can develop a lot of video coding algorithms there.


How can I plot frequency vs magnitude of wavelet transform? I am running Morlet continuous wavelet transform. I have got wscalogram of signal and now I want to plot freq-magnitude like the following picture., but I don't know how to do it:

I've used scal2freq MATLAB function to convert scales to pseudo frequencies. Also I have some frequencies in my signal that have a large damping ratio(4%) so they are not well visible in plot. How can I exaggerate these highly damped modes?

I am using MATLAB, here is my codes:

% Import the text4.txt to matlab workspace. and save it under name "data"
t=linspace(0,30,301);
Fs=ceil(inv(t(2)-t(1)));
x=data(:,4); % use x=data(:,3),x=data(:,5) too. first column is time,second is refrence
wname = 'morl';
scales = 1:1:256;
coefs = cwt(x,scales,wname,'lvlabs');
freq = scal2frq(scales,wname,1/Fs);
surf(t,freq,abs(coefs));shading('interp');
axis tight; xlabel('Seconds'); ylabel('Pseudo-Frequency (Hz)');
axis([0 30 0 1 0 60])
xlabel('Time'); ylabel('Frequency');
figure;
sc=wscalogram('image',coefs,'scales',freq,'ydata',x);
hold on

My current plot looks like this:

Continuous wavelet transform (CWT) is a method for time-scale analysis. Yes, you read it correctly, scale, not frequency. However, it is possible to map the scales to frequencies, and even quite easily. Since you are a MATLAB user, you will probably want to use this function, which does the following:

  F = scal2frq(A,'wname',DELTA) returns the pseudo-frequencies corresponding to the scales given by A, the wavelet function 'wname' (see wavefun for more information) and the sampling period DELTA.

As you can read, the frequencies are not the real frequencies in the signal, but approximations. If you want the real frequencies, you should be using the short-time Fourier transform (STFT).

So, when should one use the CWT? As a general principle, wavelets work well when the signals are transient (i.e. have quick changes).

EDIT: haha, I don't know what I'm doing here. sorry :-D completely missed the point of your question. Just read the title....


Design one of 4 type of FIR filter by window function? There're 4 type of FIR to chose from

symetry ,even
symetry ,odd
antisymetry, even
antysymetry, odd

Also, there are many type of window like rectangular, Hanning, Kaiser,etc So how to use a window function to design the N-length filter and what's the point of doing it? That's what my teacher ask me to do but I don't understand the reason

What I confuse here is that for example I want the ideal filter $H(\omega)$ that have passband in $\omega_a <\omega<\omega_b$ then what I need to do is just use the inverse Fourier transform to the frequency response of the ideal filter to find $h[n]$ that is infinite then use the window function to make it finite. Then what's the point of make it look like the form of one of those 4 type? Moreover, what I need to do after doing inverse Fourier Transform and window function to make it in the form of FIR digital filter with N-length?

The topic your referring to goes by various names such as "FIR Filter Design Using Windows", "Windowing Method of FIR Filter Design", "Window Designed FIR Filters", etc. I'll bet your class textbook covers this topic (most DSP books do). Meaningful answers to your questions require more space than is available here. I suggest you have a look at the following web page (and its following three pages) for a decent introductory discussion of that topic:

http://www.bores.com/courses/intro/filters/4_fir.htm

Then have a look at Section 2.2 of the following web page:

http://www.mikroe.com/chapters/view/72/chapter-2-fir-filters/


The notation to express a datum filtered by arbitray number of one-dimensional filters? I'm a pure novice in signal processing but I'm writing a paper that needs to be expressed in the SP way. The question is, I want to express an N-dimensional datum goes through N one-dimensional filters, each along a different dimension of the datum. The simplest example is applying Gaussian blur filter on a 2-D image. You can separate the 2-D Gaussian filter into 2 1-D Gaussian filter and apply it along each dimension of the 2-D image one after the other. But my question is, is there an elegant notation to express such operation?

...my question is, is there an elegant notation to express such operation?

If you know that these operations are to act on orthogonal dimensions, you can define their 1D version and then use indexing to express the order by which the function is applied to the multidimensional signal.

But if you are working with multidimensional signals over an orthogonal vector basis throughout, you can simply specify this in the beginning and use (for example) $\mathcal{F}(x)$ throughout, to denote a multidimensional discrete Fourier transform.

Hope this helps (?)


Can you present the convolution of sinusoidal with itself? Ladies, Gentlemen,
Because I am homeless (in France), and get internet access only in public libraries with many restrictions in timing etc, I can not write down even a simple convolution. So I ask you post here the convolution of  some sinusoidal function with itshelf.

For example: 0.951056516295, 0.587785252292, -0.587785252292, -0.951056516295, 0.

Regards


What is perfect and imperfect channel? Often in estimation of the input signal to an unknown channel, the bit error rate is plotted as the performance measure. The curve is compared with an ideal curve labelled as 'Perfect channel'. For example in this paper, http://www.iaeng.org/publication/WCECS2015/WCECS2015_pp668-673.pdf
 although the bit error rate curve is not present, in the introduction the Authors do mention the term perfect channel. One such method of estimation is the Least Squares.

What is the meaning of a perfect channel? Does a perfect channel mean that there is no noise or does it mean that the pilot symbols are used to estimate the channel?

It's not the channel that is perfect; it is the estimation. So, "perfect channel estimation" means that the receiver knows the fading coefficients perfectly. This is of course only true in a simulation; any actual system  performs imperfect estimation and the estimated channel coefficients are only approximately equal to the actual coefficients.


Understanding a simple case of fft2 in python I have 2D arrays representing surfaces. I have created them by scanning some theoretically flat objects and I ended up with what can be called "imperfections", or else, deviation from the perfectly flat surface.

I hope to get information on the dominant wavelengths of imperfections and their amplitude, so I employed numpy's fft2. Since I am a first timer, I did some digging and tried to play around with some basic examples to get a grip.

Here is one:

import numpy as np
import matplotlib.pyplot as plt

"""
fft2 playground.
"""

# Initialise an empty array
field = np.empty([29, 481])

# Set the wave amplitude
amp = 2.5

# Create the synthetic sinusoidal field
for n, i in enumerate(np.linspace(0, 2 * np.pi, field.shape[0])):
    for m, j in enumerate(np.linspace(0, 2 * np.pi, field.shape[1])):
        field[n, m] = amp * np.sin(j)

# Perform 2D fourier and shift the result to centre
f = np.fft.fft2(field)
fshift = np.fft.fftshift(f)

# Calculate the magnitude and phase spectra
magnitude_spectrum = 20*np.log(np.abs(fshift))
phase_spectrum = np.angle(fshift)

# Reconstruct the initial field
f_ishift = np.fft.ifftshift(fshift)
re_field = np.abs(np.fft.ifft2(f_ishift))

# Plot
fig = plt.figure()

fig.add_subplot(411)
plt.imshow(field, cmap='gray')
plt.title('Field'), plt.xticks([]), plt.yticks([])
plt.colorbar()

fig.add_subplot(412)
plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('Magnitude spectrum'), plt.xticks([]), plt.yticks([])

fig.add_subplot(413)
plt.imshow(phase_spectrum, cmap='gray')
plt.title('Phase spectrum'), plt.xticks([]), plt.yticks([])

fig.add_subplot(414)
plt.imshow(re_field, cmap='gray')
plt.title('Reconstructed field'), plt.xticks([]), plt.yticks([])
plt.colorbar()

plt.show()

Which gives:

Obviously the returned field is not the same as the input. The amplitude is correct but it is mirrored and 90 degrees off-phase. I feel I am missing something but I cannot nail it.

Also, I would expect that the magnitude spectrum would be something like 2 dots, one on either side of the central pixel and nothing else. What I understand is that the central pixel is the zero-frequency and moving to the edges are the higher frequencies, up to the Nyquist.

Worth noting, if I change the dimensions of the input field I get different results on the magnitude and phase spectra.

Could someone point out what am I missing please? How would I get back the initial field from the fft and why the magnitude spectrum is the way it is?

Thanx!

When you construct re_field, ditch the np.abs. It turns all negative values of the sine wave positive. Instead, do a np.realon the result of the ifft2 to get rid of the imaginary parts caused by rounding errors. Then you will get a result perfectly matching the original.

In the magnitude spectrum you cannot see the sine wave you inserted, because its frequency is to low. The vertical black line is the DC part in x direction, which is zero, the horizontal white line is the DC part in y direction, which has full ampitude because there is only DC in y direction.

Try some higher frequencies, then you will get a better idea.


Total Harmonic Distortion calculation and its origins I ve been researching about harmonic distortion and THD. From what i ve seen on the internet, HD is caused by non-linear signal processing. I want to calculate it using an digital signal processor. For this, i suppose i need to compare the distorted signal with the original one.

My question is, is it possible to generate THD in a digital signal only by digital processing? Or i need to cause it in an analogic enviroment then convert to digital for detection?

Cause i ve only found examples of it being caused by analogic processing like amplifiers.

You can create non-linear digital systems (an example would be a system that finds the absolute value of the input). You can also simulate an analog non-linear system using DSP. The easiest way is to use a power series. Assuming that the non-linear system you're modeling is time-invariant and memoryless, you can approximate it as follows:

y(t) = a0 + a1*x(t) + a2*x(t).*x(t) + a3*x(t).*x(t).*x(t) + ....

where y(t) is the system output, x(t) is the input, and an are coefficients that depend on the system you're modeling.

As an example, you can do a simple simulation of mild harmonic distortion by setting x(t) equal to the a cosine signal, and setting a1=2, a3=-0.1, and all other coefficients equal to 0. a1 is the system gain, and a3 is the non-linear gain. By setting a3 to larger negative values, you increase the nonlinear effect -- you can think of it as an amplifier getting closer to saturation.

If you plot the spectra of x(t) and y(t), you'll clearly see harmonic distortion occuring.


Computation of parameter filter to match a given frequency response I'm looking for the practical way to compute the parameters ($a_1$, $a_2$) of a digital filter to match a certain frequency response. I'm studying the 12-poles filter of the vintage component SP0256 which is used in speech synthesis, using bandpass filters to create frequency response of allophones.
The 12-poles filter is composed of six 2-poles filters :

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

I came accros a github of the FPGA implementation of the SP06256 which provide a matlab example of this filter for the 'EH' allophone. I ploted the 6 frequency responses in blue and the final one in red with matlab :

I could find the formant frequency for the other allophones, but I can't figure out how to compute the filter parameters to match the desired frequency response.

TL/DR:

For a 2nd order transfer function of a real digital filter given as:

$$H(z) = \frac{1}{1-2A_1z^{-1}-A_2z^{-2}}$$

The filter coefficient $A_1$ and $A_2$ are determined from the gain and resonant frequency to be:

$$A_1 = Kcos(2\pi f_r/f_s) = Kcos(\omega_r)$$

$$A_2 = -(K^2)$$

Where:

$f_r$: Resonant frequency in Hz

$f_s$: Sampling frequency in Hz

$\omega_r$: Normalized angular frequency ($2\pi f_r/f_s$)

$K$: Magnitude of each pole

$G$: Gain, with the gain in dB as $20Log_{10}(G)$

K and G are related with the following formula:

$$G = \frac{1}{(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)}}$$

I was not able to solve for $K$ as a function of $G$ and $\omega_r$ in closed form (however see reasonable approximation in Update section below); For a given $G$ and $\omega_r$ I would determine $K$ iteratively, or solve for the root of the following equation:

$$(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)} - \frac{1}{G} = 0$$

Where the solution for K given G is the real positive root between 0 and 1 of

$K^4 + a K^3 + b K^2 + a + c = 0$

with

$a = -2(1+ cos(2\omega_r))$

$b = 4cos(2\omega_r)+2 = -2(1 + a)$

$c = 1- \frac{1}{G^2}$

Thus the solution for K from G can be found from Matlab

roots([1 a b a c])

or Python (NumPy) using:

np.roots([1, a, b, a, c])

Note as Olli pointed out in the comments, the actual resonant frequency will vary slightly from $f_r$; for more on that see his detailed analysis at Precise Centre frequency of an All-pole digital filter

UPDATE

A reasonable approximation for $K$ given $\omega_r$ and $G$ that applies for large gains was determined by @Travis on the mathematics site at this link: https://math.stackexchange.com/questions/3035346/solving-for-ellipse-parameters-given-a-radius-and-angle-challenge-2/3035437#3035437

Resulting in the following for $\omega_r$ not approaching $0, \pi$:

$$K \approx 1 - \frac{1}{Gsin(\omega_r)}$$

With $\omega_r$ = $0, \pi$

$$K = 1- \frac{1}{\sqrt{G}}$$

DETAILS:

The basic 2nd order digital bandpass filters can be constructed as follows:

Position the poles on the z-plane to be close to the unit circle at the desired frequency location for resonance. For a real signal the poles will come in complex conjugate pairs, hence the 2nd order sections. There will also be two trivial zeros at the origin. This is demonstrated in the graphic below showing the pole locations on the z-plane. For the z-plane, the unit circle (drawn) represents the frequency axis, starting with DC at z = 1, and rotating to half the sampling rate at z = -1.

The two poles shown are located at at 30 degree angle or $2\pi \frac{30}{360}$ in radians. So if the sampling rate was 10 KHz (for example) this would represent a resonant frequency of $10e3 * \frac{30}{360}$ since the 360° around the complete circle represents the sampling rate. The closer the pole is placed to the unit circle, the tighter the resonance and the higher the gain. The gain is determined using the following relationship once doing the math starting with the pole locations in this case of complex conjugage poles:

$$H(z) = \frac{z^2}{(z- p_1)(z-p_2)}$$

$z_1$ and $z_2$ are complex conjugate poles, placed to be close to the unit circle at the angle that represents the resonant frequency:

$p_1 = Ke^{j\omega_r}$

$p_2 = Ke^{-j\omega_r}$

Substituting these in the formula for H(z) along with use of Euler's Identity ($2cos\omega_r = e^{+j\omega_r}+ e^{-j\omega_r}$) results in:

$$H(z) = \frac{1}{1 - 2(Kcos\omega_r) z^{-1} + K^2 z^{-2}}$$

where $K$ is the magnitude of the pole (close to but less than 1), and $\omega_r$ is the angle of the pole in radians, representing the frequency of the resonance as given by the normalized angular frequency.

From this it is clear that:

$$A_1 = Kcos(\omega_r)$$
$$A_2 = -K^2$$

The overall gain for the 2nd order filter is related to the inverse of the distance from the unit circle when z is at resonance to the pole. The exact relationship is complicated given the two pole locations, resulting in a gain that varies depending on the resonant frequency for a given pole magnitude as given by:

$$G = |H(z=e^{j\omega_r}| = \frac{1}{(1-K)|(1-K)cos(\omega_r) + j(1+K)sin(\omega_r)|}$$

Note, in case it helps illuminate a closed form solution for K, that the second part of the denominator has a magnitude that starts at $(1-K)$ for $\omega_r = 0$ and grows as an ellipse to $(1+K)$ at $\omega_r = \pi/2$

The simplest I could get this to is as follows:

$$G = \frac{1}{(1-K)\sqrt{(1-2K cos(2\omega_r)+K^2)}}$$

When the resonant frequency is on the real axis the gain is maximum for a given K:

$$G(\omega_r = 0, \pi) = \frac{1}{(1-K)^2}$$

And when the resonant frequency is on the imaginary axis the gain is at a minimum for a given K:

$$G(\omega_r = \pi/2) = \frac{1}{(1-K)(1+K)}$$

Where K is the magnitude of each pole and G is the resulting gain of the second order filter. (20Log(G) in dB)

The exact solution for K given G for all cases of resonant frequencies is quite entailed, so an iterative solution may be more practical (vary K to converge on G).

For example for the first curve in the OP's plot, from the graph it appears to be resonant at 600 Hz with a gain of 33 dB with an assumed sampling rate of 10KHz. From this we would get:

The resonant normalized angular frequency in radians is:

$$\omega_r = 600/10e3 * 2\pi = 0.377$$

33 dB converted to magnitude:

$$10^{(33/20)} = 44.67 = G$$

The magnitude for each of the two poles is then found from:

$$44.67 = \frac{1}{(1-K)\sqrt{(1-2K cos(2*.377)+K^2)}}$$

$$= \frac{1}{(1-K)\sqrt{(1 - 2K(0.729) + K^2)}}$$

Solving for the root of

$$(1-K)\sqrt{(1 - 1.458 K + K^2)}-1/44.67 = 0$$

Using the approach in the introduction:

$a = -2(1 + cos(2\omega_r)) = -3.4579$
$b = -2(1 + a) = 4.9158$
$c = 1-1/G^2 = 0.9995$

>>> np.roots([1, a, b, a, c])
[0.7294172+0.68505263j 0.7294172-0.68505263j 1.02993616+0.j 0.96914237+0.j]

From which we get from the real root:

$$K = 0.9691$$

Thus $A_1$ and $A_2$ for this case are:

$$A_1 = Kcos(\omega_r) = 0.901$$
$$A_2 = -K^2 = -0.939$$

$$H(z) = \frac{1}{1-2A_1z^{-1}-A_2z^{-2}}$$

= $$\frac{1}{1-1.802z^{-1} + .939z^{-2} }$$

With the frequency response as given below:


How to Improve Least Means Squares (LMS) / NLMS Filter Performance? Are there ways to increase computational performance of a normalized least squares (NLMS) filter? Multidelay block frequency-domain (MDF) filters have been proposed to do this, but they also take away from convergence speed and accuracy, because they only update the estimated impulse response once every block, not once every sample. Are there any other methods?

If you want to implement the "standard" NLMS algorithm without cutting any corners, then you're probably not going to find a structure that is significantly more efficient. Block forms of LMS filtering aim to use fast convolution techniques (like overlap-save or overlap-add) to speed that part of the process. However, as you noted, the filter coefficients are only updated per block, as the filter must be constant over the block to use the fast convolution approach.

The highly recursive nature of NLMS is going to limit you if you want to keep the sample-by-sample update characteristic. While the filtering action is non-recursive, the filter coefficients at time instant N are a function of the coefficients at time instant N-1, which limits your ability to speed the process by using parallelism or block-oriented computation. As in most cases, there isn't a free lunch: if you want pure NLMS, then you're best off just implementing that.


Constellation diagram for QPSK modulation I need to plot (in MATLAB R2011a) constellation diagram of a QPSK modulated signal after AWGN was added to the signal.  My current MATLAB code looks like this:

b=round(rand(1, Nbits));
fill=10;
T=1/1000;
fs=120000;
fc=4000;
t=0:2/fs:length(b)*T-1/fs;

for k=1:floor(Nbits/2) %Dividing pair and odd bits
    even(k)=b(2*k);
end
for k=1:round(Nbits/2)
    odd(k)=b(2*k-1);
end
E=A*sqrt(T/2);
for k=1:length(even)  %NRZ code
    if (even(k) == 1)
        evenNRZ=[evenNRZ ones(1, fill)*E];
    else
        evenNRZ=[evenNRZ -ones(1, fill)*E];
    end
end

for k=1:length(odd)
    if (odd(k) == 1)
        oddNRZ=[oddNRZ ones(1, fill)*E];
    else
        oddNRZ=[oddNRZ -ones(1, fill)*E];
    end
end
fi1=sqrt(2/Ts).*cos(2*pi*fc*t); %Basic functions
fi2=sqrt(2/Ts).*sin(2*pi*fc*t);

up=fi1.*evenNRZ;
down=fi2.*oddNRZ;

QPSK=up+down; %QPSK signal

QPSK_n=awgn(QPSK,20,'measured'); %QPSK signal + noise

This is the method we are required to use to make QPSK modulation simulation for an assignment. Now I need to plot constellation diagram from QPSK_n signal. Anyone has any idea on how to do it?

I think you may be missing some fundamental understanding of QPSK modulation, which was likely the point of the homework exercise so let me hopefully fill in some blanks, from which you can figure out how to code up in your Matlab.

This diagram of a QPSK modulator may be helpful, taken from http://www.masoodkh.com/projects/communications/18-qpsk-modulatordemodulator:

Note the similarity to your implementation, where the data is converted to +/-1, and multiplied by the sine and cosine of the carrier.  The process of doing this takes your waveform that would otherwise be centered at frequency = 0 (DC) which is the baseband waveform and move it to a carrier frequency (where it is easier to move the signal over great distances). With such a linear modulator, everything about that signal at DC is the same at the carrier in terms of the signal's amplitude and phase, we just change the frequency the signal rides on. So we start with a "DC" carrier and move it to $f_c$ which is what the block diagram is doing and what you are doing in your code. The constellation is the baseband amplitude and phase assignments that are made at DC prior to shifting the signal up to a higher frequency carrier.

We call the baseband data IQ data as it is typically a complex data stream with real and imaginary components, "I" is in-phase and represents the real portion and "Q" is the quadrature-phase and represents the imaginary components. This is certainly the case with QPSK, we take every other bit and assign it to the real axis and the others to the imaginary axis to get I+jQ at baseband.

Your original constellation and how to plot it without the noise should now be very straight forward (plot I+jQ before it gets multiplied by the sine and cosine functions).

To plot the constellation after adding noise seems to me would be very challenging to ask for a homework question, which makes me think there is an easier way than what I proceed to describe below (and hopefully someone else can provide that), but here are my thoughts:

It is not clear if the constellation with noise should be plotted for all samples over the symbol duration, or if you are to get an average of the noise to plot the one sample that a demodulator would determine before making a decision (it would be the latter if the assignment was clear that there should be only one sample per symbol), but let me describe both cases:

Note that the "carrier" is multiplying your baseband signal by $e^{j2\pi f_c t}$ since

$$e^{j2\pi f_c t}= \cos(2\pi f_c t)+ j\sin(2\pi f_c t)$$.

We assign a $j$ to the Q input of the modulator and multiplying by this identity according to Euler matches the block diagram above. The implementation is specifically equal to taking the real part as follows:

$$Re[(I+jQ)(e^{-j2\pi f_c t})] = I\cos(2\pi f_c t)+jQ\sin(2\pi f_c t)$$

Importantly based on the 4 possible dibits on the input ({1,1},{1,-1},{-1,1},{-1,-1}) we assign 4 phases (45°, 135°, 225°,315°) to the carrier. Note that the math and seeing the result is actually much easier using the $e$ instead of sines and cosines, understanding that $Ke^{j\theta}$ is a vector on the complex plane with magnitude $K$ and angle $\theta$.

Note that multiplying a signal by $e^{j\omega_c t}$ (where $\omega_c=2\pi f_c$) simply shifts the signal in frequency in one direction, depending on the sign of $\omega_c$:

A modulated signal that is varying in amplitude and phase (and in your case of unfiltered QPSK, varying predominantly in phase only) will have the same amplitude and phase variation at baseband (DC) as at any carrier, as represented in the IQ diagrams below. It's just in the case of a carrier, the amplitude and phase is plotted with reference to the carrier. This plot shows how multiplying a signal at the carrier frequency by $e^{-j2 \pi F_o t}$ (where $F_o$ was used in the diagram but is the same as your $f_c$) will move the signal from the carrier back to baseband (DC).

Ultimately to see your constellation with the noise, you want to move your signal that is at the carrier back to DC. A further complication is that your signal as modulated is completely real (note that I took the real part after multiplying by $e^{j2 \pi f_c}$ when we did the upconversion to make it match your implementation. A complex signal can be single sided (positive frequencies only for example) as we described would happen after multiplying by $e^j2\pi f_c t$ (and note that I specifically multiplied by $e^-j2\pi f_c t$ in order to match your implementation, which means I shifted the frequency to a negative frequency carrier!). Once we take the real part of a complex frequency, what was once a single sided spectrum (if it was, doesn't have to be but will always be the case is the positive spectrum is different from the negative spectrum), will become a double sided spectrum that is complex conjugate symmetric: For all real signals the spectrum will always be complex conjugate symmetric- meaning the negative frequency spectrum will equal the positive frequency spectrum in magnitude and opposite phase.

Since your signal is real, there will be two complex conjugate symmetric spectrums, one at $+f_c$ and one at $-f_c$. (When I took the real part of what was originally just a spectrum at $-f_c$, the companion spectrum at $+f_c$ is immediately created by the operation of taking the real part).

That said, to move your signal back to DC, simply multiply by $e^{-j2\pi f_c t}$ which will shift both spectrums to the left, and you will get one at DC and another that is at $-2 f_c$. It is possible here that the spectrum will be inverted, as I did not follow the math carefully in my processing described, in that case it would be better to multiply by $e^{-j2\pi f_c t}$, or simply swap I and Q and change the sign of Q which will also invert the spectrum. You then need to filter out the higher frequency signal that will be at twice the carrier frequency in order to get your baseband constellation only. This is the part I am surprised you would have to do as a homework problem; and if you had access to the noise samples at baseband this would be much easier but given the way the noise was added it doesn't appear that you do?

If you needed to plot the constellation with noise as 1 sample per symbol (IQ dibit), then you would simply take the average over a symbol period of the signal after converting back to baseband. Actually, by taking the average, the required low pass filtering will also be done (by nature of averaging is a low pass filter). So such a result for the constellation with noise is actually easier achieved when plotting one sample per symbol. Note too that taking such an average, if synchronized with the symbol period, represents the ideal matched filter for your example of QPSK with a rectangular pulse shape. This is often depicted on diagrams of receiver implementations as an integrate and dump filter (moving average where you accumulate or integrate over the symbol period, dump the result and reset the moving average).


Why is the correlation of two random markov chains so large? A very simple python code shows that the correlation coefficient and spearmann rank between two datasets of uniform random numbers drops proportionally to the square root of the number of points in the dataset. However, when one compares Markov chains produced by accumulating those numbers, there is no such behaviour, and even for large data sets, correlation can easily be like 0.7. What is the origin of such behaviour? Are two random markov chains more "related" to each other than two random datasets? Can this effect be corrected?

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

x = np.random.normal(0, 1, 10000)
y = np.random.normal(0, 1, 10000)

# Convert random data into Markov chains
for j in range(1, 10000):
    x[j] += x[j-1]
    y[j] += y[j-1]

print("Correlation: ", np.abs(np.corrcoef(x,y)[0, 1]))
print("SpearmannRank: ", np.abs(scipy.stats.spearmanr(x,y)[0]))

plt.figure()
plt.plot(x)
plt.plot(y)
plt.show()

EDIT: Unsurprisingly, the original problem is a bit more tricky. I have (at least) two experimentally obtained data sets I know relatively little about. Formally speaking, there is nothing stopping me from calculating their correlation in order to investigate their relationship. Since I was getting garbage, I tried to investigate correlations of several semi-random datasets, until I arrived at a surprising (for me) result above. Let me refine the question: is it possible to guess from the shape of the data, whether or not correlation is an informative metric?

Hi: you can't cumulate two sets of random numbers and then calculate a correlation between them. The two samples are non-stationary ( since you are cumulating them )  and what you get is referred to as "spurious correlation" or "non-sensical correlation" in econometrics. ( some would say statistics but the econometricians really figured out what was happening ).  I wouldn't do it justice trying to explain it here but, essentially, what you are doing is calculating two series whose means are changing over time so the correlation calculation is illegal statistically. Granger and Newbold figured out the statistics behind why this is in 1976 so if you want all the details, just google for "spurious regression". The econometrics literature on this is enormous and, even though what you are doing has maybe zero to do with econometrics, you are basically generating two series that are each random walks and then correlating them which is a garbage in garbage type of scenario.  I hope this helps.


DFT sample point k < N has negative frequency (From: Schaum's DSP outline, 2nd edition, page 254, problem 6.35)

A signal $x_a(t)$ that is bandlimited to 10 kHz is sampled with a sampling frequency of $f_s = 20$ kHz.  The DFT of N=1000 samples of x[n] is then computed, that is:

$$
X[k]=\sum_{n=0}^{N-1} x[n] e^{-j(2\pi / N)nk}
$$

with N=1000.

(a) to what analog frequency does the index … k=150, and k=800 correspond?

A few equations:

$$
\omega_k = \frac{2 \pi k}{N}
$$

$$
f_k = \frac{\omega_k f_s}{2 \pi}
$$

For k=150:

$$
\omega_k = \frac{2 \pi}{1000} 150 = \frac{3}{10}\pi
$$

$$
f_k= \frac{3\pi 20000}{2 \pi 10} = 3000\ \  hz
$$

Ok no problem so far:

For k=300:

$$
\omega_k = \frac{2 \pi}{1000} 800 = \frac{4}{5}\pi
$$

Here's where i hit a problem:

book says,  for k=800, we need to be careful. Because $X(e^{j\omega})$ is periodic:

$$X(e^{j\omega}) = X(e^{j\omega + 2\pi})
$$

Therefore, k=800 corresponds to the frequency:

$$
\omega_k = \frac{2\pi}{N}k = \frac{2\pi}{N}(k - N) = -200\frac{2\pi}{N}
$$

My question is this:

what's the reason we need to use a negative frequency here?  Since k=800 < N=1000?  Wouldn't it just be:

$$
\omega_k = \frac{2 \pi 800}{1000} = \frac{8\pi}{5} = 5.02...
$$

Is this just the book being strange or am i missing some understanding of DFT?  I would have thought you wouldn't need to worry about $2\pi$ unless k > 1000?

The DTFT frequency $\omega_k$ corresponding to the $N$-point DFT index $k$ is given by (frequency sampling relation)

$$ \omega_k = \frac{2 \pi}{N} k ~~~,~~~ k = 0,1,...,N-1.$$

The confusion might arise with the interpretation that the, on the continuous-time domain these frequencies correspond to:

$$  f_k = \frac{f_s}{N} k ~~~,~~~ k = 0,1,...,N-1 ?$$

But no. Because, after about $k = N/2$ which corresponds to $\omega = \pi$ and $f = f_s/2$, you will be sampling the negative portion of the shifted spectrum: $X(e^{j (\omega-2\pi)})$, which corresponds to negative discrete-time frequencies of $X(e^{j \omega})$ due to the periodicity of DTFT $X(e^{j\omega})$.

A more relevant (eventually the same) interpretation of the CTFT, DTFT, DFT frequency relation can be given by the following.

Let $x_c(t)$ be bandlimited to $-B \leq f <B$ and sampled at $F_s = 2B$ samples per second. The base period of the spectrum of $x[n]$ is between $-\pi \leq \omega < \pi$. Note that only the left member is included and the right member is discluded.

Now assume an $N$-point DFT $X[k]$ of $x[n]$ is computed. According to the frequency sampling interpretation of DFT, samples of $X[k]$ are given by :

$$ X[k] = X(e^{j \frac{2\pi}{N} k}) ~~~~ , k = ...-3,-2,-1,0,1,2,3... $$

Note very carefully, that eventhough most of the time we restrict our attention of the DFT index $k$ into the range of $0 \leq k < N$, indeed, the relation is valid for all $k$ which is a manifestation of DFT $X[k]$ and DTFT $X(e^{j \omega})$ being periodic by $N$ and $2\pi$.

Hence, instead of taking the range of $k$ in $[0,N-1]$, we can also take it in the range: $[-N/2 , N/2-1]$ (assuming $N$ even for the moment).

So, if you take $k = -N/2,-N/2+1,...,-1,0,1,2,...,N/2-1$ as the range of $k$ for DFT $X[k]$, then the mapping into the negative frequencies become obvious.

Note that, conventional FFT software always take the range of $k$ in $[0,N-1]$. Therefore, one has to relate them too which is simple: for $k$ in $[-N/2,-1]$ it corresponds to $k'$ in $[N/2,N-1]$ in the conventional DFT range, and for $k$ in $[0,N/2-1]$ it corresponds to the $k'$ in $[0,N/2-1]$ range in the conventional range.


how to interpret the 2D FFT I know how to compute the 1D  FFT  (and interpret values from 0 to Nyq).
When computing the 2D FFT,  do we compute the FFT of row[1]
then the FFT of row[2]  then the FFT of row[3] up to the last row.
And then compute the  FFT of col[1]   col[2]  for each of columns across ?
But to report the RESULT of the 2D FFT ( do we report the complex value of each final   row from col[0] to  col[N/2]   for each row from top to bottom ?

What INFO does the  2D complex FFT result convey ?

Yes, due to the separability of the kernel :

$$e^{-j \left(\frac{2\pi}{N_1} n_1 k_1 + \frac{2\pi}{N_2} n_2 k_2 \right) } = e^{-j \frac{2\pi}{N_1} n_1 k_1} \cdot e^{-j \frac{2\pi}{N_2} n_2 k_2}$$

the 2D-DFT sum
$$X[k_1,k_2] = \sum_{n_1} \sum_{n_2} f[n_1,n_2] e^{-j \frac{2\pi}{N_1} n_1 k_1} e^{-j \frac{2\pi}{N_2} n_2 k_2} $$

can be implemented using row-column (or column-row) decompositions.

$$X[k_1,k_2] = \sum_{n_1} \left( \sum_{n_2} f[n_1,n_2] e^{-j \frac{2\pi}{N_2} n_2 k_2} \right) e^{-j \frac{2\pi}{N_1} n_1 k_1} $$

The following MATLAB/OCTAVE code shows how to apply :

N = 8;            % length of columns-rows
x = randn(N,N);   % row data

X = fft(x,N) ;    % 1D-fft along columns of x[n1,n2]
X = fft(X.',N).' ;% 1D-fft along rows of intermediate X[k1,k2]


How to handle a logarithmic term in Kalman filter? I am trying to implement a Kalman filter for an echo pulse detection application as similar to this paper. (an open source version is here (pg 16))
The measurement variable is $h(x,t)=A_0 (\dfrac{t-\tau}{T})^\alpha  \cdot \text{exp}{(\dfrac{t-\tau}{T})}$ and the state vector is $x=[A_0 \; \text{ } \alpha \; \text{ } T \text{ } \; \tau \text{ }]^T$.
In order to compute the measurement sensitivity matrix $H_k=\frac{\partial h}{\partial x}\Bigr|_{\substack{x=\hat{x}_{k^-}}}$, we need to compute the partial derivative with respect to $\alpha$. I obtain it as $\frac{\partial h}{\partial \alpha}=h \cdot log(\dfrac{t-\tau}{T})$.

$\tau$ would correspond to the start of the envelope of the pulse so it will be greater than first couple of time steps practically speaking. When I iterate the Kalman filter initialized with a real $\tau$ for example, then this partial derivative term becomes undefined for the start values of iterations because at $t=0, 1  \; time step, 2 \; time step \cdots$etc and so on are smaller than $\tau$ and the logarithm is undefined for negative values.

First of all, you're trying to evaluate the derivative of an exponential. If the base of the exponential is positive, the derivative exists. However, if the base is negative, the derivative does not exist. It would be pointless to use a Kalman filter when you have an exponential with a negative base.
Secondly, your model of the echo envelope is
$$
A(t) = A_0(\frac{t-\tau}{T})^{\alpha}e(\frac{t-\tau}{T})
$$
However, since it is an echo, it does not make sense for the echo envelope to be something else than 0 for $t < \tau$ where $\tau$ is the time of flight
You could add a unit step to the equation to make sure that A(t) = 0 for t < $\tau$
$$
A(t) = A_0(\frac{t-\tau}{T})^{\alpha} e(\frac{t-\tau}{T})  u(t-\tau)
$$


Output of Zero-forcing equalizer doesn't match theoretical calculations Let's say I have a sequence $x[n]$ (signal) and perform convolution with another sequence $h[n]$ (channel), so:
$$
y[n] = h[n]*x[n] = \sum_{k}{h[k]x[n-k]} \leftrightarrow Y(z) = H(z)X(z)
$$

I'm looking for $C(z)$ such that:
$$
C(z)= \frac{1}{H(z)}
$$

So that I can compensate the channel:
$$
c[n]*y[n] = \sum_{k}{c[k]y[n-k]} \leftrightarrow C(z)Y(z) = \frac{1}{H(z)}Y(z) =
 \frac{1}{H(z)}H(z)X(z) = X(z)
$$

Let's assume I know the channel $h[n]$ and I use that to obtain $c[n]$ using Zero-Forcing equalization.

My question is, if I do:

x = [1 2 3 4 5 6 7];
h = [1 0.9 0.3];
c = zeroforcing(h); % c=[-4.7619, 4.2857, -1.4286]

y = filter(h,1,x);
x_d = filter(c,1,y);

Shouldn't x_d be a delayed version of x ? But I don't seem to get that in my simulations.

(The way I obtain c is converting h to a toeplitz matrix and doing $c[n]=H^{-1}\delta [n]$ as explained here:

$$
\begin{bmatrix}
c_1 \\ c_2 \\ c_3
\end{bmatrix}
=
\begin{bmatrix}
h_2 & h_1 & 0 \\
h_3 & h_2 & h_1 \\
0 & h_3 & h_2
\end{bmatrix}^{-1}
\begin{bmatrix}
0 \\ 1 \\ 0
\end{bmatrix}
=
\begin{bmatrix}
0.9 & 1 & 0 \\
0.3 & 0.9 & 1 \\
0 & 0.3 & 0.9
\end{bmatrix}^{-1}
\begin{bmatrix}
0 \\ 1 \\ 0
\end{bmatrix}
$$

It seems like there might be something weird going on with that blog, but I only quickly looked. Here is the idea for the ZF equalizer:

Figure out the convolution matrix.

convMatrix = convmtx(h, length(x));

This will be a $9 \times 7$ matrix because the input vector $\mathbf{x}$ is length 7 and the output convolution will be of length $7+3-1=9$. It is not a square matrix but we can still get an inverse by computing the Moore-Penrose inverse (https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse), which in MATLAB can be done by calling:

invConvMatrix = pinv(convMatrix);.

Choose your delay. invConvMatrix is a $7 \times 9$ matrix and we are multiplying with an all zero vector except for a single one at some delay $\mathbf{\delta}$. Lets choose the delay to be 5, so we have:
$\mathbf{\delta}=[0, 0, 0, 0, 1, 0, 0, 0, 0]^T$.
Compute the ZF equalizer:

 zfEqualizer = invConvMatrix * d;

Generate the received signal, and generate signal after equalization

 y = conv(h, x);
  z = conv(zfEqualizer, y);

Now lets take a look at the resulting signals, and specifically notice where the equalized signal starts, right at the delay of 5 that we choose earlier. I hope this helps.


Coefficient Scaling of IIR filter to obtain unity gain response I am designing a 2nd order IIR digital filter :
My tf equation with coefficients is
b = [0 1.209e09]
a= [9.2175   -2.6952    1.0000]
sys=tf(b,a,0.1,'Variable','z^-1')
bode(sys)

I have couple of questions :

How to get unity scaling .. as a formula... so that i can always multiply it with the 'sys' tf.
I usually seen that coefficients are less than 1 , is there a way i can make them so ?
My coefficients are calculated from some adaptive filter optimization algorithm, so they keep on changing, is there a way i can reduce the filter to filter coefficients changes to a a minimum.?

It's important to specify at which frequency you want unity gain. But assuming you mean DC ($\omega=0$), because that filter has a low pass characteristic, the DC gain of an IIR filter is given by
$$G_{DC}=\frac{\sum_kb[k]}{\sum_ka[k]}\tag{1}$$
It's also common to normalize the denominator coefficients such that $a[0]=1$. In your example that would give
a = [1.00000  -0.29240   0.10849]
and
b = [0.00000   131163547.59967]
Finally, normalizing by the DC gain $(1)$ will give you a filter with unity gain at DC. This leaves the denominator coefficients unchanged, and the new normalized numerator coefficients are given by
b = [0.00000   0.81609]


Good Heuristic Point Where Sending FFT to GPU is Computationally Worth It? This is a rather simple question, but I haven't found much posted within the last 5 years or so regarding a good rule of thumb for computing an FFT on the GPU vs the CPU. I know that most of the overhead of GPU-calculated FFTs is that you have to first write the data to a buffer and then transmit that data to the GPU, from the CPU, in order to compute the transform.

This implies naturally that GPU calculating of the FFT is more suited for larger FFT computations where the number of writes to the GPU is relatively small compared to the number of calculations performed by the GPU. But the issue then becomes knowing at what point that the FFT performs better on the CPU vs GPU. The best I've found is on the lines of "when you're computing larger FFTs", but that's a little relativistic to be particularly meaningful guideline for practitioners, especially considering that GPU technology has been accelerating so rapidly in the past few years.

I'm just looking for a more concrete estimate of at what point it becomes efficient to compute the transform on the GPU vs CPU (perhaps there's a benchmark somewhere out there?). Thanks in advance for any help you can provide.

I'm afraid there won't ever be a rule of thumb!

The reasons are manifold, mainly, that both the systems you're considering and the problems you're trying to solve vary across a very large range.

Problem-based aspects

You say you want to do an FFT – but that's always but half of what you actually want to do!

Need the FFT converted to it's abs², then mapped to colors, then displayed on a screen? Do it in the GPU, it's right where it belongs; fosphor does that, at an easy 200MS/s on capable PC/GPU combinations:

In that case, the size of the FFT doesn't even matter. Your data is going to be processed further on by the GPU, so do the FFT there.

On the other hand, you might want to do something that depends on a lot of checks on individual elements in the FFT, in the CPU? Possibly just one FFT, and then not many more?

In that case, your theoretical throughput doesn't help you at all. Just waiting for the data to get out of your CPU cache, back into coherent memory, so it can be DMA'ed to the GPU, where you then initiate an FFT (possibly wasting an interrupt/context switch on the way), just to wait till it's done, the GPU DMA'ed the data back to your main memory and you got it into your CPU's cache won't pay, even for medium sized FFTs.

So: this whole "high-latency math accelerator business" really only ever pays off if you can do something sensible while you wait. If you can't, there's a huge latency penalty.

System-based aspects

Ok, not going into too much detail here, but:

DSP systems are either CPU or Memory-bandwidth limited
If your GPU operation helps with CPU limititation, but puts additional data movement load on the memory interface, while in fact the rest of your system is memory-bandwidth limited, you're hurting yourself.
Same applies the other way around: Maybe your algorithm (the FFT in your specific size of interest) is CPU limited, but your GPU acceleration leads to additional interrupts
What is an FFT size that your CPU can do very well? That is probably defined by the sizes of it's L1 and L2 caches. A Xeon number-crunching CPU will have dozens of Megabytes of those, whereas an ARM running in a Jetson NVidia SoC will not.
What is the FFT size your graphics card is good at? There's a humongous difference in numbers of parallel threads, their flexibility and memory bandwidth across cards.
What is a metric for "good", at all? Just some strange ratio of throughput and latency, but maybe also energy and leaving the right kind of ressources free for other work?
What is your CPU<->main memory interface? Is it a quad-channel DDR4 interface running at nearly 2GHz, or is it single-channel DDR?
What is your GPU<->GPU memory interface?
What is your GPU<->main memory interface?
How well does your CPU<->GPU communication work for your specific use case?
Is there a high load on e.g. the PCIe bus, for example, because the same PCIe switch has to handle the data flowing in and out of your high-rate system (eg. storage, but more likely, 10 Gigabit ethernet or video data)?

So, the answer probably won't be satisfactory, but it really is:

Somewhere above 64 bins, somewhere below 220 bins, for a single-precision FFT. It depends.


How can the order of a transfer function be derived from its equivalent state space representation? Suppose I have a discrete state space model:

\begin{align}
\theta[k+1] &= A \theta[k] + B u[k]\\
y[k] &= C \theta[k]
\end{align}

I know that the equivalent transfer function can be found by solving

$$
Y(z) = C (zI-A)^{-1} B U(z)
$$

But is there a quick way to determine the order of the transfer function for each input? I.e., the degree of its numerator and denominators, assuming that the dimensions of $A$ will be $n\times n$, $B$'s will be $n \times m$ and $C$'s will be $1 \times n$? (Assuming here there is only one variable being observed.)

One trick will be whether there are any pole-zero cancellations in the overall transfer function.  For that, you will probably need to check for the minimality of the system (i.e. check whether the system is both controllable and observable).

The question of "determine the order of the transfer function for each input" seems odd, to me.  If the system is minimal (it can't be represented by a smaller order system), then you need it as-is.  Pulling out one of the $m\times 1$ systems on its own may be of interest, but you can't reduce it without destroying the overall system.

Can you give more information about why you want to examine just one of the $m$ transfer functions?


How to translate "phase" into "delay" Suppose that a system induces a phase of -90 degree, how do I practically measure the delay in the signal coming out of the system? I think time lag would be a more useful measure than phase.

Clearly these have different units, what is the conversion factor that tells me how much time lag is induced due to the system?

The phase delay to time delay conversion depends on the frequency of the signal for which that phase delay is associated.

Note that unless the phase delay is given as correctly "un-wrapped", the absolute time delay could be ambiguous by integer multiples of the period of the frequency at which the phase delay was specified or measured.


MATLAB rand and C++ rand() I am trying to port a MATLAB code into C++ and found the usage of rand in MATLAB.

Is MATLAB's rand function and C++ rand() function same? If not, is there any inbuilt function for MATLAB's rand in C++ or OpenCV?

It is a programming question, not a signal processing question - keep that in mind, and next time use the StackOverflow.

MATLAB rand function returns uniformly distributed variable on interval $[0,1]$.
C++ rand function returns integer between 0 and RAND_MAX, additionally it is not very 'random' function. You would also have to normalize its output. For more info please refer here: C++ rand() between 0 and 1. Although I think that you should use the uniform_int_distribution function, providing you can take the advantage of C++11 standard.

I believe that you should use OpenCV after all, instead of playing with new standards and some tricky functions. It contains the RNG class that you can easily use for generating of uniformly distributed sequence. Here is the example:

RNG rng;
float r_num = rng.uniform(0, 1);


What does 0 mean in YUV420 I know whats the difference between various formats such as YUV444,YUV420,YUV422.
But what exactly these numbers mean? They are mentioned as chroma subsampling in horizontal and vertical directions. In YUV420 Cb and Cr are sampled at half the resolution of luma in both vertical and horizontal directions. So what does zero in YUV420 mean?

As explained on the Wikipedia Page on Croma Subsampling:

  The subsampling scheme is commonly expressed as a three part ratio ''J'':''a'':''b'' (e.g. 4:2:2) or four parts if alpha channel is present (e.g. 4:2:2:4), that describe the number of luminance and chrominance samples in a conceptual region that is ''J'' pixels wide, and 2 pixels high. The parts are (in their respective order):

  ''J'': horizontal sampling reference (width of the conceptual region). Usually, 4.
  ''a'': number of chrominance samples (Cr, Cb) in the first row of ''J'' pixels.
  ''b'': number of changes of chrominance samples (Cr, Cb) between first and second row of ''J'' pixels.

  …

In other words, your 0 says there's no vertical chroma subsampling at all. Instead, you just assume that columns within the four pixel wide subsampling region have identical chroma.


possibilty of extracting a tone from MIDI values I have a project of recording human's sound via Android's SoundRecoder API, and then displaying the extracted note on The screen. It is specific for Android platform (can't use Jmusic or Jfugue API as javaSE environment is not supported).

Now suppose I get to record a sound using this API . How can I then deploy FFT to extract notes from this sound?  Basically as it is not pure sinusoidal wave it might involve multiple frequencies and overtones. Are there any good algorithms to do so? If yes, where I can simulate and the last of all whole of the extracted tone should be rendered on musical sheet?

You correctly note that a musical note is usually not just a simple sine wave. Thus, an FFT may only be one of the steps used in a pitch detection/estimation algorithm (there are many: http://en.m.wikipedia.org/wiki/Pitch_detection).  For music, a good place to start might be the research papers published on MIREX: http://www.music-ir.org/mirex/wiki/Audio_Melody_Extraction


Transfer function determination from input and output data I have some input and output data that I believe adequately includes excitation of the important dynamics of a system. I know it is at most a 4th-order transfer function.

How can I identify the transfer function? I have Python available. (I have MATLAB but do not have access to the System ID Toolbox and need to limit my dependencies on MATLAB)

I know how to do least-squares fits, but I have no idea how to apply the idea to transfer functions, where using the absolute error in the transfer function doesn't make sense (relative error seems to be the appropriate thing to fit).

I'm vaguely familiar with ARX/ARMA/ARMAX approaches but they always seem to overemphasize the high-frequency response, because they solve for the impulse response.

if your input excitation signal (call it $x[n]$) is arbitrarily long and so also your output $y[n]$, then you should be able to take as long of block (let's call that length $N$) of each and, using the DFT and division compute
$$ H[k] = \frac{Y[k]}{X[k]} $$

in audio, we call that method of system identification the "dual-channel FFT".  the DFT of $h[n]$ is equal to the the DTFT of $h[n]$ (wrapped-around and overlapped and added every $N$ samples) which is the Z-transform of the same wrapped (and aliased) $h[n]$ sampled at $N$ equally-spaced points around the unit circle.
$$ H(e^{j 2 \pi k/N}) = H[k]  $$

so if your 4th-order IIR transfer function is
$$ G(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + b_3 z^{-3} + b_4 z^{-4}}{1 + a_1 z^{-1} + a_2 z^{-2} + a_3 z^{-3} + a_4 z^{-4}} $$

what you want to do is choose $ \{ b_0, b_1, b_2, b_3, b_4, a_1, a_2, a_3, a_4 \} $ in such a way to minimize
$$ \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H(e^{j 2 \pi k/N}) \right|^2 \ = \  \sum_{k=0}^{N-1} \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$
or, perhaps a weighted version of that,
$$ \sum_{k=0}^{N-1} W[k] \left| G(e^{j 2 \pi k/N}) - H[k] \right|^2 $$
where $W[k]$ is known.


Difference in information between IQ Sampling vs. Real Sampling at $2F_s$ As per my current understanding we use IQ sampling as way to increase the useful bandwidth to Fs rather than $F_s/2$ when viewing the FFT because the IQ sampling removes the negative frequencies and the entire $0-F_s$ band contains useful frequencies and not images of first half spectrum.
My questions are:

Is there any difference between this and sampling at twice the original sample rate?
And if I am sampling at $2F_s$ can that data be used to create the same IQ data as sampling I and Q components individually at $F_s$ or would there be a difference?

The bandwidths might be the same.  I agree with Jason’s assessment.

Is there a difference?

Yes.  Let’s say we construct both signals such that their frequency domain’s are equivalent from $0-F_s$.  For the real signal the frequency domain would be limited to $X(F_s-f) = X^{*}(F_s+f)$.  For the complex signal it would be $X(f) = X(F_s+f)$.

Can the IQ signal be recreated from the real signal?

Yes.  Again, we assume the frequency domains are equal from $0-F_s$.  First, we filter the real signal to remove the frequency content from $F_s-2F_s$, which would result in an interpolated version of the IQ signal. Then we downsample by throwing away every other sample, yielding the IQ signal.  There are practical considerations, but that would be one way to go about it.


Difference in Interpretation: $ω$ (rads/s) vs. $ω$ (rads) and $X(ω)$ vs. $X(e^{jω})$ The fourier transform of a continuous time signal $x(t)$ is $X(ω)$ where the unit of $ω$ is radians/second. And for a discrete signal $x(n)$, the DTFT is $X(e^{jω})$ where the unit of $ω$ is radians. Clearly the $ω$ for FT and DTFT shouldn't mean the same thing.

Question

What is the difference in interpretation of $ω$ (rads/s) and $ω$ (rads)?

If we interprete $X(ω_i)$ as the ampitude or $|X(ω_i)|^2$ as the energy associated with the frequency component $ω_i$, then how should we interprete $X(e^{jω_i})$ or $|X(e^{jω_i})|^2$?

Thanks in advance.

Most people learn (or try) DSP from this book

  Oppenheim, and Schaefer, (1999). Discrete-time signal processing. Pearson Education.

and there $\Omega$ is used for continuous time and $\omega$ is for discrete time, which is how they keep their notation consistent. Usually both are continuous, which can lead to some confusion when you consider the DFT and/or Fourier series.

In the Olden Days, Yea DSP was primarily practiced by people who had a background in continuous time signals derived from some physical problem, and notational conventions for $\omega$ were typically $2\pi f$ (in yonder rads/sec).

Most books back then, tended to start with some general Linear Systems theory, and then go to sampling.  I think the most recent version of Oppenheim and Schaefer wait until Chapter 4 to cover sampling.

People who taught and wrote about DSP, recognized that including a sample rate wasn't really necessary.  My first DSP class was taught by a guy who did adaptive beamforming and he introduced the concept of unity sampling, $f_s=1 \text{ sample/second}$ so all our filter frequency responses were plotted between -1/2 to 1/2 Hz.  You just needed to scale everything by your actual sample rate when you got back to your physical problem.

It became clearer over time that there were time series, like Stock Market time series that were inherently discrete in time and a lot of DSP didn't really require reference to time, and things like image processing really didn't either, so to make things more general, $\omega$ as radians seems to have become the most common notation. You can also sell your books to finance majors now.

The notation $X(e^{j\omega})$ is handy when you plug $e^{j\omega}$ into a Z transform (if the unit circle is in the region of convergence), and reinforces the concept that a function or composition of a periodic function $e^{j\omega}$, is periodic.

The notation $X(\omega)$ needs to be interpreted in the context of where it is, could be discrete time, could be continuous time.

So, the answer to your question is: Context.

As an American English speaker, it would be natural for me to assume that this is the complete story but someone who learned from the Russian (or French, or ....) DSP literature would most likely have a different explanation.