What filter to use in audio analysis filterbank instead of FFT? Standard bandpass filters can make super precise analysis filterbanks with 1024 to 4096 filters, on reaktor4. I tried in code to used cookbook BandPass and the result was aweful.

Does someone know a precise BandPass Filter that transmits narrow bands of an intended frequency without noise and irregularity? resonance is an advantage because it takes the detected frequency's value and adds a large pure amount of that frequency as a sine wave, which is cool because ideally i would be letting through individual sine waves of given frequency. i want something of that kind. what kind of filter should be used on filterbanks? do you have an example in code?

This has already been addressed in depth here:
https://stackoverflow.com/questions/5901483/simple-audio-filter-bank

and i found some c# code on the subject here:
https://waveletstudio.codeplex.com/

Copied answer from SO here

Using FFT to split an Audio signal into few bands is overkill.

What you need is one or two Linkwitz-Riley filters. These filters split a signal into a high and low frequency part.

A nice property of this filter is, that if you add the low and high frequency parts you get almost the original signal back. There will be a little bit of phase-shift but the ear will not be able to hear this.

If you need more than two bands you can chain the filters. For example if you want to separate the signal at 100 and 2000Hz it would in pseudo-code somewhat like this:

low  = linkwitz-riley-low (100, input-samples)
temp = linkwitz-riley-high (100, input-samples)

mids = linkwitz-riley-low (2000, temp)
highs = linkwitz-riley-high (2000, temp);

and so on..

After splitting the signal you can for example amplifiy the three output bands: low, mids and highs and later add them together to get your processed signal.

The filter sections itself can be implemented using IIR filters. A google search for "Linkwitz-Riley digital IIR" should give lots of good hits.

http://en.wikipedia.org/wiki/Linkwitz-Riley_filter


What is the meaning of the delay in a FIR filter I'm using a C# library that has the FIR filter implemented and I'd like to know what's the meaning of the ${z^{-1}}$ operation (delay) mathematically.

The triangle is multiplication, the circle is a sum, but what does the delay means in terms of mathematical operations?

In mathematical terms, time discrete systems are are most often expressed in terms of difference equations. A delay is simply $y[n]=x[n-1]$ so the difference equation for your FIR filter would be $$y[n]=\sum_{i=0}^{M}b_{i}\cdot x[n-i]$$.

See https://ccrma.stanford.edu/~jos/filters/ for a good primer


What is a correct way to downsample an impulse train So I have a signal which is caused by an impulse train. I can determine the impulse train with 500hz but sample my output signal with only 10hz. Then I want to calculate an impulse response function. So I have a binary "signal" for the impulses with 50 times more samplepoints then my actual signal. If I resample my binary signal using fourier which I would do with an ordinary signal, I get Gibbs ringing and I doubt my impulse response function would still make sense. I can choose the nearest neighbor to be 1 in my "downsampled" impulse train but then the timing is less accurate, I could divide the impulse over the two neirest neighbors bu then no single point will have the value of 1. What would be a correct way to go about this?

The usual way to do this is with a low-pass decimating FIR filter.  If you are going to decimate by a factor of 50 then you will probably want to break it up into multiple LPF decimation filters.  I would probably do it as /10 followed by /5.


Why do we need to normalize a WAV file before calculating the FFT? I am trying to calculate the FFT of a signal stored in a WAV file using SciPy. I found a solution here, but it seems like we need to perform this step before the FFT:

b=[(ele/2**8.)*2-1 for ele in a] # this is 8-bit track, b is now normalized on [-1,1)

where the signal samples are stored in array a. Can someone explain the significance of this step?

Thanks.

It seems that the WAV file used in that example has samples with values between 0 and 255 (likely stored as unsigned chars). The normalization step just changes the samples to floating point values in the range [-1,1). This is, strictly speaking, not necessary to perform the FFT, but it is a good idea. For instance, it removes the large DC component that results from interpreting the unsigned char values as actual signal amplitudes.


What are the fundamental differences between speech and music in frequency domain? What are the fundamental differences between speech and music in the frequency domain?
I can think mainly of features in the time domain that can differentiate speech from music signals, e.g. short silences [Wolfe2002Speech].
Maybe fundamental frequency could be used for distinguishing speech from music? If yes then what would be typical ranges of $F_0$ for speech and music?

One way is to look at the energy below 80 Hz. Unless you have a recording of James Earl Jones, there will be almost none for speech and any music with bass and/or drums is going to have a lot down there.
Note that there is huge variability in music recording and a lot in speech. Creating an algorithm that can reliably distinguish between the two for the majority all recordings is a fairly complicated. Looking at the "average spectrum" isn't that useful since the any individual instance can (and does) look completely different than the "average spectrum". This is actually a good candidate for machine learning although many deterministic algorithms do exist as well.


Weighting function for output error IIR filter design I noticed, when I try to fit filter coefficients to a given complex transfer function with the output error method, implemented for example in the MATLAB function invfreqz.m, I get way better results when I use a weighting function that drastically emphasizes the lower frequencies, for example $1/(2\pi f)$.

Can anybody explain or give me a hint in the right direction, why emphasizing the low frequencies in the weighting function gives so much better fitting results?

My frequency vector is logarithmic, in case that makes a difference and I try to get the best fit for "relatively simple" transfer functions with a very low filter order.

The reason for the worse fit at low frequencies (without additional weighting) is the low density of frequency points at low frequencies due your choice of a logarithmic frequency vector. I'm pretty sure that you will get a better fit at low frequencies if you choose a (more) linear frequency grid.


What does "Instantaneous Angle" means My teacher said to me that Instantaneous Angle means Argument of Complex Function thats all. Someone please tell me what significance does this "instantaneous" Word has here ?

Instantaneous usually means the value of something right at this time. I think what this means in your case is that the function $f(t)$ is complex and you are taking the argument (finding the angle):
$$
\theta(t) = \arg\left( f(t) \right) = {\rm atan}(\Im[f(t)],\Re[f(t)])
$$
where $\Im[ \cdot ]$ indicates taking the imaginary part and $\Re[ \cdot ]$ indicates taking the real part.


Estimate a microphone's transfer function I am trying to compare the transfer functions of k different microphones and need some advice on my current approach. As of now, I have attempted the following:

Play a pink noise test tone $(X)$ from a speaker and record the tone on all $k$ target microphones. Let's call the outputs as $Y_{i}$.
Compute the inverse pink noise tone and convolve it with the recorded microphone tones $Y_{i}$ separately. The output of this convolution operation should be the estimated transfer function ($H_{i}$) of each microphone (it will not be 100% accurate because there is an effect of the room conditions, but let's ignore it for now).

Now to validate that the estimated transfer functions of all microphones are correct, I plan to compute the inverse transfer function $(H^{-1})$, multiply it with each microphone output $Y_i$ in the frequency domain and recover the input signal X. Finally, if all goes well, I expect the recovered X to be the same for all microphones.

$$ Y_i = H_i . X  $$

$$X = H_i^{-1} . Y_i $$

Are there any other approaches to do this? Any pointers for Python implementations of it?

Your $H^{-1}$ is effectively what is called a zero-forcing equalizer for the channel represented by $H$. That channel incorporates DAC, amplifier, speaker, your room, the microphone, amplifier, ADC and several filters necessary to fulfill sampling conditions.

Problem is that for any frequency where $H$ has a notch, $H^{-1}$ has to have a very strong amplification! Sadly, that means that the frequencies where your estimate is based on the least received energy have the largest amplitude in the resulting inverted channel model. That leads to amplified noise.

There's a lot of ways to estimate a channel inverse; the MMSE estimator is pretty popular in communications.

Because you ask for a Python implementation: GNU Radio is usually used from Python, and there's multiple equalizers that come with it (try first!); and, there's modules like gr-adapt which allow for calculating (and using!) channel-reversing filters adaptively. What I especially like about gr-adapt is that its documentation and comments are plentiful and come with references.

For example, have a look at https://github.com/karel/gr-adapt/tree/master/examples :

Here, gr-adapt's LMS Filter is used to estimate (and invert) the channel (which is simulated using the Decimating FIR Filter block).


Time Setting of $z$ and Laplace Transforms I'm aware that the z-transform and the Laplace Transform have an analogous relationship but I want to be doubly-sure that the z-transform only works in discrete-time and that the Laplace transform only works in the continuous time settings?

Thanks.

In typical (or all) applications, the Laplace Transform is used for continuous time systems and the z-Transform for discrete systems. However, the Laplace Transform for Discrete Time Systems certainly exists but would be more complicated than necessary to solve. The z transform exists as a mathematical simplification of the Laplace Transform that can be applied if a consistent sampling rate is used. This may be clearer with the formulas given below:

The (one-sided, for causal systems) Laplace Transform is given as:

$$X(s) = \int_{t=0}^\infty x(t) e^{-st}dt$$

Thus the Laplace Transform for a discrete time system is given simply by setting t= nT for integer n (where the integral then becomes a summation):

$$X(s) = \sum_{n=0}^\infty x(nT)e^{-snT}$$

At this point we can make a substitution to get rid of the more complicated exponential (which is one possible mapping from s to z):

$$z = e^{sT}$$

resulting in the z-Transform for a discrete-time sequence:

$$X(z) =  \sum_{n=0}^\infty x(nT)z^{-n}$$

Equivalently the z-Transform for a discrete sequence (no units of time, simply sample number) is:

$$X(z) =  \sum_{n=0}^\infty x(n)z^{-n}$$

Thus we see that the Laplace Transform of a discrete time system exists (but we have no need to solve that typically as working with the z transform once discrete is much simpler). I am not aware of a reverse case although we may be able to create one by taking the limit $T\rightarrow 0$ of the discrete time z-Transform. I haven't worked through this but we may find the same situation (we can force it to be mathematically equivalent but the math involved would be more complicated than necessary vs taking the Laplace Transform directly).


Manipulating line level signals I need to manage several line level signals (physically, not in software), and need to know what to expect. I was expecting AC signals, but am not quite sure. Are line level audio signals AC or DC? What Vpp can I expect?

Audio signals on line level outputs are AC signals with no DC bias, so they run positive and negative around the 0V level. Many devices support 2VRMS output for full scale signal, so that is about 5.66 Vpp.


Matlab: xcorr: Having trouble when I try to calculate delay between two sinusoids I have a problem when I try to use Matlab built-in function for cross-correlation 'xcorr'.
Background: I created two sinusoids, one of them delayed in time. I choose center frequencies of these signals slightly different than another. In this way, I am trying to analyze the effects of frequency mismatch on cross-correlation results. Problem: The problem is when I try to calculate delay using xcorr function in MatLab, I could not find the exact lag. The result is quite far away from my offset-time. I try to increase my signal length or sampling frequency. But nothing changes.
Do you have any idea? Here is my simple code, below.
clc
clear all
close all

Fs=48e3;
dt=1/Fs;
t=0:dt:(4*Fs-1)*dt;
CarrierFreq=5e3;
Dev=50;
CarrierFreqDev=CarrierFreq +Dev*(2*randi(2,1)-3);

Offset=10000;
OffsetTime=Offset*dt;
NoisePower=-95
DesiredSNR=20;

% Signal Power
ReceivedSignalPower=NoisePower+DesiredSNR;
ReceivedSignalPowerWatt=10^(ReceivedSignalPower/10);
NoisePowerWatt=10^(NoisePower/10);

% Generate Noise @ specified power
Noise1 = wgn(1,length(t),NoisePowerWatt,'linear');
Noise2 = wgn(1, length(t),NoisePowerWatt,'linear');

% Generate Signal1 @ specified power
Signal1=cos(2*pi*CarrierFreq*t);
pSignal1 = sum(Signal1.^2) / length(Signal1);
Signal1 =  Signal1/sqrt(pSignal1);
Signal1_Norm=sqrt(ReceivedSignalPowerWatt)*Signal1;
Signal1=Signal1_Norm+Noise1;

% Generate Delayed Signal @ specified power
Signal2=cos(2*pi*CarrierFreqDev*t);
pSignal2=sum(Signal2.^2) / length(Signal2);
Signal2_Norm=sqrt(ReceivedSignalPowerWatt)*Signal2;
Signal_Delay=circshift(Signal2_Norm, -Offset);
Signal2=Signal_Delay+Noise2;

[Coeff, Delay_Time] = xcorr(Signal1,Signal2); %% with dev.
[ ~,Indx]  = max(abs(Coeff));
Delay_Time = Delay_Time(Indx)*dt;

A sine wave is cyclo-stationary so is not the best waveform to use for estimating delay between transmit and received signals. It will correlate whenever the delay causes the signal to be an additional $2\pi$ in phase offset. Consider using waveforms with better autocorrelation properties for this purpose such as pseudo-random noise sequences or other generated white noise sources.
Regardless, the effect of frequency offset on cross-correlation is that the magnitude rolls-off as a Sinc function with frequency offset with the first null at $1/T$ where $T$ is the total time duration of the signal, and the phase decreases linearly from $0$ at $f=0$ to $-\pi$ at $f=1/T$. At a fixed frequency offset between the signals, a fixed phase offset would result, and for a narrow band signal this can be approximated as a fixed delay according to the carrier frequency of the signal.
Consider that correlation is the sum of the element by element product between two signals (the dot product), and observe how the product of two signals will be at the frequency difference between them (the beat note). Thus if your total time duration is exactly 1 complete cycle off in frequency, the result will average to zero (showing how the first null is at $1/T$).
The following code demonstrates this experimentally with a random white noise waveform of total time duration 0.1 seconds:
T = .1;
fs = 10e3;
N = T*fs;
t = [0:N-1]/fs;
x = randn(1,N);

index = 0;
foffset = 0:.01:20;
for f = foffset
    index++;
    y = x.*exp(j*2*pi*f.*t);
    # dot product is the correlation, and normalized
    corr(index) = x*y'/(N*std(x).*std(y));
end

figure
plot(foffset,abs(corr))
figure
plot(corr)

This creates the two plots shown below with the first plot of the magnitude of the normalized correlation versus the static frequency offset between the two signals, and the second plot as a polar plot showing both magnitude and phase components.


could anybody explains me why IIR filters do not have linear phase? Please tell me. Also I dont know why phase is linear with FIR filters. I would like quantitative analysis. And why linear phase is not achieved by IIR filters ?

For digital filters, linear phase places the following requirement on the transfer function:

$$ H(z) = H(z^{-1}).$$

That restriction implies a linear phase IIR filter would need to have poles both inside and outside the unit circle, making it unstable. Similar arguments apply for analog filters.

That being said, there are any number of approximations which may be "close enough" to linear phase, depending on the application--especially if the causality of filter is sacrificed. For a review of techniques, see the introduction to this paper:

  S.R. Powell, P.M.Chau, A Technique for Realizing Linear-Phase IIR
  Filters, IEEE Trans. Signal Processing, Vol 39, No 11, Nov 1991, pp
  2425-2435.

The algorithm in that paper achieves linear phase with acausal block processing, rather than the usual offline "Forward-Backward" zero-phase approach.


How to make pop-free, crisp two-beep ringtone? I hope this ultra-noob question is not entirely out of place here.

I'm trying to make a super simple ringtone: 0.5s silence + 1s "pure tone" + 0.5s silence + 1s "pure tone" + 0.5s silence.  By "pure tone" I mean something like a sinusoidal or square wave with a fixed frequency (let's say 350Hz for concreteness).  This profile results in a very standard beep-beep ringtone.

With Audacity, it's relatively simple to generate this profile (with Generate > Silence... and Generate > Tone...), and one even has a choice of wave shapes (sinusoidal, square, etc.).

If I pick the sinusoidal option when constructing the profile above, and export it as a WAV file, when I play the ringtone in my phone, I can hear a slight "pop" at the beginning of each beep.

I've tried to get rid of those pops, with little success.  I figure that the pops are an artifact of the discontinuity of the wave at the beginning of each beep (though, curiously, I don't hear a pop at the end of each beep).  So I tried fading the sound in and out, which does get rid of the pop, but, as one would expect, makes the beeps sound less crisp.

I've heard many ringtones of this sort that are crisp, pleasant-sounding, and pop-free, so there must be a way to get rid of these pops.  Any suggestions would be appreciated!

P.S.  If I pick the square wave + no aliasing option (which, to my ear sounds a little "cleaner" than the option with aliasing, for some reason), the resulting ringtone does not have the "pop" at the beginning of each beep, but it sounds very harsh and loud.  I tried to adjust the loudness of the beep by lowering the amplitude of the wave generated by Audacity, and in fact, when Audacity plays this reduced-amplitude wave it indeed sounds quieter, but when it is saved as a ringtone in my phone, the phone plays both versions of the square-wave ringtone equally loudly.  I figure that the phone probably normalizes the ringtone's loudness somehow.  What's surprising to me is that the ringtone made with the sinusoidal wave (the one with the slight "pops" at the beginning of each beep) sounds considerably less loud than the ringtones made with the square wave.  In any case, even if I could find a way to reduce the volume of the square-wave ringtone that the phone did not "normalize away", I still would prefer a pop-free version of the sinusoidal-wave ringtone.

The "click" is caused by the discontinuity in the waveform and its derivatives - even if the waveform stops at a zero-crossing you might still hear a pop if there's a discontinuity in the higher order derivatives! Its loudness depends on the amplitude at the discontinuity - and it is thus influenced by the length/frequency/phase of the sinusoidal tone. Note also that if the tone is contained within a clip in an audio editor, a short fade-in/fade-out might be added automatically by it, attenuating the problem. I know a few audio editors that dynamically apply a short fade (say 32 samples) at the boundary of each audio clip during mixing during playback.
You hear the click clearly with the sine wave because the tone has very few harmonics (ideally: none). You might not hear it with the other waveforms (square) because they already have high frequency content that will cover part of it (frequency masking phenomenon). For example, you won't hear a click when a burst of white noise stops.
The solution is to apply a short fade-in / fade-out. 1ms is a good start. Above 5ms it'll start feeling "slow". There is a trade-off to find between "clicky and snappy" and "slow" - this is a problem even well-established synth manufacturers have to educate the public about. With a sine-wave, an envelope curve like a raised cosine function ensures that the resulting signal has continuous derivatives throughout the fade.
Aliasing is the bad thing to get rid off, so it is normal that the "no-aliasing" option produces a better quality output.
A sine wave with an amplitude of 1V has twice as less energy as a square wave with the same amplitude. Assuming that the phone normalizes amplitude, it is normal that with a given dynamic range/amplitude budget, the sine wave will sound quieter. Energy/RMS is a rather bad proxy for the perceptual concept of loudness, but it's much more accurate at this task than peak-to-peak amplitude!


What is a frequency bin? I'm doing a research on the FFT method, and a term that always comes up is "frequency bin". From what I understand, this has something to do with the band created around the frequency of a given sinusoid, but I can't really figure out how. I also figured out how to go from a given bin to its related frequency, but still no intuition on what a frequency bin is.

It's simpler than you think. When we discretize frequencies, we get frequency bins. So, when you discretize your Fourier Transform: $$e^{-j\omega} \rightarrow e^{-j{2\pi k}/{N}}$$ Our continuous frequencies become $N$ discrete bins.

This is exactly why the following is true: $$n^{th}bin = n*\dfrac{sampleFreq}{num(DFT points)}$$ Note that the FFT represents frequencies 0 to sampleFreq Hz.


Controller in the feedback path? What happens if the controller of a closed-loop system is placed in the feedback path before the summation point, instead of between the summation point and the plant to control? Are we allowed to do so? If so, what are the main differences?
Thanks!

Assuming $P$ is the transfer function of your process and $C$ the transfer function of your controller.
The closed-loop transfer function of a standard control loop, with the controller in the normal path, is
$$ G(s) = \frac{C(s)P(s)}{1+C(s)P(s)}$$
while the transfer function of your alternative controller is
$$G_{alt}(s) = \frac{P(s)}{1+C_{alt}(s)P(s)}
$$
If you set your new alternative to controller to
$$
C_{alt}(s) = \frac{1+P(s)C(s)-C(s)}{P(s)C(s)}
$$
You will get the same closed-loop transfer function (G(s))
That being said, not sure it's a good idea for 2 reasons
1 - To get to same result as the original controller, you need to cancel poles and zeros in the alternative closed-looped transfer function. If you cancel a RHP zero with a RHP pole, your final transfer function will seem stable however your controller will not be, the output of your alternative controller could become unbounded.
2 - Even if you don't cancel RHP zeros with RHP poles, this method is way more complicated.


Mathematical expression of Multicarrier frequency equalizer I have a signal $X$ with length of $N$, multiplied with any unitary matrix, i.e  the transpose of DCT matrix as:
$x = D' X$
where $D'$ is the transpose of the of DCT matrix. Then let's add $CP$ guard interval into the signal $x$. The resultant signal after adding the cp into the head of the signal is $x_2$.
We transmit $x_2$ through a channel $h$ leading to have the signal $y$ as follows:
$y = h*x_2$   ......  where $*$ is the convolution operation.
Hint: the convolution is not circular in that case as we used DCT instead of DFT.
At the receiver, after removing the delay of the channel and the CP guard interval, I used the following steps:
$x_3 = ifft(fft(y)./H)$ .....  I used ZF equalizer, $H$ is the frequency-domain channel found by $H = fft(h,N)$
Then, the signal is recovered  by multiplying by the DCT matrix
$x_4 = D * x_3$
My question,  how does it work the step when I get $x_3$?? is there mathematical expression for that,  , I think it's called frequency domain equalizer, but I didn't get its mathematical expression. Most importantly, why does it work however we don't have circular convolution when getting $y$? !!
thank you

Hint: the convolution is not circular in that case as we used DCT instead of DFT.

Luckily, that's not true! At the receiver, you're doing a DFT, and you've prepared the transmit signal such that a convolution with the channel looks like a circular one by prepending the CP.
The fact that the data inside has something to do with the DCT is irrelevant there.

Then, the signal is recovered by multiplying by the DCT matrix

Yep, because that inverts $D'$.

I  think it's called frequency domain equalizer,

No, the frequency-domain equalizer is $\text{IFFT}(\text{FFT}/H)$.


How to do addition of sound pressure RMS in a DFT Based on this answer. Once you have your DFT in the unit sound pressure RMS, what is the correct method for getting the sound pressure RMS sum over multiple bins (if for example you're trying to get the total RMS over a frequency band)?

The power spectrum of the DFT is the square of the magnitudes and is related to the power spectrum of the discrete time signal by Parsevel's theorem:

$$E_x=\sum_{k=0}^{N-1}\lvert x[k]\rvert ^2=\frac{1}{N}\sum_{k}^{N-1}\lvert X[k]\rvert^2.$$

The RMS is just the square root of the power. (I assume you've already taken care of the 1/N since you mentioned your DFT is already in the correct units).

$$\textrm{RMS}_{x}=\sqrt{\sum_{k=0}^{N-1}\lvert X[k]\rvert^2}.$$

You can compute the RMS of a range of bins $[i,j]$ by constraining $k$.

$$\textrm{RMS}_{x_{ij}}=\sqrt{\sum_{k=i}^{j}\lvert X[k]\rvert^2}.$$


Is there any way to get intensity value of a certain wavelength for given pixel of a normal  image? I understand I can get a color from wavelength but it's not true for the opposite direction.
In case this picture would help.

To make this question simple, let's lighting and calibration aside to an optimized condition.

I wonder if there is any way to get intensity value of a certain spectral wavelength for given pixel of a common image. That is, Wavelength and rgb of a pixel are both known. Based on that, I want to compute intensity of the wavelength for the pixel.

If I understand correctly - you want to know what wavelength a given RGB value will show?

Only some RGB values represent a pure wavelength, those on the edge of CIE diagram (there is no wavelength of brown for example) see http://www.fourmilab.ch/documents/specrend/ for details and code


Period of a continuous signal So I have
$$x_1=2 \cos(.6\sqrt\pi x+\pi/6)$$ and $$x_2= \sin(1.2\sqrt\pi x-\pi/3)$$ and need to find the period of $(x_1+x_2)^2$.

Let $$a=.6\sqrt\pi x+\pi/6 ~~~\text{and} ~~~ b=1.2\sqrt\pi x-\pi/3$$ so that
$$x_1=2 \cos(a) ~~ \text{and} ~~ x_2=\sin(b).$$

And finally we get:
$$(x_1+x_2)^2=1.5+2\cos(2a)+2\sin(a+b)+2\sin(a-b)-.5\cos(2b)$$

Now, the how do I calculate the time period of it?

Edit:

$$\text{period of the individul sinusoids in the square}~~T_1=5\sqrt\pi/3 ~~ T_2=10\sqrt\pi/9~~T_3=10\sqrt\pi/3 ~~ T_4=5\sqrt\pi/6 $$

$$\text{The ratios are}~~T_4/T_3=1/4~~T_4/T_2=3/4~~T_4/T_1=1/2$$

$$\text{LCM of denominators is 4 } $$

$$\text{Hence period }~~T=4T_4=10\sqrt\pi/3$$
is that correct?

Let $T_0$ be the period of the signal  $s(t) = x_1(t) + x_2(t)$.
$$s(t) = s(t+T_0) \tag{1}$$

Then the period for $s^2(t) = ( x_1 + x_2)^2$ will be either $T_0$ or $T_0/2$ :

First, it's obvious that $T_0$ is also a period for $s(t)^2$ as it satisfies $$s^2(t) = s^2(t+T_0) \tag{2}$$

However, $T_0 / 2$ may also be period for $s^2(t)$. Since we are looking for the smallest number to satisfy eq(2).

Hence do the following:

Find the period $T_0$ of $s(t) = x_1(t)+x_2(t)$.
Check if $T_0/2$ satisfies $s^2(t) = s^2(t+T_0)$.
If yes; then period of $s^2(t)$ is $T_0/2$; otherwise $T_0$.


Alias frequency Formula I'm taking a multimedia systems class in my MSc Computer Science, and I'm having some trouble understanding the formula for the alias frequency - this could stem from my misunderstanding of the alias signal.

My understanding of an alias signal is that if you undersample your input signal (i.e. sample at a rate which is less than twice the maximum frequency) then we can get aliasing because we're not sampling frequently enough to capture the high frequency details.  The aliasing signal is the result of taking these sample values and joining them with a smooth curve.

Therefore, the resulting signal has a frequency of half the sampling frequency, since a pure sinusoid will need two samples per oscillation (1 for each turning point) - this would mean that the alias frequency should just be a function of the sampling frequency.

The formula for the alias frequency is the absolute difference of the signal frequency and the closest integer multiple of the sampling frequency - can someone explain this to me?  Thanks in advance!

Suppose that the sampling is done at a rate of $1000$ Hz, one sample every millisecond. Suppose also that the signal being sampled is at $3200$ Hz,
the first sample is at the peak of the sinusoid. The next sample,
will be taken one millisecond later during which time the sinusoid will
have gone through $3.2$ periods, and so the next sample will have the
same value as if the sinusoid had gone through $0.2$ periods, not $3.2$
periods. The one after that will be $0.4$ periods away from the peak, and
so on.  This is exactly the same set of samples that we would gave gotten
if we had been sampling a $200$ Hz sinusoid.  In one millisecond it would have
progressed through $0.2$ of its period of $5$ milliseconds and so on.
In other word, just by looking at the samples alone we cannot tell whether
the samples came from a $3200$ HZ signal or from a $200$ Hz signal.

If the signal being sampled was at $2800$ Hz, then we would get samples
corresponding to $0$, $-0.2$ of the period, $-0.4$ of the period and
so on.  But because sinusoids look the same in either direction in time,
these samples also look like they are the result of sampling a $200$
Hz signal.  This is the reason why the formula you are given, viz.

  Aliased frequency is the absolute difference between the actual signal frequency
  and the nearest integer multiple of the sampling frequency.

works to give you the right answer.


Approximating $ {L}_{0} $ Norm Minimization with Non Linear Convex Inequality Constraints using Reweighted $ {L}_{1} $ Minimization I have an optimization problem consisting of the $ {\ell}_{0} $ norm as the objective and a nonlinear (convex) constraint as well as a linear constraint. I am wondering if the reweighted $ {\ell}_{1} $ norm minimization algorithm would be applicable to this problem. The original paper "Enhancing Sparsity by Reweighted $ {\ell}_{1} $ Minimization" assumes linear constraints. References to relevant papers would be appreciated.

Edit:
This is the problem:

$$\begin{align}
\arg \min_{x} \quad & {\left\| x \right\|}_{0} \\
\text{subject to} \quad & f \left( x \right) \leq d \\
& 0 \leq {x}_{i} \leq p
\end{align}$$

where $ x \in \mathbb{R}^{n} $, $ d, p \in \mathbb{R} $ are constants, and $ f \left( x \right) $ is a convex function. The notation $ {\left\| \cdot \right\|}_{0} $ refers to the $ {\ell}_{0} $ norm, which is the number of nonzero entries the vector for which it is being computed.

I wish to relax the objective using reweighted $ {\ell}_{1} $ minimization, but am wondering if the constraints must be linear for this heuristic to apply.

Also this is not compressive sensing. I just wish to make use of a technique in compressive sensing.

I think testing it practically rather that its theoretical inspection might help more. Therefore, I tried to solve the above mentioned problem using CVX toolbox in Matlab.

cvx_begin;
   variable x(1000);
   minimise(norm(x,1));
   subject to
      x.^2<=150; % f(x) here, is x^2
      0<=x<10;
cvx_end

I think the answer to your question depends on your function f(x), and sometimes even if there is a solution based on convexification of the stated problem, I do not see any guarantees that the answer is must be right. I mean using l1 norm as proxy for l0 norm (as what CS suggests) in this case, might be wrong and bring no fruit. Be careful with blind application of l1 norm minimization method. This great lecture from Prof. Baranuik called "Compressive nonSensing" might be help you more on this regard:
You might find it here:
http://www.norbertwiener.umd.edu/FFT/2015/


white noise filtering If I have a white noise fed into a filter what would be the output of the filter - what would you expect to see?
I know that white noise has a unit power spectral density (PSD), how these thing relate to each other?

For a linear time-invariant filter, the output of the filter would be coloured noise.

Again assuming linear time-invariant filter, the effect of the filter on the output can be further characterized in terms of its power spectrum density (PSD).
Specifically, for an input $X(f)$ with power spectrum density $S_x(f)$ to a linear time-invariant filter with frequency response $H(f)$, the power-spectrum density $S_y(f)$ of the output of the filter is given by:

$$
  S_y(f) = |H(f)|^2 S_x(f)
$$

So, for a white noise input with unit power spectrum density $S_x(f)=1$, the corresponding PSD of the output would be
$$
  S_y(f) = |H(f)|^2
$$


Determining instruments in audio data What exactly is one looking for in audio data to find different instruments or energy sources?

So the FFT can extract frequencies, but how can e.g. percussion be detected alongside e.g. a bass line?

Would the frequency bins be checked into a frequency mapping of instruments? This would lead to collisions in possible instruments.

Are the frequencies of a FFT output directly mappable to pitch frequency charts or should there be a conversion first?

There seem to be four questions here. I can't give complete answers for each of them but let me address them in turn.

  What exactly is one looking for in audio data to find different instruments or energy sources?

There's not really a single answer to this. There's lots of research into many closely related but ultimately different problems. In general instruments vary in the way their spectral properties vary over time. Tracking the changes of all the harmonic and inharmonic components of a note amongst other sounds is no small task.

  So the FFT can extract frequencies, but how can e.g. percussion be detected alongside e.g. a bass line?

Generally you'd use a Short-time Fourier transform so that there's time data as well as the spectral information. From here you've got a few problems to think about:

Source separation: look up techniques such as ICA, NMF and CASA.This will divide up the energy between the detected sources allowing you to detect the separate instruments.

Fundamental frequency detection: I recommend summary auto-correlation of frequency bands rather than purely time or frequency domain methods.

also Source identification (which I don't know much about).

  Would the frequency bins be checked into a frequency mapping of instruments? This would lead to collisions in possible instruments.

You need to apportion the frequency information not just assign it all to one instrument.

  Are the frequencies of a FFT output directly mappable to pitch frequency charts or should there be a conversion first?

The fundamental frequency, which is given on a pitch-frequency chart, may not actually present in a note ('the missing fundamental problem'). So a fundamental frequency detection algorithm such as the one I suggest above will be needed before you can use the chart.


What does boundary discontinuity in DFT imply? I have been studying about boundary discontinuity in DFT and why it is not used as transform in image compression. What i understand is that DFT is N periodic and that causes a discontinuity at the boundary. This means that more number of terms are needed to represent the function and the rate of convergence is slow.

Can somebody explain as to what does it mean when they say that "more number of terms are needed to represent the function?" . And how does it help to use DCT which is 2N periodic?

Reference article is at : http://en.wikipedia.org/wiki/Discrete_cosine_transform

Thanks!

Because the left (top) edge of an image is unlikely to be a reflection of its right (bottom) edge, there are discontinuities all along the edges of an image when it is viewed as N-periodic. These discontinuities are represented in the frequency domain by high-frequency coefficients. From an image compression point of view, using the 2-D DFT as a transform would require to store these high frequency coefficients (or non-zero terms in the decomposition of the image as a sum of elementary components), which do not describe anything meaningful in the image itself, but are a nasty side-effect of an inappropriate choice of transform.


FFT of a Square wave, what are frequency of each Harmonic here After getting suggestion to first study FFT on square wave to understand FFT of discrete signal, I'm trying to understand the FFT of square wave.

I've coded a program, here is the details,
Frequency of Wave, 100KHz

f= 100e3;

Sampling Frequency, 1MHz

Fs = 1e6;
Ts = 1/Fs;   %Sampling Rate
% Ts = 0.000001

t = 1/f;       %Time period of 1 Oscillation = 1/f

% t = 0.00001
n = 0.000001:Ts:t;  %Generating Samples

% t/Ts = 10

i.e 10 Samples will be generated

% 1st  Sample at 0.000001s
% 2nd  Sample at 0.000002s
% 3rd  Sample at 0.000003s
% 4th  Sample at 0.000004s
% 5th  Sample at 0.000005s
% 6th  Sample at 0.000006s
% 7th  Sample at 0.000007s
% 8th  Sample at 0.000008s
% 9th  Sample at 0.000009s
% 10th Sample at 0.00001s

defining and plotting the square wave,
Since Amplitude is not defined it will be 1 by default.

x=square(2*pi*f*n);
subplot(2,1 ,1)
stem(n,x)

Output wave,

Now plotting its FFT

subplot(2,1 ,2)
F=fft(x)
plot(2*abs(F)./(t/Ts));

Output wave,

For the First point,

Fs / # of Samples

1000,000/10 = 100KHz <- First point = 100KHz

Number of  Harmonics should be = 10 but What are the Frequency of each Harmonic ?
also, why the amplitude has changed to 1.2 which is supposed to be 1

From the Wikipedia article on the Discrete Fourier Transform:

  The sequence of $N$ complex numbers $x_0, ..., x_{N−1}$ is transformed into an $N$-periodic sequence of complex numbers according to the DFT formula:
  $$ X_k=\sum_{n=0}^{N-1} x_n e^{-2\pi ikn/N}.$$

What you have done is taken $N=10$ integers, $n=\{1,2,...,N\}$, and turned them into $10$ time samples $t=\{T_s,2T_s,...,NT_s\}$, or $t_n=nT_s$. (To clarify my notation change, this second set is what you called $n$.)

Then you created the signal $x_n=\text{square}(2\pi ft_n)$ with period $f=(NT_s)^{-1}$and fed it into MATLAB's FFT algorithm, attempting to take the DFT of it (as the FFT is just a fast DFT). However, this is NOT what the FFT expects to see!

Your signal starts with $x_1$ and ends with $x_N$. However, if you see the definition of the DFT I gave above, it expects the signal to start with $x_0$ and end with $x_{N-1}$. Ordinarily, this would not be much of a problem, because if your signal $x$ is actually $N$-periodic, then $x_0=x_N$.
However, as noted in the comments by Jim Clay and Jason R above, the signal you start with is not actually a square wave. As you can see in your screenshot, there are six "1" values and only four "-1" values. The square wave should have an equal number of "1"s and "-1"s. I do not know why the values you put in are not a proper square wave, and I suspect there is some odd detail in how MATLAB implimented the function $\text{square}$. To create a square wave, you should change the line

n = 0.000001:Ts:t;  %Generating Samples

to

n = 0:Ts:t-Ts;  %Generating Samples

or, even better, to

N=t/Ts;
n=(0:N-1)*Ts;

which makes it clear that you are sampling at integer multiples of your sampling time. The signal $x$ you generate in this way will be equivalent to what Jim Clay has generated in his answer.

As to why your signal has magnitude $\approx 1.2$ instead of $1$, you need to remember how the square wave is defined. From the Wikipedia article on the square wave:

  $$x_{\mathrm{square}}(t) =\frac{4}{\pi}\left (\sin(2\pi ft) + {1\over3}\sin(6\pi ft) + {1\over5}\sin(10\pi ft) + \cdots\right ).$$

The first term of this function has frequency $f$ and magnitude $\frac{4}{\pi}\approx 1.27$. If you look at Jim Clay's plot this is exactly the magnitude in bin 2 of the function he has plotted.
Up to bin $N/2+1$, the value that will be plotted in bin $k$ is the coefficient of the term in the square wave with frequency $(k-1)f$. (The $-1$ comes from MATLAB indexing beginning with one rather than zero).


How to determine the power of a finite signal? I am trying to write mathematically how to calculate the power of a finite signal. I know the general formula $P_{\mathrm{x}}=\lim_{T\to\infty} \frac{1}{T} \int_{t=0}^{T} x^{2}(t) \mathrm{d} t$ used to calculate the power of a signal which has infinite Energy and finite power. But how do you express it when the does not go to infinity, is not period and does not have infinite energy?

I know its a basic but i have searched all over and have not been able to find a proper answer.

Thanks in advance.

First of all, note that your formula for signal power is only valid for real-valued $x(t)$ that satisfy $x(t)=0$ for $t<0$.

A more general formula is

$$P_x=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T}|x(t)|^2dt\tag{1}$$

For signals with finite energy

$$E_x=\int_{-\infty}^{\infty}|x(t)|^2dt\lt\infty\tag{2}$$

the power $P_x$ defined by $(1)$ is clearly zero.


State space physical meaning Let $T(s)$ be a transfer function that describes a mechanical system, where the input is force and the output is position.
And let $[A,B,C,D]$ be the equivalent state-space representation of $T(s)$, where:
$$
    \dot{x}(t) = Ax(t) + Bu(t) \\
    y(t) = Cx(t) + Du(t)
$$
and let $[A_d,B_d,C_d,D_d]$ be the discretized model of $T(s)$, using zero-order hold on the inputs:
$$
    x_d(k+1) = A_dx_d(k) + B_du(k) \\
    y(k) = C_dx_d(k) + D_du(k)
$$
I know there are infinite possible realizations of $[A,B,C,D]$ that represent the same transfer function $T(s)$.

How can I force the state vector, $x(t)$, to have a simple physical meaning? in this example, since it is a mechanical system, the elements of the state vector $x(t)$ would probably be position, velocity and acceleration. I can guess the matrix $A$ must have this form:
$$
A = \left[\begin{array}{ccc} * & * & * \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{array}\right]
$$
to express the derivative relationship between the state vector elements.
Same question, but for the discrete time case: How can I force the state vector, $x_d(k)$, to have a simple physical meaning? In this case, I'm not sure how "velocity" and "acceleration" can be expressed.

As you pointed out, there are many state-space realizations of one particular transfer function. The reason is that a transfer function only represents the input-output behavior of a system (observable and controllable dynamics) and not the internal states.
That being said, you can directly write state-space realizations from a transfer function with the so-called vertical and horizontal companion forms (equivalently, the observable and controllable canonical forms). The problem is that the physical meaning of the states is not usually guaranteed.
So unless you have a model of your system in terms of differential equations (first principle modeling like Newton's laws or Lagrangian mechanics for mechanical systems), the physical meaning of the states will usually not be accessible. The same holds in the discrete-time domain.


Yes, this is always true if the input to the DFT is real valued. It's called the "conjugate complex symmetry", because
$$
X_{N-n} = {X_n}^*
$$
where $X_n$ is the DFT output and $()^*$ denotes the conjugate. It can be proven by inserting the property into the transformation formula of time domain sequence $x_k$:
$$
X_n = \sum_{k=0}^{N-1}x_ke^{-j\frac{2\pi k n}{N}}\\
X_{N-n} = \sum_{k=0}^{N-1}x_ke^{-j\frac{2\pi k (N-n)}{N}}\\
=\sum_{k=0}^{N-1}x_k e^{-j 2\pi k}e^{j\frac{2\pi k n}{N}}
$$
Using $exp(-j2\pi k) = 1 \:\: \forall \: k$ we find
$$
X_{N-n} = \sum_{k=0}^{N-1}x_k e^{j\frac{2\pi k n}{N}}
$$
Now we exploit that $x_k$ is real for all $k$ to rewrite the above:
$$
X_{N-n} = \sum_{k=0}^{N-1}\left(x_k e^{-j\frac{2\pi k n}{N}}\right)^* = X_n^*
$$
It's now straightforward to show that the same is valid for the IDFT. (real spectrum -> complex conjugate symmetric time domain sequence). Moreover, the opposite statement is also true: if the input to the DFT (IDFT) is complex conjugate symmetric its output is real valued.


Showing specific values in an STFT plot I have a short term fourier transform plot that I plot using matplotlib's pcolormesh() function:

If I zoom in the plot, I see this:

Now, I want the plot to just show the zoomed-in range on the y-axis - till 4 or 3 kHz. Is there any way to do this without having to manually zoom in? I want to do a comparison for multiple waveforms in the frequency range as this is my passband in the passband filter I apply on the data.
Here is the code I use to compute and plot the stft:
sampleFreq, segmentTimes, stftX = sp.stft(filteredY, fs)
segmentTimes += (1.5 / 1000)
plt.pcolormesh(segmentTimes, sampleFreq, np.abs(stftX), shading='gouraud')
plt.ylabel("Frequency (Hz)")
plt.xlabel("Time (s)")
plt.title("Sort Term FFT")
plt.show()

I add 1.5 to the segmentTimes because I trim my original data from 1.5 to 2.5 ms, and for some reason, if I don't add 1.5 here, the stft starts from 0.0 s. And the division by 1000 is for conversion to seconds from milliseconds.
Also, sampleFreq and segmentTimes are simple lists but stftX is a list of lists.

Use matplotlib.pyplot's axis() function as follows:
# show section based on bandpass filter cut-off frequencies
plt.title("Zoomed Short Term Fourier Transform")
plt.xlabel("Time (s)")
plt.ylabel("Frequency Magnitude")
plt.pcolormesh(self.segment_times, self.sample_freq, self.abs_stft, shading='gouraud')
plt.axis([self.segment_times[0], self.segment_times[-1], high_pass, low_pass])

plt.show()

where high_pass and low_pass are the frequency cut-offs within which you want to display the y-axis data.


How to process video with DFT/FFT - getting the frequency of green color changes? Goal: I have a video stream from a webcam. I need get frequency of the green color changes on the video.

Video color space is RGB
The signal looks sinusoidal with a frequency range 0.5 - 2.4 Hz

I performed the following steps:

Capture Video frames into collection.
Processing frames by Gaussian pyramid down level 4. (frame size now is ~ 37x32)
use DFT/FFT with sample rate 30
get spectrum
calculate frequency

on DFT I put data from each frame:
video frames quantity = 60

Method 1:
I collect green channel from frame and join with the green channel from a next frame, in result I have a big data array.

37*32=1184

1184*60=71040 the size of the sample buffer

one long array with 71040 data point of all pixels in the green channel of each frame.
(all green pixels from one frame + from next frame)

Method 2:
I collect the green channel from an each pixel on the frame and join with the same green channel pixel from a next frame, as a result I have a big data array.

1184 arrays with 60 data points
(each green pixel from one frame + from next frame)

Which method is the correct one?

If you need to track changes in time, then the most natural and obvious way is to form one temporal (green data) vector per pixel that you want to observe, then take its DFT (it's way 2, but with more than 2 frames).

The actual vector length and DFT function that you should use will depend on the actual case you are working on.

Suppose you work from a recorder video. Then, you can take the 1D DFT of the vector formed with the value of one pixel of interest in the whole sequence.
On the other hand, if you work from a live feed, then you need to use a FIFO data vector and probably apply a windowed 1D Fourier transform to avoid artifacts, see short-time Fourier transform.


$Q$ factor refers to what on a lowpass filter? (EE kind of definition) For what I'm learning, the bandwidth (of a BPF, for example) refers to the $\Delta$ between $f_L$ and $f_H$ at (usually) $-3\textrm{ db}$ from the $f_0$ (from wiki). For a peaking filter instead, $f_L$ and $f_H$ are situated (usually) at the peak $\textrm{dB Gain}/2$.

Q Factor become $f_0 /\Delta$. It's clear till here.

Now, if I take a lowpass filter, how is calculated $Q$ Factor?
I have $f_0$, but how would I calculate the $\Delta$ of a lowpass?

There is no a "peak" as reference, so $-3\textrm{ dB}$ (or $\textrm{dB Gain/2}$) make no sense as references points for the "created" bandwidth.

For lowpass and highpass filters, Q factor is not well defined. What is defined is the cutoff frequency and bandwidth.

Q factor is useful for bandpass or bandstop (notch) filters as it shows the selectivity of the BPS. Quantitatively Qf is the ratio of the center frequency of the pass-band to the bandwidth of the pass-band. Assuming the system has only a single pass band.

For a highpass filter, since the bandwidth is infinite, this Q will be zero.
And for a lowpass filter it will always be 1/2.
As a result, Q is not used for HighPass or Lowpass analog filters.


FIR-Decimation and Low-pass filter (taps vs number of input points vs number of decimation stages) I am new to digital filters and trying to understand parameters behind them.

Scenario

I am trying to get DC value from signals sampled at 100kHz sampling rate. such high frequency is used due to system hardware limitations. A very sharp low-pass filter is to be used with cut-off frequency at 1Hz. The plan is to capture 10000 ponints first and then FIR-Decimate the income signal with several stages:

stage 1: decimation factor 10, FIR cutoff frequency 1/20 (128 taps), 100kHz->10kHz

income 10000 points, result 1000 points

stage 2: decimation factor 10, FIR cutoff frequency 1/20 (128 taps), 10kHz->1kHz

income 1000 points, result 100 points

stage 3: decimation factor 10, FIR cutoff frequency 1/20 (128 taps), 1kHz->100Hz

income 100 points, result 10 points

stage 4: decimation factor 10, FIR cutoff frequency 1/20 (128 taps), 100Hz->10Hz

income 10 points, result 1 point

Then a FIR with cut off frequency will be implemented with cutoff frequency 1/10.

Questions

In the FIR decimation, the FIR cut-off frequency is always stricter than the decimation factor to avoid aliasing, is this design necessary?
What are the tradeoffs of choosing number of decimation stages and decimation factors?
At stage 4, there are only 10 income points. However, FIR has 64 taps... Would the one point result still be valid?

I'll answer your questions first:

You should look at the final decimation factor when designing the earlier filters. Hence, for the first stage you can have a cutoff frequency of 1Hz and stopband frequency of 100000/10 - 1 = 9999 Hz, since this is what is folded/aliased into your band of interest. You would actually prefer to have several stopbands 9999-10001, 19999-20001, 29999-30001, 39999-40001, and 49999-50000 Hz /you'll have to normalize yourself...). This would require a rather low order filters, which is nice since it runs at the highest sample rate. For the next stage, 999-1001, 1999-2001, and so on. Higher order (about ten times), but ten times lower sample rate. And so on.
Roughly more stages => fewer multiplications per input sample (given that you do a proper design as above). This holds up to some limit which can be quite straightforwardly determined by searching some combinations.
After the first stage you will have $1000 + 2 \times 127$ samples (with 128 taps) and so on.

If you want to get the DC value you can just average the sequence. This is, apart from the possibly obvious observation that the DC level should be the average "offset" of the otherwise AC signal, also what the DC component of the DFT evaluates to.


I often see here formula expressed in term of $z$. But what is $z$? While searching resources for generating pink noise (and with your help in the comments and answers of other questions), I came to such kind of formula:

  $$ H(z) = { .041 - .096z^{-1} + .051z^{-2} - .004z^{-3} \over{}  -2.495z^{-1} + 2.017z^{-2} - .522z^{-3}} $$

  (from https://dsp.stackexchange.com/a/376/46389)

What is $z$?

I know it should be obvious for anyone familiar with digital signal processing but for people coming from a different field, it is not obvious to know what this represents and how it relates to the signal amplitude or frequency or any other "tangible" parameter.

Edit: amusingly enough, the site suggested me the z-tranform tag for this question. I suspect I should consider that as a clue.

$z$ is interpreted as the time advance operator. $z^{-1}$ is the time delay operator.

So for a difference equation like
$$ y[n]=a x[n]+b x[n-1] $$
in the $z$ domain
$$ Y(z)=(a+ b z^{-1}) X(z)$$


Window functions with rippleless spectra On Wikipedia I found the Hann-Poisson window, and the article claims the spectrum is smooth, but it isn't theoretically smooth, as it turns out. In practice you achieve partial smoothness by jacking up its $\alpha$ parameter.

Can I find what I'm looking for in a function with finite support? Is it mathematically possible?

That's an interesting question, I personally don't know if, in theory, that's possible or not (would be interesting to investigate). But as a practical matter perhaps you could design a Chebyshev Type II filter, take the IFFT, and use a sufficiently large number of taps. Chebyshev Type II filters have ripple in the stopband but not the passband (although that is monotonically decreasing). I'm not sure how much the truncation would influence the ripple, or if you can handle that large of a FIR. But it might get you close enough depending on your application.

This page details some of the results.
http://iowahills.com/B2PolynomialFIRFilters.html


Physically understanding autoregressive model (really basic question) I would like to basically understand what an autoregressive model is used for (so I don't really attempt maths in the answers).

I just started a signal processing course and the model was introduced.

This is what I understood but I don't think I understood it well :

We have $n-1$ samples of a signal and we want to guess what the $n$'th will be.

With Young Walker, we have an approximation of the $n$'th sample that is : $\tilde{x}[n]=\sum a_k x[n-k]$

But $\tilde{x}[n]$ is just an approximation of $x[n]$.

So we have : $x[n]=\tilde{x}[n]+\epsilon$ where $\epsilon$ is the "error" I make when I use $\tilde{x}[n]$.

The autoregressive models assumes that $\epsilon$ is a white noise with average zero (and this makes sense because there is no reason in our model to have more often $\tilde{x}[n]>x[n]$ than $\tilde{x}[n]<x[n]$ ).

Did I understand the physical meaning of the model well ?

Thanks !!

Your understanding is right. Just a few points...

I think you meant Yule-Walker equations.
It is better to distinguish $N$ from $n$. In an AR($N$) model we consider last $N$ terms to approximate the next term.
An important interpretation that you missed was the fact that when we consider
$$\epsilon=x[n]-\tilde{x}[n]=x[n]-\sum_ka_kx[n-k],$$ or in a more standard way (representing the process by variable $y$ instead of $x$
$$y[n]=\sum_ka_ky[n-k]+\epsilon[n],$$
it can be seen as an input-output relationship of an LTI filter, $\epsilon[n]$ being input and $y[n]$ the output. More specifically,

   AR model can be seen as the output of an all-pole IIR filter which is excited by a white noise input.


Demodulating upper sideband (USB) signals A modulating signal $m(t)= \sin(2 \pi \nu_m t)$ is transmitted via a carrier of frequency $\nu_c$ using upper sideband (USB) modulation. Show that the USB modulated signal can be demodulated using a reconstructed carrier signal and a low pass filter.

Attempt:

USB (SSB signals in general) can be demodulated by multiplying them with a reconstructed carrier $\cos(2\pi \nu_c t).$ So I first tried to write an expression for the transmitted USB modulated signal:

$$f(t)=A \Big[ m(t) \cos(2 \pi \nu_c t) - \hat{m}(t) \sin(2\pi \nu_c t) \Big]$$

$$f(t)=A \Big[ \sin(2 \pi \nu_m t) \cos(2 \pi \nu_c t) - \hat{m}(t) \sin(2\pi \nu_c t) \Big] \tag{1}$$

Where $A$ is some amplitude, and $\nu_c >> \nu_m$. So to proceed further, I think we need to calculate the Hilbert transform $\hat{m}(t).$ It is possible to do this by finding the imaginary part of the analytic signal:

$$m_a (t) = m(t) + j \hat{m} (t)= 2 \int^\infty_0 M(\nu) \exp(j 2 \pi \nu t) d\nu \tag{2}$$

But the Fourier transform of $m(t)$ is tedious, so I am wondering if there is a simpler way of doing this?

Apart from scale factors, your USB signal is

$$f(t)=m(t)\cos(2\pi\nu_ct)-\hat{m}(t)\sin(2\pi\nu_ct)\tag{1}$$

(as you've correctly stated in your question). Now if you multiply with a (coherent) carrier $2\cos(2\pi\nu_ct)$, you get

$$\begin{align}g(t)&=2\big[m(t)\cos(2\pi\nu_ct)-\hat{m}(t)\sin(2\pi\nu_ct)\big]\cos(2\pi\nu_ct)\\
&=2m(t)\cos^2(2\pi\nu_ct)-2\hat{m}(t)\sin(2\pi\nu_ct)\cos(2\pi\nu_ct)\\&=m(t) (1+\cos(4\pi\nu_ct))-\hat{m}(t)\sin(4\pi\nu_ct)\tag{2}\end{align}$$

Applying a low pass filter to $(2)$ will filter out all terms at $2\nu_c$, leaving you with the message signal $m(t)$.


Interpolation based on sinc function I implemented an interpolation method in C++ based on this equation

y(t) = x[n]*sinc((t-nT)/T)

where:
n - index of input signal,
x - input signal buffer,
t - is an actual time,
nT - time of input samples,
T - period of input signal.

I calculated it on 33 samples:1+2*16. I tested it on Matlab and C++ on VS 2010. It works great on close frequencies (eg. 24kHz to 20kHz) but when I tested situation where the input signal is 24kHz and has a peak on 8kHz and then I wanted to interpolated it to 700Hz it gave me a peak at 300Hz so it is totally wrong. I expected that the main couse of this situation is that I did not use the cutoff frequency for low pass filter in this equation above. Seriously, I do not know how to add low pass filtration to my method.

Thanks for all your suggestions or advise.

I'm still not quite sure I understand what you're trying to achieve, but here's an example I cooked up in scilab.

The input is 33 points of a sine wave:

fs = 100; // 10kHz sampling rate
Tin = 1/fs; // Input sampling period
N = 16;
tin = [-N:N]; // Times of samples
f0 = 5;

x = sin(2*%pi*f0/fs*tin + 1.123791823)

clf
subplot(211)
plot(tin,x)

The sampling frequency is 100Hz, the frequency of the sine wave is 5 Hz.

Then, I try to interpolate this up to a sampling frequency of 700Hz using the sinc formula:

Tout = 1/700;
tout = -N:Tout:N;
y = zeros(1,length(tout))
for n = 1:length(tin)
    y =  y + x(n)*sinc(%pi*(tout-tin(n)))
end

plot(tout,y,'r:')
title('Input (blue line) and interpolated (red dotted line)')

subplot(212)
dx = splin(tin,x);
x_interp = interp(tout,tin,x,dx);
plot(tout,y-x_interp)
title('Error')

The results are as plotted below.


What is the relation of the transition band's width and the filter order for the FIR windowing method When designing an FIR filter using the windowing method, how can one estimate the filter order ? It's obvious that the type of window and the transitions band's width has some effect on the order.

There are only heuristic formulas for estimating the filter order. For a Kaiser window (which is probably the most frequently used window for filter design) the required filter order can be estimated from [1]

$$M=\frac{A-8}{2.285\,\Delta\omega}\tag{1}$$

where $A=-20\log\delta$ ($\delta$ is the maximum deviation from the desired response), and $\Delta\omega$ is the (smallest) transition bandwidth. This formula is of course only valid for the approximation of ideal frequency selective filters (low pass, high pass, etc.). This formula is implemented in Matlab's kaiserord.m function.

[1] Oppenheim, Schafer, Buck, Discrete-Time Signal Processing, 2nd ed., p. 476


How can this equation hold $$\sum_{m=0}^{N-1} \sum_{n=0}^{N-1} |u(m,n)|^2 = \sum_{k=0}^{N-1}\sum_{l=0}^{N-1} |v(k,l)|^2$$? I am trying to learn digital Image processing by myself and now stuck at a problem in the two dimensional unitary transformations. It states that let $U$ be the input image and $V$ be the transformed image and $A$ is the unitary transformation matrix.
To prove above equation I am rather trying to prove $$U^2=V^2$$ which I think is equivalent to proving equation in question.We have $$ V=A U A^T $$ now to prove the above relation: $$V^2=V^{*T} V$$ where
 $$ V=A U A^T \mbox{ and } V^{*T}=A^*[A^*U^*]^T$$ which on solving gives me $$ A^* U^{*T} U A^T $$ which gives $$ A^* U^2 A^T $$but it should be equal to $$U^2$$, as I understand. Though $$ A^T A^*=I $$ as $A$ is Unitary matrix, but as much as I know we can't change order.

OK, so let's have a look at whether

$$\sum_{m=0}^{N-1} \sum_{n=0}^{N-1} |u(m,n)|^2 = \sum_{k=0}^{N-1}\sum_{l=0}^{N-1} |v(k,l)|^2\tag{1}$$

is equivalent to

$$U^2=V^2\tag{2}$$

Let's have a look at a simple $2 \times 2$ example:
$$
U = \left [
\begin{array}{cc}
u_a & u_b\\
u_c & u_d
\end{array}
\right ] \\
V = \left [
\begin{array}{cc}
v_a & v_b\\
v_c & v_d
\end{array}
\right ]
$$
so that
$$
U^2 = \left [
\begin{array}{cc}
u_a & u_b\\
u_c & u_d
\end{array}
\right ]
\left [
\begin{array}{cc}
u_a & u_b\\
u_c & u_d
\end{array}
\right ] = \left [
\begin{array}{cc}
u_a^2 + u_b u_c & u_au_b + u_bu_d\\
u_cu_a + u_d u_c & u_c u_b + u_d^2
\end{array}
\right ]
$$
and
$$
V^2 = \left [
\begin{array}{cc}
v_a^2 + v_b v_c & v_av_b + v_bv_d\\
v_cv_a + v_d v_c & v_c v_b + v_d^2
\end{array}
\right ]
$$

Then (1) becomes
$$
|u_a|^2 + |u_b|^2 + |u_c|^2 + |u_d|^2 =
|v_a|^2 + |v_b|^2 + |v_c|^2 + |v_d|^2
$$

The two are completely different.  Why do you think they are the same?


GNU Radio QPSK 4 Symbols I'd like to generate a simple QPSK signal in GNU radio, but so far I haven't been able to get it to perform the way I'd like and I'm not sure if this has to do with some misunderstanding I have about how QPSK works or if I just am not familiar enough with GNU Radio.

Essentially, my conception of QPSK modulation is this - 2 bits are fed into the modulator (4 possible symbols - 00, 01, 10, 11) and, based on which combination of those 2 bits was received, one of 4 possible waveforms is generated (corresponding to +/-.707 +/- .707j). Thus, when using a QPSK modulator in GNU Radio I would expect to see, in the constellation diagram, a square with 4 points at each corner, but instead I see diamonds with thousands of points, even in the absence of any introduced noise. When I print the output of the QPSK modulator to a file sink, I also see more than 4 distinct points.

Is there a way to take a random input stream of 0's and 1's and produce a modulated signal from it, like QPSK, that has only one of four values? It is important for the simulation that I'm running that there are only 4 possible modulation values before noise is added to the system. Also, I will be sending these signals between USRP's, so if there is some way to generate a QPSK signal and transmit it on a carrier frequency between USRP's that doesn't rely on the PSK Mod block (I'm not too familiar with how to implement a raised cosine filter) that would be appreciated.

It is important for the simulation that I'm running that there are only 4 possible modulation values before noise is added to the system.

So, you made a mistake when building your simulation, or your simulation already includes a frequency offset. The way you ask this, you probably did the former, but since we have zero knowledge of what exactly you've built: Can't help you with that, make sure your receiver simulation does exactly the same as what your transmitter simulation does. ("fix your program!")

I'll answer as if this was a system that's generally correctly designed:

  Thus, when using a QPSK modulator in GNU Radio I would expect to see, in the constellation diagram, a square with 4 points at each corner, but instead I see diamonds with thousands of points, even in the absence of any introduced noise.

That's not surprising! No two oscillators are exactly the same, and so your receiver and transmitter have a frequency offset.

A frequency is the linear change of phase over time. That means that your symbol phase at the receiver has a phase offset that grows constantly every time you look. Hence, the "diamond".

You need what you've been taught in your digital communications course: a frequency synchronization!

You can do that with a second order costas loop.


How would a Linear Gray Level Transform affect an image? There are three basic gray level transformation.

  Linear
  Logarithmic
  Power – law

  The overall graph of these transitions has been shown below.

I am assuming a Linear Gray Level Transform is somehow related to the Identity line from the graph.

  Identity transition is shown by a straight line. In this transition, each value of the input image is directly mapped to each other value of output image. That results in the same input image and output image. And hence is called identity transformation.

From this statement I understand that in order to convert a colored image to a gray scale one (which I think Linear Gray Level Transform should do - converting to gray scale) I have to take each value of the input image and make a function that does a kind of conversion for each value.

I surfed the internet, but didn't managed to find a Linear Gray Level Transform example.

My question is if this is the right way (some kind of a small algorithm, if these are the steps) for a Linear Gray Level Transformation or if there is an example of images or an explanation of how does this Linear Gray Level Transform works.

The graph you've shown above describes a pixel-wise operator, basically a function which maps input pixel values $r$ to a new value $s$ in the output image. i.e. each of the functions in the figure you provided describes a transformation $s = T(r)$.

These transformations are typically used to improve image contrast through histogram equalization, basically spreading out the pixel values to make it easier for us to distinguish objects in the image.

To implement something like this, you have to iterate through every pixel, get its grey level $r$, and replace this value with the output $s$ of the function $T(r)$.

For converting color to greyscale, it's the same idea, although this time your function $s = T(r, g, b)$ has three inputs (r = red, g = green, b = blue), since each pixel in a color image will return these three values. Implementation is the same, just iterate through every pixel and replace these three values with just one value $s$ based on the function.

While the identity transform above (given by $s = T(r) = r$) is linear, linear transformations do not necessarily mean identity.

An example of a linear rgb color to greyscale transformation (the one used by MATLAB) is $s = T(r, g, b) = 0.2989 * r + 0.5870 * g + 0.1140 * b$. However, as long as you keep the function in this form (and just change the coefficients to constant values), your transformation will remain linear; this does not guarantee that the transformation will produce good results, though.

EDIT: What I described is only a subset of linear transformations. A linear transformation is any function $T(r)$ where $T(A * r1 + B * r2) = A * T(r1) + B * T(r2)$


Why is discrete cosine preferred to FFT in neuroimaging GLM High-pass filtering is often used in neuroimaging data analysis. Commonly, whenever a general linear model is fitted to the data (as for instance in statistic parametric mapping) a number of columns of the design matrix are assigned to discrete cosines (DC). See the figure below for a graphical example.

I am curious why this is preferred to taking the fourier transform (FT) of the data, and applying a filter or cutoff. Arguments I think would speak for the FT alternative are:

You would not have to re-do the high-pass filtering upon every application of a GLM
Since DCT is discrete I am guessing a FT would be more accurate in extracting all high frequeny components (looking at the figure above, I find it hard to believe that those few cosines cover all the low-frequency spectrum in between them)

Reasons I think DCT might be the preferred method:

Maybe it is faster than FFT.
Maybe there is a reason why it is desirable to perform the filtering concomitantly with your GLM fitting.

Could you tell me what you think about this?

Why is DCT preferred?
Would FFT not be better?

Say we have a signal,

The DFT assumes the signal is periodic, e.g.,

While the typical DCT assumes the signal is even at the boundaries, e.g.,

In the DFT example you can see there is a big jump from one side of the signal to the other. Since it assumes a periodic signal, this big jump introduces a lot of high frequency components. Thus the energy of the signal is spread out over more frequencies. For example, the log of the absolute values of the DFT is

Because the DCT assumes the signal is even at the boundaries, generally there is no big jump and thus no extra high frequency components. If this assumption is correct, more of the energy is concentrated in the lower frequencies. For example, the log of the absolute values of the DCT is

Actually, point at which the DCT assumes the signal is even can be moved. It doesn't have to be at the boundaries. However, if your signal generally has some large ramp component, then it makes sense to use the boundaries.

For the edge condition the DCT is equivalent (with multiplication by a scalar) to the DFT of the mirrored signal,

The log of the absolute values of the DFT of this signal is

which is the same as the DCT result above.

Note

The DCT doesn't automatically reduce boundary effects. Consider the extreme example of a sine signal, its DFT and DCT:

Because the even assumption didn't hold, the DCT actually resulted in many more high frequency components.

Why use it

I'm guessing because of drift in the nueroimaging signal, it is not periodic, and thus there would be a big jump from start to end. This would introduce lots of high frequency components with the DFT. Using the DCT assumes the signal is even and reduces the boundary effects.


Best way to resample a signal (image)? I have an original image $\ I(x_i,y_j)$ where

$\ x_i = i*\Delta x, 0<= i < n_i $

and

$\ y_j = j*\Delta y, 0<= j < n_j$

I am trying to implement a function that resamples the image to:

$\ I(x_k,y_l)$ where

$\ x_k = k*\Delta xx, 0<= k < n_k $

and

$\ y_l = l*\Delta yy, 0<= l < n_l$

The case when $\ \Delta xx <= \Delta x $ and $\ \Delta yy <= \Delta y $ : upsampling,
is quite easy (?).

Linear interpolation between nearest neighbours should do a decent job?

Downsampling on the other hand is more tricky.

To make sure Nyquist is followed so that no aliasing occurs; I must lowpass filter before decimating.

How do I do this an good yet simple manner?

I was thinking maybe using the FFT to do both upsampling and downsampling:

take the FFT of the image
zero pad or truncate the FFT
inverse FFT

In the case of zero padding this amount to sinc interpolation (which is better than linear interpolation)?

In the case of truncation this amount to (ideal) lowpass filtering?

So this should be a both simple to implement and very good solution?

I guess one thing to look out for is edge effects.

How do I best avoid this?

I was thinking repeating the image twice in x and y and then circulary shifting it half the width and height of the image.
That should make the padded image cyclical periodic.

An FFT approach to filtering/resampling should work fine, but you don't want to just truncate or pad with zeros.  That will add effects to the image.  See here for more information about FFT resampling.

You could also just do a FIR filter.  If you upsample by a ratio of $\frac{p}{q}$ then you would want your anti-aliasing filter to have a normalized bandwidth of $\frac{q}{p}$.  You could make the bandwidth a little smaller to get rid of more alias energy, as desired.  If you downsample by a ratio of $\frac{p}{q}$ then you would want the filter bandwidth to also be about $\frac{p}{q}$.  Again, you could make it a little smaller to get rid of more alias energy.

Once you have your filter you just need to do a two-dimensional convolution and trim the extra resulting samples on each edge.  Many signal processing libraries have built-in two-dimensional filtering routines that can do this for you.


Numerically Controlled Oscillator (NCO) for phasor implementation? If I want to simulate a digital oscillator (or phasor) to modulate an arbitrarily length signal:
$$y(t) = \cos(2\pi f_ct)\quad\text{where}\quad t = \frac{n}{f_s}\quad\text{for sample}\quad n$$

What is the advantage of using a numerically controlled oscillator (Answered here) versus simply incrementing $n$ over the length of the input?

main():
  n = 0                 # total samples processed
  N = len(input)        # num samples in input
  blksz = N / 16        # don't create phasor all in one call
  # calculate phasor over multiple passes
  while ( n < N )
      update_phasor(blksz)
# calculate phasor
update_phasor(blksz):
  for (i = 0; i < blksz; i++)
      y[i] = exp(j*2*pi*fc*n/fs)
      n++

The NCO referenced above and the pseudocode implementations give different results and I am trying to understand which makes more sense.

Also, what is the best way to synthesize the phasor if the desired frequency needs to change on the fly?

The NCO is a cyclical counter that can go on indefinitely but is otherwise similar to what you suggest in that you are increment n to set the output rate.  It basically is a look up table of all the values in one complete cycle, and "wraps" on an overflow so that it will output continuous cycles with no discontinuity.

I think the NCO is ideal for what you are trying to do; given its simplicity, it's ability to run indefinitely keeping track of its own position in time, its quantified noise levels that you can set to what you need, it's fixed-point implementation (with no multipliers) and its ability to change the frequency of your phasor "on the fly" as you require.  I think this bottom lines the difference with the alternate approach you describe, which would not be as efficient (given the multiplications required and possible implementation in floating-point without careful mapping and scaling in which case you may as well go down the NCO path).

A little more theory may help you to see all the advantages (and simplicity) of an NCO.

First, referring to the diagrams below for the basic NCO Architecture:  A digital Frequency Control Word (FCW) sets the count rate to an extended precision accumulator (counter), optionally a Phase Control Word (PCW) for phase modulation can then be added to the output of the accumulator. The most signficant bits of this summation are then used as the address pointer to a Look-Up-Table (LUT) which holds the values for one complete cycle of a sine wave (you could also imagine having two pointers for sine and cosine to enable a complex I and Q output).

Now see the same block diagram with a mathematical view below that helps give further technical insight into the operation of the NCO. An accumulator (counter) is the digital counterpart to an integrator (if you don't see that right away, imagine sending all 1's into a counter, the output would be a ramp: 1,2,3,4,... just as you would expect an integrator output to be with a constant level at its input).   The input FCW which is just a digital signal, which can change will time (as you were looking for), is a waveform that represents frequency vs time. What value of frequency for each digital word you can input I will elaborate on later, but for now know that its value at any given time is directly proportional to the output frequency. The integral of frequency is phase (and if you are less familiar with that it may be easy to see that frequency is a change in phase versus a change in time, therefore $f=d\phi/ dt$; frequency is the derivative of phase, and therefore phase is the integral of frequency.) Since our FCW input to an accumulator is the digital representation of a frequency quantity, and the accumulator is a digital integrator, then the value at the output of the accumulator represents phase versus time (which is why we can add a phase offset with PCW at this point if desired) and the accumulator counts from 0 to 2$\pi$, rolling over upon overflow.

Since the accumulator output represents phase that is changing with time, and we want to generate a sinusoid output ($sin(\theta)$), then we can simply use a LUT to perform the trigonometric function. (Note: If you have plenty of extra cycles but no memory, other techniques to calculate the sine of an angle can be done, notably the CORDIC algorithm).  Beautiful right? So now how do we decide on specifics to design our NCO, and what happens when we lose all the least signficant bits in our phase word? Read on!

First, the accumulator sets the frequency resolution, and usually an extended precision accumulator is used with either 24, 32 or 48 bits typically used depending on the application. This is easy to see, imagine first FCW =1: The accumulator will step through every value, meaning the address pointer to the LUT will also step through every value in the stored sine wave, so the sine wave output will be at the slowest rate, and that rate will be as given by the "step size" in the formula below. Why step size? Because then imagine setting FCW =2, and the counter will now count by 2's and therfore go twice as fast before rolling over (and also upon rollover the counter must continue to count and that is why the NCO will continue to output the desired sine wave indefinately), put in FCW = 3 and it will count 3 times as fast etc.  Therefore,

$$F_{out}= FCW\frac{f_{clock}}{2^{accum_size}}$$

So regardless of how many bits we decide to use for the LUT, the output frequency is strictly set by this formula and nothing else.

Now to briefly explain phase truncation and primary considerations: Phase truncation is when we decide to only use the most significant bits of the accumulator output to send to the LUT, and in doing this we are truncating the phase word (rounding down). To understand the implications of this, first consider the diagram below of what would happen if we did not have any phase truncation (meaning a very very large look up table, or very course frequency step size if the accumulator is small). What this picture is showing is that an impossible implementation containing a perfectly sinusoidal analog source at the specific frequency shown (with no phase noise), sampled with a perfect 100 MHz clock with a perfect 12 bit A/D converter will produce the IDENTICAL results of the NCO output, using a noisy 100 MHz clock. In fact for the NCO with no phase truncation, all of the output frequencies in multiples of $f_{step}$ as provided by the formula above will be this precise, with quantization noise at the output being the only noise source (which you can control by setting the output word width). You can see if you imagine different cases that without phase truncation, the look up table will provide the exact output at any given point in time that is required (limited to the quantized frequency choices with FCW, but that step size can be very small with large accumulators). The waveform will be very smooth without any skips or hiccups, in other words, pure.

So this is great, consider the example with a 32 bit accumulator, and 12 bit output of the LUT; providing very fine frequency resolution with great spectral purity (and 6 dB better for every additional bit you add to the output width)...until you get to the memory requirements! An herein lies the motivation to consider phase truncation.

Phase Truncation

With phase truncation, memory requirements are significantly reduced at the expense of an additional noise source (phase noise). What we will see if that the noise will be well understood and can be planned to be well below any given requirements (as a trade with memory required).

Also to mention for memory optimization only a quarter cycle is needed since the remaining portions of the cycle can be derived from the first quarter cycle. There are many other memory optimizations as well such as interpolation between values (most common), and to mention without explanation, Hutchison Algorithm and Sunderland Algorithm, as well as the Cordic Rotator previously mentioned.

The phase noise pattern itself from phase truncation will be a sawtooth function of phase versus time, representing the truncated values that are missing.  From this, the useful relationship of SNR to phase truncation is given as in the picture below.  Here SNR is the power of the desired sine wave output relative to the power of all spurious output due to the phse truncation. This formula applies when the small angle criteria applies (when $sin(\theta) \approx \theta$) and comes from the rms value for a sawtooth function (or equivalently the standard deviation for a uniform distribution) which is $\frac{D}{\sqrt{12}}$ where D is the peak to peak height of the ramp or width of the distribution.   This formula is combined with the quantization noise contributions from the LUT for a digitized signal (using a similar formula: 6 dB/bit + 1.76 dB also derived from our $\sqrt(12)$ factor since quantization noise can be modeled as a uniform distribution!), to address all noise sources in the NCO.  To use this formula, the number of bits in the formula is the number of bits sent to the NCO (number of phase bits not truncated).

Finally we may be interested in the spurious free dynamic range (SFDR) which would be the power level of the strongest spur relative to our output signal (as opposed to the sum power of all spurs in SNR).  The power of the strongest spur due to phase truncation is simply 6.02 dB/bit, where again bit is the number of bits sent to the LUT. (This can be derived by taking the Fourier Transform of the ramp pattern which represents our phase error, again applicable to small angle approximations). All the spurs are integer harmonics of the fundamental output frequency, many of which will have digitally folded in the first Nyquist zone of our implementation, as suggested in the diagram below. Unlike the diagram, the 2nd harmonic is not necessarily the strongest spur, but helps to give context to the idea of spurs and SFDR.

Dithering

Dithering is the process of adding a small amount of noise (for example, use a LFSR generator as the PCW input) which will improve the SFDR at the expense of SNR. The overall noise power increases (due to our additional additive noise), however the spurious levels can be substantially reduced in the process.


Whitened Matched Filter I am seeking for an advice on whitened matched filtering technique. I have looked into the literature and I do understand its the purpose and how to select the filter in order to achieve the desired response.

However, what I don't understand is that if the signal is matched to the channel filtered signal plus noise using the matched filter, and then if we apply the whitening filter, we get the inverse again which should actually be the channel filtered signal and noise. So in essence we end up where we started from.

I am sure I am missing something but any comments will be highly appreciated.

Thanks
Milos

In the ideal AWGN channel we have the received signal is $r(t)=s(t)+n(t)$, where $s(t)$ is the transmitted signal and $n(t)$ is white Gaussian noise. In this case, the transmitted symbols can be estimated using a matched filter whose output is sampled at the symbol rate. Note that in general the noise at the output of the matched filter is correlated, and no longer white; however, at the sampling times the noise is uncorrelated.

In the ISI channel we have $r(t)=c(t) \ast s(t) + n(t)$, where $c(t)$ is the channel response. We can think of this system as an AWGN channel where the transmitted signal is $g(t)=c(t) \ast s(t)$, and then we can use $g^*(-t)$ as a matched filter. However, in this case we no longer have uncorrelated noise samples. Correlated noise is more harmful, and thus this situation is undesirable.

A whitening filter after the matched filter decorrelates the noise samples and improves the system's error rate performance.

Note that the whitening filter does not revert what the matched filter did. The reason is that the purpose of the filter is not to turn the noise back into white noise; its purpose is to decorrelate the noise at the sampling instants. If the transmitted symbols are $a_k$ for integer $k$, and the (discrete) whitening filter has taps $f_n,\,n=0,1,\ldots,L$, then the output of the whitening filter is $$v_k=\sum_{n=0}^L f_n a_{k-n} + w_k,$$ where the noise samples $w_k$ are uncorrelated. The symbols $a_k$ can then be optimally obtained from $v_k$ by the Viterbi algorithm, or (perhaps sub-optimally, but easily) from another type equalizer (ZF, LS, etc.).

See also How to 'whiten' a time domain signal? and Decorrelating Stationary Colored Gaussian Noise -- Effect On The Desired Signal for more details.


How can you intutively figure what the impulse response will behave like? I was trying to understand convolution from this book and came across this figure:

Here, the input signal is a combination of a sine wave and a ramp.
In fig a, the input signal is convolved with a special type of impulse response that is said to be a 'low pass filter' and the output is just the ramp.

And in fig b, impulse response is actually a high pass filter and the output is just the sine wave.

What I'm trying to do is, just by looking at the wave-forms for the impulse response functions, trying to find an intuitive way of knowing what those impulse responses will do when convoluted with a given input signal. I can't 'see' how the first waveform will result only in the ramp as the output signal when given the mentioned input.

Can you guys figure this intuitively in your head? If so, then how?

Here's an idea...

In the first case, the impulse response has only positive coefficients. Groups of nearby samples will be averaged to yield the output signal. Thus, short term bumps/variations will be "averaged out", so the signal will be smoother. Furthermore, it looks like all coefficients sum up to 1 (there about 30 of them, in a triangle-like shape, the maximum is about 0.06, I suspect they sum to 1), so the overall amplitude / magnitude of the output waveform will be preserved. Last bit of information: the impulse response has its center of gravity at $n=15$, so the output signal will be delayed by 15 samples.

The second impulse response looks like a kronecker delta (at $n=15$) minus the impulse response from the previous example. Because convolution is a linear operation, the result will be the original signal shifted by 15 samples, minus the smoothed signal obtained in the previous example. In other words, the result will be whatever has been smoothed out in the previous example because of averaging - only the fast wiggles will remain.

Summing the output signals from example 1 and example 2 will yield the input signal delayed by 15 samples.


Finding the dominant tone in a signal I am given 256 samples at a sampling rate of 48k and am asked to find the frequency of dominate tones and the amplitude of the largest tone. The data looks like this when plotted:

I assume I must use a FFT transform to do this. I've tried this in matlab:

Fs = 48000;
x = dataSet; % 256 samples data set
xdft = fft(x);
maxAmp = max(abs(xdft));
%Not sure how I can grab the tones here
freq = 0:Fs / length(x):Fs/2;

Not really sure how to proceed, Any ideas?

I think you are having trouble constructing your frequency axis properly. Once you do this, you can do a simple peak pick. I have re-written the code for you:

Fs = 48000;
x = dataSet;            % 256 samples data set
fftLength = length(x);  % Always make sure to be at least as long as your data.
xdft = fft(x,fftLength);
maxAmp = max(abs(xdft));

freq = [0:fftLength-1].*(Fs/fftLength); % This is your total freq-axis
freqsYouCareAbout = freq(freq < Fs/2);  % You only care about either the pos or neg
                                        % frequencies, since they are redundant for
                                        % a real signal.

xdftYouCareAbout = abs(xdft(1:round(fftLength/2))); % Take the absolute magnitude.

[maxVal, index] = max(xdftYouCareAbout); % maxVal is your (un-normalized) maximum amplitude

maxFreq = freqsYouCareAbout(index); % This is the frequency of your dominant signal.


What's the ideal FFT window for measuring a group of signals of differing amplitudes but close in frequency? In my system, I need to measure the amplitudes and phases of four lasers, each of a different, known wavelength, chopped at a 50% duty cycle via TTL at frequencies spaced evenly (or not) around some fundamental frequency f0. The lasers are attenuated by some analyte medium, and the resulting laser power (which is measured by a single detector) tells me what I need to know about the medium.
The signals are recovered by performing I/Q demodulation. The application is similar to the one in this question. The signals' individual frequencies can be controlled by me, but should not be more than 6-7 Hz away from f0.
In single-signal applications of my setup, I've used the Exact Blackmann window, which is known to be good for single-frequency measurements, at least according to National Instruments. My question is, is there a more robust choice for recovering four extremely close signals?
Other details: sample rate is 8192 or 16384 Hz for two seconds, f0 is in the neighborhood of 1230 Hz. Typically I have the frequencies evenly spaced, at f0±6 and f0±2, so they're all 4 Hz apart. The input signal is demodulated and windowed four times, once at each of the frequencies of the lasers. Followup question: could spacing them unevenly help?

As for ideal windows for FFT frequency resolution and dynamic range, many assume that the Gaussian window would be best since it has the minimum time-bandwidth product, and due to this it is the filter most often used in spectrum analyzers. However, this is only true when the time domain extends to infinity. For finite time durations the window that actually has the minimum time-bandwidth product is the Discrete Prolate Spheroidal Sequence (DPSS, or Slepian) window which maximizes the energy concentration of the main lobe. The Kaiser window very closely approaches this optimum and both have an controlling parameter to trade frequency resolution (ability to discern two closely spaced frequencies) and dynamic range (ability to discern two frequency tones of different amplitudes). For further details on the comparison of the two, see this link by Julius Smith: https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html
That said, the signal after IQ demodulation could be windowed and then FFT'd to provide the phase and amplitude (relative to the demodulation signal used) of closely spaced modulation frequencies. For this application I would use a Kaiser window with Beta = 8 which provides 80 dB dynamic range of signals captured over a 2 second duration while being able to resolve two signals that are 4 Hz apart.


Is possible to perform algebric operation with bitwise operations in PDM domain? I have to compute the some operations (say sum and subtraction) between 2 signals.
I have data in PDM domain:
x1_pdm(t)=[0101011....];
and
x2_pdm(t)=[0101011....];
What I do now is:
convert x1_pdm and x2_pdm in PCM domain (low pass filter and decimation factor M=64) and that compute the sum and difference (or any other operation). Is it possible to perform this sum and sub in PDM domani with bitwise operation ? and then compute one pdm to pcm conversion...
thank you
mattia

Yes, this is sort of possible.   While PDM signals are full of high-frequency noise, at lower frequencies the "density" domain is limited to signals representing values from -1.0 thru +1.0.   This makes it difficult to implement true addition/subtraction since 1+1 exceeds the range of possible outputs, however if the output signal is scaled back into the allowable range, so averages like 1/2(A+B) and 1/2(A-B) can be calculated by multiplexing the two bitstreams into a single output stream (at twice the bitrate).

The best explanation of this I've seen was in a series of five articles titled "Signal Processing in the Density Domain" published in 2011 by Dave Van Ess in Electronic Design magazine.   The web copies of these articles exist but the equations and figures from the original versions seem corrupted by the web's evolution over time.   I'll include this link to the finale, Part V which wraps up how addition/subtraction and multiplication work.

While not easy, the missing figures can be retrieved from the internet wayback machine.


Confusion with ADC conversion in MCP3424 I am learning about ADC MCP3422 Chip, because I want to use ADC Pi chip with my Raspberry Pi which is using MCP3424 ADCs, but I am having trouble understanding how the ADC calculation is done in the datasheet.

(Data sheet links are available in the end)

From what I know of ADCs,

An ADC has a bit size.
and ADC tends to generate the total number of level based on the Reference Voltage.

for example,

If there is a 3-bit ADC then and the Reference Voltage is 2volts, then,

2^3 = 8-bits

Now 2volts will be divided in 8 levels.

8/2 = 0.25volts (each level will represent 0.25 volts

and,

000 = 0 -> 0.25 volts
001 = 0.25 -> 0.50
010 = 0.50 -> 0.75
011 = 0.75 -> 1.00
100 = 1.00 -> 1.25
101 = 1.25 -> 1.50
110 = 1.50 -> 1.75
111 = 1.75 -> 2.00

Now, in case of MCP3424, I have so much confusions

The bit size is 18.
Reference voltage is 2.048volts

2^18 = 262144

Does that mean there will be 262144 levels ?

So, 2.048/262144 = 0.000007812volts
Does that mean each level will represent, 0.000007812 volts ?

Also, the data sheet says,

The data rate for analogue to digital conversions is 3.75 (17 bit), 15 (15
bit), 60 (13 bit) or 240 (11 bit) samples per second. Data rate and
resolution can be configured within software using the I2C interface.

How the data rate is calculated here, and how come reducing the bit is increasing the Data Rate ?

MCP3424 Datasheet: www.pridopia.co.uk/pi-doc/mcp3424.pdf

ADC Pi Datasheet: http://www.abelectronics.co.uk/docs/stock/raspberrypi/adcpi2/Datasheet-ADCPiV2.pdf

If the conversion is 18-bit, there will be 262144 levels.

However, it seems to me that on the ADC Pi board, one of the differential inputs of each MCP3424 converter channel is grounded ; so this actually halve the range, and the first bit of a word read from the ADC will always be zero. That's why the ADC Pi datasheet reports half the resolutions compared to those published on the MCP3424 datasheet (17-bit, 15-bit, 13-bit, 11-bit).

There is a relationship between data rate and bit rate because the converter has a delta-sigma architecture. There is an internal conversion/integration clocked at a fixed 983040 Hz, and this conversion "budget" can be used at your whim to achieve different data rate / bit-rate ratios with:

$$s_r = \frac{983040}{2^{b_r}}$$

Intuitively (and this is a gross mis-simplification), reducing the sample rate gives the ADC more time to extract the finer details of each sample.


Compensate for FIR filter attenuation Is there a simple way to compensate for the reduction in amplitude of a low pass filter of a noisy sine wave signal that has been smoothed using a FIR filter so that the amplitude of the smoothed sine wave matches that of the approximately known amplitude of the noisy original signal? The approximate amplitude of the original data could be estimated empirically from the RMS value over one full (estimated) period. The FIR filter coefficients I am using are

h = [ 1 2 3 3 2 1 ] / 12

FIR definition:

$$ y[n] = \sum_{k=0}^{N} { b_k x[n-k] }  $$

Sinusoid signal definition:

$$ x[n] = M \cos( \alpha n + \phi ) $$

A whole bunch of math:

$$ y[n] = M \sum_{k=0}^{N} { b_k \cos( \alpha (n-k) + \phi ) }  $$

$$ y[n] = M \sum_{k=0}^{N} { b_k [  \cos( \alpha n + \phi ) \cos( \alpha k ) + \sin( \alpha n + \phi ) \sin( \alpha k ) ] }  $$

$$ y[n] = M \cos( \alpha n + \phi ) \sum_{k=0}^{N} { b_k  \cos( \alpha k ) } + M \sin( \alpha n + \phi ) \sum_{k=0}^{N} { b_k \sin( \alpha k ) }  $$

$$ y[n] = A \cos( \alpha n + \phi ) + B  \sin( \alpha n + \phi )  $$

$$ A = M  \sum_{k=0}^{N} { b_k  \cos( \alpha k ) } $$

$$ B = M  \sum_{k=0}^{N} { b_k  \sin( \alpha k ) } $$

$$ y[n] = M_2 \cos( \alpha n + \phi + \theta ) $$

$$ y[n] = M_2 \cos( \theta ) \cos( \alpha n + \phi ) - M_2 \sin( \theta ) \sin( \alpha n  + \phi ) $$

$$ A = M_2 \cos( \theta ) $$

$$ B = -M_2 \sin( \theta ) $$

$$ A^2 + B^2 = M_2^2 = M^2 \left[ \left( \sum_{k=0}^{N} { b_k \cos( \alpha k ) } \right)^2 + \left( \sum_{k=0}^{N} { b_k \sin( \alpha k ) } \right)^2 \right] $$

Your desired equation:

$$ \frac{M}{M_2} = \frac{ 1 }{ \sqrt{ \left( \sum_{k=0}^{N} { b_k \cos( \alpha k ) } \right)^2 + \left( \sum_{k=0}^{N} { b_k \sin( \alpha k ) } \right)^2 } } $$

I think I've done the math right.  The $b_k$s are your FIR coefficients and $\alpha$ is your frequency in radians per sample.  $M_2$ is the amplitude of your smoothed sinusoid and $M$ is the original amplitude.

So you want multiply your smoothed results by $ \frac{M}{M_2} $.

Hope this helps.  I just did this and haven't tested it.

Ced