Why negative harmonics copies in sampling i have a signal of 120Mhz and sampling at 800Mhz. we look at the region of fs/2=400MHz so obvios we have harmonics of the fundamental which is 240 360. But in the document i see that there is another harmonic -4*120M+800M I cant see why we have negative harmonics and their copies. Upon what theory its based?
they say the -120M is the only one copied to the 0-400M region but i thought we could take also
-6*120M+fs=-720M+800M=80MHz
I cant see what why we have negative harmonics and which one is copied to the 0-fs/2 region.
Thanks.

First, about "why we have negative frequencies and which one is copied"...
In a real signal, there is no difference between a negative frequency and a positive frequency in this sampling context. Consider that you have captured 1000 samples that you know contains a single sine wave. You analyze it and find it contains exactly ten cycles of a sine wave. Knowing the sample rate, you determine that it's a 100 Hz sine wave. Now analyze the capture from the other direction—end to start. You would also find that it's 100 Hz, but since we moved backwards through time, it would be -100 Hz.
That is, if you choose to plot negative frequencies, they will always be a mirror of the positive frequencies. You will never find -80 Hz without 80 Hz, for instance. There is never a choice of "which on is copied". If you find somehow that there is a -80 Hz component, you can be sure there is an 80 Hz component.
Sampling is a type of amplitude modulation (AM). Pulse Amplitude Modulation (PAM) in the analog domain, Pulse Code Modulation (PCM) in the digital domain. In the analog domain, it's a multiplication of the signal with an impulse train at the sampling frequency. Looking at just the first harmonic of the impulse train, it's a sinusoid (cosine in particular) at the sampling frequency. In your case this means you are multiplying your signal by 800 MHz (as well as its pulse harmonics, but we need only consider the closest one here).
Amplitude modulation produces sum and difference frequencies, but for this problem we only care about the difference frequencies at the moment. That means if your signal has a 120 MHz first harmonic, sampling will also produce 800-120 HMz = 780 Mhz. And the fourth harmonic of your signal, 480 MHzHz, would produce a new component of 800-480 = 320 MHz. Note that this is exactly the same answer as your "-4*120M+800M". Except that I came up with my answer without considering negative frequencies (800M minus 4x120M).
To summarize, negative frequencies always have positive counterparts. And there is no particular need to consider negative frequencies with what you are looking at—it seems it just adds confusion, but may make more sense in context of the document you refer to.


preprocessing raw ECG data I need to preprocess raw ecg data in R, here is a sample already standardized. I'm not an expert in signal processing nor experienced in working with medical data, so I need concrete answers ideally with code.

Below is a plot of how the data looks like, so perhaps this can give an idea of the problems present.

Question 1)
What strategy is needed in terms of subsampling, filtering to preprocess the data to clean get a clean ECG signal (google how it should look).

Question 2)
Would this strategy work for EMG data also?

I need to preprocess raw ecg data in R, here is a sample already standardized. I'm not an expert in signal processing nor experienced in working with medical data,...

Not being an expert on how the heart works and how its phases manifest themselves on the ElectroCardioGram (ECG) is not a problem. But it would help immensely if you mentioned what is the purpose of this, i.e what do you want to do, what are your specifications and constraints?

For example "We put the driver in the car, they go through the track a few times and then we want to know what was their heartbeat at a specific point in the track".

Excellent, then, do not even use ElectroCardioGraphy (ECG) because:

The electrical environment of the racing car is incredibly noisy already. For comparison, you are trying to measure a signal that has a max amplitude of a few mV through long leads that have to run from the chest of the driver to the ECG box. Just a short distance away from this long electrical loop, you have sparks going on and off, relays turning big loads on and off, electrical motors or stepper motors turning. All of these things create fields that the loop will pick up and they will contaminate your signal.
ECG measures the electrical activity of the heart via contact electrodes. Which means that the skin of the driver would have to be prepared, ideally shaved, cleaned with alcohol, then have a conductive gel applied to it and finally apply the electrodes and fix them with medical tape. Throughout the track, I would not expect the driver to pay attention to the electrodes or, to not sweat. Even in controlled conditions, if you tap the electrodes with your hand, you can observe changes in the signal.

Instead, use Pulse Oximetry. This is the way that most gadgets (e.g. smart watches, fitness wrist bands, etc) obtain Beats per Minute from the heart today. Pulse oximetry needs to have contact with the skin, but not because it measures electricity. Very briefly: It measures the attenuation of light that is proportional to the blood flow below it. Because the heart is pumping blood in a rythmic way, the flow is rythmic as well and this signal can be used to derive a relatively accurate BPM count.

It will also be much easier to cross reference this with the rest of the telemetry because it will simply be a time series of BPM measurements. Pulse oxymeters for people to manage their heart conditions are more accurate than gadget type instruments. But then again, they might get a bit of getting used to by the driver. (For an example, see this)

  Below is a plot of how the data looks like,

No way.

Here (and here ) is what typical ECG signals look like.

To extract BPM, you need to "detect" those spikes which correspond to heart contractions and then measure the timing between them (also known as R-R interval, because the main spike of the ECG is the R peak).

Very briefly: QRS Detectors work by Pre-filtering -> Rectification -> Integration -> Thresholding -> QRS Complex pulse sequence -> R interval extraction -> Beat Detection. For much more detailed information about this, have a look at Moody's work.

What is depicted in the plot is saturation and possibly lots of noise from the contact. There is no trace of "typical" ECG complexes to work with.

  Question 1) What strategy is needed in terms of subsampling, filtering to preprocess the data to clean get a clean ECG signal (google how it should look).

Why subsample? The signal is already in a bad state, why make it even worse? You could try a bandpass filter between 0.1 - 100Hz. This would reduce DC and cut out any high frequency noise but you have a lot of spikes which means that most of time times you will be "hearing" the ringing of the filter.

  Question 2) Would this strategy work for EMG data also?

EMG is even worse because it picks up signals at higher frequencies. What is the application?

If you are trying to see the timing of the driver as to WHEN do they start to accelerate (i.e. leg muscles start working) versus when the car actually accelerates or brakes (i.e. break signal from telemetry), then maybe you can get away with very simple processing using thresholds. But if you are trying to make more fine measurements, you might be back into the same problem of lots of noise from various sources. In the case of EMG you will have to work around those, which means, twisted pair wires, shielding, as short lengths as possible, very good contact with the driver's skin, very good preparation of the driver's skin, perhaps even slightly different electrodes that attach on to the skin rather than simply be in contact with it.

Hope this helps.

EDIT:

I am sorry, this is an unsalvageable situation. I have a couple of comments only,but I do not believe they will be very much of help.

So, here is the full signal you shared (thank you, by the way):

As you can see, there is nothing usable here. My perception is that the ECG amplifier is either saturating (because of induced voltage) or the contact with all three electrodes is not stable. If you notice a Holter monitor, which is used to "catch" abnormal heart functioning wherever it may happen, you will see that the leads are as short as possible with the device being mounted ON the chest. If we were to use the same device in an application similar to yours, we would have to use shielded leads and make sure the pads are firmly attached on the body.

Going a bit closer, it almost looks like "reversed bits"

What I mean by this is: Say for instance, you created your own ECG circuit. So, take an off the shelf instrumentation amplifier, set the gain and then feed that to an "Arduino"ish type of board with some 10bit Analog To Digital (ADC) converter that establishes a serial like channel to a pc or stores values internally. Because you have 10bits (or more) and you would have to send two bytes over the wire (or store two bytes) you would have to decide if they are sent/stored/received/retrieved Big Endian or Little Endian.

So, when you read a binary file from the disk and you get the endianess wrong, your integer values, the levels recorded by the ADC, would look a little like this second figure right here. There are glimpses of what looks like a signal inter dispersed with wild swings from the minimum to the maximum value. Also, your "zero" is not at 0, but maybe you used just two ECG leads.

I am just leaving this here as a comment in case it is what is going on, though I doubt it.

If you look more closely, the signal is just all over the place:

In fact, if you have telemetry from the rest of the car, try to plot together the RPM of the engine and the ECG signal that you have. What I would expect to see is the "density" of disturbances in the signal to correlate with the RPM of the engine, assuming that most of this interference is coming from the electrical system of the car (ignition, injection, throttle stepper motors, etc).

Hope this helps.


Resampling a digital sound signal As a non maths-pro, I'm looking for some pointers.

I am rewriting the audio core in my emulator to improve accuracy, and am getting a bit stuck on the specifics of the technique.

I am trying to resample a sample buffer of 4,466 samples, at a sample rate of 223,807.5 Hz (from a simulated Yamaha SN76489), to a sample buffer of 880 samples at a sample rate of 44,100 Hz, which can then be output to a sound card.

Now from what I understand, as the ratio between 44,100 and 223,807.5 Hz is 5.075, I need to turn this into a fraction, which works out as 203/40, and then interpolate by the denominator (leaving 0 samples inbetween), run a low-pass filter, and then decimate by the numerator, giving me my 880 samples of 44,100 Hz output.

(1) Am I along the right lines here? Before I was just downsampling by a straight 5 (i.e. missing out 4 of every 5 samples) and whilst it worked well enough, some notes are slightly out of key in some games, and there is aliasing. Is this new approach correct?

(2) Also, I am not sure what would do the trick in terms of a low-pass filter. I found the following pseudocode on wikipedia that seems to fit the bill:

    // Return RC low-pass filter output samples, given input samples,
    // time interval dt, and time constant RC
    function lowpass(real[0..n] x, real dt, real RC)
       var real[0..n] y
       var real α := dt / (RC + dt)
       y[0] := x[0]
       for i from 1 to n
           y[i] := α * x[i] + (1-α) * y[i-1]
       return y

But I have no idea how to get RC and dt - what are these in relation to the filter I wish to achieve? As I don't know how to select suitable values here, the output signal is heavily garbled at the moment. If a better/more understandable formula is available then by all means suggest it. Sorry for the newbish question and I hope this is the right place for it. Thanks :-)

you can use this approach provided it does not yield efficiency burdens.

% Note that matlab arrays start from index 1, not 0.
F1 = 223807.5;     % Sampling rate of Yamaha chip
F2 = 44100;        % New sampling rate of standard Audio-CD
M = 203;           % Necessary Decimation ratio
N = 40;            % necessary Expansion ratio
Lx = 4466;         % Length of Audio Input Buffer at F1 rate
Ly = Lx*N/M;       % Length of Resampled Output Buffer at F2 rate

x = sin(2*pi*1234*[0:Lx-1]/F1);  % Just a test sine wave
y = zeros(1,Ly);                 % allocate output signal memory

for i=0:Ly-1              % run per each output sample
   n = i*(M/N);          % compute exact sampling position inside the long buffer
   ind = floor(n);       % convert it into largest integer smaller than itself for array indexing
   d = n-ind;            % take the difference for proportional weighting
   y(i+1) = (1-d)*x(ind+1) + d*x(ind+2); % sum the weighted mixture.
end

figure,plot(x);            %plot each waveforms...
figure,plot(y);


Blind Estimation of Signal Parameter and Noise Variance Let $y[n]= h*x[n] + w[n]$, where $h$ is an unknown but deterministic parameter, $x[n]$ is a BPSK random variable with equal probability of +1 and -1, $w[n]$ are i.i.d. Gaussian with zero mean and unknown variance . Given $N$ observations find the estimate of $h$, and variance of noise.

$x[n]$ and $w[n]$ are independent.

Squaring Method

The input signal $x[n]$ can take two values: $+1$ or $-1$. After being multiplied by $h$, the signal becomes $|h|e^{j\angle h}$ (if $x[n]=+1$) or $|h|e^{j(\pi + \angle h)}$ (if $x[n]=-1$). By squaring this we get, $|h|^2e^{j2\angle h}$. Both points get mapped to this. Now we can pick out the magnitude and phase of $h$. In summary, the steps are:

Square $y[n]$ to get $z[n]=y^2[n]$.
Get the magnitude, $\hat{|h|}=\sqrt{\bigg|\frac{1}{N}\sum_n z[n]\bigg|}$. It was pointed out in the comments that this is a biased estimate with bias $\sigma^2$, so it is not favorable for a large $\sigma$. I describe another way to estimate $\sigma^2$ (without using $\hat{h}$) so that we can use this $\hat{\sigma}^2$ to remove the bias from $\hat{|h|}$.
Get the phase, $\angle h=\frac{1}{2N} \sum_n \angle z[n]$
Use $\hat{h}$ to get the noise variance, $\hat{\sigma}^2=\bigg(\frac{1}{N}\sum_n z[n]\bigg)-\hat{h}^2$. This comes from the fact that $E\big[z[n]\big]=h^2+\sigma^2$ (also see the clustering method is a better way to go about this as it does not depend on any other estimates).

Caveat: This method is valid only for phase offsets within $\pm \frac{\pi}{2}$. Anything beyond that, the received signal becomes indistinguishable and there are multiple phase shifts which could give the same signal.

Clustering Method

If $h$, can be correctly determined, then you can take hard decisions and calculate the noise variance using the distance of the received symbols from the decision you made.

If $h$ can't be determined right away, then this clustering based estimate of $|h|$ and $\sigma^2$ can be used: We know there should be two clusters ($\pm 1$) so we can use a simple method like k-means (https://en.wikipedia.org/wiki/K-means_clustering) with $k=2$. The algorithm is:

Run k-means on the received samples and get a list of labels back labeling each sample to either cluster $1$ or $2$. Let $S_1$ be the set of sample indices in cluster $1$ and $S_2$ be the set of sample indices in cluster $2$.
Now we are going to take all the points in each cluster, shift them be zero mean and take the variance (this will give us a noise variance estimate for each cluster). Then we will average over all of the cluster noise variance estimates to form the final estimate. $\hat{\sigma}^2=\frac{\bigg( \text{var}\big(y[n]-E[y[n]] \big)\big|_{n \in S_1} + \text{var}\big(y[n]-E[y[n]] \big)\big|_{n \in S_2} \bigg)}{2}$
Now we can take the estimate from step $2$ and subtract off the bias term which we estimated here as $\hat{\sigma}^2$ to get the refined estimate $\hat{|h|}_0=\hat{|h|}-\hat{\sigma}^2$. You could also get $|h|$ using the magnitude of the cluster means, $|\hat{h}|=\frac{\text{E}\big(y[n]\big|_{n \in S_1}\big)+\text{E}\big(y[n]\big|_{n \in S_2}\big)}{2}$.

There is also a way get $\angle h$ from the clusters but it too has the criteria that $\angle h$ must be less than $\pm \frac{\pi}{2}$ so I will not detail it here.

Long story short: $|h|$ and $\sigma^2$ are doable, but $\angle h$ is only doable for some cases.

Code: https://github.com/B-William/DSPSE/blob/master/blindEstimationOfSignalParameterAndNoiseVariance.m


Confusion about speech being piecewise stationary From wikipedia:
"Modern general-purpose speech recognition systems are based on Hidden Markov Models. These are statistical models that output a sequence of symbols or quantities. HMMs are used in speech recognition because a speech signal can be viewed as a piecewise stationary signal or a short-time stationary signal. In a short time-scale (e.g., 10 milliseconds), speech can be approximated as a stationary process. Speech can be thought of as a Markov model for many stochastic purposes."

I get that if you break a speech signal into small enough parts, then a every part will make a very specific, basic sound, which I presume, has the same (or very similar) type of spectral envelope, as many other parts around it/close to it.. but I don't quite understand what the stationary process is, especially not in this context.

To be exact, I don't understand what exactly the stochastic process/it's random variables are, or even the motivation for describing a piece of signal as a stochastic process.
I mean the actual input is not a stochastic process, it's a concrete input, it's an observation - is the quote saying that basically any small enough piece of the input can be considered as a certain sample of some stationary process?

Do we call it stationary, because while we know that the small part has some sort of quality of being very repetitive, we don't know exactly where we've seperated it/where it starts? Similar to something like a stationary process of sin(t+Y), Y being a random variable.

I'm pretty confused about this whole thing- if I see a specific sample of white noise, is it correct to call it 'stationary'? I mean, the original distribution that generated this sample, that is a stationary process, but does that mean you call the sample/signal itself 'stationary'?

How exactly does the fact that sound is piecewise stationary imply it a good idea to use HMM? I simply get stuck at that - maybe a slow, thorough explanation of this implication (and what exactly it means that it's piecewise stationary) would be the most helpful to me.

Sorry if the questions are not understandable or not concrete enough. I'm simply very confused and don't know enough about this.

In this context, stationary usually means that the spectral density is (theoretically) the same throughout the observation window.
So, piecewise stationary means that if you divide the signal into short enough time frames, the signal in each frame can be described, accurately enough, as sum of sinusoids in presence of white noise (which is stationary as well).


Different number of modes while performing the EEMD (Ensembled Empirical Mode Decomposition) I am performing an Ensemble Empirical Mode Decomposition (EEMD).
I wrote a MATLAB routine to perform the EMD. I noticed that when adding white Gaussian noise to the input signal, sometimes I obtain a different number of modes, even if the noise variance is not so high. Clearly, if this happens, I cannot sum compute consistent ensemble-averages for each mode.
Do you think this depends more on the nature of the input signal or on my EMD algorithm/routine?

It depends on both on the signal content and on the EMD implementation. EMD can be sensitive to time shifts, especially with strong transients, to noise power. Robust or constrained EMD on single signal or multivariate signals have been developed.
Your question relates to signal morphology, processing purpose, and quality metrics. I doubt they can be answered globally.


What does the phase discriminator portion of the Costas Receiver do mathematically? What does the phase discriminator portion of the Costas Receiver do mathematically?

The output of the I-channel is $ \frac{1}{2}A_C \cos \phi \, m(t) $. Which means for small deviation of phase $ \phi $ , $ \frac{1}{2}A_C \cos \phi \, m(t) \approx \frac{1}{2}A_C \, m(t) $. Or it can be said that $  \frac{1}{2}A_C \, m(t) $ would be attenuated by a small amount. However, we have to keep in mind as the text says (as it not might be as innocuous as it seems):

The output of the Q-channel is $ \frac{1}{2}A_C \sin \phi \, m(t) $. Which means for small deviation of phase $ \phi $ , $ \frac{1}{2}A_C \sin \phi \, m(t) \approx \frac{1}{2}A_C \phi \,  m(t) $.
As per the text, the phase discriminator consists of a multiplier followed by a low pass filter.
Which means multiplication of $ \frac{1}{2}A_C \cos \phi \, m(t) $ and
$ \frac{1}{2}A_C \phi \, m(t) $ would yield:

$ g(t)=\frac{1}{2}A_C \cos \phi \, m(t) * \frac{1}{2}A_C \phi \,  m(t) $

$ g(t)= \frac{1}{4}A^2_C \phi \cos \phi \, m^2(t)$

Now, let's say $ m(t) $ is band limited to $ M $. then the term $ m^2(t) $ in the frequency domain would spread across  $ -2M $ to $ +2M $ centered around $ f=0 $ and the term $ \frac{1}{4}A^2_C \phi \cos \phi $ is a constant.
Now, if $ g(t) $ is subjected to a LPF, then the term $ \frac{1}{4}A^2_C \phi \cos \phi \,  m^2(t) $ would be retained as it is(if the cutoff frequency of the LPF is slightly > $ 2M $).
So, what does phase discriminator in Costas receiver do mathematically?

Also, if the phase error is significant (that is we have no idea at all about the probable phase of the carrier), then how does the analysis change? Because this is a possibility that we might have no idea about the phase of the carrier.

Text Used: Communication Systems By Simon Haykin

This pdf has the answer for you:
The short answer is that you have conflicting assumptions in number 4 above. You use the small angle approximation for sine but not cosine. If $\phi$ is small then $cos(\phi) \approx 1$. So you would end up with only $\phi$ in that term times the amplitude and $m^2(t)$. Once you integrate, or low pass filter, the signal component $m^2(t)$ for long enough that becomes basically constant and you're left with the result of something that is proportional to $\phi$. That makes the VCO advance or retard the frequency momentarily based on the phase error.
I should add that the small angle approximation generally is a good model even for when there are large errors. Generally when the error is large the error term and VCO push things in the right direction, even if the precise amount isn't perfect. Keep iterating on that and you eventually get close enough for the small angle approximation to be valid.


How should someone start to learn machine learning and computer vision for mobile applications? I am very interested in learning Machine Learning and Computer Vision for mobile apps, for iOS and Android.
But I have no clue or even a slight idea on how to even begin. I know this is heavy on math and some theory I am not familiar with.

Should I use a library for this kind of thing, when I build a mobile app? Or should I do native coding? For example, to create  photo effect from a camera?
What are all things I need to know to be able to learn this?

For machine learning there is a very good free MOOC on Coursera that can give you the basics of machine learning, with a lot of examples. Prerequisites are basic maths and computer programming comprehension.

For computer vision I would advise you to first learn about DSP in general.
This book has had great reviews everywhere and although I haven't read it yet, I'm pretty sure it's a good start in the vast area of signal processing. There is also an online book, The Scientist & Engineer's Guide to Digital Signal Processing, with plenty of information on digital filtering.

Once you have understood the fundamental theory of Machine Learning and DSP, you can start learning programming in C++. One very famous library for computer vision is OpenCV. You will be able to find many tutorials online.


EM algorithm and Kalman filter I am new to the subject of Kalman filtering and therefore my question might seem trivial.

I see that there is a tight connection between Kalman filter and EM algorithm when one wants to predict the next state of a given signal. My question is what is the connection and how we can combine them in the prediction process?

Thanks

EM-algorithm jointly estimates the state space model parameters as well as state estimates. The E-step is a Kalman filter, which uses the current estimates to predict the new states. The M-step uses this result in an MLE procedure to obtain the parameter estimates.

I would suggest that you read the following paper, which explains a generative model, from which all of the tools such as HMM, Kalman Filter, VQ and etc. are derived:

  Read Roweis, Ghahramani, "A Unifying Review of Linear
  Gaussian Models", Neural Computation, Vol. 11, No. 2, 1999.


What are the units of amplitude of my ultrasonic recording equipment I'm not familiar with the topic of sound. I have ultrasonic audio recordings sampled at 250kHz (16bit) in WAV format. The equipment used to record them is linked here:
http://www.avisoft.com/usg/usg116Hn.htm

I've use the matlab programming language to plot the audio signal. While the X axis is just going to be the sample number (divide this by the sampling rate to get time), I don't quite know what the units for the amplitude are.

Typical values (which can be positive or negative) have very high precision and are of the order of 10^-4 or 10^-5. An example amplitude value is:

-3.051757812500000e-05 (at least according to Matlab)

Anyone have any idea of what units I am dealing with? I did not use the equipment myself, my supervisor did. When I asked him what units the amplitude were, he advised I look at the link, but I am not too certain I understand.

Given that the gain of the system and/or ADC seems to be unknown, the scale of the sample values can only be determined from samples of a calibrated source of known value at a known distance from the transducer.  With a variable enough gain knob (or AGC) the scale of the samples could be anything.


why is moving towards higher frequencies so important in frequency spectrum? As we know, as frequencies becomes higher, bandwidth becomes higher.And, according to channel capacity theorem, channel capacity increases with higher bandwidth.

Also,  energy is directly proportional to frequency(E=hf). So, if frequency increases, signals possesses higher energy and can travel far.

Which one of them would be more relevant to prove our need of higher frequency? Why not the other one?

So, if frequency increases, signals possesses higher energy and can
  travel far.

Not really, it's more the other way round:

When transmitting signals with radiowaves in free-space, most telecommunications engineers use the Friis equation which gives the amount of received power $(P_r)$ according to emitted power $(P_t)$, distance between emitter and receiver ($R$) and wavelength ($\lambda$) :
$$
\frac{P_r}{P_t} = G_tG_r\left(\frac{\lambda}{4\pi R}\right)^2
$$
$G_t$ and $G_r$ being the transmitter and receiver's antenna gain.
We see with this equation that $\frac{P_r}{P_t}$ is proportionnal to $\frac{1}{f^2}$ (because $\lambda = \frac{c}{f}$) which means that the higher the frequency, the higher the loss (and ^2 higher !).
The equation $E = h\nu$ is more related to quantum physics and has no practical use in telecommunications. Maybe it means that higher frequency signals are "more energetic" that low frequency ones, I don't really know honestly, maybe someone else* can answer better.

So the need to increase frequencies is mostly due to :

Saturation of the currently used spectrum : if you read international and national regulations and standards, you'll see that being able to emit on the RF spectrum is very expensive, because it's actually saturated, every frequency is reserved for a certain application (e.g. GSM signals is my country are transmitted in 200 KHz wide canals around 900 Mhz-ish frequencies). Still in my country, a new mobile phone operator had to pay billions of € to the governement to access some mobile comm. bands. Increasing the frequencies will give everyone more "room" to transmit radiowaves.
Increasing data rates, because having a higher carrier frequency also allows to widen the communication channel
Miniaturization (higher frequencies = lower wavelengths = smaller antennas and systems)

*edit :

  your claim on signal energy and frequency relation, which is based on photo-electricity (E=hv) effect is improper. First of all, the physical nature of waves has nothing to do in developing the mathematical notion of waves and their spectral representation by periodic signals. And even though it's true that the magnitude-squarred |x(t)|2 signal power relation is basicly adopted from power calculations in physics, that power (and energy) is implied to be the energy carried out by the waveform of a classical physics kind.

from Fat32's comment above


Butterworth filtering behaving unexpectedly -MATLAB I'm new to signal processing but I'm trying to understand how the Butterworth filter works. To do so, I did the following: I assumed a sampling frequency of 100Hz and wanted to attentuate everything after half the Nyquist frequency. The code I used is the following

[b,a] = butter(2,0.5);
fs = 100; %Sampling frequency
t = 0:1/fs:1-1/fs;
y=0.123*sin(2*pi*20*t+0.234)+0.123*sin(2*pi*10*t); %Pair of sine curves
plot(abs(fft(y))) %This behaves as expected, with two peaks at 10 and 20
dataOut = filter(b, a, y);
figure,
plot(abs(fft(dataOut)))

I expected the second plot to be identical to the first because there are no frequencies above the Butterworth filter's cutoff anyway. However, the signal is artifically inflated (by a small amount) for all components below 25Hz and then goes to zero. See the original signal and the filtered signal below.

Is this normal behavior? If so, it is effectively adding noise in my passband, is it not? Also, this effect is a lot worse if I use a higher order Butterworth filter.

The effect you observed is mainly due to spectral leakage. It is a property of the DFT and it is only indirectly related to the Butterworth filter. Note that there are an integer number of periods of both sinusoidal components of y inside the DFT window. That's why you see two nice peaks in the frequency domain without any leakage. If you truncate y only by a few samples, you'd get spectral leakage even without the Butterworth filter, because now the period of the signal doesn't line up anymore with the DFT window:

plot(abs(fft(y(1:97))))

The output signal of the Butterworth filter shows some transient before it converges to its steady state. So it is not purely sinusoidal, and - more importantly - the number of periods of the strong sinusoidal components are not aligned with the length of the DFT window. If you remove the initial transient and if you make sure the DFT window aligns with the periods of the sinusoids you get two perfect spectral peaks:

plot(abs(fft(dataOut(11:100))))


Edge detection and area measurement I am a novice in image processing. As a part of my school project, I have to do the following steps in Matlab on the below image.

Detect the edges of a leaf.
Calculate the area of a leaf (excluding the stem portion)

My friend suggested me an idea of calculating the total number of pixels in the leaf (excluding stem portion) divided by the number of pixels per cm2.

But I don't know how to do these in Matlab. Someone, please guide me.

Here's my suggestion:

Convert the image to grayscale
Define a threshold and binarize image (only black and white)
To remove the stem apply an erode operator (the size will depend on your image)

Function syntax: https://www.mathworks.com/help/images/ref/imerode.html

apply a dilation operator (same size as in previous step)

Function "imopen" combines both steps: https://www.mathworks.com/help/images/ref/imopen.html

Sum all over the columns (considering that leaf pixels have value 1 and if you are using Matlab colsum=sum(a) )
Sum the previous vector and you will have the area ( area=sum(colsum) )


Understanding the definition of mean/autocorrelation I was studying about the definitions of mean, expected value and autocorrelation. I wanted to verify my understanding the evaluation of mean, expected value and autocorrelation. At the same time to verify the definition of autocorrelation as the convolution of the signal with this delayed copy mathematically which I was not able to.

Here is my analysis:

Basically, I see mean as an average of all values in a population. For example, consider an experiment where we measure the outcome of the random variable that can go anywhere between $[1,100]$. By performing the above experiment $50$ times, we get $50$ outcomes:$a_1, a_2, \ldots, a_{50}$ which has values ranging from $1$ to $100$.

The mean of the above random variable is given by:
$$
\mu = \frac{a_1 + a_2 +\ldots + a_{50}}{50}
$$

The expected value is given by:
$$
E(x)=(1*\textrm{probability of getting 1})+(2*\textrm{probability of getting 2})+\ldots+(100*\textrm{probability of getting n})
$$

which can also be represented as:
$$
E(x)=\prod_{k=1}^{100} kP_k
$$

where $P_k$ is the probability of getting $k$. Now taking the actual definition of an expected value which is:
$$
E(X)=\int_{-\infty}^{\infty} xF_x(X)dx
$$

The above seems somewhat satisfying as we could relate the above with the example - taking $x$ as one of the outcome and $F_x(X)$ as the probabilities (actually it's a probability density function to be accurate).

For a random process, the above expands to (at time $t_1$):
$$
E(X_{t_1})=\int_{-\infty}^{\infty} x_{t_1}F_{X_{t_1}}(x_{t_1})dx
$$
Similarly, expected value can be computed at several time instances since the random variable changes every time instant.

Now, lets take a Strict-Sense Stationary Random Process where expected value (mean) is constant and autocorrelation depends up on the time difference. So, the definition of Autocorrelation is given as per wikipedia:

  $$
R_x(\tau)=\frac{E\left[\left(X(t_1) - \mu\right)\left(X(t_1+\tau) - \mu\right)\right]}{\sigma_{t_1}\sigma_{t_2}}
$$

The above is with Normalization. Taking out the normalization, we have:

  $$
R_x(\tau)=E\left[X(t_1)X(t_1+\tau)\right]
$$

Till now it's pretty clear. Now applying the definition of an expected value in the above expression:
$$
R_x(\tau)=E\left[X(t_1)X(t_1+\tau)\right]=\int_{-\infty}^{\infty} x_{t_1}F_{x_{t_1}}(X_{t_1})x_{t_2}F_{x_{t_2}}(X_{t_2})dx\quad\text{where}\quad t_2 = t_1+\tau
$$

But in Wikipedia, it's given as:

  $$
R_x(\tau)=E\left[X(t_1)X(t_1+\tau)\right]=\int_{-\infty}^{\infty} x(t)x^*(t-\tau)dt
$$

This is where I am stuck at - I am not able to match the last two expressions. I am very well aware that autocorrelation of a random process is convolution with itself at a time lag of $\tau$ but I am not able to relate the expressions mathematically. I am stuck here for several days.

The definition of the autocorrelation function $R_x(\tau)$ depends on the nature of your $x$.

If $x$ is a deterministic signal with finite energy then: $$R_x(\tau)=\int_{-\infty}^{+\infty}x(t)x^*(t-\tau)dt$$
If $x$ is a deterministic  signal with finite average power$^{(1)}$ then: $$R_x(\tau)=\lim_{T\to+\infty}\frac{1}{T}\int_{-T/2}^{+T/2}x(t)x^*(t-\tau)dt$$
If $X$ is a (real-valued) random process then: $$R_X(t_1,t_2)=E[X(t_1)X(t_2)] = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}u_1u_2p_X(u_1,u_2)du_1du_2,$$

where $X(t_1)$ and $X(t_2)$ are two random variables (the realizations of the random process $X$ at times $t_1$ and $t_2$, respectively), $u_1$, $u_2$ are their respective values and $p_X(u_1,u_2)$ is the joint probability density function of these two random variables. This statistical autocorrelation becomes only dependent on $\tau = t_1-t_2$ if the random process is stationary and we can write:
$$R_X(\tau)=E[X(t)X(t-\tau)].$$

If, furthermore, the random process is ergodic then the statistical autocorrelation equals the temporal autocorrelation, i.e:

$$E[X(t)X(t-\tau)]= \lim_{T\to+\infty}\frac{1}{T}\int_{-T/2}^{+T/2}x(t)x(t-\tau)dt^{(2)}$$

$^{(1)}$: Signals with finite energy are also known as energy signals whereas signals with finite average power are also known as power signals. Moreover, the average power of an energy signal is zero and the energy of a power signal is infinite.

$^{(2)}$A random signal is by convention considered as a power signal.


Comfort Noise Insertion in VAD algorithms I am in a intention of Porting Voice Activity Detection algorithm in real time hardware,so i was successful in classifying speech and silence/noise but now the problem is replacement of silence frames with comfort noise(preferebly any stationary noise)

so can any one suggest any opensource CNG algorithms or how a comfort noise can be of any type(pink,white etc)

to do this usually we use a rand() functions but where as these type of functions are not available with real time DSP so how to port such type of Algorithms?

You may want to try a two-pronged approach.  First, have a standard pre-packaged CN to inject (like jan and Jason R suggest) at the start of the conversation/exchange.

As the conversation progresses, gather low-volume data, perhaps as part of your noise reduction algorithm I assume you already have, by sampling the lowest volume signal coming in between VAD end and start points. Compare your low-volume samples to your standard CN, and replace (or convolve to simply dilute your synthetic CN) if it matches some threshold or characteristic you established. All you need is some rand and normalization functions, and perhaps FFT, all of which are readily available in C (and assembly).


Fit a line pattern on curve with unknown number of points I've got a sample curve which ends theoretically with decreasing exponential. The curve end falls into noise. The sample points are given in log scale.
What I want to do, is to find and fit the linear part of the log curve to retrieve the exponential factor. The trick is that I don't know the starting point neither the ending point of the linear part of the log curve.

The strategy I'm using is to fit a line at each point with at least say 20 points and until the end of the curve is reached. Then from all these regressions, I keep the one with the best determination coefficient.

I did several tests and found that the RMS error increases with the number of points systematically as the curve becomes more and more noisy, so the extracted factor is always computed on the minimum number of samples (20 in this example case).

My question is, is there a more efficient method to compute this factor ? Does increasing the minimum number of samples improve the accuracy of the fit ?

Here a sample curve I want to fit

The computed slope regarding the position in the curve

And the corresponding RMS error

A useful approach could be to take the numerical derivative of the function. Then you will get some constant line at some $y=a$, where $a$ is non-zero. This would help to eliminate all the other parts and to determine the starting and eding points of the line.

About the fit Matlab has a polyfit command. If the curve of interest is linear than basically this command with parameter $2$ as input will give you the pair $(a,b)$ s.t. the best $ax+b$ fits into the numberical numbers that you provide.


Mathematical advantages of the ZT, DTFT and DT? I apologize if this question is too general to answer concretely, but I was hoping more to perhaps be pointed towards some resources that could help more extensively. Essentially, I have a Discrete-Time Signal Processing Final exam coming up soon, and although I understand and am familiar with the mathematical descriptions of the ZT, DTFT and DFT, I was hoping to gain some intuition about which one to use in which scenario.

To me it seems that the DFT is a discretely sampled version of the DTFT, and the DTFT is the ZT specified on the unit circle. In this way, they are all essentially the same thing with increasing generality towards the ZT.

However, there are mathematical subtleties associated with each one (can Parseval's only be applied for DTFT and DFT?), even between the bilateral and unilateral ZT, for example with respect to the time shift property; when you shift a sample to the left in the time domain, then you must subtract the value of the ZT at n = 0 in the unilateral case, but no such subtraction (or addition) is required in the bilateral case. So why use the unilateral ZT, when it is simply the same as the bilateral ZT for causal sequences?

More generally, I was hoping for some insight into when each of these 4 transforms is most useful. For practical purposes I suppose the final transform, the DFT, is useful because computers process discrete information so it is useful to have samples that are discrete in both the time and the frequency domain, but this kind of answer is not what I am looking for (although I would appreciate more information about this as well). In my final I will be given a system in the form of perhaps a difference equation or a block diagram and asked to find the impulse response sequence, or perhaps given the impulse response sequence and asked to find a description of the system. Sometimes I will be asked to use a specific transform, but other times it is up to me to decide which is the quickest/easiest path to what I am looking for.

I understand that part of this kind of analysis is an art and the intuition comes with time, but maybe some of you know of some resources which can accelerate that process.

Thanks in advance for any responses.

My friend, i hope you can gain a proper command of those techniques for achieving an excelent result at your exam. Of course the question is broad and more discussionfull, but the purpose of this is "clarify" not to turn this into a "closary" of questions.

Indeed. there are a lot of variants for the transforms, but only a few of them have a clear direct meaning, and the others are generalizations seldom used.

At first, the standard, classical FT, is the main entry for your analysis.
$$\mathbf{F}=X(f)=\int_{-\infty}^{\infty}x(t)e^{-i2 \pi ft}dt$$
$$\mathbf{F}^{-1}=x(t)=\int_{-\infty}^{\infty}X(f)e^{i2 \pi ft}df$$
That is the base, plain and simple, sweet, the FT and the IFT are the same (of course exchanging the signs and variables :) ) and with the easiest expressions $\mathbf{F}\{e^{-\pi t^2}\}=e^{-\pi f^2}$, $\mathbf{F}\{1\}=\delta(f)$, $\mathbf{F}\{\cos(\pi t)\}=\frac{1}{2}\delta(f-1)+\frac{1}{2}\delta(f+1)$, $\mathbf{F}\{\Delta(t)\}=\Pi(f)$, $\mathbf{F}\{\Pi(t)\}=\Delta(f)$, $\mathbf{F}\{|||(t)\}=|||(f)$ (where $\Delta$, $\Pi$ and $|||$ are the triangle, step and comb functions:)), and so on.
You use this on the handwritting, papers, reports, analysis, etc.

We can safely ommit the DTFT and DFFT, and the tricky $i$ and $k$ symbols for discrete time and frequency and keep $t$ and $f$, discretized, with the sampling time $dt$  and the frequency periodic spacing $df$ (again with the same notation as the integral). Thus for the DTFT direct transform is:
$$\mathbf{F}_{DT}=\sum_{t=-\infty}^{\infty}x(t)e^{-i2\pi ft}dt$$
and the DFFT inverse is:
$$\mathbf{F}^{-1}_{DF}=\sum_{f=-\infty}^{\infty}X(f)e^{i2\pi ft}df$$

because the $|||(t)$ (and $|||(f)$) comb function:
$$|||(t)=\sum_{i=-\infty}^{\infty}\delta(t-i)$$

is all what we need for recovering the related expressions and cumbersome hard to learn theorems regarding sampling (DT) and perioditizicing (DF). And you recover all your symmetries just by multiplying by your combs:

$$\mathbf{F}_{DT}(x(t))=\mathbf{F}(x(t)\frac{1}{dt}|||(\frac{t}{dt}))=X(f)\frac{1}{dt}\frac{1}{F}|||(\frac{f}{F}), F=\frac{1}{dt}$$

Hence, at my second place, you have your well forgotten ZT, unilateral, again considering the $t$ as a discretized variable under a sampling time $dt$:
$$X(z)=\sum_{t=0}^\infty x(t)z^{-t}$$

You have a DAQ with a lot of sensors and controllers, and you have some fixed sampling time and you need to know the continuous parameters of the system based on the discrete behaviour, or tuning some controllers, or identifying a system. Soon or later you need to pass onto the Z domain and check if your IIR filters are stable or nor, or if your FIR smoothers are too short on delay, and comparing your ARMA models with your differential equations.


Are Fast Wavelet Transform coefficients constant if the input signal frequency coefficients are constant? I'm currently studying the Fast Wavelet Transform. As I currently understand, the Fast Wavelet Transform is implemented as a QMF filter bank where the frequency resolution decreases as the signal is low-pass filtered and sub-sampled. Pictured below:

However, I'm curious as to why this "works".
Let's say I have a true 750Hz signal that is sampled at 2kHz. The amplitude of this 750Hz signal is 2Vpp.
The Fast Wavelet Transform says that I can high-pass filter with a pass region of 500Hz - 1000Hz, decimate by a factor of 2, and it will output coefficients that correspond to the signal frequency. Because I know the signal is constant frequency 750Hz, I would hope for constant coefficients of magnitude 1 at Level 1.
But I'm confused. Because we're subsampling at 1KHz (half of the original sample rate), we're essentially aliasing the 750Hz signal into the 0Hz - 500Hz range. It seems like the coefficient value will vary depending on where we're sampling. However, because we know the true signal is a constant 750Hz signal, it would be desired that the wavelet coefficient is also constant.
How exactly do the coefficients of the Fast Wavelet Transform correspond to the true frequency content of a signal?
Edit: Are there any good resources for understanding the general properties of Wavelet coefficients? e.g. should a wavelet coefficient remain constant if the signal has constant frequency content?

[Specific answer] If you have a real sine or cosine, the discretization and finite support may entail issues in the possibility of having a constant: edge samples or periods that cannot be divided by powers-of-two can affect the  wavelet coefficients.
But more important: let us suppose the discrete sine period is well chosen. Then a low-pass or high pass filter on a sine will yield something close to a sine (because complex exponentials are invariant under linear systems). So even after downsamling, the coefficients in a wavelet subband are likely to behave as sinusoidal values.
If you expect "constant values", you are likely looking at a complexified and absolute-valued representation of the coefficients.
With a real signal and discrete wavelet, one will likely see ripples on several subbbands, instead of a constant.
[Generic answer] The purpose of most transforms (whether fast or not) is to help in handling or interpreting the content of some data, or classes of data. Invertible transformations keep all information, sometimes with redundancy.
Discrete wavelet transforms (DWT) are meant to capture information in a non-redundant way from not-so-stationary signals. They often  are really bad as capturing information from a purely periodic sine signal.
As most wavelet filter bank filters are not perfect, and the signal is of finite length, filtering and representing coefficients will suffer artifacts.  So in a DWT, a sine would spread on many subbands, and suffer from aliasing, which is well known. However, despite the above caveat, the magic of wavelets can however recover the original sine, from the spread and aliased coefficients.
So :

don't expect DWT coefficients of a sine to be constant
don't expect them correspond to the true frequency content of a signal

They are not meant for that. However, if a signal has constant frequency content, but phase shifts (thus non-stationary), then  wavelets (maybe continuous once) can be great at detecting their location or onsets.


What does it mean for audio to be 'sped up'? How can I undo it? I very simply have a .wav audio file, but there is a note with it that says "This audio has been sped up by 10x".

I am wondering, what does it mean exactly from a signal processing perspective for an audio signal to be 'sped up' by some factor, and secondly, how can I undo this? Is it possible?

It could mean:

That the metadata in the .wav file wrongly indicates that the signal was recorded at a sample rate 10x higher than the true value (for example it is advertised as having been recorded at 80kHz while it has been recorded at 8kHz). Simply fixing the metadata undoes this.
That the signal has been downsampled by a factor of 10, but is still tagged as having been recorded at the original sample rate. This can be undone by dividing by 10 the playback sample rate (by tweaking the metadata); or by upsampling the signal by a factor of 10 without changing the playback sample rate. Note that in any case, 90% of the bandwidth of the original signal (the high frequency content) has been lost and cannot be restored (or only very approximately through spectral band replication).
That a time-stretching (time-scale modification) algorithm has been applied to it. This is a lossy process, and cannot be undone without severe loss at such a high stretch factor. The "antidote" is simply to use any time-stretching algorithm with a stretch factor of 10x.


Digital demodulator (SDR) for receiving legacy analog modulated wave I understand that in SDR, after the ADC the signal processing is fundamentally based on I/Q modulated waves where any discrete value that can be transmitted corresponds to a symbol which is represented as amplitude and phase through I and Q signals.
How can a modern SDR receive legacy AM waves?
The SDR uses I,Q components whereas legacy AM doesn't have any such component applicable.
I presume SDRs are able to receive and demodulate existing legacy AM,FM transmissions - please correct/add.
Do current SDRs receiving AM waves switch to a different technique wherein it doesn't use I/Q demodulation?

Sure, an SDR can transmit and receive conventional, analog AM. I've done it many times as an introductory exercise in my courses.
A quadrature signal $s_{\text{QAM}}$ is made up of two signals, $s_I(t)$ and $s_Q(t)$, that are transmitted over the same bandwidth and modulated in quadrature (i.e. with carriers phase-shifted by 90 degrees):
$$s_{\text{QAM}}(t) = s_I(t)\cos(2\pi f_c t) + s_Q(t)\sin(2\pi f_c t)$$
There is no restriction on the shape of these signals, beyond the physical limits of the hardware or the channel (in bandwidth or in power). In particular, they can be a DSB-LC AM signal!
$$s_{\text{QAM}}(t) = (A_m + m(t))\cos(2\pi f_c t) + 0\sin(2\pi f_c t)$$
The regular downcoversion and filtering process applied to $s_{\text{QAM}}$ will produce $A_m + m(t)$. You can also apply a traditional envelope detector. That's the beauty of SDR: it does not care about the meaning or structure of the in-phase and quadrature signals; you define that in software.
To conclude, note that not only it is indeed possible to transmit/receive AM with an SDR; you can actually transmit/receive two AM signals on the same bandwidth!
$$s_{\text{QAM}}(t) = (A_{m1} + m_1(t))\cos(2\pi f_c t) + (A_{m2} + m_2(t))\sin(2\pi f_c t)$$
You select one of the IQ branches and ignore the other, and then you can listen to one "station" or the other. Of course, this assumes that you are transmitting your own quadrature AM signals (for example, over one of the ISM bands).


Phase information of a signal when creating a spectrogram Sorry for the uninformed question but I need some help understanding this. I'm trying to understand what phase information is and it would be very helpful if someone could correct my understanding of it.

Fourier transforms are used to create spectrograms (I understand the vanilla DFT). They allow us to see the different frequencies that make up a signal, and with this knowledge we can plot create an amplitude spectrogram. During this process, phase information is discarded. Correct?

Therefore we can not reconstruct the original audio signal because we are missing phase information. Now what I don't understand is what this phase information is telling us. The way I understand it, is that phase tells us how the different frequencies are aligned/arranged in the original signal? Is this correct?

Yes, you are correct. If you take a DFT of square wave and only look at the amplitudes, doing an inverse DFT but using different or random phases for the sine components, it does not look like square wave in time domain any more. But it will have a matching spectra. Kind of like two racecars on a circular track, they might always have same velocity, but without knowing the initial position, or even the direction they are going, it will be impossible to know if they side by side or on the opposite side of track.


Resources on Non-Gaussian Channels I understand why most of the classic textbooks only consider Additive White Gaussian Noise channels in digital data transmission. Unfortunately in my case noise is not Gaussian and it might even be colored (non-white).

I spent a couple of hours searching around but only found some special cases that were eventually approximated by Gaussian models. I appreciate if you could point me to some resources that non-Gaussian noise is considered, in concise and direct manners.

Some notes: I'm interested in evaluating the performance of a receiver in the presence of Gaussian and other non-Gaussian noise sources. For example, I'd like to see the effect of the ADC quantization noise passed through some filter, like an equalizer (and there is not much support to use Central Limit Theorem to make Gaussian assumption), or some residual jitter from my clock recovery loop that again doesn't fit nicely to a Gaussian model.

There are many varied impairments in a transceiver that we could call "noise", (see my graphic below), and ultimately I resort to simulations with good noise and impairment models for each effect, that is very dependent on the receiver type and operation.

This question is to broad to give a specific answer. However, let me describe my treatment of colored noise in the receiver in case this adds some insight. In most cases my approach is to first have a thorough understanding of the frequency domain characteristics of the noise process involved, and then to evaluate the equivalent frequency response of the receiver translated to baseband prior to final detection, and integrate the total noise across that bandwidth to get the total noise power in band. The signal bandwidth will have a low frequency limit (typically set by carrier and timing recovery loops, or the duration in time that a packet is observed-- you need to observe a signal for infinite duration to go down to true DC so it will always be a high pass!) and then an upper frequency limit set by the symbol rate if the receiver is optimized (typically a percentage over R/2 where R is the symbol rate and the percentage is based on the pulse shaping factor), or the actual pre-detection bandwidth if not optimized.

Here is an example with phase noise specifically. It is usually clear that the receiver will filter out the higher frequency components beyond the channel bandwidth, but often there is confusion with how to handle a low frequency limit, and knowing that the phase noise keeps increasing as we approach DC, we could get very different answers depending on what low frequency limit we choose!

If there is a carrier tracking loop involved, the low frequency limit is often established by that tracking loop, since it works by tracking out phase rotation (so slow rotations or movements of phase within the loop BW of this tracking loop would be tracked out).  Depending on the loop this could be a first order high pass response as I show in the figure below.

It is the resulting phase noise density after being multiplied by this "filter" of the receiver that I would integrate (and double as half of the band is shown) to assess the phase noise impact on my receiver SNR.

The key point in this example is that with a thorough understanding of the frequency domain characteristics of the noise (the phase noise response, and beyond the scope of this example, how it maps to my receiver channel of interest), and a clear understanding of the receiver both in the low frequency limit and the upper frequency limit of the baseband channel with all operations prior to final detection; one can integrate the resulting colored noise density over the band of interest to establish the total noise in the receiver prior to detection.  I also do time domain histograms to capture and access the effects of different noise distributions, but in most case the combination of multiple independent noise sources makes the Gaussian estimate reasonable.


Zero order hold interpolation and Nearest-neighbor interpolation? Is there any difference between Zero order hold interpolation and Nearest-neighbor interpolation
I want to perform zero order hold interpolation in MATLAB,but there isn't any information about zero order hold in the description available on MATLAB help page link below,also shown in attached snapshot
So ,in MATLAB,if we want to implement Zero order hold interpolation,can we use the method"Nearest-neighbor"?
https://www.mathworks.com/help/matlab/ref/interp1.html#btwp6lt-1-method

Zero-order hold and nearest-neighbor interpolation are conceptually the same thing; when extending the original dataset to cover new data points, you are just picking a value in the same area of the signal to use as the interpolant, instead of something more sophisticated.
The main difference I could see between the two is that zero-order hold is usually used in the context of interpolation of a 1-D signal, particularly to describe the action of a digital-to-analog converter. In that model, the zero-order hold is usually rendered as starting at the corresponding sample time $n$ and extending until the next sample instant $n + 1$.
Nearest-neighbor interpolation is usually drawn somewhat differently, where the region of influence of a particular sample is centered on the sample position, so instead of the value being held from sample times $n$ to $n + 1$, that sample value would instead influence the interpolated output from $n - 0.5$ to $n + 0.5$.


Two real time signals convolving This might be a stupid question but is it possible to convolve two real-time signals together? I know that generally for running convolution you have the IR and the block of the real time signal and then you convolve them both, making sure you have zero-padded, and then you overlap and add these, but would it be possible, theoretically, to convolve two real-time signals? Like, have a block of one and the block of the other and convolve them real-time, and if so, would the process be fast enough to create something musically dignified?

Thanks

Gabe

Convolution is what happens when a signal passes through an LTI system.
$$y(t)=x(t)\star h(t)=\int_{-\infty}^{\color{red}\infty}x(\tau)h(t-\tau)d\tau$$
The keyword in this question is real-time. The associated keyword in the answer would be weather the LTI system is causal or not.

If the system is causal, then $h(t)=0, \ \forall t<0$, and therefore $h(t-\tau)=0,\ \forall \tau>t$ which means
$$y(t)=\int_{-\infty}^{\color{red}t}x(\tau)h(t-\tau)d\tau$$
It means we only need the inputs only from $-\infty$ up to the present time $t$ which is all available for a real-time implementation.

So the convolution can be implemented in real-time if the system we implement is causal.

To conclude and generalize this idea for two signals, if at least one of the signals are right-sided (which is practically the case), then the real-time implementation is possible.

Regarding your other question, it really depends on the delay performance of your system.


Splitting a magnitude across multiple bins in an IFFT By Georg-Johann own work, CC BY-SA 3.0.

I’m trying to synthesise a sine wave using an IFFT at a frequency between  bin centre frequencies.

Here are the details

-Samplerate : 44100
-FFT Size   : 512
-SignalFreq : 110 hz
-FreqPerBin : 86.1328125
-Bin1Freq   : 86.1328125
-Bin2Freq   : 172.265625

This desired Signal Frequency means that the signal is at midway between $\textrm{bin}_1$ and $\textrm{bin}_2$ (magnitudes) therefore I am wanting to split the synthesis magnitude across 3 bins: $\textrm{bin}_0 , \textrm{bin}_1$ and $\textrm{bin}_2$.

What am I doing wrong?  Here’s how I’ve approached my calculation.

First I calculated the synthesis frequencies percentage between the bins.
0.554195011337869  = (SignalFreq − Bin1Freq) ÷ (FreqPerBin ÷ 2)
Then I multiplied the percentage by PI ($\pi$) to convert it to radians (note I have assumed that bin centre frequency to the bin frequency midway between $\textrm{bin}_1$ and $\textrm{bin}_2 = \pi$1.74105497627516 = (0.554195011337869 * PI)
Then process it through the $\mathrm{sinc}$ function 0.566059652593099 = Sin(1.74105497627516)/1.74105497627516
Which give me the scaling value of $\textrm{bin}_1$
To calculate bin 2, I take the radian value of $\textrm{bin}_1$ and add PI to it (because this moves the location into the next bin for the $\mathrm{sinc}$ function )4.88264762986495 =  (bin1 radian value) 1.74105497627516 + PI
Process it through the $\mathrm{sinc}$ function-0.201845607081636 = sin(4.88264762986495)/4.88264762986495
Which gives me the a value for $\textrm{bin}_2$.  (I’m assuming this should be the absolute value)
To calculate bin 3, I take the radian value of $\textrm{bin}_1$ and subtract PI from it-1.40053767731463 =  (bin1 radian value) 1.74105497627516 - PI
Process it through the sinc function0.703687584403628 = sin(-1.40053767731463)/-1.40053767731463

My initial magnitude for $\textrm{bin}_1$ was $1.0$ so this means that splitting $\textrm{bin}_1$ across $\textrm{bin}_0, $\textrm{bin}_1 and $\textrm{bin}_2$ looks like this

bin0 = 0.703687584403628
bin1 = 0.566059652593099
bin2 = 0.201845607081636

This doesn’t seem right to me.  It seems like bin0 (apart from this bin wanting to have a value of 0.0) has a value higher than the bin matching the frequency of the signal.

What am I doing wrong?

Your frequency in term of bin number is approximately:

SignalF / binDelta = 110 / 86.1 = 1.28

So, for the positive frequency image :

w0 = pi * (0 - 1.28)
w1 = pi * (1 - 1.28)
w2 = pi * (2 - 1.28) = pi * 0.72
wx = pi * (x - 1.28)
binmag(x) = sin(wx) / wx

Also, one needs more than 3 points to avoid windowing artifacts.

But the big problem with Sinc interpolation for any f near bin 0 (DC) is that the signal plus its complex conjugate mirror image at the nearby negative frequency -f both wrap across bin 0, and thus constructively or destructively interfere within all those nearby bins, depending on the phase of f measured relative to the center of the window.  See http://www.nicholson.com/rhn/dsp.html#4 for a more complete Sinc (or Dirichlet) interpolation formula.  A similar problem occurs for frequencies near bin N/2 for an IFFT of length N. This creates a problem with IFFT (re)synthesis of any fractional frequencies near DC or N/2.


generating spectrograms in python with less noise I'm trying to do something similar to Shazam (paper here) in Python. Basically the idea is to use FFT with sliding window to transform both WAV files and recorded audio to the same spectral representations for recognition purposes.

The only problem is the spectrograms I extract from certain songs just seem to generate lots of useless peaks that become false positives when heard in other songs. They look a lot like noise. Examples below.

Spectrogram of SONG_1 [channel 1]:

Spectrogram of SONG_2 (dominates and generates WAY more peaks) [channel 1]:

Spectrogram of SONG_3 [channel 1]:

Spectrogram of SONG_1, recorded through my laptop mic (looks like noise) [channel 1]:

Here's the code I use to generate the spectrograms:

import pylab
from scipy.io import wavfile

fs, frames = wavfile.read("song1.wav")

channels = [
    np.array(frames[:, 0]),
    np.array(frames[:, 1])
]

# generate specgram
Pxx, freqs, t, plot = pylab.specgram(
    channels[0],
    NFFT=4096,
    Fs=44100,
    detrend=pylab.detrend_none,
    window=pylab.window_hanning,
    noverlap=int(4096 * 0.5))

I've implemented the hashing, storage, recognition, etc as well but they are rather useless unless the extracted peaks are distinct to the song heard.

I should mention my peak extraction algorithm is not just a threshold - I find points in Pxx that are the highest value in a neighborhood of around 15-20 adjacent points.

What am I doing wrong? How can I get more distinctive spectrograms out of audio sample channels?

None of these spectrograms show music. I would validate your algorithms by substituting a known signal, e.g. a 440 Hz pure tone. That should be a single line in your spectrogram. Audacity can generate that as a .wav file, but you can also numerically generate the signal in Python. Next, mix two known tones and check if your spectrogram is lineair (i.e. showing two lines).

Audacity can also show you what the spectrogram should be for your inputs.


How to find the impulse response of filter that outputs the autocorrelation of the input? I am trying to solve question 9.61 in the attached image.

I guess, by evaluating the autocorrelation for a specific instance,  say $ t = 0 $ , that the impulse response, $h(t)$, should be : $$ h (t) =  x(-t) $$.

But can I get this from equating the convolution integral to the given output integral, i.e.

$$y(t) = \int_{-\infty}^{\infty} x(\tau)x(t + \tau) d\tau = \int_{-\infty}^{\infty} x(\tau)h(t - \tau) d\tau $$ .

Is it then right to equate as below?

$$ x(t + \tau) = h (t - \tau) $$

Thank you!

HINT:

Note that

$$\int_{-\infty}^{\infty}x(\tau)x(t+\tau)d\tau=\int_{-\infty}^{\infty}x(\tau-t)x(\tau)d\tau$$


In SDR could you isolate I and Q I'm just curious if this could be a way of transmitting stereo?

Yes I and Q are in quadrature and therefore orthogonal so they are naturally isolated. For example you could consider QPSK as two independent BPSK channels transmitted on the same carrier frequency (or you could consider it having twice the spectral efficiency as BPSK since you can send twice as much data at once within the same bandwidth, and at the same SNR level).

To maintain isolation in the SDR you need to have good IQ balance or correct for IQ quadrature and amplitude error. Correcting for analog errors with software as depicted in the diagram below is easy to do: if the data content is equally distributed in the long term then the data can be used directly, and this would be a reason to "whiten" the actual data with scrambling prior to transmission. Otherwise training sequences could be used for such calibration of the receiver hardware- to the extent that is needed to be sufficient for audio purposes. Amplitude is adjusted by comparing the rms magnitudes of I and Q, and phase is adjusted by correlating I and Q (average of the product of I and Q over all samples) and scaling by the amplitude (to get phase).


Using a Butterworth high-pass filter together with Fast Fourier transform I have a signal and I'm interested in identifying the dominant frequencies greater than 1/3600 Hz. To do this, I use Fast Fourier Transforms and examine the locations of the peaks in the frequency domain, which reveals that the largest peaks are below this frequency. As a noob, it seems like there are two potential options:

Simply ignore the FFT results outside my desired range and focus on my region of interest.
Filter the signal (e.g., with a Butterworth high-pass filter) before I perform the FFT.

What are the pros and cons of these approaches? Is there a third (or fourth, fifth, sixth...) option that is better than both of these methods at achieving my goal?

Edit:

As the comments say: it really depends on your signals. If the energy in the "don't care" region is significantly larger than the energy at the frequency of interest, it's helpful to pre-filter it first. Otherwise you may see "spectral leakage" of the high energy components into your target area.

On the downside: the pre-filter will affect the group delay the transient response of your overall system.

This being said: if the filter is not too steep and the cutoff frequency not too close to your target frequencies, pre-filtering is generally a good idea. The pre-filtering is much cheaper than the FFT, so it have almost no impact on CPU or memory requirements.


Different results for seperable filtering v/s direct filtering on 2-D image in MATLAB I compared the difference between seperable and direct filtering for MATLAB's cameraman image. The difference is huge. I would like to understand the reason for this. Code below.

I = imread('cameraman.tif');

% Do seperable 1-D filtering
H = 11;
sigma = 2;
g1 = fspecial('gaussian',[1 H],sigma);

I2 = double(I);
I_seperable = double(I);
for i = 1:size(I,1)
   I2(i,:) = imfilter(double(I(i,:)),g1,'circular','same','conv');
end
for j = 1:size(I,2)
  I_seperable(:,j) = imfilter(double(I2(:,j)),g1,'circular','same','conv');
end

% Do direct 2-D filtering
g2 = fspecial('gaussian',[H H],sigma);

% How much do the results differ ?
I_direct = imfilter(double(I),g2,'circular','same','conv');
delta_I = double(I_direct) - double(I_seperable);
norm(delta_I,'fro')

The result of last computation

ans =

    2.1879e+03

Turns out that when doing filtering by columns, i.e. for the following section

for j = 1:size(I,2)
   I_seperable(:,j) = imfilter(double(I2(:,j)),g1,'circular','same','conv');
end

we need to make sure the filter is aligned along columns. Changing g1 to g1' makes the difference almost zero as seen below.

ans =

   1.9782e-11


Windowing after filtering a discrete signal to get it back to its original length What happens if you filter a discrete bandpass signal of length $N$ with a bandpass filter (Hann truncated sinc) of length $L$ and apply a window of length $N$ to the filtered signal of length $N+L-1$ to get back the signal to its original length $N$?

You cannot "really" apply a window whose length is different from the size of the signal, since it resorts to a sample-wise product. However, you can for instance apply a composite window, made of a standard window $w$ for the first $N$ samples, and filled with zeros for the remaining samples (to the right). If $w$ is a rectangular window, you get with conv.m what you would have with filter.m, plus the flat tail, as explained by @MattL; the red lines (filter.m) and blue crosses (conv.m) are matching perfectly:

Of course,  you could use a non-rectangular window (here the Hann), and you get a nicer apodization at the edges (the signal is flattened to $0$), but you could as well window the output of filter.m


Can we calculate moments (mean, var, skewness, kurtosis) of a signal in frequency domain? The computation of the mean, variance, skewness and kurtosis of a time series in the time domain is straightforward from their definition formulae.

I need these statistic values in different frequency bands. Currently I apply different bandpass filters to the signal and repeat to compute the statistics. I'm curious if it is possible to compute them from the spectrum in the frequency domain?

If your time series $x(t)=s(t)+n(t)$ where $s(t)$ is a deterministic signal,and $n(t)$ is independent Gaussian noise,  it has a pdf of the form
$$
p(x(t))=\frac{1}{\sqrt{2 \pi} \sigma} \exp\left( -\frac{1}{2} \left(\frac{x(t) - s(t)}{\sigma}\right)^2\right)
$$
the first moment is
$$
E\left\{ x(t) \right\} = s(t)
$$
which is NOT the DC component of the Fourier Transform of a sample. The first moment is not a constant and the time series is not WSS, and the first moment is not equal to the time average.


Lesser bins histogram looks more equalised than higher bins histogram after histogram equalisation? I tried to do a histogram equalisation on an image, and then display its histogram. I realise that when I show the histogram with 10 bins, it looks more "equalised" than when I show the histogram with 256 bins. Here's is an image of the 10-bin and 256-bin histograms of the same image after running a histogram equalisation:

So again, in the image above, clearly, the histogram with 10 bins looks more equalised than the histogram with 256 bins. Since both the histograms reflect the colour pixels for the same image, is there reason as to why the 10-bin histogram looks more equalised than then 256-bin histogram?

I'm not sure if this is the best answer to your question out there, but I would disagree that your 256-bin histogram isn't equalized. Think of it this way.

Each bin represents the count of pixels of a certain intensity. When you're performing histogram equalization, you're moving that bin to a new intensity value, according to the empirical distribution (i.e. histogram of the original image). You're shifting all of those pixels to a new value, which means that if 30% of all your pixels had a certain value, all 30% will again be in the same bin, just with a different value.
Histogram equalization tries to equalize the density of pixel values within any neighborhood. For example, since that huge bin with 30% of all pixels cannot be made smaller (all values of those pixels change to the same new value), equalization will try to move all other bins away from it, so that that neighborhood around the big bin frees up a bit. This is exactly what happens in the picture on the right. Notice that all tall bins have large spacing between them, while all short bins are stuck close together. This ensures that on average you have the same amount of pixels "per unit length of histogram's x-axis" if you will. This brings us to the last point.
When you're using less bins, (10 in your case), you're basically summing or averaging bins of your high-resolution histogram into a low-resolution histogram. Remember that the amount of pixels in a certain intensity neighborhood is what's really equalized, and that really manifests itself when you're summing or averaging components in one of those neighborhoods. Because taller bins have more space between them, and shorter bins have less spacing, you're getting approximately equal height bins when you sum or average. This shows that histogram equalization works as intended.

Hope this helps.


How to compare many bandpass filters (for their signal propogation delay)? Let's say we have an audio signal sampled at 96 kHz, and we want to compare several bandpass filters to find the one with the lowest rectangular-envelope-propagation-delay (see graph below).

The band is interest is $[f_0-30\,\text{Hz}, f_0+30\,\text{Hz}]$, with $f_0 = 1000 \,\text{Hz}$, and we aim for an attenuation of, say $-50\,\text{dB}$ at $1300\,\text{Hz}$.

Is there a rule of thumb or a ready-to-use calculator (either as online website, or in Python) to find the best filter with the lowest propagation delay, for example in the case of a signal modulated by a rectangular envelope?

Of course I can run tests in Python, with many filters, and many variations of these parameters:

FIR filters (how many taps?)
Butterworth filters (which order, ...)
Elliptic, etc.

but it's a bit long, and a bit random.

Example: here with a Butterworth of order 2, we have a delay of about ~ 25ms for the propagation of the envelope (the x-axis is in ms indeed).

PS: is this more or less close to group delay?

For a linear phase FIR filter, the group delay will be about the same delay as what you have defined to be a rectengular envelope delay.

For a nonlinear phase IIR or FIR filter (like the Butterworth), the envelope delay is associated with the dominant time constant of its impulse respone $h[n]$.

Indeed what's happening in your plot is that the filter reaches its steady-state response, after the transients decay. And that would happen after a few time constants. For the relation between group delay of the filter and the time constant, it can be said roughly that, short transition bandwidths in the frequence response (sharp filter or narrow band filter with large time constants) yield more aggressive phase changes which also yields larger group delays at the passband frequencies.

A minimum-phase filter would yield the minimum delay (among the same frequency response magnitudes). A minimum-phase filter's impulse response has the fastest decay and begins more aggresively than the other filters; i.e., it's energy is concentrated at the very beginning of its impulse response samples.


Extracting different frequency bands (alpha, beta, gamma) from brain activity data I am new to signal processing and despite googling about my question, I still couldn't find the solution.

I have an Electrocorticography (ECoG) dataset. This consists of recording brain activity from many channels. The data is stored in a matrix. Data in each row relates to the brain activity recorded from each channel. Brain activity data was recorded with 10kHz frequency (Therefore, the continues 10,000 elements of each row is the brain activity data of 1s)

Similar to many analysis such as, I need to extract powers of different frequency bands (1-8 hz, etc) from this data and run some regression anlaysis.

For example, as a sentence from a previous study:

  A fast Fourier transformation (FFT; EEGLAB v5.03) was performed for each 1-secondsignal to obtain the power of each of the 3 frequency bands(2–8, 8–25, and 80–150Hz) for each electrode

I am still confused how should I apply the fft function of matlab to correctly obtain power of different frequency bands.

Considering I only have this matrix of dataset (each row: brain activity recorded from each channel) and also the fact that I am not very familiar with methods in signal processing, I appreciate if someone can kindly help me how should I do such an analysis.

MATLAB's fft returns the FFT of the columns of a matrix X, so by transposing your input matrix before feeding it to the algorithm should work fine: F = fft(X.'); and then you will have the FFT for each channel of yours.

The FFT can be used to obtain the power spectrum. Here is a link with examples and explanations.

Once you have your power spectrum, then there are plenty of ways to extract different frequency bands (mostly based on filtering).


Recovering a signal after nonuniform sampling Let $x(t)$ be a bandlimited signal such that $X(j\omega) =0 $ when $|\omega|>M$. Also $p(t) = p_1(t) - p_1(t-\Delta)$ is a nonuniformly spaced periodic pulse train where $$p_1(t) = \sum_{k = -\infty}^{+\infty}\delta\left(t - \frac{2\pi k}{M}\right), \quad \text{with}\quad \Delta= \frac{\pi}{2M}$$

Let $x_p(t) = x(t)p(t)$ and apply an ideal low-pass filter with cutoff frequency $\omega_c = M$ to $X_p(j\omega)$. The result is $z(t)$. Design a system which recovers $x(t)$ from $z(t)$.

My try:
It's easy to see that
$$P(j\omega) = \left(1 - e^{-j\Delta \omega}\right)\left(M\sum_{k = -\infty}^{+\infty}\delta(\omega - kM)\right) = M\sum_{k = -\infty}^{+\infty}\left(1 - e^{-j\Delta kM}\right)\delta(\omega - kM)
$$
Also we have
$$
X_p(j\omega) = \frac{1}{2\pi}X(j\omega)\star P(j\omega) = \frac{M}{2\pi}\sum_{k = -\infty}^{+\infty}\left(1 - e^{-j\Delta kM}\right)X\big(j(\omega - kM)\big)
$$
Since $\Delta= \frac{\pi}{2M}$ we have
$$
X_p(j\omega) = \frac{M}{2\pi}\sum_{k = -\infty}^{+\infty}\left(1 - e^{\frac{-jk\pi}{2}}\right)X\big(j(\omega - kM)\big)
$$
After low-pass filter we get
$$
Z(j\omega) = \begin{cases}
\frac{M}{2\pi}\bigg[(1+j)X\big(j(\omega - M)\big) + (1-j)X\big(j(\omega + M)\big)\bigg],  & \lvert \omega\rvert < M\\
0, & \text{O.W}
\end{cases}
$$
I've got stuck here. How can we recover $x(t)$? I tried to use phase shifter and $z(t)\cos(Mt)$ but it didn't work.

I agree with your result for $Z(j\omega)$. Apart from scaling, what you have is the original spectrum with positive and negative frequencies swapped and multiplied with factors $1-j$ and $1+j$, respectively. In order to restore the original signal, we need to add right-shifted and left-shifted versions of the spectrum, while getting rid of the complex factors. Final lowpass filtering removes the redundant components in the frequency range $[M,2M]$.
The necessary shifting and complex scaling is achieved by multiplication with
$$e^{jMt}(1+j)+e^{-jMt}(1-j)=2\cos(Mt)-2\sin(Mt)\tag{1}$$
Modulation with $(1)$, lowpass filtering with cut-off frequency $M$, and appropriate scaling restores the original signal.


Matlab: Estimating power spectral density of an experimental data? How do I estimate PSD in matlab ? There are lots of methods to estimate, depending on application.

I am using an experimental data to identify faults ( e.g. interturn fault of electric motor), which gives periodic/regular waveforms (pre-fault) but erratically different (or could have some periodicity) waveforms (post-fault).
In this application, my voltage signal, v(t) is sampled and collected using DAQ systems.

Questions:

Can v(t) be considered a wide-sense-stationary random process (WSS-RP) ?

I understand that there are multiple ways to get PSD.
 What is the best suggestion ? Is it ...

By fourier transform the autocorrelation, if only v(t) is WSS-RP. [as described by Wiener-Khintichine Theorem]
By

splitting the signal,v(t) into multiple frames( or segments)
applying window function : multiply each segment by a hamming window (i.e. convolution in freq. domain )
FFT each of the windowed signal
either multiply its conjugate or absolute square, to get the spectral density of the segment
resulting spectra are averaged to get the ultimate averaged spectrum
Just like the explanation here Fig 9.1, Sect.11: Averaging Periodograms, pg14

What is the difference in the mathworks' application note in (www.mathworks.com/help/signal/ref/spectrum.html)

Can someone please guide in the codings ?

Feel grateful in the midst of dsp gurus !
kit

MATLAB has built-in functions taking care of the steps you mention in 2. You can check out the pwelch function (here) which uses Welch's method for PSD estimation.(here)
The choice of segment length, number, overlap and windowing function presents a trade-off between bias and variance.


Standard deviation of spectral flatness -- so what am I measuring, conceptually? In my never-ending quest to identify snores, I've found that "spectral flatness" seems to be a fair measure of signal "quality".

I'm computing spectral flatness as the geometric mean of the FFT power $(R*2 + I*2)$ data points divided by the arithmetic mean of the same points.

I then (a little twist here) am computing the running (over 50 frames) arithmetic mean and standard deviation of the spectral flatness and computing a "normalized" standard deviation as the running standard deviation divided by the running mean.

For my samples I find that this metric is greater than about $0.2$ (ranging up to $0.5$ or so) when audio is "good" (ie, I have reliable tracking of the breathing/snoring sounds of a sleeping subject) and it generally slips down below $0.2$ when audio is "in the mud".  (I can improve on this discrimination somewhat by using a threshold that moves with other factors, but that's presumably a different topic.)  I also observe that the measure goes over $1.0$ when there's substantial background noise (eg, someone enters the room and rustles about).

So, my basic question is:  Is there a name (beyond "normalized standard deviation of spectral flatness") for what I'm measuring, and can anyone offer a conceptual explanation of what the metric "means"?

(I've tried a dozen other metrics for signal "quality", and this one seems to be the best to date.)

Added:  I probably should admit that I don't have a particularly good conceptual handle on what simple spectral flatness is measuring (just the Wikipedia article), so any further explanation of that would be appreciated.

Since you are interested in the "flatness" of your spectrum, in fact, you are interested in how close your signal is to a white noise (that by definition has a flat spectrum + random phases). If you step back, one measure would be the "distance" of your observation to the white noise reference.

The obvious measure in terms of information theory is Kullback-Leibler divergence. You do not need to understand every part of it, but it measures in bits (if you use log base 2) the distance between both distributions.

The good thing in your case is that your reference is flat, so that what remains is the entropy of your spectrum. There are many existing implementation (for instance in scipy).

Note that you are still on the safe side: if your distribution is approximately gaussian, both measures (entropy and std) will be proportional. The entropy is however more general and more principled. As an extension, you will be able to generalize to other types of noises (1/f for instance).


Extracting data from matlab graph I do not seem to know how to extract from a Matlab plot. For example in the plot given below, I have a simulated $P_b$ for Rayleigh fading and theoretical AWGN channel Pb plot against SNR ($\large \frac{E_b}{N_o}$).

My question is how do I get the value of $\large \frac{E_b}{N_o}$ for a particular value of $P_b$. For example, if I want the read the $\large \frac{E_b}{N_o}$ from the graph to get a probability bit error for the Rayleigh channel, what code can I use?

The data that is used to create the curves in the plot are stored in the axes object of the figure that contains the plot. You can obtain the underlying data by retrieving the axes' xdata and ydata properties. For example:

h_axes = gca;
x_data = get(h_axes, 'xdata');
y_data = get(h_axes, 'ydata');

After extracting the data you can use interpolation (see e.g. interp1())to find the value that you are looking for.

Is this what you meant? If not, please let me know and I'll try to explain better...


What does the phase of a FFT bin of an audio signal represent? I currently run real time audio through an FFT and the library (vDSP) returns magnitude and phase values (polar coordinates) for data in the FFT bins.

Does the bin with the highest magnitude always represent the fundamental frequency?

What exactly can the phase of the frequency bins be used for in regards to audio? What does it represent? Would that be used for audio identification, e.g. matching amplitude and phase pairs within the same bins of different audio segments?

The frequency bin with the highest magnitude is very often not the fundamental pitch frequency, especially for lower (male) voices and larger stringed instruments.  It could instead represent a strong overtone or harmonic.

The FFT phase result is very important in representing any transient waveform shape or timing, for instance whether a note onset or percussive hit occurred toward the beginning or end of the FFT frame.  It's thus not very useful in matching sounds unless the FFT frame timing is somehow very precisely matched to the sound's timing.

For a very narrow-band signals that are stationary across multiple FFT frames with known time offsets, the FFT phase results can also potentially be used for frequency estimation finer than the FFT bin spacing resolution.


Simple question about Sinc() Suppose I use the ideal lowpass filter eq to attempt creating a IIR filter using it
:
$$
h[n] = \frac{\sin(\omega_c \, n)}{\pi \, n}, \qquad -\infty < n < \infty
$$

Obviously, I can't really (unless I'm wrong) realize a filter with this impulse response because I don't have the finite difference equation to make the recursion happen.

But I'm told that it's also an unstable system (textbook solutions).
My question is.. why?

The stability condition states that a system is stable if its impulse response satisfies:
$$
\sum_{n=-\infty}^{\infty} \Big|h[n] \Big| < \infty
$$

From what I see, if I look at the condition blindly, it does satisfy the stability criterion. Doesn't it? Then why it's unstable? Is it because of something like a discontinuity at $n=0$ (even if technically its equal to one)?

Thanks!

No, the signal
$$
h[n] = \frac{\sin(\omega_c \, n)}{\pi \, n}, \qquad -\infty < n < \infty
$$

is not absolutely summable, as the sum

$$
\sum_{n=-\infty}^{\infty} \Big|h[n]\Big| = \sum_{n=-\infty}^{\infty} \left |\frac{\sin(\omega_c \, n)}{\pi \, n} \right|
$$

diverges to infinity. Hence the absolute sum does not converge and therefore it's considered as an unstable system's impulse response. Yet, its square is absolutely summable:

$$
\sum_{n=-\infty}^{\infty} \Big|h[n]\Big|^2 = \sum_{n=-\infty}^{\infty} \left |\frac{\sin(\omega_c \, n)}{\pi \, n} \right|^2
$$
does not diverge. This kind of condition is stated as mean square convergence; i.e., a typical condition for discontinuous frequency responses that has finite jumps in its frequency response which should be continuous for all $\omega$ according to the Fourier theorem that would require uniform convergence at all frequency points $\omega$.


How to prove that Dilution of a time domain signal results in duplication in Frequncy domain? While reading a description on FFT from this book, I found the following statement:

  ...diluting the time domain with zeros corresponds to a duplication of
  the frequency spectrum.

How can I prove this intuitively or otherwise? I have read some properties of Fourier Transform, such as a  circular shift in time domain and frequency domain etc, but I am unable to deduce how duplication in frequency domain will happen. Following illustrate this:

Time domain samples

x(0)=a
x(1)=b
x(2)=c
x(3)=d

Frequency domain samples

X(0)=A
X(1)=B
X(2)=C
X(3)=D

Dilution with 0 in time domain

x(0)=a
x(1)=0
x(2)=b
x(3)=0
x(4)=c
x(5)=0
x(6)=d
x(7)=0

Frequency domain samples having duplicate values

X(0)=A
X(1)=B
X(2)=C
X(3)=D
X(4)=A
X(5)=B
X(6)=C
X(7)=D

Frequency domain samples have duplicate values. Why?

Unfortunately I can't come up with an intuitive explanation for this, but lemme try to prove this. During reading please be very careful with indexes as they can be confusing! I use zero-based indexing. And please correct me if I'm wrong.

Let  be our undiluted time domain of N elements.

So, diluted vector looks like this:

Frequency spectrum of x is given by .

Frequency spectrum of y, teh diluted x is . Note the length of Y equals 2N.

Now I'd like to prove that the right half of the diluted frequency domain is "copied" from the left half, that is . I guess it's sufficent to prove only the last equation, but the first two are good for grasping the idea.

By expanding a sum for Y[0] and Y[N] we get:

The fact that Y[0] and Y[N] are equal can easily be seen by noticing that  because of it's periodic nature- see Euler's relation: . Later we will exploit this property of exp(j*phi) to "move" exponent back or forward by 2nPi.

The next (unnecessary, but good to examine) step is to prove that Y[1] = Y[N+1]. Let's expand:

Because Y[1] = Y[N+1] must hold regardless of x values, the only way that it can do so is when exponents next to proper xn's are equal. Let's prove that . If we can do this, then when a=0, we prove that exponents next to x0 in Y[1] and Y[N+1] are equal. When a=1, exponents near x1 are proven to be equal and so on, until a=N-1.

As you can see, it was not hard. It's more like showing than proving. If you are wondering why did 2pi*a dissapear, look again at the Euler's relation and maybe rewrite the middle step of the proof in cos - isin form. It should be clearer now that you can add or subtract 2pi from argument of the sine or cosine function without affecting it's value.

Okay, let's move on to the next, (the only necessary) part. Y[k]=Y[k+N].

Like the last time, let's show equality of exponentials near xn's:

The idea is the same, only the last point has the least numbers and the most variables, making it the most general case. It's general enough that I think it's the end of the proof.

Hope it helps :)


Can LMS filter converge to sparse solution? Consider received signal $r=h*s$, where $h$ is the channel and $s$ is transmitted signal. Now from $r$ ( $s$ is known) can we predict $h$ using LMS algorithm under all conditions?  I was under impression that both $s$ and $h$ needs to be stationary to estimate $h$ properly. In case, $h$ is sparse (few non zero only) LMS cannot converge to sparse $h$. Is it true or false? Can someone explain my doubt pleas?

...can we predict $h$ using LMS algorithm under all conditions?

LMS converges under most "reasonable" conditions.

The weights of the filter are adapted by a "signal" that is derived by the mean square error of the difference between two signals, which is convex. That is, it has a single minimum (when it does) and it will trickle towards it. If you notice the formula by which each weight is adapted, it does not include a cumulative term. It is just adapting the weights by the magnitude of the error. Therefore, it is supposed to "track" the situation smoothly until it eliminates the difference that is driving it. But if you now crank the $\mu$ up, it will be overshooting that single "minimum" in either side and it will go into "oscillations". Not in the "system is unstable and it is now whistling" sense but in the way that the logistic map goes into "oscillations"

  I was under impression that both s and h needs to be stationary to estimate h properly.

They would have to be stationary for the weight adjustment process to "converge". Otherwise, it will keep "hunting".

  In case, h is sparse (few non zero only) LMS cannot converge to sparse h. Is it true or false?

It is false...most of the times.

Given a reasonable choice for the learning parameter $\mu$ and enough time, LMS will converge to something. But, there are two "problems" that will affect your $h$ estimate: The fact that there is noise and the weights will be trying to minimise even that and the fact that you estimate $h$ from a small window in time. Therefore, the signal(s) better have exhausted the full repertoir of their variance within that window, otherwise, LMS will be "hunting" (smoothly, or wildely) around the minimum again.

Here is an example: $h_{channel}$ is innocent, maybe $h_{channel} = [1,0,0,0,0,0,0,0,0]$. You start with some random $h_{predicted}$ and feed through a sinusoid. The filter will now start adapting. Through the first few iterations, you get totally random "guesses". But there will be one guess, purely because of chance, that a coefficient will land higher than the rest. LMS will pick this up immediately and start improving on it. In the end, it will converge to something that looks like $h_{channel}$, if only the noise perturbations were not there. So, your $h_{predicted}$ is not going to look exactly like $h_{channel}$ with a few strong returns to show reflections and zeros in between, it will still have some values interdispersed with the strong reflections.

Try to think "vertically" into the coefficients of $h_{predicted}$. Because the LMS is like an Automatic Gain Control "Bank". LMS doesn't optimise the "filter" as a whole (i.e. looking at it "horizontally", as one $h$). It uses the error signal to tune each filter coefficient separately and through convolution it applies the effect of each one on the signal. Try to work out LMS, by hand, with a two coefficient $h$ to see what I am trying to say with this.

Hope this helps.

EDIT:

I am adding this small "post script" here, given the additional information that was provided in the comments.

It seems that what is sought after is the impulse response of a "corridor" over the medium of Radio Frequencies.

Depending on the room size, its construction and frequency of operation, the system will not even have the chance to "resolve" between returns.

A nanosecond of propagation at the speed of light is about 30cm of "distance". The distinct pulses that are expected here, with large intervals of "silence" between them, will pile up into "averages", creating a "complex" $h$. That $h$, by the way, is not going to be stationary, depending on the building construction, because you cannot control what the RF is bouncing off of OUTSIDE the limits of the corridor.

You might be better off with a slow chirp to obtain the frequency response of the channel (assuming that this is what you are after).

Hope this helps


Find the dominant low frequency I have this signal, with an arbitrary x-axis, and I want to find the (visually) dominant frequency of the signal. There may be a high frequency component, but I'm interested in the component with a period of around 120 pixels.

First I implemented the Pisarenko Harmonic Decomposition, but only got a very high frequency out (0.04 Hz) (one bin, one second). The I tried with the Fourier transform, and when I plotted the frequency amplitude spectrum I also only got a range of very high frequencies.

Can someone help me solve this problem? What am I doing wrong? It feels like such a simple problem, but for some reason I just can't solve it.

The raw dataset is located here: http://pastebin.com/5ghKbBRm

Have a look at the Python code below. It calculates the auto-correlation of the data, after the DC part of the signal was removed. The autocorrelation shows nice peaks with much less noise than the original signal. Then, we detect the local maxima and measure the distance between them. This distance is then the estimated period of your signal:

import csv
import requests
import numpy as np
from io import BytesIO
from scipy import signal

url = 'http://pastebin.com/raw/5ghKbBRm'
response = requests.get(url)
A = np.genfromtxt(BytesIO(response.content)) # read the CSV file

# begin signal processing

A = A - np.mean(A)  # remove the DC part of the signal
plt.subplot(121)
plt.plot(A)

plt.subplot(122)
C = np.correlate(A, A, 'same')  # calculate autocorrelation to smooth out the noise
plt.plot(C)

plt.subplot(122)
P = signal.argrelextrema(C, np.greater)[0]  # caluclate local maxima (this step can be improved)
P = P[C[P]>15000] # filter maxima, where the autocorrelation values is too low
plt.stem(P, 30000*np.ones_like(P))  # plot where the detected peaks are

print ("Peak distances: ", np.diff(P))  # calculate distance between peaks and the mean of the distances
print ("Estimated period: ", np.mean(np.diff(P)))

Program output:

Peak distances:  [117 118 117 117 117 117 118 117]
Estimated period:  117.25

There is room for improvement. For example, the detection of local maxima is not very robust and the threshold was set manually, but this should enable you to find out the rest on your own.