Amplitude modulation in pure sine wave I am generating sine wave in MATLAB with following code

fs=2000e6;%sampling frequency
t=0:1/fs:(1e-6)-1/fs; %1us seconds duration
fc=501e6;%frequency of the cosine wave
x=cos(2*pi*fc*t);
plot(x)

The wave it generates is somehing like this

My question is: Why do we see amplitude modulation in the wave. Shouldn't amplitudes be 1 at all peaks of the wave? like this wave

To explain what you are seeing (the envelope you see specifically is not a plotting artifact, but a fundamental DSP concept worth explaining) consider if you instead sampled a 500e6 Hz cosine wave at your sampling rate of 2000e6 Hz, which is exactly 4x the sampling rate. This would result in samples 1, 0, -1, 0, 1.... and since Matlab as you plotted would connect the adjacent samples, it would look like the bottom picture. (And if you used plot(x,'o') you would see just three lines).

Now consider your case where you have offset the frequency by 1 Hz, with this the sample location will "roll" along the sine wave, one cycle in the 1 us duration you are using.

This is exactly what you are seeing. The magnitudes 1 and -1 given in the first case are now moving slowly as the sample locations change slowly along the plot due to the slightly higher frequency, and if you use plot(x, 'o') for that case you will more clearly see the sample locations I am speaking of rather than the lines connecting the adjacent samples that are positive and negative (you will see what appears at 4 shapes varying sinusoidally and in quadrature. The envelope is just the artifact of the relatively low sampling rate (you have slightly less than 4 samples per cycle)- not how Matlab plots it.


Why is a wavelet transform implemented as a filter bank? The mother wavelet function $\psi(t)$ must satisfy the following:

$$\int\limits_{-\infty}^{+\infty}  \frac{|\psi(\omega)|^2}{\omega} d \omega < +\infty,$$
$$\psi ( \omega ) \bigg|_{ \omega =0} =0,$$
and
$$ \int\limits_{-\infty}^{+\infty}  \psi(t) \ dt = 0$$

To serve as the wavelet basis for wavelet transform
$$ \gamma (s, \tau ) = \int\limits_{-\infty}^{+\infty} f(t) \ \psi_{s, \tau }(t) \ dt$$

where $ \psi_{s, \tau }(t) \triangleq \psi\left( \frac{t-\tau}{s} \right)$.

While I understand that the wavelet must be an oscillatory function having no frequency component at $\omega=0$ and effectively have a band pass filter like spectrum, from the equation of wavelet series or wavelet transform can you tell me why is it that the wavelet transform is implemented as a filter bank? What is the intuition behind it? What makes it possible?

I am asking this question since the fact that practically the DWT is implemented as a filter bank means that it is not a DWT anymore, it is just a set of low pass and high pass filters. It is mind bogling.

First of all the basic idea of wavelet transforms lies in multi-resolution analysis. What this means is that the signal is looked at from different scales.

It is probably easier to understand this with images (which are 2D signals). The idea of multi-resolution is like zooming in and out of a reference signal and compared with windows of the image according to the amount of zoom.

Similarly for 1D you are trying to see how the mother wavelet at different scales (think of it as compressing and expanding in the 1D) compares with the signal at different delays (at different points on the discrete 1D).

Now when you look at a given set of discrete points you want to transform based on a mother wavelet it easy to see them as filters (based on mother wavelets) that are applied at different scales.

For a better understanding of wavelets and DWT I suggest you read Polikar's tutorial : Wavelet Tutorial


Filtered backprojection implementatoin I am wondering the detail of the filtered backprojection.

Suppose I have a projection file:

$p(s,\theta)$=[1 2 3 1 2]

Then to get its filtered backprojection version using ramp filter, I should

apply FFT to it, to get its Fourier version $P(\omega,\theta)$
To use |$\omega$| to multiply $P(\omega,\theta)$
Backprojection

I am a little confused at second step here.

Using the numeric example above, to apply FFT to [1 2 3 1 2], we will get a complex sequence

   9.0000 + 0.0000i
  -1.0000 - 1.1756i
  -1.0000 + 1.9021i
  -1.0000 - 1.9021i
  -1.0000 + 1.1756i

Should I apply |$\omega$| to both real and complex part  of the sequence above, then do inverse FFT? Or should I only take the magnitude of above sequence, apply |$\omega$| to the magnitude, then do inverse FFT?

Either way, it seems I will get a complex sequence as the final results, but shouldn't the filtered projection be a real sequence? Anyone give me some pointers here? Thanks

I did it by creating a convolution kernel, then FFT the kernel, FFT the projection file, then multiply these two, then do IFFT to the product. Then the result of the IFFT is the filtered backprojection.


Variance of White Gaussian Noise It could seem an easy question and without any doubts it is but I'm trying to calculate the variance of white Gaussian noise without any result.

The power spectral density (PSD) of additive white Gaussian noise (AWGN) is $\frac{N_0}{2}$ while the autocorrelation is $\frac{N_0}{2}\delta(\tau)$, so variance is infinite?

White Gaussian noise in the continuous-time case is not what is called a second-order process (meaning $E[X^2(t)]$ is finite) and so, yes, the variance is infinite. Fortunately, we can never observe a white noise process (whether
Gaussian or not) in nature; it is only observable through some kind of device,
e.g. a (BIBO-stable) linear filter with transfer function $H(f)$ in which case
what you get is a stationary Gaussian process with power spectral density $\frac{N_0}{2}|H(f)|^2$ and finite variance
$$\sigma^2  = \int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,\mathrm df.$$

More than what you probably want to know about white Gaussian noise
can be found in the Appendix of this lecture
note
of mine.


Detect marker orientation inside an image I've been reading how augmented reality works and how markers are detected. Anyway I can't find any resource that would help me in detecting an arbitrary complex image as, for instance, a book's cover.

I don't think that falls into the "template matching" field since I would also need its orientation to perform AR rendering on it.

What algorithms are used to detect an arbitrarily complex marker into an image and have its transformation extracted as well?

In essence, you are looking for a way to do planar object detection.

Well, the most basic approach is to find a correspondence between your planar template and the image using a feature extractor, such as SIFT. You could then cast the detection problem as finding the subset of points, which have similar descriptors. OpenCV has that implemented already:
http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_feature_homography/py_feature_homography.html

http://docs.opencv.org/trunk/dc/dc3/tutorial_py_matcher.html

The next thing is to use template matching kind of approaches. For severe viewpoint this might be hard to achieve, but there are ways, such as this one.

If you can also design a marker in the image space, then there are approaches which let you to read a randomly generated marker. For example:

  T. Birdal, I. Dobryden, S. Ilic  X-Tag: A Fiducial Tag for Flexible
  and Accurate Bundle Adjustment  International Conference on 3DVision
  (3DV), Stanford University, California, USA, October 2016
  http://campar.in.tum.de/pub/tbirdal20163dv/tbirdal20163dv.pdf

Take a look at the related work section of the references I gave here. This could broaden your knowledge in the field. You might also like to read a survey in the field, such as:
http://www.irisa.fr/lagadic/pdf/2012_fcv_uchiyama.pdf


Convolution output signal is very large, do I need to scale inputs or outputs? I am writing some simple matlab code to refresh my convolution knowledge (it has been awhile) in regards to DSP.
I have an input signal x and a filter -
x = [1 3 6 2 8];
h = [2 9 5 3];
y=conv(x,h);

which results in
y = [2 15 44 76 73 100 46 24];

I understand these values are large because convolution has lots of multiplication and addition, however
if this was an audio signal I was trying to filter, these values seem very large.
Do I need to scale the input filter or output signal some how?  If so which one and by how much?

In practice you would make sure that your filters don't have very large gains, at least not over a large frequency range. The filter you're using in the given example has a maximum gain of almost $26$ dB (at DC), and the gain is no less than $12$ dB for any frequency. I.e., you're applying a flat gain of $12$ dB, and on top of that you boost a relatively wide band around DC by another $14$ dB. This is not what a typical filter for audio signals would do, unless you make sure that the input signal is small enough and the additional gain is really needed.


Phase spectrum of the DTFT of a rect pulse using FFT I'm trying to use the FFT to plot the phase spectrum of the DTFT of a sampled rect impulse function.
Here's the Matlab code I've written.

fs=10; % sampling frequency Hz
Ts=1/fs;
T=4; % length of time vector in seconds
t=-T/2:Ts:T/2-Ts; %vector of 40 time samples -2:0.1:1.9
y=rectpuls(t); %unit length rect impulse
z=fft(y,10000).*Ts;
%fft multiplied by Ts to get rid of the 1/Ts factor in the spectral replicas
%of the DTFT of the signal
f=linspace(0, (length(z)-1)/length(z),length(z));
figure
plot(f,abs(z));
grid on

%plot of the phase of z as it is, so with the time sequence
%considered to be shifted by length(y)/2 samples by fft.
figure
plot(f,unwrap(angle(z)),'r');
grid on

%I try to get rid of the linear phase delay component to get the phase
%spectrum of the DTFT of the signal in its frequency period.
%I divide each fft sample by its own exponential delay factor.
delay=length(y)/(2*length(z));
expd=zeros(1,length(z));
for index=1:1:length(z)
expd(index)=exp(-1i*2*pi*(index-1)*delay);
end
z2=z./expd;
figure
plot(f,unwrap(angle(z2))./pi,'r');
grid on

Now the problem is that this phase plot is wrong, the DTFT is the sum of real sinc functions so its phase spectrum con only assume k*pi values.
There's still a linear component in the phase plot and the slope is positive so the exponentials that I'm using are over compensanting the delay.
I have found out that the problem is that the fft function introduces a linear phase delay as if the sequence was shifted by 19.5 samples instead of 20 samples as it should be in my example.
Infact if I set in my code delay=19.5/10000 I get the correct phase plot (or at least one that makes sense :D):

So what's going on here, is it a precision problem?
Are there mistakes in my Matlab code and/or my dsp theory formulas?
Thank you!!

Unfortunately I can't link more than two images in the post :/

In order for the result of an even number N length FFT to be strictly real, the input vector data vector has to be exactly symmetric around sample N/2, which is the sample point before t=0, given your sample window, which is offset slightly to the left.


Kalman filter for random observation matrix $G_t$ I have a problem that is similar to the state space model in Kalman filter but the observation matrix $G_t$ of $$y_t=G_tx_t +w_t,$$ is random. The elements of $G_t$ are i.i.d. random variables with a given distribution. Is it possible to use the Kalman filter to this problem? If it has already been used please provide a reference.

Above $w_t$ is noise vector of Gaussian distribution with known covariance matrix, $x_t$ is the state of the system at time $t$ and $y_t$ is the observation at time $t$. The state tarnsision is given by $x_{t+1}=F_tx_t+B_tu_t+v_t$ where $F_t$ is a known matrix and $v_t$ is Gaussian noise with known covariance, $B_t$ is also a known matrix and $u_t$ is Gaussian distributed random input.

Thanks a lot.

Very interesting question!

I do not know of any reason why you can't do this, but I have had a hard time finding an example.  I finally found this paper, which does an analysis of what happens when the measurements are randomly lost due to transport delays / errors in communication systems.

The main problem that they report is that the error covariance matrix iteration is now stochastic rather than deterministic (as in the usual Kalman filter formulation).


DCT step suppress high frequencies in MFCC extraction process I'm trying to get MFCC from speech audio data. I'm using 26 Mel Filter on 0-24kHz spectrum. My MFCCs doesn't look correct to me, cause always first bins have high values and values getting smaller through the high frequency bins, even with high frequency audio data.

Then I realized that the cos((PI/N)*(n+0.5)*k) part in the DCT Type II formula, causes this result. Here is how the value of that cos() part changes by increasing the number of mel filter bin:
http://img402.imageshack.us/img402/5994/dctresult.png

When I multiply the mel filter bank energies with that cos() part, the MFCCs have only low frequencies and high frequencies are getting suppressed. MFCC result of a "s" sound has a frequency of 5khz-10khz:
http://img9.imageshack.us/img9/4279/mfccresult.png

Here is my code of DCT(Discrete Cosine Transform) step:

numofMelFilters = 26;
numofMFCCs = 26;    // This ll be 12..
for(int k=0; k<numofMFCC; k++){
   for(int n=0; n<numofMelFilters; n++){
   MFCC[k] += melEnergies[n] * cos((PI / (double)numofMelFilters) * (n + 0.5) * k);
   }
}

I can't find out what I am doing wrong, any help will be appreciated. Thanks in advance..

I think you misinterpret the output coefficients of the DCT. If you want to look at the frequency content of your input signal, you need to look at the output of the triangular filter bank, before the DCT. So if there is a fricative, like an "s" sound, the outputs of the high frequency filter bank kernels will have a higher energy than when the input is a vowel. However, after the DCT you do not have this type of frequency information, because you take the DCT of a frequency-domain signal, not of a time-domain signal. So the "high-frequency" bins of the DCT output correspond to fast variations in the shape of the spectrum, not to high frequencies in the input (time-domain) signal. The very effect you see, is exactly what the DCT is supposed to do. It compresses the input data, in the sense that almost all information is contained in the lower DCT bins. Usually the higher DCT bins are not used, which in fact means that you are doing lossy compression of your data. Luckily it turns out the the information that is being thrown away is actually irrelevant for the performance of the speech recognizer, and neglecting this information makes the recognition even more robust.

Having said that, of course there can still be a bug in your code, but the effect you describe is actually a desired result and nothing to worry about.


Decimation with IIR filter altering amplitude of signal I'm downsampling a signal through decimate (scipy implementation), but the result has a different amplitude than the original signal, which if I understood correctly, shouldn't happen. The original signal is sampled at 128Hz, and I want to downsample it to 1Hz. Decimation code is:

    sig_d = signal.decimate(sig, 8)
    sig_d = signal.decimate(sig_d, 8)
    sig_d = signal.decimate(sig_d, 2)

I then checked in the first 1 minute window of the original and decimated signals:

    df = pd.DataFrame(sig[0:7680])
    print(df.describe())
    df = pd.DataFrame(sig_d[0:60])
    print(df.describe())

and this is the output statistics:

          Original
count  7680.000000
mean     94.191307
std       0.551882
min      92.580971
25%      94.028488
50%      94.086906
75%      94.794744
max      94.969666

       Decimated
count  60.000000
mean   91.011716
std     0.516294
min    89.658212
25%    90.847055
50%    90.886592
75%    91.589404
max    91.689528

As you can see, there is a huge decrease in amplitude for the same time window in the decimated signal. Per the statistics, it doesn't seem to be related to peak removal, since this window in the original signal has very small variability, and not really any pronounced peaks. I did plot the two signals to check if I wasn't messing up the alignment, and they were aligned perfectly, their shape in this time window is the same (besides sampling rate), but the decimated signal seemed to be "pushed down" in amplitude. Unfortunately, I cannot post the plots (for proxy reasons).

This seems to happen only with IIR filter(though FIR has other issues). Is this normal behavior of decimation? Or am I doing something wrong?

First, 0.3dB is generally considered small potatoes in signal processing.

Second, the documentation for scipy.signal.decimate states that it uses an 8th-order Chebychev filter by default (without specifying the ripple, dangit!).  Chebychev filters have some ripple in the passband; if the bulk of your signal's energy is falling into the troughs of the ripple, then you should expect some attenuation.  In particular, it appears that your filter has a lot of DC content, and even-order Chebychev filters tend to have a gain of less than unity at DC.

So I suspect that it's doing exactly as expected.


At what stage do we compute the approximations and details while performing a DWT? $$
\int_{-\infty}^{+\infty}{f(t)\psi_{j,k}^\ast(t)dt}\ \textrm{with}\ \psi_{j,k}(t)\ =\ a_0^{-j/2}\psi(a_0^{-j}t\ -\ b_0k)
$$
If this is the expression for the wavelet transform, how does this lead to wavelet decomposition and at what stage do we compute approximations and details?

This expression is more a discretization of a continuous wavelet transform than an actual DWT (discrete wavelet transform), provided $\psi$ is a genuine wavelet. It only computes the wavelet coefficient $c_{j,k}$ associated to a specific shifted and dilated continuous wavelet $\psi_{j,k}(t)$. This yields a frame-like wavelet decomposition, if you picture all possible $c_{j,k}$. A frame is a kind of overcomplete basis satisfying some bounds; in other words, more projection vectors than the dimension of the space, still ensuring a certain stability.

It is not a discrete wavelet transform per se, since time is still continuous. This is a bit like Fourier series. Approximations and details for DWT are more classically obtained from a direct discrete scheme. From a discrete signal $x[n]$, at one single level, it is possible to build two downsampled sequences:

$a(n) = (h_0\ast x)\downarrow(2)[n]$ (approximation)
$d(n) = (h_1\ast x)\downarrow(2)[n]$ (details)

with $h_0$  and $h_1$ complementary (low-pass and high-pass) filters, such that one can recover $x$ from half-length sequences  $a(n)$ and  $b(n)$ only. Coefficients for $h_0$  and $h_1$  can be obtained from certain wavelet shapes only. Passing the signal in a filter before downsampling reduces aliasing artifacts, while preserving invertibility.

Approximations and details (discrete wavelets and wavelet packets) can be obtained by iterating the above at different levels.
Stationary wavelets are a version of DWT without subsampling (but generally the same wavelets). They often behave better in case of higher noise,for denoising, or model mismatch in adaptive filtering (some details in A Primal-Dual Proximal Algorithm for Sparse Template-Based Adaptive Filtering: Application to Seismic Multiple Removal).


Extrapolation after taking inverse Fourier transform Suppose I have 2000 point(time series) signal and I take 512 point FFT, after this some operation is done in frequency domain. Now if I take IFFT to get 512 points in time domain, Is there a way to map the output 512 IFFT points to the original time domain 2000 points.

Inverse  of this problem would be 400 point(time series) data and 512 point fft, and finally 512 IFFT points, here we can use Overlap and Add method to get 400 point output, is there any similar method for above problem.

An FFT possesses linear features. Let us use a linear analogy. If one tells you that the sum of two numbers is $a+b=3$, what can you say about these numbers? Not much. More generally, $n$ linear equations on $m$ variables is underdetermined when $n<m$: you do not have enough information to retrieve the variables without more conditions. This can be experimented with a very simple DFT case.

Let us assume the vector $x = [1,-1]^T$. A simple two-point DFT is defined by the matrix:

$$ \begin{bmatrix}
 +1 & 1\\ 1 & -1
\end{bmatrix}$$

Its products with $x$ yields either $0$ (for the DC) or $2$ (up to a  $\sqrt{2}$ Fourier definition scaling) at Nyquist. Can you recover the two elements from $x$ from only a single coefficient $0$ or (xor) $2$?

Probably not, without further assumptions, or on very specific signals.
This happens whenever you have less DFT points than samples. There is no bijection between two finite sets of different sizes. This is indeed one definition of infinity: a set is infinite if it can be put in bijection with a proper subset.

However, using sparsity, positivity properties, there are frameworks to retrieve (at least approximately) the data.


Why do we use discrete-time sampling? I have been reading Signals and Systems (2nd ed.) recently. Chapter 7 is about sampling. Oppenheim uses discrete-time signals and continuous-time signals to explain sampling.

But why do we need discrete-time signals? We use the A/D to receive and convert the continuous-time signal, so the signal is already discrete. In other words, the discrete signals are the sampled version of continuous signals.
Is there any practical meaning to sample discrete-time signals?

There is a difference between a quantized signal, i.e. taking only integer values, and a sampled signal, i.e. the value of which is only looked at on spaced instants.

A physical signal is usually continuous in values and in time.

An A/D converter turns a physical signal in a quantized one; if the converter is asynchronous, it generates a time-continuous signal.

Discrete-time/continuous-value signals are created by sample & hold circuits.

And digital signals processed by computers are usually both quantized and time-discrete.


How to discard deviating values from standard deviation and running average I'm building a inductive vehicle loop detector. By measuring the frequency of a LC-oscillator where the cars will pass over the inductor of the LC-Oscillator.

If the cars pass over the inductor the frequency will change over time according to one of the following profiles.

The frequency will also slightly change / drift based on environmental variables: like temperature, humidity or other variances in the environment of the loop.
To filter these last variations out, I used running average and std_variation based on: This post/comment on stack exchange. And was planning to detect vehicles when the frequency would go too many standard deviations away from the average.

The measurements made when there is a vehicle present shouldn't be added to the average calculation. This is easy to do when a certain threshold is passed.

But the problem is when a car is approaching the frequency will start to change but hasn't yet passed the threshold so these samples are added to the statistics. When a car approaches very slowly the standard deviation will be changed so much that the car will never be detected.

What solutions / alternatives / suggestions do you have to detect approaching cars and discard the values from the statistics for the reference empty loop.

Thanks in advance

Update: because I didn't get any response in the DSP section, I asked the question in the electronics section.

Let's ignore the issue that you already have a discrimination problem in your hands (bus/car/truck/van) there is an initial discrimination problem that has to be solved for the data to be considered valid for any moving object.

The issue is one of rate of change, so the problem should probably be moved into derivative space. That is, ignore the mean frequency and only use the derivative of the curve. Extracting derivatives of real data is a tricky issue as a true derivative amplifies the noise, so you end up with a bandpass filter with a single zero at DC that carries out a first derivative only for the bandwidth of interest and attenuates the noise. (That's a field all in itself, but this should get you started).

You are left with a curve that represents the rate of change.

Ignoring the sign and the additional gains in information that you could extract by using the statistics of the signal, the problem is basically reduced to that of a binary classifier, which in many fields is described by the receiver operating characteristics curve (ROC). You can see these in the following figure. Depending on the location of the threshold you get a different mix of false positives and negatives.

Clearly in this image from wikipedia the signal and noise are represented by Normal distributions, but that is seldom the case. In your case it might be a safe assumption for the noise, but not for the signal.

Given that you have additional statistics that you can use for your vehicle classifier, you can allow for a higher proportion of false positives to get through.

A very similar problem to this is in the biological field with the detection of spikes in neural action potentials. The way this is commonly solved is by using the standard deviation (i.e., power) of a section of the input that is "known" to be free of signal and is therefore noise and using that standard deviation to set the threshold (in these applications the voltage threshold is set between 4 to 6 standard deviations of the noise). Of course, this has not solved the problem, as you still have to figure out how to get the signal out of your noise calculations.

A clever way to do this, is to rely on the statistics of the noise being Gaussian, and only use the tail of the total input distribution where you are almost assured that all you will see is noise. And infer the noise standard deviation only from that tail. A family of algorithms can be derived from this, in neural spike detection these are known as the ADAband family of algorithms.

The idea is:

calculate the power (squared magnitude) of the input
calculate a short-term estimate of the standard deviation using a small number of samples.
store a number of estimates, let's say 100 (a longer time segment) in an array.
sort the stored values (this is in a sense similar to a median filter).
use only the bottom bins (e.g., the first 15 represent the bottom 15% of the corrected standard deviation estimates.
weigh each bin appropriately according to the chi-squared distribution it comes from and combine them to produce a better standard deviation estimate.
appropriately filter/combine successive estimates to reduce its variability.

Do note that this is just to obtain a reasonable estimate for the threshold. Any classification algorithm (including CUMSUM) would require a threshold to be defined. The advantage here is that the threshold adapts to the conditions of the input.


Formant analysis LPC (Matlab) I wanted to do formant analysis using MATLAB and did everything according to the following tutorial:
https://de.mathworks.com/help/signal/ug/formant-estimation-with-lpc-coefficients.html
I have speech signals with a sampling frequency of 44.1kHz, therefore I changed the order of the LPC estimation to 46.
But I have the following problem: the algorithm sometimes works well, but most of the time the formants are completely off. Especially for female voice. I tried to play around with the required bandwidth and lower limit for the formant frequency, but nothing yielded sufficiently good results for various voices.
I hope, that someone might have some tips for me to improve the code, such that it works more generally.

According to the documentation and my past experience with the ETSI GSM 06.10 "GSM Full Rate Speech Transcoding" standard, you should be careful in the requirements of each block of the system and the length of the segments after the pre-emphasis or segmentation blocks.
Probably you made it, but check the length of the signal before running the LPC function. And remember that we can obtain p/2 formants of the voice. GSM uses p=8 in the standard.

Cheers.


How to convert time-domain signal to complex envelope? Matlab and Simulink Communications Toolbox digital demodulators are defined to only work on the complex envelope representation of a baseband signal.

To obtain the time-domain representation of this signal, I believe one takes the real part of the complex envelope.

Given just this real part, how does one convert the time-domain representation of modulated data back into a complex envelope such that a Matlab demodulator will demodulate it?

Application: I'm trying to simulate a simple audio FSK system in Simulink.

These are the main ideas:

Consider a receiver that picks up a signal $r(t)$. This signal has bandwidth $W$ and is centered on carrier frequency $f_c$.
Using the Hilbert transform, eliminate the negative frequencies of $r(t)$. The resulting signal, $r_+(t)$, is called an analytic signal.
Now, downconvert the analytic signal using a complex exponential, so that the spectrum is centered around 0. This signal, called $\tilde r(t)$, is the complex envelope of $r(t)$. Its spectrum goes from $-W/2$ to $W/2$, and its bandwidth is $B=W/2$.
The complex envelope is almost always complex, but if the spectrum of $r(t)$ is symmetrical around $f_c$, the CE is real. Examples of this are AM DSB and BPSK.

In the transmitter you would do the opposite:

Start with a complex envelope $\tilde s(t)$ that represents the information you want to transmit. Its spectrum should go from $-W/2$ to $W/2$. An example would be a QAM signal with complex symbols.
Upconvert the CE to frequency $f_c$, multiplying by a complex exponential. The result is an analytic signal with only positive frequencies. The signal will be centered around $f_c$ and cover the frequency range from $f_c-W/2$ to $f_c+W/2$.
Take the real part of the analytic signal to convert it to a real signal that can be physically generated and transmitted.

All of these steps can be accomplished in Matlab. Be sure to read the documentation for the hilbert command. I can't comment on the corresponding Simulink blocks, since I don't use them.

Also note that it is very likely that you don't have to do all this. If you can perform all your simulations using the complex envelope, you'll save a lot of computing time. The reason is that the transmitted and received signals need to be sampled at least at a rate $2(f_c+W/2)$, whereas the complex envelope can be sampled at $W$ real samples per second -- a huge difference for large $f_c$.


Designing A Filter Using SPTool In Matlab I am currently working on learning DSP and for my course we use Matlab. I have been tasked with designing a low-pass filter using SPTool with the following specification:

  A Low-Pass Kaiser Window FIR filter with a cut-off at 200hz.

Now, when i open the SPTool and choose "new" i get the filter design and analysis window. I use the following settings

Response Type = "Lowpass"
Design Method = "FIR" and "Window"
Window = "Kaiser"

I would like to know which values i should put into the "Fs" and "FPass" fields. The signal i am applying this filter to is a periodic 100hz wave.

Thanks for any feedback.

John Smith,

You are designing a low pass filter, which means that your filter will be applied, (convolved) against a signal, whereby it will attenuate certain frequencies, and leave other frequencies more or less untouched.

So in a low pass filter, generally that means that you will filter out (heavily attenuate) frequencies that are 'high'. Your $f_{pass}$ just quantifies this for you. How high is high? Well, if you want a 'cutoff' at 200 Hz, then that generally means you want to keep frequencies at roughly 200 Hz and below, so pick $f_{pass}$ to be 200 Hz.

But bear in mind, you usually have a specification about transition bandwidth. If your filter is infinitely sharp, you can say that you accept all frequencies less than 200, and attenuate all frequencies above 200. Realistically however you might want to accept all frequencies less than say, 130, and then from 130 to 200 you are transitioning, and then heavily attentuate all frequencies greater than 200.

Regarding the sampling rate, that is really dependent on your simulation. If you pick 10,000 Hz as your sampling rate, then you should make sure that when you use your filter on a signal, that it too was sampled at 10,000 Hz.

It is good that you are using MATLAB. It is a very powerful tool in understanding DSP concepts.


Robust phase extraction of STFT bins I am interested in tracking how the phase of a signal at a particular frequency changes over time.

The method I am using calculates a series of STFTs and takes the arctangent of the imaginary and real components of the desired bin for each FFT, however, due to noise in the system I am having issues with phase unwrapping.

Are there any alternate and perhaps more robust methods of extracting phase from a signal in this manner?

Thank you

like hotpaw, unless you know something prior about the signal (like it's periodic with period $P$), i'm thinking that you're windowing in the neighborhood of samples $m \le n < m+N$.

STFT: $$ X_m[k] \triangleq  (-1)^k \, e^{-j\pi k/N} \ \sum\limits_{n=0}^{N-1} x[n+m] w[n] \, e^{-j2\pi nk/N} $$

$w[n]$ is a decent window centered about $n \approx \tfrac{N-1}{2}$, specifically $$w[n]=w[N-1-n] \ .$$

if $x[n]$ and $w[n]$ are both real, then $X[k]$ is Hermitian symmetric: $X[N-k]=\overline{X[k]}$.  and $X[0]$ must be real.  so the complex angle of $X[0]$ is either $0$ or $\pm \pi$ depending on if $X[0]$ is positive or negative.

$$ \arg\{ X[0] \} = \begin{cases}
0       & \text{for } \Re\{X[0]\} \ge 0 \\
\pm \pi & \text{for } \Re\{X[0]\} < 0 \\
\end{cases}$$

because $ \Im\{ X[0] \} = 0 $.

then you can calculate phase increments:

$$\begin{align}
 \arg\{ X[k] \} -  \arg\{ X[k-1] \}  &= \arg\left\{ \frac{X[k]}{X[k-1]} \right\} \\
\\
                                     &= \arg\left\{ \frac{X[k]\overline{X[k-1]}}{\Big|X[k-1]\Big|^2} \right\} \\
\\
                                     &= \arg\left\{ X[k]\overline{X[k-1]} \right\} \\
\\
                                     &= \arg\left\{ \big(\Re\{X[k]\}+j\Im\{X[k]\}\big)\big(\Re\{X[k-1]\}-j\Im\{X[k-1]\}\big) \right\} \\
\\
                                     &= \arg\left\{ \Re\{X[k]\} \Re\{X[k-1]\} + \Im\{X[k]\} \Im\{X[k-1]\} + j\big( \Im\{X[k]\} \Re\{X[k-1]\} - \Re\{X[k]\} \Im\{X[k-1]\} \big) \right\} \\
\\
                                     &= \arctan\left(\frac{ \Im\{X[k]\} \Re\{X[k-1]\} - \Re\{X[k]\} \Im\{X[k-1]\}} {\Re\{X[k]\} \Re\{X[k-1]\} + \Im\{X[k]\} \Im\{X[k-1]\}} \right) \\
\end{align}$$

or recursively

$$ \arg\{ X[k] \} = \arg\{ X[k-1] \} + \arctan\left(\tfrac{ \Im\{X[k]\} \Re\{X[k-1]\} - \Re\{X[k]\} \Im\{X[k-1]\}} {\Re\{X[k]\} \Re\{X[k-1]\} + \Im\{X[k]\} \Im\{X[k-1]\}} \right) $$

for $ 1 \le k \le \tfrac{N}{2} $.

this is assuming that every little phase increment is smaller in magnitude than $\tfrac{\pi}{2}$.  this is how you deal with unwrapping phase naturally.


Approximating a Gaussian Process Suppose that $\theta_t$ is an exogenous variable with known Gaussian process. Next, suppose that for any index $i\in [0,1]$,
$$ a_{i,t} = (1-\beta)\mathbb E[\theta_t|\mathcal I_{i,t}]+\beta \mathbb E[a_t|\mathcal I_{i,t}],$$
where $a_t=\int_0^1 a_{i,t} di$, and
$$\mathcal I_{i,t} = \{\dots,\theta_{t-1}+\epsilon_{i,t-1},a_{t-2},\theta_{t}+\epsilon_{i,t},a_{t-1} \}.$$
Finally, for all $i$ and $t$, $\epsilon_{i,t} \sim \mathcal N(0,\sigma^2)$.

I am trying to find an algorithm to numerically approximate the stochastic process for $a_t$. I have tried many different ways of applying the Kalman filter to do this but none of them converge.

  Is there an algorithm that is likely (or better yet guaranteed) to converge?

You can use an ARMA model for $\alpha_t$. But, there is an equivalence between ARMA models and Kalman filters.

Therefore, what I am thinking is that maybe the formulation of the problem is a little bit off and instead of having to produce one estimate, you really have to estimate two. Here is why:

You present an $\alpha_{i,t} = (1-\beta) \mathbb{E}[\theta_t|\mathcal{I}_{i,t}] + \beta \mathbb{E}[\alpha_t|\mathcal{I}_{i,t}]$. What does this say? It says that $\alpha_{i,t}$ is a "mixture" of the average value of two processes:

The expectation of $\theta$ given the current $\mathcal{I}$
The expectation of $\alpha$ given the current $\mathcal{I}$

With $\beta$ controlling the mix.

It is then added (as a constraint (?)) that the value of $\alpha_t$ is known to be the sum of all $\alpha_i$ (here, two of them) at the current time instance. So, is $\mathbb{E}[a_t|\mathcal{I}_{i,t}]$ different than the constraint on $\alpha_t$ or the same thing? Because, both seem to determine $\alpha_t$ (?).

Then comes the definition of $\mathcal{I}_{i,t}$ and I think that this one should really be "split" into two different things. In its current form it is specified as $\mathcal I_{i,t} = \{\dots,\theta_{t-1}+\epsilon_{i,t-1},a_{t-2},\theta_{t}+\epsilon_{i,t},a_{t-1} \}$. Notice here that $\mathcal{I}$ contains both $\theta$ and $\alpha$ elements.

This is the bit that I don't get. Do $\theta, \alpha$ have the same units / dynamic range? Are they "similar"? If they are, why don't we call them all the same name / type? Is $\mathcal{I}$ meant to say that $\theta,\alpha$ are considered together for the same time instance? That's "fine" to a point but it would not work from the point of view of estimation. Imagine, for example, if you had to get a derivative of $\mathcal{I}$. It would work across the $t$s, but not across the $i$s because you would be differentiating different things, in different "references" ($\frac{\alpha - \theta}{(what?)}$)

Therefore, given what is presented in the question, my recommendation(s) would be that:

$\mathcal{I}$ is split into two parts, one for the $\theta$s and one of the $\alpha$s.
These are then plugged back into the estimation of $\alpha_{i,t}$ but essentially what we get here is two things that need to be estimated prior to obtaining the final value for $\alpha_t$.
The constraint (?) on $\alpha_t$ is clarified. If $\alpha_t$ is known to be this integral, then there is no point in also deriving the expectation of $\alpha_t$ from $\mathcal{I}$. It is either one or the other. Except if it is a constraint that means to say "whatever you work out your $\alpha$s to be, their sum will have to be unity"...but at the moment it does not look like this (?).

Finally, if you do have some vectors of the quantities you are trying to estimate, it would be good to provide one or two plots, just to see what you are dealing with and what this "non-convergence" looks like.

Hope this helps.


What is QEEG technically, and how is it computed? I've been working with EEG signals for a while. I wonder how to compute QEEG (or Brain Map) representation. I think it is extracted from raw EEG and I want to know what is technical method for QEEg extraction? Is it simply a Fourier transform or something complicated is utilized?

In a nut shell:
 - What is qEEG exactly and how to compute it technically to come up with something like in figure above?

I feel it should be related to short time-Fourier transform.

qEEG or quantitative EEG is an umbrella term that encompasses many different signal processing methods that compute a number/statistic/quantitative biomarker from EEG signal traces. Algorithms usually boil down to some form of spectral analysis either using Fourier transforms, or wavelets.

The images in shown in the question are most likely generated through spectral analysis of individual EEG channels. In your figure, it seems like the author chose Delta to mean 1-3.5Hz. Some authors use Delta as all frequencies <4Hz. (https://en.wikipedia.org/wiki/Electroencephalography)

To compute the contribution of each band, one resorts to standard power spectrum analysis techniques. For example, for a given EEG channel $i$, the Delta qEEG biomarker value can be computed as:

$$
\mbox{Delta}_i = \frac{\int_{1}^{3.5} P_{i}(f) df}{\int_{0}^{\infty}P_i(f)}
$$

where $P_i(f)$ is the power spectral density of the $i^{th}$ EEG channel signal. In practice, you can use, say, Welch method to estimate $P_i(f)$ from sampled EEG channel traces, the integrals will get replaced by sum over frequency bins, and you'd only go up to Nyquist.

Repeat this process for each frequency band, for each EEG channel to get sequences $\{\mbox{Delta}_i\}_{i=1}^N$,$\{\mbox{Theta}_i\}_{i=1}^N$, $\{\mbox{Alpha}_i\}_{i=1}^N, \cdots$, where $N$ is the number of EEG channels.

The final step is to map these numbers on the image of the scalp. This is a visualization step which involves interpolation. You can use a tool like this: https://martinos.org/mne/stable/generated/mne.viz.plot_topomap.html#mne.viz.plot_topomap or https://www.nbtwiki.net/doku.php?id=tutorial:plot_a_time-frequency_representation_and_power_spectrum_of_a_single_channel#.WJtLMxC2pJI.


Is a good practice to use one PID in rise time and another PID for steady time in order to eliminate oscillations during steady time? I'm trying to control temperature with PID(Kp,Ki, Kd). During the rise time the system use one PID and when target is reached the system switch to another PID with other parameter in order to eliminate the oscillation. I know that I can remove oscillations increasing Kd but in this case it does not work.

Results of real time system

Control temperature with constant PID:

Control temperature using PID1 during rise time (blue color) and when temperature reach the target, the controller is changed to PID2 (Steady state error: black color):

void setTarget(...)
{
    mPid = pid1;
}

void steadyStateError()
{
    float e =  fabs( temperature - target);
    if(e<=0)
       mPID = pid2;
}

Switching between 2 PID controller is called gain-scheduling. There are various ways to implement gain scheduling. But basically, the idea is that the gain should change smoothly. You should not change suddenly. You could have one set of gains if T < T1, another set if T > T2 and use interpolation between the 2 sets if T1 < T  < T2.

You can also use 2 PID controller with the unused PID set in tracking mode while the other controls the actuator also. However, it will be harder to have a smooth transition between the 2 controllers.

Secondly, before implementing gain scheduling, make sure your issue is really due to non-optimal gains. Typically, there are 2 problems with control loop for temperatures

1 - Actuator saturation which causes integrator windup. This happens often in temperature control when the difference between the setpoint and the measured temperature is significant. This can be fixed with an appropriate anti-windup scheme.

2 - Asymetry between heating and cooling. Typically, if you have a heater, you can only "heat" you cannot cool your system. This makes the choice of gain a bit tricky as the cooldown rate is not controllable by your controller. If you have some kind of Peltier cooler, you can both cool down and heat up, but the cooling rate is usually slower.


Normalizing Cepstrum by Energy I have separately computed the cepstrum and the energy of a finite discrete-time signal. If the energy of the signal is high, we would expect (on average) larger cepstral values than if the energy of the signal is low. I am therefore interested in normalizing the cepstral values by the energy in some fashion to adjust for this, because I want to be able to compare the cepstral values from different signals in an energy-invariant manner. However while I know the units of the signal's energy $(\text{volts}^2\cdot\text{seconds})$ I have no idea what the units of the cepstral values are, and therefore do not know how to transform them in order to ensure that the normalization is performed appropriately. What would be the most unit-consistent way to perform this normalization?

Nevermind, I figured it out. For my purposes, all I really want is the cepstrum of the energy-normalized signal, which after working through the algebra is equivalent to adding a delta function to the original cepstrum. (See below.)

For a real, discrete, finite-time signal $\vec{s}$, we can compute its energy as $E(\vec{s})=\vec{s}\cdot\vec{s}$. To energy-normalize the signal (i.e. ensure that $\text{energy}=1$), we multiply $\vec{s}$ by the scalar $E(\vec{s})^{-1/2}$.

The cepstrum of the the signal $\vec{s}$ is defined as $\text{cep}(\vec{s}) = \mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big)\big)$.
Therefore the cepstrum of the energy-normalized signal can be computed as: $\text{cep}\big(E(\vec{s})^{-1/2}\cdot\vec{s}\big)\\ = \mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}\big(E(\vec{s})^{-1/2}\cdot\vec{s}\big)\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(\big|E(\vec{s})^{-1/2}\cdot\mathfrak{F}(\vec{s})\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(E(\vec{s})^{-1}\cdot\big|\mathfrak{F}(\vec{s})\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(E(\vec{s})^{-1}\big)\cdot\mathbf{1} + \text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big)\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big) - \text{log}\big(E(\vec{s})\big)\cdot\mathbf{1}\big)\\ =
\mathfrak{F}^{-1}\big(\text{log}\big(\big|\mathfrak{F}(\vec{s})\big|^2\big)\big) - \mathfrak{F}^{-1}\big(\text{log}\big(E(\vec{s})\big)\cdot\mathbf{1}\big)\\ =
\text{cep}(\vec{s}) - \big[\text{log}\big(E(\vec{s})\big)\;\; 0\;\; 0\;\; \dots\;\; 0\big]$

Where $\mathbf{1}$ denotes the vector of all 1's.


shannon interpolation in image processing I tried to implement Shannon interpolation on a 2D array.

First, implemented it on a 1D signal, just for sanity-check:

t = -7:7;
t2= -7:0.01:7;
figure
plot(t, sin(pi*t / 3))
upSamp = zeros(size(t2));
for iT = 1:length(t2)
    upSamp(iT) = sum(sinc(t2(iT) - t) .* sin(pi*t / 3));
end
hold
plot(t2, upSamp)

Looks ok.

Then, implementation on a 2D array. The reference array is:

figure
surf(sin(pi*T1 / 3).*sin(pi*T2 / 3))

Then the interpolation is done by a convolution with a Sombrero function:

[T12, T22] = meshgrid(t2, t2);
val = zeros(size(T12));
for iInterp=1:numel(T12)
    xdist = (T12(iInterp) - T1);
    ydist = (T22(iInterp) - T2);
    val(iInterp) = sum(sinc(sqrt(xdist.^2 + ydist.^2)) .* sin(pi*T1 / 3).*sin(pi*T2 / 3), 'all');
end
figure;surf(val, 'EdgeColor', 'none');

seems to work good only on the edges :(

Any suggestion?

The error source is the kernel function, which should be a multiplication of two 1D sincs (without rotational symmetry) instead of the Sombrero function, which is a sinc function with rotational symmetry.


Fourier transform of the magnitude of the fourier transformed signal I've come to realize that the Fourier transform of an already Fourier transformed signal gives the time-reversal signal.
$$\mathcal F(\mathcal F(x(t)))=x(−t)$$
ref 1,ref 2
However, my question is, if Fourier transformation is applied to the magnitude of the Fourier transformed signal (excluding phase information), is there a way to interpret the returned signal?
$$\mathcal F\{|\mathcal F\{x\}|\}=?$$

If it were the squared magnitude, the result would be the auto-correlation function, apart from a time inversion and possibly a scaling factor. If it's just the magnitude, then the result can be characterized as the only even function that has the same Fourier transform magnitude as the original function.


Why do I need inverse FFT? I have magnetometer time-series measurements in mV. I am trying to understand what is theory behind software that the company has provided.

Let $x$ be your time series vector, which is an array of double, and $W$ your
window length. $x$ is sliced into $L$ windows, each window is detrended and
(if not rectangular window function is chosen) a Hanning window is
applied. Then the FFTW is used to obtain an array of complex numbers $X$.
From $W$ real input values we get $P = W/2$ complex numbers in $X$. $X$ is calibrated with

$$
    X'= \frac{X}{\sqrt{f_s*P}}
$$
Transfer function is given by
$$
   H(\omega) =0.8V/nT*F_m(\omega)*F_t(\omega)*V
$$
where
$$
F_m(\omega)=\frac{j\omega}{j\omega+2\pi 4}\\
F_t(\omega)=\frac{2\pi 1591}{2 \pi 1591 + j\omega}
$$

If a transfer function $H$ is given we obtain

$$
    X'' = X'/H
$$

Finally if you use a windowing function the result has to be multiplied with the inverse of the integral of that window.

My question is why do I need the inverse of the integral of that window?

Are my units going to be $nT$?

I am giving part of the calibration file

 FREQUENCY    MAGNITUDE      PHASE
    Hz                     deg
Chopper Off
+4.0000E-01  +1.0714E-02  +1.4808E+02
+5.5580E-01  +1.3325E-02  +1.3911E+02
+7.7230E-01  +1.5775E-02  +1.2956E+02
+1.0731E+00  +1.7732E-02  +1.2025E+02
+1.4910E+00  +1.9094E-02  +1.1174E+02
+2.0718E+00  +1.9931E-02  +1.0442E+02
+2.8787E+00  +2.0392E-02  +9.8155E+01
+4.0000E+00  +2.0603E-02  +9.2520E+01
+5.5580E+00  +2.0611E-02  +8.7089E+01
+7.7230E+00  +2.0408E-02  +8.1384E+01
+1.0731E+01  +1.9949E-02  +7.4971E+01
+1.4910E+01  +1.9006E-02  +6.7315E+01
+2.0718E+01  +1.7761E-02  +5.8714E+01
+2.8788E+01  +1.5629E-02  +4.9042E+01
+4.0000E+01  +1.3116E-02  +3.9247E+01
+5.5580E+01  +1.0350E-02  +3.0040E+01
+7.7230E+01  +8.0188E-03  +2.2525E+01
+1.0731E+02  +5.9868E-03  +1.6408E+01

It seems that people do not understand my questions.Manuel can dowloaded from here

https://geo-metronix.de/ocl/public.php?service=files&t=9b54b82730e4b8b869e083574923d99a

You can go to tsmp and download tsmp.pdf.Calibration is on page 40.

I'm still not sure that I understand what your question is, but this specific one is easy;

  My question is why do I need the inverse of the integral of that window?

The issue is that you'll apply a window to the data $x$:
$$
x_2(t) = x(t) \cdot w(t)
$$
that means that the sum of $x_2$ is going to be different from the sum of the (rectangular windowed) $x$.  All the inverse of the integral is doing is rescaling things to attempt to make the sums equal.

As far as your other question goes:

  Are my units going to be nT?

That's all down to calibration. You have millvolts.  Your expression for $H(\omega)$ has a calibration factor $0.8V/nT$. To me, that means that 0.8V (or 800 mV) is 1 nT.

Tesla are one kilogram per second squared per ampere in base SI units.  Usually, with any measurements calibration is what equates the voltage (or current, usually) measured to the units of you're actually interested in.

All $H(\omega)$ seems to be doing is applying a bandpass filter to your data.


Difference between jpg compression and mp3 compression My understanding of jpg compression is that, if done repeatedly, any image will eventually be reduced to a single flat color, due to the way it deals with neighboring pixels. However, years ago I was told that mp3 does not have this problem, a 192kbps song will stay at the same audio quality with little file size variation no matter how many times it is compressed at the same setting.

Why is that? What's the differentiating factor?

Thanks

Actually, it's kind of the other way around. If you reuse the same JPEG encoder at the same quality level (without any smoothing steps as built-in prepcosessing) and a decoder which faithfully decompresses the images, I expect the image quality not to degrade from generation to generation. This is because quantization (the lossy part) is done the same way during the 2nd compression. You can imagine scalar quantization like this: given some value v and a quantizer step s, the quantized value is round(v/s)*s. Example:

v = 6.3
s = 2
q = round(v/s)*s = 6

So, the first "compression" includes an error of 0.3 in this case. But if you feed this into the encoder a 2nd time

v = 6
s = 2
q = round(v/s)*s = 6

it does not change anymore. It stays 6. If you change the quality level from one compression to the next, you also change the quantizer step s. This would add a new rounding error.

Keep in mind that this is just a simple example. JPEG Quantization may include a "dead zone quantizer" which changes the threshold a bit near 0. For example, you might want to round any value v to zero in case |v|<0.7. Something like this might give you a slightly higher quality-per-bit ratio because zeros can be much more compactly represented. But this should not change what I described above.

With MP3 the "quantization steps" are computed dynamically and depend on the actual signal. This is what the "psychoacoustic model" is for. It analyzes the audio data and tells us how strongly we can quantize certain portions of the signal without us noticing a difference. But this quantization adds "quantization noise" which could make the psychoacoustic model compute other quanzitation steps during the 2nd compression even in the case where you told the encoder to use the same "quality level". In the case the 2nd compression uses different quantization steps, you'll introduce new rounding errors.

Actually, something like this has been done before. Some artist did such an experiment where he compressed some audio file over and over again (like 1000 times) and the result had little to do with the original. Unfortunately, I don't remember any names or references.

There are other issues which make lossy audio compression more complicated. For example, a decoder might add or remove some leading samples (encoder+decoder delay mismatch) which would make the encoder work on different blocks for the 2nd pass. This would also add new rounding errors. But if I remember correctly, such a "block position mismatch" was shown to be beneficial w.r.t. generation loss by some experiments done by the Hydrogenaudio crowd a couple of years ago. So, while there was still a generation loss, the loss was less annoying than with matching block positions.


Why does Pre-warping work with a notch filter but not a low pass filter? Hey guys a little bit of a continuation from a previous post I made here, I found the solution as why the pre-warping was not working at the time,
However I was just playing around and found out pre warpings were working with notch filters pretty good oppose to the low pass filters, now I am more then confused. I wonder if anyone can explain this to me?
Ill give two case sceneries where the Critical frequency is close to and not close to the Nyquist frequencies.
Diagrams
Test 1:

Critical Frequency = 40kHz
Sampling Frequency = 96kHz
Ratio = 40kHz/48kHz = 83.33% pretty close.

Test 2:

Critical Frequency = 20kHz
Sampling Frequency = 96kHz
Ratio = 20kHz/48kHz = 42% halfish close.

Test 3:

Critical Frequency = 10kHz
Sampling Frequency = 96kHz
Ratio = 10kHz/48kHz = 21% closeish.

What the heck is going on, Why is it working for notch and not a low pass?!

The bilinear transform compresses the analog frequency in the s-domain associated with  $f=\infty$ to the point $z=-1$ in the z-domain (which is where the frequency is $f_s/2$ where $f_s$ is the sampling rate. Consider what the magnitude of your notch filter is at $f=\infty$: It has no effect on the magnitude or phase of your signal; there is no filtering impact there whatsoever.
Compare in contrast to this your low pass filter (and how much of your "desired" shape is already beyond the Nyquist boundary which no pre-warping would be able to reconstruct). As you approach the Nyquist boundary the ability to compensate with pre-warping becomes increasingly difficult (you will see the same with the notch if it is brought close enough such that there is continued attenuation beyond the boundary).
The degree of this frequency warping and how it pronounced is depicted in the chart example that I pasted below (not the OP's sampling rates) which shows which frequency on the vertical axis for a given sampling rate of 20 rad/sec (125.6 Hz) would map to the digital frequency domain extending from $\pi$ to $\pi$ rad/sample.
With prewarping, especially at high frequencies, it is very difficult to match an analog shape over a wider frequency range given this mapping is according to the following transcendental equation (would require Taylor series expansion out to an increasingly higher order as we near the Nyquist boundary to match effectively, but is easily matched with the first or second terms only at the lower frequencies where the slope becomes approximately linear):
$$\omega_a = \frac{2}{T}tan(\omega_d T/2)$$
(as referenced here:  https://en.wikipedia.org/wiki/Bilinear_transform)
But for a single point, or narrow range as typical for a notch filter, this repositioning of a single notch zero would be trivial to do and with high accuracy for the repositioned notch location (but distortion right at the Nyquist boundary if the target filter has residual insertion loss beyond Nyquist).


Determining how similar audio is to human speech While looking for an answer to this problem, I found this board so decided to cross post this question of mine from Stack Overflow.

I am searching for a method of determining the similarity between an audio segment and a human voice, which is expressed numerically.

I've searched quite a bit, but what I've found so far (detailed below) doesn't really fit what I need:

One method is to use speech recognition software to obtain words from an audio segment. However, this method is unable to come up with how "similar" audio is to human speech; it can often tell whether or not there are words in the audio, but if there are no definite words, it can't tell close the audio is to having such words.Examples: CMU Sphinx, Dragonfly, SHoUT
The more promising method is referred to as Voice Activity Detection (VAD). However, this tends to have the same problems: the algorithms / programs using VAD tend to just return whether or not the activity threshold has been reached, and no "similarity" value before or after such threshold. Alternatively, many just look for volume, not similarity to human speech.Examples: Speex, Listener, FreeSWITCH

Any ideas?

This kind of problem is usually solved using machine learning techniques.

Break down the signal into a sequence of 20ms or 50ms frames. Extract features on each frame. MFCC are generally good for this kind of application, though there are feature more specific to voice detection (4 Hz modulation energy - which is roughly the rate at which people speak ; zero-crossing rate).

Then, using a training set of audio you have manually labelled as being speech / not speech, train a classifier (Gaussian mixture models, SVM...) on the frames features.

This will allow you to classify unlabelled frames into speech/non-speech classes. The last step consists in smoothing the decisions (a frame classified as non-speech surrounded by hundreds of speech frame is likely to be classification error), for example using HMMs, or just a median filter.

A few references:

Robust speech / music classification in audio documents (Pinquier & al)
Speech / music discrimination for multimedia applications (El-Maleh & al)
A comparison of features for speech / music discrimination (Carey & al)

Note that the features and classification techniques they describe are also relevant for the 1-class problem of detecting speech (instead of discriminating speech vs something else). In this case, you can use 1-class modeling techniques such as 1-class SVM, or just take the likelihood score out of a GMM trained on a speech data as a "speechiness" measure.

If, on the other hand, your problem is really discriminating speech vs something else (say music), you could also very well use unsupervised approaches which are focused on detecting the boundaries between similar audio content - rather than identifying this content itself.


DT unit impulse function properties proof How to prove the following properties of DT unit impulse function. If anyone got link to a proof please mention it. I have search on the web. But only found the properties, not a proper method of proof.

Thank you.

Hint :
Use definition of $\sigma[n]$

consider
 2 cases
$$m=n$$
$$m\ne n$$

also, use definition of convolution

as an aside, it took me a while to do this but a good habit is to use $[n]$ (square brackets) for discrete time index, and $(t)$ (curved brackets) for continuous time


What is the spectrogram of an impulse response used for? I saw that in acoustics an estimated impulse response of a room can be transformed using STFT/wavelets to a spectrogram as an analysis tool. However, I did not understand how the spectrogram is interpreted (in contrast to the spectrum of the impulse response which gives us the frequency response for example).

Also, are there other cases in which the spectrogram of impulse responses is used, and what are the intrepretations in those cases?

Thanks!

I think you are thinking of waterfall plots and/or cumulative spectral decay plots.  They are spectograms in the widest sense of the word since their windowing function tends to be exponential.  They are popular for plotting the response of speaker systems as they show group delay for various frequencies a lot more, well, graphically than the impulse response itself or a magnitude/phase plot of the transfer function would.  Resonators can be used for straightening out an otherwise wobbly (for whatever unrelated reason) frequency response magnitude at the cost of introducing delays.

The resulting lack of "impulse accuracy" which leads to problems with sound clarity and spatial perception is not easily apparent with other representations.


What are the definitions of symbols and samples related to PCM? What are the definitions of symbols and samples related to PCM? For example in binary encoding $\{0,1\}$ with 8 levels of quantization, is one symbol either of $ \ 0 \  or \  1\ $? or is it the collection of  $(x,y,z)$ where $x,y,z \ \epsilon \ \{0,1\}$.

Similary for ternary encoding $\{0,1,-1\}$,suppose with 9 levels of quantization, is a symbol either of $-1,1,or \ 0 \ $ or is it complete $(x,y)$ where $x,y \ \epsilon
 \ \{-1,0,1\} \ $ ?

Quantization and encoding are largely independent. "Symbols" is another word for "pulses", and the line encoding can also play a role in how the information is transmitted.

Say you quantize one sample to $256 = 2^8$ levels, or 8 bits/sample. In order to transmit those 8 bits, you can, among other options:

use binary encoding, which requires transmitting eight symbols, each with amplitude in  the set $\lbrace +, - \rbrace$, where $+$ and $-$ correspond to some predefined voltage levels.
Use ternary encoding with AMI (alternate mark inversion), which also requires eight symbols, but now with amplitudes in the set $\lbrace +, -, 0 \rbrace$.
Use quaternary (or 2B1Q) encoding, which requires four symbols.

Other examples may include Manchester encoding, CMI, block line codes.... but in all cases, one symbol correponds to one transmitted pulse.


ICA - Exercise Problem from book - Matlab Code Im trying to do a problem of ICA from a Image Processing book, but when I finish I got the Matrix C different from the one in the exercise.
This is the Exercise:

This is my Matlab Code:

t=1:1000;
f=sin(t*pi/19);
g=sin(pi*t/5)+sin(pi*t/31);
h=mod(.3*t,1.5);
m1=.3*f+.4*g+.3*h;
m2=.5*f+.2*g+.3*h;
m3=.1*f+.1*g+.8*h;
mn=[m1;m2;m3];
m1p=m1-mean2(mn);
m2p=m2-mean2(mn);
m3p=m3-mean2(mn);
mnp=[m1p;m2p;m3p;];
C=(mnp*mnp')/3;
[V,D]=eig(C);
eigenvalues=diag(D);
eigenvector=V;

My result in the Matrix C is:

C =

   77.2656   60.7483   26.2664
   60.7483   64.0480   23.0485
   26.2664   23.0485   56.4664

This are the steps of the algorithm:

Then there is some info about applying to Signal Processing:

My problem is I don't get the same Matrix C from the exercise.

In the line

C=(mnp*mnp')/3

you're dividing by 3 when you should really be dividing by the number of observations, which is 1000 in your case.


Pixel time series trend: What are the causes? I have taken an video of the following static scene for 10 mins (3 frames/Sec); totally 1800 frames for 10 mins.

I also show here the pixel time series, that is, a set of values of a particular pixel say (1,1) for 1800 frames..

I could see a linear trend in the series. For some pixels the trend is increasing and for some the trend is decreasing. Due to trend, though the scene is static, there appears to be a small movement in the object. It is possible to de-trend it.

My Question is: what are the causes of this trend? I use a standard camera fixed on a immovable plane and I illuminate with light in order to photograph it.

You are very likely to observe trends in either your camera or your light source.

Sources from the camera can include:

thermal effects inside the sensor (basically, the more you use it, the more current will traverse it in various parts and it will get warmer)
drift in the ADC of your sensor (probably linked to the aforementioned).

Also, light sources can also drift with time: thermal effects will influence light bulbs.

You can try to mitigate these effects by:

changing your light source (I would assume LED lamps are more stable, but this needs to be checked)
waiting long enough before starting the measurements in order to avoid transitory phenomenon
carefully playing with the settings of your camera. Typically, I would fix manually all the exposure-related settings (exposure time, gain, framerate...) so that the autocontrols do not make vary the exposure values when an object enters the scene. However, this will prevent the camera from tracking natural scene trends.
simply fitting a linear regression to your signal in order to remove the trend while processing.


Why is the FFT "mirrored"? If you do an FFT plot of a simple signal, like:

t = 0:0.01:1 ;
N = max(size(t));
x = 1 + sin( 2*pi*t ) ;
y = abs( fft( x ) ) ;
stem( N*t, y )

1Hz sinusoid + DC

FFT of above

I understand that the number in the first bin is "how much DC" there is in the signal.

y(1)  %DC
  > 101.0000

The number in the second bin should be "how much 1-cycle over the whole signal" there is:

y(2)  %1 cycle in the N samples
  > 50.6665

But it's not 101!  It's about 50.5.

There's another entry at the end of the fft signal, equal in magnitude:

y(101)
  > 50.2971

So 50.5 again.

My question is, why is the FFT mirrored like this?  Why isn't it just a 101 in y(2) (which would of course mean, all 101 bins of your signal have a 1 Hz sinusoid in it?)

Would it be accurate to do:

mid = round( N/2 ) ;

% Prepend y(1), then add y(2:middle) with the mirror FLIPPED vector
% from y(middle+1:end)
z = [ y(1), y( 2:mid ) + fliplr( y(mid+1:end) ) ];

stem( z )

Flip and add-in the second half of the FFT vector

I thought now, the mirrored part on the right hand side is added in correctly, giving me the desired "all 101 bins of the FFT contain a 1Hz sinusoid"

>> z(2)

ans =

  100.5943

Real signals are "mirrored" in the real and negative halves of the Fourier transform because of the nature of the Fourier transform.  The Fourier transform is defined as the following-

$H(f) = \int h(t)e^{-j2\pi ft}dt$

Basically it correlates the signal with a bunch of complex sinusoids, each with its own frequency.  So what do those complex sinusoids look like?  The picture below illustrates one complex sinusoid.

The "corkscrew" is the rotating complex sinusoid in time, while the two sinusoids that follow it are the extracted real and imaginary components of the complex sinusoid.  The astute reader will note that the real and imaginary components are the exact same, only they are out of phase with each other by 90 degrees ($\frac{\pi}{2}$).  Because they are 90 degrees out of phase they are orthogonal and can "catch" any component of the signal at that frequency.

The relationship between the exponential and the cosine/sine is given by Euler's formula-

$e^{jx} = cos(x) + j*sin(x)$

This allows us to modify the Fourier transform as follows-
$$
H(f) = \int h(t)e^{-j2\pi ft}dt \\
= \int h(t)(cos(2\pi ft) - j*sin(2\pi ft))dt
$$

At the negative frequencies the Fourier transform becomes the following-
$$
H(-f) = \int h(t)(cos(2\pi (-f)t) - j*sin(2\pi (-f)t))dt \\
= \int h(t)(cos(2\pi ft) + j*sin(2\pi ft))dt
$$

Comparing the negative frequency version with the positive frequency version shows that the cosine is the same while the sine is inverted.  They are still 90 degrees out of phase with each other, though, allowing them to catch any signal component at that (negative) frequency.

Because both the positive and negative frequency sinusoids are 90 degrees out of phase and have the same magnitude, they will both respond to real signals in the same way.  Or rather, the magnitude of their response will be the same, but the correlation phase will be different.

EDIT: Specifically, the negative frequency correlation is the conjugate of the positive frequency correlation (due to the inverted imaginary sine component) for real signals.  In mathematical terms, this is, as Dilip pointed out, the following-

$H(-f) = [H(f)]^*$

Another way to think about it:

Imaginary components are just that..Imaginary! They are a tool, which allows the employ of an extra plane to view things on and makes much of digital (and analog) signal processing possible, if not much easier than using differential equations!

But we can't break the logical laws of nature, we can't do anything 'real' with the imaginary content$^\dagger$ and so it must effectively cancel itself out before returning to reality. How does this look in the Fourier Transform of a time based signal(complex frequency domain)? If we add/sum the positive and negative frequency components of the signal the imaginary parts cancel, this is what we mean by saying the positive and negative elements are conjugate to each-other. Notice that when an FT is taken of a time-signal there exists these conjugate signals, with the 'real' part of each sharing the magnitude, half in the positive domain, half in the negative, so in effect adding the conjugates together removes the imaginary content and provides the real content only.

$^\dagger$ Meaning we can't create a voltage that is $5i$ volts. Obviously, we can use imaginary numbers to represent real-world signals that are two-vector-valued,  such as circularly polarized EM waves.


Using Spectrogram to check whether aliasing will occur I was wondering if the spectrogram can be used as a MATLAB tool to check whether or not aliasing will occur when the signal is downsampled from let's say, 44.1kHz to 8kHz? And if yes, how do I go about it?

If the spectrogram displayed is sensitive enough to show spectral energy from half the old sample rate to half the new sample rate (e.g. anything at or between 4 kHz to 22.05 kHz), that energy will be aliased unless low-pass filtered out before or during downsampling (to 8KHz).  If there is nothing there, then there is nothing that will be aliased by properly interpolated (Sinc kernel, et.al) downsampling without any further low-pass filtering (such as a wider windowed Sinc FIR kernel, etc.)


What does kchip/s datarate mean? I am working in w mbus and came across this term kchip/s in datarate like 37.768 kchip/s what does the unit specify?

Protocols like w-mbus or Zigbee typically use some form of frequency-spreading like DSSS. In DSSS we do not send symbols directly using some pulse shape and encoding the information in the sign but instead we use a rapidly alternating sequence of $\pm 1$ as a pulse. This leads to an increase in the occupied bandwidth ("spreading"). The receiver multiplies the incoming signal with the same sequence, coherently combining back all the information ("de-spreading"), which overall increases the robustness to fading and reduces the peak power spectral density. See the link for more details.

What it means practically is that each symbol is further sub-divided in time domain into smaller units called "chips". The rate at which these chips are sent is called chip rate which itself is a measure for the occupied bandwidth. For instance, a chip rate of 500 kchips/s would mean that 500'000 chips per second are sent, which requires something like 500 kHz (practically more, due to the required guard bands).

If you want to know your data rate, you need to know the code (i.e., how many chips per symbol are used) and the modulation and coding scheme (i.e., how many bits per symbol are sent).

For example, this TI reference design states 100kcps with 11-bit DSSS coding and 2-FSK modulation. This gives 100kchips/s / 11chips/symbol * 1bit/symbol ~ 9kbit/s (minus packet overhead).


minimum frequency in periodic components Suppose we have the following signal model

$y(t) = A_1\sin(\omega_1 t+\phi_1) + A_2\sin(\omega_2 t+\phi_2) + ... + A_p\sin(\omega_p t+\phi_p) + z(t)$

When we are sampling at sampling frequency $f_\mathrm s$ we know that the sampling frequency must be at least 2 times the maximum frequency, but what about a condition for the minimum frequency (largest period)? Is there any reasonable limit below which it is useless existence of frequency? Suppose that the sampling frequency is $100\,\text{Hz}$, then the maximum frequency can be $50\,\text{Hz}$, but what about the minimum frequency?

No, there is no lower bound to the sampling theorem. Any signal $x(t)$ sampled at rate $f_\mathrm s$ can be reconstructed from its sampled version if it contains only spectral components lower than (1) $f_\mathrm{s}/2$. I.e. the frequencies contained in $x(t)$ are allowed to be as low as 0 Hz. Note however, that perfect reconstruction of $x(t)$ is only possible if an infinite number of samples $N$ has been drawn from $x(t)$.

When it comes to spectral analysis of $x(t)$ using its sampled version $x_n = x(n/f_\mathrm s)$ the number of acquired samples $N$ plays an important role because it determines the frequency resolution of the DFT. In order to be still distinguishable the lowest frequency cannot be arbitrarily low. This is intuitively understandable as the lower the frequency the longer the period and thus the signal should be measured for a sufficently long time. Here is some further reading about DFT frequency resolution and spectral leakage

Number of FFT points required for a specific frequency resolution for an oversampled signal

FFT of sine wave not coming as expected i.e single point

(1) Actually $x(t)$ must not contain a sine wave with exact frequency $f_\mathrm s$. For an interesting discussion here on dsp.se see Nyquist Frequency Phase Shift


Is it possible to design a digital filter that rejects 60Hz noise while keeping step response in the time domain below one mili-second Have been trying to design IIR and FIR 60Hz notch filters but all have step response way above 10ms.

You can trade off steepness of the notch against time domain ringing. That's a fundamental trade off and there is no way around it. When using an IIR notch you can simply adjust the "Q" (quality factor). A Q of 3 will give a step response that's roughly 10ms long.


Maximum Likelihood Derivation as you can see on page 314 the derivation of the ln[f etc.] equation is done via the matched filter that i have also provided below. i would like to know how this equation is used to derive the ln[f etc.] from its previous probability functions, and how it then finally becomes in the form of ||a-b||<||a-b||

any help would be greatly appreciated. thanks

The formula you need, which I don't see in the scanned pages, is the N-dimensional Gaussian distribution of the received vector given that the signal $s_m$ was transmitted:

$$f(\vec{r}|s_m)=\frac{1}{(\pi N_0)^{N/2}}\exp\left[-\frac{\sum_{j=1}^N(r_j-x_{mj})^2}{N_0}\right]\tag{1}$$

The maximum likelihood receiver seeks to maximize this conditional probability by choosing the appropriate $s_m$. Since the logarithm is a monotonic function, you can as well maximize the logarithm of (1):

$$\log(f(\vec{r}|s_m))=-\frac{N}{2}\log(\pi N_0)-\frac{1}{N_0}\sum_{j=1}^N(r_j-x_{mj})^2=\\=-\frac{N}{2}\log(\pi N_0)-\frac{1}{N_0}||\vec{r}-s_m||^2\tag{2}$$

Since the first term in (2) is constant (i.e. independent of $s_m$), maximizing (2) is equivalent to minimizing the distance $||\vec{r}-s_m||^2$. If $||\vec{r}-s_m||^2$ is minimal, then obviously also $||\vec{r}-s_m||$ is minimal. So we arrive at the ML decision rule given in your book, i.e. choose the signal $s_m$ for which the distance to the received signal is a minimum.


Matlab: Help in understanding if the example for Maximum Likelihood Equalizer is properly functioning I am stuck in the part where I need to apply the MLSE Equalizer with Viterbi  code. The equalizer is an optimal. I am using the Communication Toolbox  http://www.mathworks.com/help/comm/ref/comm.mlseequalizer-class.html#properties

and cannot understand pretty much if the simulation is working properly or not. Basically, the equalizer is an FIR model and the channel is a multipath (not Rayleigh).  Following the example in the documentation, I cannot understand

(1) What this statement does? Why is it calling the ML Equalizer here?

  modSignal = step(hMod, data);

(2) In the example code, there is no statement for additive Additive White Gaussian noise. The example given in the documentation does not clearly mention if the additive white Gaussian noise has been incorporated into the transmitted signal. It only says Introduce channel distortion . What kind of distortion is unclear.

(3) Is the example properly functioning because the simulation BER = 0 for all SNR levels. I am sceptical if I am doing it properly.

I shall be really obliged if some help is provided in how I should properly use the MLSE equalizer. Thank you.

UPDATE : The code throws error :

Error using comm.MLSEEqualizer/step
Complexity mismatch with input 1; expected complex, got real.

Error in MLSEEqualization_test (line 47)
      eqSignal = step(hMLSEE,chanOutput);

I wanted to see if I can apply MLSE for real signal also such as BPSK and to minimize the use of built-in functions for data generation, modulation and demodulation. This would help me to understand how the theoretical concepts can be put together in the code. That is why I am not using the communication toolbox for data generation and other steps; I am using only the MLSE EQualizer module provided in the toolbox. How can I go about this error message? The channel coefficients are real, noise is real and the data is also real. Thank you.

N = 10^3 % number of bits or symbols
rand('state',100); % initializing the rand() function
randn('state',200); % initializing the randn() function

n = 1/sqrt(2)*[randn(1,N) + j*randn(1,N)]; % white gaussian noise, 0dB variance
Eb_N0_dB = [0:10]; % multiple Eb/N0 values
h=[1 0.6 0.3]';

L = length(h);
for aa = 1:length(Eb_N0_dB)

    % Transmitter
ip = rand(1,N)>0.0; % generating 0,1 with equal probability
data = 2*ip-1; % BPSK modulation 0 -> -1; 1 -> 1

   chanOutput = filter(h,1,data');

    %% ----- ADDED
    EbN0 = Eb_N0_dB(aa);
    SNR_linear  = 10^(EbN0/10)*log2(M);
    SigPWR      = 1; %unit circle symbols
    NoisePWR    = SigPWR/SNR_linear;
    NoiseFactor = sqrt(NoisePWR/2); %half power complex half imaginary

    hMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
  %  hDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
    const = step(hMod,(0:M-1)'); % BPSK constellation

            %% ----- ADD NOISE !!
  %    noiseVec   = NoiseFactor*randn(N,1)+1j*NoiseFactor*randn(N,1);
   noiseVec   = NoiseFactor*randn(N,1);
      chanOutput = chanOutput+noiseVec;

    hMLSEE = comm.MLSEEqualizer('TracebackDepth',10,...
                      'Channel',h, 'Constellation',const);

    hError = comm.ErrorRate;

 % Equalize the channel output and demodulate
      eqSignal = step(hMLSEE,chanOutput);
      demodData = real(eqSignal) > 0;
   %   demodData = step(hDemod,eqSignal);
 % Compute BER
      errorStats = step(hError, data, demodData);

    fprintf('Error rate = %f\nNumber of errors = %d\n', ...
      errorStats(1), errorStats(2))

    %assume 1 symbol error = 1 bit error (needs Gray coding in real system for this to be true)
    BER(aa) = errorStats(1); % SER ~= BER with undelrying Gray coding
end

%% changed, this is already divided by N*N2
simBer = BER; % simulated ber

theoryBer = 0.5*erfc(sqrt(10.^(Eb_N0_dB/10))); % theoretical ber
close all
figure
semilogy(Eb_N0_dB,theoryBer,'b.-');
hold on
semilogy(Eb_N0_dB,simBer,'mx-');
axis([0 10 10^-5 0.5])
grid on
legend('theory', 'simulation');
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for BPSK modulation');

I had a look at your code, and I know why you get 0 BER.

The thing is that you do not add any noise, you should use Eb/N0 to calculate the noise variance, and add a complex noise on your channel output.

The distortion added by "chanOutput = filter(chCoeffs,1,modSignal);" merely adds linear ISI, by convolving the channel impulse response with the input signal. This linear distortion is exactly known to the equalizer since you pass it the exact channel coefficients. This is never the case in real applications.

Even assuming you do know the channel coefficients, you still need to add noise immediately after the channel distortion and before the equalizer, to get a meaningful result of BER over Eb/N0.

Noise can be generated with randn having unit variance, then you need to multiply this noise by a constant.

I did not spend too much time understanding what you want to verify, but i did correct the script to include noise.

clear
N = 10^4; % number of symbols
N2= 10;   %repeat with same Eb/N0

% Eb_N0_dB = [0:25]; % multiple Eb/N0 values
Eb_N0_dB = (0:10); % multiple Eb/N0 values
nTx = 1;  %SISO
nRx = 1;

M = 4; %QPSK

for ii = 1:length(Eb_N0_dB)

    %% ----- ADDED
    EbN0 = Eb_N0_dB(ii);
    SNR_linear  = 10^(EbN0/10)*log2(M);
    SigPWR      = 1; %unit circle symbols
    NoisePWR    = SigPWR/SNR_linear;
    NoiseFactor = sqrt(NoisePWR/2); %half power complex half imaginary

    hMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
    hDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
 % Channel coefficients
    chCoeffs = [.986; .845; .237; .12345+.31i];
    hMLSEE = comm.MLSEEqualizer('TracebackDepth',10,...
                      'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
 % Create an error rate calculator
    hError = comm.ErrorRate;
    for n = 1:N2
      data= randi([0 3],N,1);
      modSignal = step(hMod, data);

 % Introduce channel distortion.
      chanOutput = filter(chCoeffs,1,modSignal);

      %% ----- ADD NOISE !!
      noiseVec   = NoiseFactor*randn(N,1)+1j*NoiseFactor*randn(N,1);
      chanOutput = chanOutput+noiseVec;

 % Equalize the channel output and demodulate
      eqSignal = step(hMLSEE,chanOutput);
      demodData = step(hDemod,eqSignal);
 % Compute BER
      errorStats = step(hError, data, demodData);
%       nErr=errorStats(1) ;
    end
    fprintf('Error rate = %f\nNumber of errors = %d\n', ...
      errorStats(1), errorStats(2))

    %assume 1 symbol error = 1 bit error (needs Gray coding in real system for this to be true)
    SER(ii) = errorStats(1); % SER ~= BER with undelrying Gray coding
end

%% changed, this is already divided by N*N2
simBer = SER; % simulated ber

EbN0Lin = 10.^(Eb_N0_dB/10);
theoryBer_nRx1 = 0.5.*(1-1*(1+1./EbN0Lin).^(-0.5));
p = 1/2 - 1/2*(1+1./EbN0Lin).^(-1/2);
theoryBerMRC_nRx2 = p.^2.*(1+2*(1-p));

close all
figure
semilogy(Eb_N0_dB,theoryBer_nRx1,'bp-','LineWidth',2);
hold on
semilogy(Eb_N0_dB,theoryBerMRC_nRx2,'kd-','LineWidth',2);
semilogy(Eb_N0_dB,simBer,'mo-','LineWidth',2);
axis([0 25 10^-5 0.5])
grid on
legend('theory (nTx=1,nRx=1)', 'theory (nTx=1,nRx=2, MRC)', 'sim (nTx=2, nRx=2, ML)');
xlabel('Average Eb/No,dB');
ylabel('Bit Error Rate');
title('BER for BPSK modulation with ML equalizer');


How can I calculate the Expectation for a particular known vector? Let us say I have a row vector $X = [x_1, x_2, x_3, \ldots, x_n]$ and another row vector $Y = [y_1, y_2, y_3, \ldots, y_n]$. I want to check whether the two vectors are statistically independent or not.

Now two vectors are said to be statistically independent if $E(X Y) = E(X)  E(Y)$, where $E$ is the expectation operator.

So how would I calculate the expectation value for each of the concerned vectors?

Thank you for your suggestions - I'm really new to statistics.

First off, I would read Examples of Independent and uncorrelated data in real-life, and ways to measure/detect them.
It will likely prove useful to you.
If you are looking for something quick, easy, and useful, I would simply cross-correlate the two vectors and normalize the result.  A small result (less that 0.05?) indicates, but does not prove, independence.  A large result (more than 0.2?) would indicate some amount of dependence.  If the vectors have a non-zero mean, you may need to subtract out the mean before doing the correlation, depending on why they are non-zero mean.
You could also go a more academic route and use the method pichenettes outlines (estimate the probability distribution using a kernel density estimator) in the thread that I linked to at the beginning of the post.


Why discrete time signals are defined as sequence of numbers? When a discrete-time signal is obtained by sampling a continuous-time signal, why we are removing the information about sampling interval and only deal it as a sequence of numbers?

Samples of (regularly sampled) discrete time signals, like pixels of a digital images, are sequences: they are unitless on the ordinal axis (albeit time or space), and unitless in amplitude. A sequence is

an enumerated collection of objects in which repetitions
are allowed and order matters

This is close to what people call raw data.
However, to make a genuine "signal" or "image" object (meaningful in the real world or for devices) that can be better interpreted or manipulated, it is important to keep ancillary information, like: sampling rate, quantization, time-space coordinates, size, channels, amplitude units, calibration parameters, etc. Those additional parameters are found for instance  in the header of most digital sound or image raster formats, and then followed by the unitless digital sequence.
Below is an example from the WAVE PCM soundfile format (converted in png for display). The Discrete signal is typically in the bottom "data" sub-chunk, yet requires the above information to be read, listened to, analyzed.


How should you decide on the number of sweeps to be used in FMCW radar? I am designing an FMCW radar in MATLAB and I would like to know on how should I decide on the number of linear sweeps in FMCW .Does the number of sweeps affect the performance?

The basic performance requirement of FMCW radar is to measure the radial distance and velocity.

The number of sweeps or linear sweeps relates with the velocity measurement, thus if for a constant integration (processing) period you have more number of chirps (which can be achieved by decreasing the sweep duration) the max measurable velocity is increased.

For the case when the number of sweeps are increased by increasing the processing (integration) time the velocity resolution is enhanced.

In order to read more detail about the basic of the FMCW parameters, I would recommend to go through this paper FMCW Signal Processing


change FFT bins position I'm sampling a signal at 44100kHz, and I have 1024 samples on which I perform an FFT.
My bin resolution is therefore around 43 hz/bin, so the frequencies I get energy for are multiples of 43hz: 43Hz,86Hz,129Hz etc...

The problem is that I the frequency I need to calculate the energy for lies between those bins! Specifically, I have:
BIN[389]= 16750 Hz
BIN[390]= 16793 Hz
but I wish to calculate the energy for 16781Hz, 16824Hz,...

Is there a mathematical manipulation I can pull off to get the info I need? zero padding? Or perhaps multiplying my sample by an exponent before performing an FFT?

Thanks for the help

It's actually possible to pre-multiply the signal with a complex exponential and thus shift the center frequencies of the FFT bins. Particular useful is multiplication that shifts by half a bin since it maintains complex conjugate symmetry for real valued time domain signals and you end up with N/2 complex at N/2 frequencies instead with N/2+1 frequencies with N/2-1 complex numbers and two real numbers (Nyquist and DC).

However, it really depends what you are trying to do. The interpretation of the FFT magnitude as a "energy" is a bit problematic. You are NOT just getting energy at the center frequency but a weighted sum of of the energies over a bunch of frequencies, where the weighting function has roughly a sin(x)/x shape. The frequency selectivity of an FFT is actually not particularly good and it's not a great tool to actually calculate the energy in a specific frequency band.

Zero padding can also be used to increase the frequency resolution, however you need to keep in mind that the FFT assumes discrete and PERIODIC signals so zero padding can only be done on signals that are trailing off sufficiently (like the impulse response of a filter) or when windowing is applied to mitigate the transition between signals and zeros.


Right algorithm for fourier transform on physical heights I have data from a LIDAR unit that I would like to get the spectral density of. Unfortunately, the only thing I remember from my Fourier analysis class are the methods that I know will not work.

The data comes from a 1D LIDAR scan of a (mostly flat) surface, which returns radial distance at evenly spaced $d\theta$'s. I can convert this data to x-y data which looks distinct depending on the type of surface (e.g. grass looks a lot choppier than concrete). I think Fourier analysis would be a good way to distinguish between the different types. Unfortunately the data has some problems that makes it difficult to analyze:

It is unevenly spaced in x.
It has some missing values.
It is not strictly a function in x (if there is some overhang, an earlier part of the sweep might hit under the overhang to get a distant x value, and a later part could strike the overhang, to get a closer value.)

It has occurred to me to do a Fourier analysis on the original $r-\theta$ data, but since I am concerned with the deviations around perfectly flat ground, I am not sure how I could use it.

Your question seems to have a couple of nested issues: First off, the computation of the Power-Spectral-Density is as straight forward as the computation of the signals' Discrete Fourier Transform, (DFT), followed by its absolute magnitude squared. The computation of the $O(N^2)$ DFT is accomplished very efficiently using the $O(NlogN)$ FFT algorithm, that comes in canned form from many different libraries.

Thus, if $x[n]$ is your original domain signal, (sampled in space, time, or whatever other quantity), the PSD is given by:

$$
PSD[k] = |\sum_{n=o}^{N-1} x[n] \ e^{\frac{-j \ 2 \pi n k}{N}}|^2 = |X[k]|^2
$$

(Where $X[k]$ is the Discrete Fourier Transform (DFT) of $x[n]$). (Also note, the above formula is for same-length DFTs).

The second aspect of your quandary seems to rest on:

What domain to best represent your orignal LIDAR measurements in.
How to deal with non-evenly spaced samples.
How to deal with missing samples.

Regarding (1): I do not see why you cannot retain the data in its original $r - d \theta$ form. You are after all after radial distance per spatial sample, and this is exactly what LIDARs do for you. You want to then use the DFT as a transformation that will ostensibly magnify and separate different surfaces for you, by virtue of the repetitiveness in the spatial domain, corresponding to particular deltas in the Fourier domain. (Grass with high spatial repetitivity, VS concrete with low repetitivity).

Regarding (2) & (3): Both those problems can easily be solved using linear regression, and/or interpolation, whichever one suits your fancy. That is, you would collect all your $r- d \theta$ data, replete with missing samples and non-even ones. What you seek however is a uniformly spaced sampling on a $d \theta$ grid. You can do a simple polynomial fit via Least-Squares-Estimation, (LSE). That is, your objective will be to decipher the co-efficients of a best-fit polynomial (of degree constrained by the LIDAR physics and spatial bandwidth of the target surface), that best fits your data.

For example, let $d \theta$ be the LIDAR independent variable. You have elected to use a degree-2 polynomial ($D = 2$) to best fit your data in the Least Squared (LS) sense. You have collected a bunch of points from your LIDAR, at some given $d \theta$'s and corresponding radial distances $r$. You would like to find the co-efficients $p_i$, that give you the line of best fit. That is, you would like to find the $p_i$'s in:

$$
r = p_0 + p_1 \ (d \theta) + p_2 \ (d \theta)^2
$$

Thus let $\bf{r}$ be a $1$ x $N$ vector of recorded radial distances to the scanned surface. Let $\boldsymbol{d \Theta}$ be a $(D+1)$ x $N$ data matrix, such that the first row is all $1$s, the second row is composed of each $d \theta$ for each value of $n$, the third row is composed of each $(d \theta)^2$ for each value of $n$, etc. $N$ of course, is the total number of data points from the LIDAR scan of the surface. Finally, let $\bf{p}$ be the $(D+1)$ x $1$ vector of co-efficients to be solved for. Then, the LSE solution is:

$$
\bf{p} = (\boldsymbol{d \Theta} \ \boldsymbol{d \Theta}^T)^{-1} \boldsymbol{d \Theta} \ \bf{r}^T
$$

With the coefficient vector $\bf{p}$ in hand, you can now go back and solve for $r$'s at any $d \theta$ grid. That is, simply construct a uniformly spaced $d \theta$ grid, and solve for corresponding $r$'s. Then, you will be ready to transform this regressed result into the Fourier domain, and look for characteristic spatial periodicity corresponding to different materials.


How to choose a fixed adaptation step for decision feedback equalizer How do I choose $\mu$ - a value of a step size for adaptation in a decision feedback equalizer (DFE) with adaptive reference control (ARC)? For a regular adaptive FIR filter an adaptation step depends on a variance of input, but DFE gets symbol values as input. Another problem is that DFE employs a sign-error LMS adaptation algorithm, while ARC - a simple LMS. So do $\mu$ value for $ARC$ and $h_i(k+1)$ can be the same?
The algorithm for the DFE and ARC is the following (based on this M.Sc. thesis by Mirna Hage as a reference).
With $x^{in}_k$ being a $k^{th}$ sample:
\begin{align}
x_k &= x^{in}_k - \sum_{i=1}^N h_i(k) \cdot a_{k-i}\\
a_k &=
\begin{cases}
+3 & \text{if} & x_k > 2 \cdot ARC_{k-1}\\
+1 & \text{if} & 0 \leq x_k \leq 2 \cdot ARC_{k-1}\\
-1 & \text{if} & -2 \cdot ARC_{k-1} \leq x_k < 0\\
-3 & \text{if} & x_k < 2 \cdot ARC_{k-1}
\end{cases}\\
e_k &= x_k - ARC_{k-1} \cdot a_k\\
ARC_k &= ARC_{k-1} + \mu\cdot e_k \cdot a_k\\
h_i(k+1) &= h_i(k) + \mu\cdot \operatorname{sgn}(e_k) \cdot a_{k-i}
\end{align}
where

$N$: DFE length
$x_k$: DFE output
$ARC_k$: Signal level estimate
$a_k$: Symbol estimate for a 4-level PA
$e_k$: estimate error
$h(k)$: DFE weights.

I tried $\mu=0,125 * 2/N$, but it seems to be wrong in some cases.

I'm leaving the answer here, if somebody ever stumbles upon the same question.
About LMS: Both DFE and ARC employ Least Mean Squares (LMS) adaptive algorithms: DFE is an adaptive filter and ARC can be considered as an additional tap of that filter. Hence the adaptation step size should be selected accordingly.
The general update equation for the LMS algorithm is the following:
$$
\mathbf{w}(n+1) = \mathbf{w}(n) + \mu\cdot e(n) \cdot \mathbf{x}(n),
$$
where $\mathbf{w}(n)$, $e(n)$ and $\mathbf{x}(n)$ are respectively filter weights, an error and tap values at $n^{th}$ sample and $\mu$ can be a constant scalar.
With $p$ being the filter order, the bounds for $\mu$ can be defined as
$$
0 < \mu < \frac{2}{(p+1)E\{|x(n)|^2\}},
$$
where $E\{|x(n)|^2\}$ represents the power in $x(n)$ and can be estimated using an average:
$$
\hat{E}\{|x(n)|^2\} = \frac{1}{N}\sum_{k=0}^{N-1} |x(n-k)|^2.
$$
Then $\mu$ value can be chosen as a fraction of $\mu_{max}$, e.g. $0.125 \cdot \mu_{max}$.
However, in the sign-error LMS, which is an LMS algorithm with reduced complexity with the following update equation:
$$
\mathbf{w}(n+1) = \mathbf{w}(n) + \mu\cdot sign(e(n)) \cdot \mathbf{x}(n),
$$
replacing an error value with its sign changes the magnitude of the correction. If $e(n)>1$, then the magnitude of the correction is smaller compared to the standard LMS. Hence the $\mu$ value can be chosen to be larger.
To the questions:

For example in a question, the input signal is a sequence of evenly (presumably) distributed PAM4 symbols, i.e. "-3", "-1", "1" and "3", which gives a power equal 5. This gives $\mu_{max} = \frac{2}{(p+1)\cdot 5}$.

DFE and ARC employ different LMS algorithms: sign-error LMS and standard LMS respectively. This calls for separate $\mu$ values for DFE and ARC. Hence the recommended adjustment would be to use sign-error LMS for both functions. Not only does this simplifies the $\mu$ value selection, it also simplifies the further implementation.

Considering the need to simplify the hardware implementation of the algorithms, $\mu$ can be chosen to be a power of two. For this particular example, $\mu$ equal to $\mu_{max}$ rounded to one significant binary digit happened to be a good choice.
References: A good reference on the LMS algorithms and step size selection is Chapter 9 "Adaptive Filtering" in Statistical Digital Signal Processing and Modeling by M.H. Hayes. Another good option is a Least_mean_squares_filter entry on Wikipedia.