defining the frequency range with fast fourier transform In the following example I have used various examples found on the web to understand fft and detecting signals with matlab.

fs = 1000;
t = 0:1/fs:1.5-1/fs;

x = 3*cos(2*pi*20*t + 0.2) + 1*cos(2*pi*30*t + 0.3) + 2*cos(2*pi*40*t + 2.4);
N = length(x);

% define the frequency components of the series
freq = 0:N-1;
freq = freq*fs/N;

% normalise the data and compute the fourier transform
X = fft(x)./N;
cutOff = ceil(N/2);
X = X(1:cutOff);
freq = freq(1:cutOff);

figure;
stem(freq,abs(X));
xlabel('frequency');
ylabel('amplitude');

Resulting in:

Why is it that the amplitude of the final plot are 1.5, 0.5, and 1 respectively? The initial amplitudes were 3, 1, and 2 so why is the amplitude of the signal half of that?

Also, could someone explain the section

% define the frequency components of the series
freq = 0:N-1;
freq = freq*fs/N;

I don't quite understand how this works.

There are a few things going on here.  First, one of the fundamental properties of the Fourier transform is that it conserves energy, i.e. the energy of the time domain is exactly equal to the energy of its Fourier transform.  What is the power of a sinusoid?  Assuming we restrict the calculation to integer multiples of cycles it is as follows (I suggest working out this equation for yourself if you are not sure about it):

$$
P(Asin(ft)) = \frac{A^2}{2}
$$

The energy of a sinusoid is thus $\frac{A^2}{2}T$, or, in the case of discrete data, $\frac{A^2}{2}N$.  The purpose of the following line of code-

X = fft(x)./N;

is to change the FFT output from an energy-based measure to a power-based measure ("based" because you have to square the transform results, just like you would in the time domain, to get the actual energy or power).  Power is usually more interesting since it is an inherent characteristic of the signal, whereas energy is dependent on the number of samples that you are using.

So, the power of your sinusoids is, according to the equation, 4.5, 2, and 0.5.  You got 1.5, 1, and 0.5, which, if you square them, become 2.25, 1, and 0.25.  You will notice that those results are exactly one half your expected power.  The reason for the halving is that you are only showing the positive frequencies from your Fourier transform.  That is what the following lines of code are doing-

cutOff = ceil(N/2);
X = X(1:cutOff);

They are cutting off the negative frequency part of the spectrum.  If you plotted the entire set of results it would look like the following-

The spikes on the left are the positive frequencies that you saw before.  The spikes on the right are the negative frequency components of the sinusoids.  Note that they have the exact same magnitude as their positive frequency partners, and thus comprise the missing factor of two in the power.

Lastly, you ask the following-

  Also, could someone explain the section

% define the frequency components of the series
freq = 0:N-1;
freq = freq*fs/N;

Samples sampled at frequency $f_s$ can represent frequencies from $-\frac{f_s}{2}$ to $\frac{f_s}{2}$ (see Nyquist Sampling Theorem), or, equivalently, frequencies from $0$ to $f_s$ (see aliasing).  The code implicitly assumes that the samples represent from $0$ to $f_s$.


Normalized Cross Correlation Operation According to my question about separability of Gabor filters in this link, I want now to convolve my image with this separable filter by using the normalized cross correlation operation. Assume my Gabor filter is G, my image is I. My Gabor is separated into Low-Pass gaussian filter f(x) and Band-Pass gaussian filter g(y). Therefore the image is convolved with the Gabor using the following equation:

I(x,y)*G(x,y) = (I(x,y)*f(x))*g(y).

But I want to achieve this separable convolution using the normalized cross-correlation operation described below:

Where ^G is the zero mean, unit normal version of the filter and H(x,y) represents a filter with all ones and the same size of the Gabor filter.

1) I didn't understand what is ^G. What should be its value? what it differs from G ?

2) How the normalized cross-correlation is computed for the separable Gabor? I don't know if I use correctly the formula: I(x,y)*f(x)*g(y) / I^2(x,y)*H(x,y) . I don't think that it's true. because I didn't understand what should be the value of zero mean, unit normal version of the Gabor.

From what the paper described, I think $\hat G$ is achieved by standardized normalization on $u$ and $v$ before you implement the formula (1) on $G$:

x = -filtSizeL : filtSizeR;
y = x;
u = x * cos(theta) + y * sin(theta);
v = -x * sin(theta) + y * cos(theta);
u = (u - mean(u))/std(u);
v = (v - mean(v))/std(v);

The filter is still separable after this normalization step, but I don't think you need to implement that, because it is just equivalent to select a different $\sigma$ value in the filter.


How to average MFCCs from different speakers saying the same word? I am trying to create a simple recognition system. i have recordings of several speakers saying the same word, let's say "west". I compute the MFCC feature vectors of word "west" from each speaker. If i have 10 speakers, i then have 10 matrices, each of them of dimensions N_frames * 13 (13 coeff and N_frames varies from speaker to speaker). I would now like to store a final matrix of dimension
M * 13 which will serve as a reference pattern for word "west". My question is, how to "average" or combine MFCCs from different speakers? (i call this process some kind of training)

What you can try is the DTW (Dynamic Time Warping) distance (https://en.wikipedia.org/wiki/Dynamic_time_warping).

for each instance of the spoken word "west", you compute MFCC's per frame and then calculate the DTW distance between each of the pre-recorded instances of word. this gives the rough estimate of the valid DTW distance. now when a new instance of the word comes we again calculate DTW with one of the earlier instance and if the DTW distance lies roughly in the range calculated before, it is the same word.


Sinc function's energy Given this signal:
  $$4000\,\mbox{sinc}(4000t)+D1000\,\mbox{sinc}(1000t)$$
  with $ D=6000$, find its energy.

To find the energy of sinc I use this transform:

$$A\,\mbox{sinc}(Bt) \rightarrow A\,\mbox{rect}(f/B)/B$$

and compute the energy $E$ using:
$$E=\frac{A^2}{B}$$

In this case what is the right answer?

A) $4000^2/4000+D^21000^2/10000$ or
B) $1/4000+A^2/1000$

Neither A nor B.

Hint: the Fourier transforms of the two sinc function can be added in the frequency domain due to the Fourier transform linearity property, but the energies do not since the spectrums overlap.

   Note that the energy formula comes from computing the product of the square of the frequency domain amplitude by the frequency range of a rectangular pulse. Applying your formula, the amplitude for a single $\mbox{sinc}$ is given by $A/B$ and covers a range of $B$, yielding an energy of $(A/B)^2 B = A^2/B$.
 The amplitude of the first rectangular function (corresponding to the $4000\mbox{sinc}(4000t)$ function in the time-domain) in the frequency-domain is $1$, and covers frequencies in the range $[-2000,+2000]$.
 Similarly, the amplitude of the second rectangular function (corresponding to the $D 1000\mbox{sinc}(1000t)$ function in the time-domain) in the frequency-domain is $D$ and covers frequencies in the range $[-500,+500]$.
 Now, computing the energy of more than one $\mbox{sinc}$ is a little trickier... You should realize that there is an overlap for frequencies in $[-500,+500]$, where the total amplitude is $1+D$.
 The energy in that range is thus $1000 (1+D)^2$.
 For the rest of the ${[-2000,-500], [+500,+2000]}$ frequencies, there is no overlap and the amplitude is that of the first rectangle function, i.e. $1$. The corresponding energy there is thus $3000$.
 This then gives you a total energy of $3000 + 1000(1+D)^2$. Substituting $D=6000$ gives $36,012,004,000$.


Carrier Frequency Offset Correction Limitations I know that carrier frequency offset is typically detected and corrected in OFDM systems. I was wondering if there is a way to determine the theoretical limit to how much offset can be successfully "compensated for" and still utilize the data successfully (specific to OFDM, if that matters). I'm thinking in the simplest terms where the transmitter and receiver are trying to transmit and receive at the exact same rate (for example, the receiver isn't simply oversampling). I know it has to do with sampling rate and the overall bandwidth of the OFDM packet...

edit: maybe I should have specified, but I'm referring to CFO correction done in the digital domain, post ADC (think DSP) instead of a voltage-controller oscillator type of control loop.

There is no fundamental limit but based on what frequency discrimination approach you are using but you will be trading sensitivity for range. To improve sensitivity you can use a cascade of discriminators to "tune in" on your signal of interest.  A trivial example of this is a case with a carrier at 2.4 GHz; regardless of what carrier frequency we start with as our initial "guess", we could band-bass filter the signal of interest with a fixed filter at 2.4 GHz (to eliminate content at other frequencies that would disrupt the discriminator) and then use a succession of frequency discriminators (which simply provide an output voltage proportional to the frequency of the input) to provide higher and higher accuracies of what the carrier is, which can ultimately be within the capture range of whatever you are using for OFDM. For example, this link provides a typical discriminator for OFDM:

Carrier frequency offset (CFO) estimation using cyclic prefix

And this link provides further details on carrier discrimination techniques in general:

What exactly is a 90 degree phase shift of a digital signal in FM demodulation appraoches?


How to normalize frequency response charts I want write a simplistic program that can transform a high-quality voice recording in a way that it sounds like it was recorded using a specific low-quality microphone.

I have obtained the mic's FR using a program that uses sine sweeps in a sort-of-anechoic room since there is no official spec available. Sine sweep was played by a good loudspeaker of which I have the FR specs.

Now these FR charts have different "dB ranges":

Mic FR: -90dB to -10dB
Speaker spec: starts at -45dB, moves to 85dB and stays mostly linear.

The goal is to normalize the mic FR data w.r.t to the speaker's FR spec and then use that result to create an (inverse) FIR filter.

Question: How to read, interpret and normalize these charts so I can do the math on them?

My guess would be to choose mic amplification during the sine sweep so that the max amplitude observed is 0dB (is this the same as the mic's clipping threshold?), and then divide by the min amplitude observed. And for the speaker, normalize in the 45–85dB range.

Example speaker frequency response chart:

Example microphone frequency response chart:

First of all, I wouldn't worry too much about the speaker response since it is relatively flat and the microphone has a much bigger roll-off.

Since you've captured the frequency response using sweep, why not to skip the whole part of designing the filter that mimics the frequency response and use the original impulse response? I don't know what kind of software you used for doing your measurements, but a tool such as REW can export the IR wav file. Then it's just a matter of convolving it with your input signal. If your device has some significant amount of self-noise, you can record it in a quiet environment and add it post convolution.

Regarding the decibel scaling, an easy thing would be to play a signal at a known level (for example 94 dB SPL) and record it with the microphone. Then check what is the dBFS level of that recording. This will allow you to map the dB SPL to dB FS. Once it's done, you can convolve your signals and analyse the spectra to see what level they simulate. If any gain adjustment is necessary then you can do it afterwards. This, of course, assumes that you know the original level of those recordings, which is generally not true. If you had the same sort of measurement for the device that was used to capture those recordings, then it's easy. You know what is the sensitivity of that device, you know how much it differs from your microphone, so you can correct for that either by adjusting the level of impulse response or applying gain post convolution. Finally, you can simply ignore all of that and set the gain of your impulse response so the frequency response peaks at 0 dBFS. This means no boosting - only attenuation.


MRI reconstruction using windowing based apodization I am trying to apply windowing to do MRI reconstruction. I have a 256 point one dimensional k-space samples which look as follows:

Here, the real part(green) is flat around 0 and the blue curve is the imaginary part. Now, what I would like to do is perform some apodization. For this, I construct a hanning window as follows:

from scipy import signal

win = signal.hann(256)
# sig is the
sig = signal.convolve(sig, win, mode='same') / sum(win)

The hanning window is the same size as the k-space vector (256). Next, I convolved the k-space signal with the hanning window, which looks like:

Next, I performed the inverse fft to reconstruct the image (which should have near uniform intensity everywhere) using:

from scipy import fftpack

rec = np.abs(fftpack.ifftshift(fftpack.ifftn(sig)))

This results in the recon image as:

I feel I have done something wrong as I was expecting to have a fairly uniform image rather than a spike. Do I need to do something post reconstruction to account for the fact that I have performed some windowing before the reconstruction?

You say:

I have a 128 point one dimensional k-space samples...
The hanning window is the same size as the k-space vector (256)...

Make sure that you have the appropriate sizes in your algorithm.

Next, I convolved the k-space signal with the hanning window...

Windowing is applied in $k$-space by multiplication - you simply have to multiply your window and the original $k$-space line. In image space, this is a convolution, hence you could take the Fourier Transform of the window function and convolve the image with it.

Another view:
Apodization is applied to suppresses certain frequencies. For example, the highest frequencies are simply skipped by setting them to 0 due to the multiplication of the window function and the $k$-space data. If you were to convolve the two, you would "spill" the information of the higher frequencies all over $k$-space - which certainly is not what you want to achieve by apodization.


When does the laugh track kick in? One of my weekend projects has brought me into the deep waters of signal processing. As with all my code projects that require some heavy-duty math, I'm more than happy to tinker my way to a solution despite a lack of theoretical grounding, but in this case I have none, and would love some advice on my problem, namely: I'm trying to figure out exactly when the live audience laughs during a TV show.

I spent quite a bit of time reading up on machine learning approaches for detecting laughter, but realized that's more to do with detecting individual laughter. Two hundred people laughing at once will have much different acoustical properties, and my intuition is that they should be distinguishable through much cruder techniques than a neural network. I may be completely wrong, though! Would appreciate thoughts on the matter.

Here's what I've tried so far: I chopped up a five minute excerpt from a recent episode of Saturday Night Live into two second clips. I then labeled these "laughs" or "no-laughs". Using Librosa's MFCC feature extractor, I then ran a K-Means clustering on the data, and got good results -- the two clusters mapped very neatly to my labels. But when I tried to iterate through the longer file the predictions didn't hold water.

What I'm going to try now: I'm going to be more precise about creating these laughter clips. Rather than do a blind split and sort, I'm going to manually extract them, so that no dialogue is polluting the signal. Then I'll split them into quarter second clips, calculate the MFCC's of these, and use them to train a SVM.

My questions at this point:

Is any of this making sense?
Can statistics help here? I've been scrolling around in Audacity's spectrogram view mode and I can see pretty clearly where laughs occur. In a log power spectrogram, speech has a very distinctive, "furrowed" appearance. In contrast, laughter covers a broad spectrum of frequency quite evenly, almost like a normal distribution. It's even possible to visually distinguish applause from laughter by the more limited set of frequencies represented in applause. That makes me think of standard deviations. I see there's something called the Kolmogorov–Smirnov test, might that be helpful here?

(You can see the laugh in the above image as a wall of orange hitting at 45% of the way in.)
The linear spectrogram seems to show that the laughter is more energetic in lower frequencies and fades out towards the higher frequencies -- does this mean it qualifies as pink noise? If so, could that be a foothold on the problem?

I apologize if I misused any jargon, I've been on Wikipedia quite a bit for this one and wouldn't be surprised if I got some jumbled.

Based on your observation, that spectrum of the signal is sufficiently distinguishable, you can use this as a feature to classify laughter from speech.

There are many ways you can look at the problem.

Approach #1

In once case, you can just look at vector of the MFCC. and apply this to any classifier. Since you have many co-efficient in frequency domain, you might want to look at Cascade Classifiers structure with boosting algorithms such as Adaboost based on this, you can compare between speech class vs. laugh class.

Approach #2

You realise that you speech is essentially a time varying signal. So one of the effective ways to do it is to look at time variation of the signal itself. For this, you can split signals in batches of samples, and look at the spectrum for that time. Now, you may realise that laugh might have more repetitive pattern for a stipulated duration where as speech inherently posses more information and hence the spectrum variation would be rather bigger. You can apply this to HMM type of model to see if you consistently remain within the same state for some frequency spectrum or you continuously keep changing. Here, even if occasionally the spectrum of speech resembles to that of laughter it will be more time changing.

Approach #3

Force to apply LPC/CELP type of coding on the signal and observe the residue. CELP Coding makes a very accurate model of speech production.

From the reference here: THEORY OF CELP CODING

  The redundancies in the speech signal are almost removed after the
  short term prediction and long term prediction of the speech signal
  and the residual has very little correlation left in it. Then an
  excitation is searched which synthesizes the speech and the codebook
  index and gain are searched from the fixed codebook. The optimum
  codebook index selection criterion is based on MMSE between the
  locally synthesized speech and the original speech signal.

To put it simply, after all the speech which is predicted from analyser is removed - what is left is the residue which is transmitted to recreate exact waveform.

How does that help with your problem? Basically, if you apply CELP coding, speech in the signal is mostly removed, what remains is residue. In case of laughter a majority of the signal might be retained because CELP will fail to predict such a signal with vocal tract modelling, where as individual speech will have very little residue. You can also analyse this residue back in frequency domain, to see if it is laughter or speech.


Filter bank for complex spectrum I'm using a complex-valued baseband signal as my input for subband adaptive filtering. Therefore i have a complex input spectrum.

My problem is now, how to split my complex spectrum into $N$ equally spaced subbands?

I know i could use a DFT filter bank with impulse responses of

$$\displaystyle h_i(n) = p(n)e^{(j2\pi i/N)n} ,\quad i = 0,1,\ldots,N-1$$

where $p(n)$ is the impulse response of my lowpass prototype filter.
But what bothers me is that $H_0(z)$ will be real valued and so $|H_0(e^{j\omega})| = |H_0(-e^{j\omega})| $ holds, and hence the complex spectrum isn't split right into $N$ subbands.

Am I missing something there, or is it not possible to use DFT filter banks for complex-valued input signals?
Are EMFBs (exponentially modulated filter bank) a workaround to solve this problem?

Typically what you’ll see in an application like this is a polyphase filter bank, which is designed with a single low pass filter prototype. This keeps your prototype filter real, and they can be quite efficient in comparison to directly designing an entire bank of individual bandpass filters.

The short of it is this is a multi-rate technique which aims to simplify the computational complexity while being mathematically equivalent to a direct implementation. By having only one prototype filter, the coefficients can be kept real and the filter is essentially frequency shifted N times in the backward DFT application after the polyphase filtering. This deals with your issue of complex conjugate symmetry in the frequency domain that you noted before.

This previous post had some questions regarding polyphase filter banks, and has a link to a really nice MATLAB implementation.

I’ve used filter banks like this in practice (my field is radar signal processing), where we almost exclusively work with complex IQ signals, so I can vouch for their viability in this application.


Fourier domain: temporal versus spatial I am little bit confused over these two. As per my understanding, spatial domain is the usual method we did in matlab by using the function fft2, which will return a complex value for each pixel. What about temporal? If it is a sequence of images, how do I obtain the 'temporal Fourier domain'?

If you looking at the image, is it true that the temporal Fourier domain shall have the same shape as the spatial domain image (shown in figure is the amplitude of Fourier domain)?

I have done very little work using 2-dimensional FFTs.  Most of my work deals with speech, so the short-term Fourier transform applied to a 1-dimensional speech signal is what I understand.

I do not know of a "canned" function, particularly in Matlab, that would realize the time-varying "temporal Fourier domain" characterization that you require.  However, I would suspect you could apply a 2-dimensional FFT to a set of time-slices extracted from the original video sequence that you are working with.  The time resolution of when the slices are extracted would depend on the nature of the events captured by the video.  The magnitude and phase response of each 2-dimensional FFT computation could then be determined and examined to see how the Fourier domain response changes across time.

From what I understand, viewing images in the Fourier domain is not readily intuitive.  I believe the Fourier domain is used to realize efficient image filtering.

I hope this helps.


Calculating the system output using frequency response Given an input signal $$x(n)=\cos(6\pi n +\frac{\pi}{6})$$ and system $$y(n)=0.5x(n)-0.1x(n-1)$$. In this case, the coefficients of the difference equation are $a_0=1$, $b_0=0.5$, and $b_1=$. The frequency response to this system is $$H(\omega)= b_0 + b_1e^{-j\omega}= 0.5 - 0.1 \cos(\omega)+0.1 j \sin(\omega)$$

Since the input $x(n)$ is a sinusoid of the form $$x(n)=A \cos(\omega n +\theta_0)$$, we can express the system's output as

$$y(n)=A|H(\omega_0)|\cos\left(\omega_0 n+\theta_0+angle(H(\omega_0))\right)$$

I can see that $A=1$, $\omega_0=6\pi$, and $\theta_0=\frac{\pi}{6}$ in this formula.  Hence, I have deduced that

$$|H(\omega_0)|= |0.5 - 0.1 \cos(6\pi)+0.1j \sin(6\pi)| = 0.4$$

Furthermore, since $H(\omega_0) = 0.4$ (a positive real number), I conjecture that $$angle(H(\omega_0))=0$$

Have I made the correct assumptions/conclusions?

Yes, you are correct.  The system input is rather strange, though, since $cos(6\pi n + \frac{\pi}{6})$ is equal to $cos(\frac{\pi}{6})$.


Phase difference between two sinuisoidal signal Currently, I have two signals, the main components of both signals are 60Hz, but both also have weaker response at 180Hz + small amount of noise. As shown in the photo below, I want to find the phase difference between the 60Hz component of first signal, and the 180Hz component of the second signal.

As you can see the second waveform of the figure, the phase lagg is of the 180Hz signal is 1.39ms, which coresponds to a 90 degree phase shift.

I am wondering, if I want to implment this phase measurement, could I simply take the fft of two signals, and then find the phase at 60Hz of signal 1, and subtract the phase at 180Hz of signal 2?

Or if there is better way to do so?

You can find phase between 2 sines, if you write it as a complex signal where 1st sine is a real part and second one is imaginary. But to do so it's necessary to rescale one component to fit another in its amplitude variation. Or another approach is to normalize both sines in amplitude value. There are some explanations:

We have signals $s_1(t)$ and $s_2(t)$, which are assumed to be sine waves. Before the measurement of their phase differences, we have to normalize there signals. The normalization process is:

$s^{norm}_i(t) = \frac {s_i(t)} {1.414 \cdot \mathrm{E}[\mid s_i \mid]}$

where $\mathrm{E}[\mid s \mid ]$ is mean value of module of $s(t)$. To perform normalization in real time it's necessary to track $\mathrm{E}[\mid s \mid ]$ constantly since the signal $s(t)$ can be non-stationary. Recurrent moving average filter is the case:

$s^{mean}_i(k) = a \cdot \mid s_i(k) \mid + (1-a) \cdot s^{mean}_i(k-1)$

where $a$ is inverse of integration time (or equivalently filter's length $K$) $a = \frac {1} {K}$ and $k$ sample index in discrete time representation.

Once we normalize both signals it is possible to estimate phase shift between them as

$s(t) = s^{norm}_1(t) + j \cdot s^{norm}_2(t)$

$\phi(t) = arg(s(t))$

$\phi(t) = tan^{-1}(\frac {\Im{s(t)}}{\Re{s(t)}})$

If it's real time implementation good practice is to use CORDIC to calculate the phase instead of $atan$ function.


Find minimum sampling rate of non-ideal notch filter Old test problem:

"You have been asked to design a real-time digital filtering system that eliminates a band of frequencies between 20 and 30 MHz but preserves everything else.  The system uses an analog-to-digital converter (ADC) to digitize the analog input, a digital signal processing chip (DSP) to filter the digitized signal, and a digital-to-analog (DAC) converter to convert the signal back to analog.  Assume that ideal sampling is used as well as ideal anti-aliasing and reconstruction filters.  The maximum frequency that needs to be preserved in the reconstructed analog signal is 90 MHz.

What is the lowest sampling rate allowable in this system given the above constraints and giving that the anti-aliasing and reconstruction filters are not quite ideal:  they have 6 MHz transition bands?"

If it was non-ideal it would be the Nyquist equation:

fs > 2*fmax

where fmax = 90 MHz, so,

fs > 2*90MHz

fs > 180 MHz

In the non-ideal situation, my thought would be to add the 6 MHz transition bands to fmax.  Since it is a notch filter there are two transition bands.  Meaning it would be

90 MHz + 6 MHz + 6 MHz = 102 MHz

I was told that the correct set-up to solve this problem is as follows:

fs > 2*(90MHz + (6MHz/2))

fs > 2*(93)

fs > 186

He mentioned something about not caring whether the transition bands overlap.

My question is why is the 6 MHz transition band divided by 2?  Is it because there are two transition bands?  And if it were a highpass or lowpass filter, where there is only one transition band, would it not be necessary to divide by two?

If you sample at 2*93 = 186 MHz, then the frequencies from 93 to 96 MHz will alias into the range from 90 to 93 MHz. However, since the problem requires that only frequencies from 0 to 90 MHz are preserved, you are not concerned with aliasing occuring after 90 MHz.


Implementing Welch's method for Power Spectral Density I want to implement Welch's method for PSD calculation in MATLAB. I do not want to use built in MATLAB cpsd function, such that I can change the FFT implementation in cpsd(for certain purpose).

I already tried to replicate the method based on Welch paper and explanation from the page: Nonparametric Methods. Below is the code that I already implemented with comments about steps to achieve the result. (it is used for accelerometer real data, so numaxis represent number of axis, $xyz$. Each column represent each axis, the time series data is arranged into row):

function psd = pspectra(data1,data2, window, overlap, nfft, fs)

%checking the size
if(size(data1,2) ~= size(data2,2))
    error('size of column is not equal');
elseif(size(data1,1) ~= size(data2,1))
    error('size of row is not equal');
else
    numaxis = size(data1,2);
end

%calculating window sliding step for iteration
winsize = length(window);
step = winsize - overlap;
iter = 1 + (size(data1,1) - winsize)/step;

%start and end index of first window/segment
istart = 1;
iend = istart + winsize - 1;

%start calculating fft for each window
for i=1:iter
    for j=1:numaxis
        %apply window
        data1(istart:iend,j) = data1(istart:iend,j).*window;
        data2(istart:iend,j) = data2(istart:iend,j).*window;

        %calculate fft
        fft1(:,j,i) = fft(data1(istart:iend,j),nfft);
        fft2(:,j,i) = fft(data2(istart:iend,j),nfft);
    end

    %move to next window segment
    istart = istart + step;
    iend = iend + step;
end

%obtain scale to create modified periodogram
scale = 1/(fs.*sum(window.*window));

%averaging window result and apply the scaling
psd = zeros((nfft/2)+1,size(fft1,2));

for i=1:iter
    psd = psd + fft1(1:(nfft/2)+1,:,i).*conj(fft2(1:(nfft/2)+1,:,i));
end

psd = psd.*scale./iter;

%multiply by 2 except dc and nyquist component (1 and 51)
psd(2:50,:) = psd(2:50,:).*2;

end

I have tried using multiple segment with non overlapping window/segment, and got exactly the same result as cpsd() function from MATLAB. For example, using this call:

pspectra(data1,data1,hamming(50),0,100,100);
pspectra(data1,data1,hamming(25),0,100,100);

However, when I tried overlapping window like this

pspectra(data1,data1,hamming(50),25,100,100);
pspectra(data1,data1,hamming(60),20,100,100);

The resulted power spectra is totally different. I have tried to find if my implementation is wrong (in theory or code), and unfortunately no luck yet. Do I miss something in here?

I really appreciate for any help.

I believe these lines are you problem:

    data1(istart:iend,j) = data1(istart:iend,j).*window;
    data2(istart:iend,j) = data2(istart:iend,j).*window;

you are modifying the incoming data.  When you don't have overlapping windows, you don't notice it. However, when you have overlapping windows, the next window will include data that has already had a window applied... which will not be what you want.


Estimate camera generalized velocity from consecutive frames Let's suppose I have a camera which can move freely in 3D space and has 6 degrees of freedom. Is it possible to estimate the linear and angular velocity of the camera from 2 or more consecutive frames? If yes, are there any environment or camera constraints that must hold?

It's possible, but it takes some sophisticated processing to do it effectively in a general environment.

First of all, you can only detect motion relative to what you can see in the field of view. This requires good foreground/background separation to eliminate objects that might themselves be moving relative to the background (which is presumably what you want to measure the camera motion against).

Second, the background must have enough features (edges) in it to make all six degrees of freedom unambiguously observable.

It boils down to figuring out what the trajectories through the viewspace are of all the features in the scene, and then doing a cluster analysis of those trajectories in order to assign them to distinct objects. Finally, pick the cluster that most probably and most consistently represents the background. The inverse of this trajectory is the camera's trajectory.

Once you have the geometry of the situation squared away, you still need information about the size of objects in the field of view in order to scale the results correctly.


Is it possible that the value of a continuous mother wavelet at origin is zero, i.e. $\psi(t=0)=0$? According to Fourier transform, a continuous wavelet could be written as
$$
\psi(t)=\frac{1}{2\pi}\int\hat\psi(k)\text{e}^{-ikt}\text{d}k
$$
From the equation above, we know that $\psi(t=0)$ is
$$
\psi(0)=\frac{1}{2\pi}\int\hat\psi(k)\text{d}k
$$
Is it possible that  $\psi(0)=0$?

Yes: a sufficient criterion for a valid wavelet is being zero mean in time domain, i.e. $\hat \psi (0) = 0$. For $\psi(0) = 0$ we require that $\hat \psi$ sums to zero; one example is higher order Generalized Morse Wavelets - from Olhede & Walden,

(Though above is only approximately zero-sum). Such wavelets can be admissible and analytic, enabling CWT inversion. Trivial Python example below; GMWs are also implemented in ssqueezepy.

Code
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import ifft, ifftshift

t = np.linspace(0, 1, 256, endpoint=False)
wf = np.exp(-(t - .08)**2 * 4096)
wf -= np.roll(wf, 25)    # shift by 25 samples and subtract
w = ifftshift(ifft(wf))  # take to time domain and center

plt.plot(wf)
plt.show()

plt.plot(w.real)
plt.plot(w.imag)
plt.show()


Significance of modular arithmetic in DFT? In what ways does modular arithmetic plays a part in DFT?
Why is it a so integral part of DFT?

This answer to a differently worded question answers your question also.

And this answer is relevant.

And this answer and this answer.

I confess that I take great exception to anyone who denies that the Discrete Fourier Transform (and inverse DFT) fundamentally periodically extend the finite set of data presented to it.  I'm a nazi about it.  The DFT assumes that the $N$ samples of $x[n]$ passed to it are one period of a periodic sequence where $x[n+N]=x[n]$ for all integer $n$. And the inverse DFT does the same in that the $N$ samples of $X[k]$ passed to it are one period of a periodic sequence where $X[k+N]=X[k]$ for all integer $k$.

The modulo arithmetic done to the indices of $x[n]$ and $X[k]$ simply guarantee this periodicity while also guaranteeing that neither index, $n$ or $k$, go outside of the bounds of the finite sized arrays in memory representing $x[n]$ and $X[k]$.  That is

$$ 0 \le n = \operatorname{mod}_N(n+mN) < N $$
and
$$ 0 \le k = \operatorname{mod}_N(k+mN) < N $$

for any integer $m$.


EEG signal processing with wavelet or fft? I have confusion about signal processing related with EEG signal. I have done some of my research and that made me more confused about processing and filtering the signal.
Let me jump into the problem -

I have found that people are using FFT for EEG signal processing which I don't understand. Why would you do that? FFT is mainly made for stationary waves and we know that EEG signals are non-stationary waves hence FFT is not so useful for EEG signal processing. Well that can be solved with the SFFT but then again, this can't be done with real-time signal EEG. (I might be wrong, please help)

Wavelet transform is another way to process EEG signals because, first, it preserves time and freq whereas FFT looses time resolution. Also, wavelet can be applied on non-stationary signals. (I might be wrong here too, so please go ahead help me with this too)

As in the figure below - Neurosky shows real-time EEG power bands fluctuations, which I believe is FFT on raw signal and some kind of mathematical operation (maybe frequency averaging) applied on a range of frequencies gives those bands (alpha,beta,gamma,delta,theta), Am I right? if yes, then how can you apply FFT on real time signal? if no, then what is the best way to get those bands?

I have muse, emotiv, neurosky and openbci hardware which I borrowed from my community friends, I have been playing around with those hardware so I would do something with the raw signal so I started learning about these things but as deep as I go in that rabbit hole, more I get confuse. I was gonna apply ML eventually but I seriously got stuck on my first step.

Here are links to the coding.  I am doing something else, so I have not dug into the code to see how the algorithms work.  If you need help analyzing the code set up a private "room" here and we can discuss it there.  I am sure that it requires more chit-chat than the hall monitors like here :)
a simple Processing application
ARDUINO BRAIN LIBRARY
Your picture
Thanks for the lead!  I had an application for treating certain dementia problems using something like this but haven't followed up.


Matlab FFT keeps plotting mirrored values My Matlab FFT plot keeps showing the mirrored values of my frequency of interest. In this example I am looking for 250Hz with a 3200Hz sampling freq.  The FFT plots my value at 266Hz as well as 1330 Hz.  I assume my sampling freq is a little off which is why it isn't exactly at 250Hz.  Matlab code is below, any help is appreciated.

% Tmax = (N-1)*Ts
% largest freq = Fs/2
% smalles = -fs/2

fid = fopen('250Hz.bin');
AA = fread(fid, 'short');
length(AA);
flenbytes = ftell(fid);
flenshorts = floor(flenbytes/2);
FF = floor(flenshorts/16)*16;
fseek(fid,0, 'bof');
BB = fread(fid, [8 (FF/8)], 'short');
CC = BB';
DataSet = CC ./2^3;
Z1Axis = DataSet(:,4);

N = length(DataSet)
Fs = 3200;
ts = 1/Fs;
tmax = (N-1)*ts;
t = 0:ts:tmax;  %start at 0 increment by ts up to tmax
plot(t,Z1Axis)

f = (0:N-1)*(Fs/N)/2; % start with -Fs/2, increment by Fs/(N-1), goes to Fs/2
z = fft(Z1Axis);
%zreal =
%plot(f,z)
plot(f,abs(z));

It is a mathematical property of the DFT that its output is conjugate symmetric if its input is real-valued. That means that for the magnitude of the DFT output signal $X_k$ of length $N$ it holds
$$
|X_{N-k}| = |X_k|.
$$
So what you see is to be expected. (In your example $3200\,\mathrm{Hz}/2 - 266\,\mathrm{Hz}\approx1330\,\mathrm{Hz})$ If you're only interested in the positive frequencies, just plot the first $N/2$ elements of the DFT output, the rest is redundant.


Adding/subtracting sinusoids I'm trying to calculate resultant function from adding two sinusoids:

$9\sin(\omega t + \tfrac{\pi}{3})$ and $-7\sin(\omega t - \tfrac{3\pi}{8})$

The correct answer is $14.38\sin(\omega t + 1.444)$, but I get $14.38\sin(\omega t + 2.745)$.

My calculations are (first using cosine rule to obtain resultant $v$ as):
$\sqrt{9^2 + (-7)^2 - (2 \cdot 9 \cdot (-7) \cdot \cos(\pi - \tfrac{\pi}{3} + \tfrac{3\pi}{8}))} = 14.38$

And the angle (using the sine rule):
$\pi - \arcsin(|-7| \sin(\pi - \tfrac{\pi}{3} + \tfrac{3\pi}{8}) / 14.38) = 157$ ° or $2.745$ radians.

The easiest way (to my mind) to solve the problem is to

Use the identity $\sin(A\pm B) = \sin A \cos B \pm \cos A \sin B$, substituting the known numerical values of $\cos B$ and $\sin B$,
Gathering the results to express your sum of sinusoids in the form of $C \sin A + D \cos A$,
Expressing the resulting function as $\sqrt{C^2+D^2} \sin\left(\omega t + \theta\right)$


Fast algorithm for the MATLAB's $\tt impz$ function? Does anybody here know what is the code algorithm behind the impz function?
I know it uses the recursive algorithm but how is such implemented in impz that it can get the impulse response of very high $N$ and $M$ orders (numerator and denominator coefficients)?

I can make a script that gets the impulse response of a second order LCCDE (linear constant coefficient differential equation) but I want to know a GENERAL code that can get the impulse response of any LCCDEs regardless of their order.

the impz function of matlab, simply computes N samples of the impulse response of a given LTI system from its LCCDE description coefficients a nd b.
As far as I know its efficiency is a secondary issue, since computing an impulse response in real time is not a frequently encountered task.

One equivalent code is simply this:

  % given a,b vectors and number of samples as N
  x = [1 zeros(1,N-1)];
  h = filter(b,a,x);

I guess you know how the filter function is implemented and what its output means. In the simplest terms, the filter() matlab function computes the output vector y for an input vector x of an LTI system: $ \sum_{k=0}^{k=N}a_k y[n-k] = \sum_{k=0}^{k=M} b_k x[n-k]$ , by the recursive approach:

$$ y[n] = - \sum_{k=0}^{k=N} \frac{a_k}{a_0} y[n-k] + \sum_{k=0}^{k=M} \frac{b_k}{a_0} x[n-k] $$

Where the initial conditions of the system is assumed to be zero, inline with impulse response definition and LTI system with initial rest conditions.

Note that in the impulse response computation, the signal $x[n]$ has only a single nonzero entry at $n=0$ and the rest of the computation is a recursion on the signal $y[n]$ alone.


Scaling Function and Wavelet Function I don't understand the meaning of scaling function when I used CWT/DWT.
For me, I understand that the wavelet function (mother wavelet) is the function that we used to decompose our signal in subsignal corresponding to the frequency of the mother wavelet (for example db8). (Maybe DWT is different from wavelet decomposition ?)
For example, if I decompose my signal into 5 levels (with wavedec in MATLAB), I have :
D1, D2, D3, D4, D5 and A5, with D the detail's signal and A the approximation's signal.
But here I never used the scaling function to make the decomposition of my signal.
Can you help me on how to understand the scaling function ?

In my experience, scaling functions are not so central in the CWT. However, there are quite important on the DWT. Scaling functions act like low-pass operators, and wavelets like band-pass ones. Indeed, in the DWT, one seldom uses th wavelet and the scaling function directly, since they are hidden in the construction of the low-pass and high-pass filters that are combined in cascading filterbanks. The actual wavelet (and scaling function) are only obtained as a limiting process involving an infinite cascade of discrete filters.
A complementary discussion: Scaling Function in Wavelet Analysis


Method of reconstructing a band-limited signal from discrete samples I understand current cell phones use digital communications. Given that this industry brings in billions of dollars each year, there is much incentive to get the best performance possible. So what method of reconstructing an audio signal from discrete samples is normally used in a good cell phone these days? It seems there must be methods that work better than the ones I read about in text/reference books.

Some recent cell phone models use something like a Cirrus Logic CS42xx series audio IO chip, which seems to use a digital polyphase interpolation filter, a sigma delta modulator, followed by a switched capacitor DAC and low-pass filter.

Sinc interpolation (or, given finite hardware, a polyphase FIR kernel similar to a windowed Sinc) is one high quality method of reconstruction from bandlimited samples.

But reconstruction isn't the main issue with respect to audio quality.  Quality or fidelity is more likely constrained by the audio data compression algorithm used (or required by standardization), and how much spectrum bandwidth is allocated by the comunications standard and regulations.


How to generate random data with a specific bandwidth I generate some random data (bits, actually, using MATLAB's randi) and then I convert such bits to BPSK (1 and -1). After that, I use a RRC pulse (roll-off factor of 0.35, 8 symbols and 4 samples per symbol) to create the waveform.

The bandwidth of the pulse-shaped BPSK then would be ($\beta$ is the roll-off factor):

$$
BW = R_S(1+\beta)
$$

My problem here is how to determine $R_s$ (data rate). I mean, I'm using randi, not an actual data source so initially I tried to make that value up. I suppose $R_s=1MHz$, after pulse shape the sampling frequency is $f_s = sps*R_S = 4MHz$ ($sps=4$ because of the 4 samples per symbol stated above), then the bandwidth should be something like $BW=(1+0.35)1MHz=1.35MHz$ and tried to plot the DFT of the pulse-shaped sequence for a final sampling rate of $4 MHz$:

As you can see, the bandwidth is nowhere near $1.35 MHz$.

So these are my questions:

How can I manipulate/control $R_S$
What's wrong with my calculations?

Thanks!

In the case of BPSK, the data rate equals the symbol rate (one bit per symbol). One way to control $R_S$ is to change the modulation. For example, you can increase the number of bits per symbol by choosing a higher order modulation like 16-QAM which has 4 bits per symbol.
There is nothing wrong with your calculation. The plot is a single sided spectrum so you only see half the bandwidth. It roughly dies off around half of 1.35 MHz so assuming your code is correct, I see nothing particularly wrong with the plot or calculation.


Find a small image pattern in another similar image I have this pattern in the first image which I want to find in the second image. However the two images are totally different so I just need to find the closest set of indices where these two images match. So how to go about it? Are there existing algorithms doing the same?

The simplest algorithm that is mathematically tractable is the matched filter. It is designed to find both the presence of (yes/no) and the location on a "target" in a "scene". However, the matched filter does not work incredibly well if there are rotations or scaling of the "target" in the scene relative to the target's image.

Assuming the "target" is in the "scene", the matched filter gives incredibly good results if the only thing corrupting the "scene" is noise. To deal with rotations, one can try "matching" rotated and scaled versions of the target image in the scene. There are many variations of the matched filter for many different signal processing applications, but in image processing, it is sometimes called "cross-correlation".

The modern approach to this "pattern recognition" problem is machine learning, using Neural Networks or Support Vector Machines. These concepts are much more difficult to learn.

MATLAB has a lot of useful tools for accomplishing these tasks, and their online documentation has a lot of examples.

A variation of the matched filter that is more robust to different types of noise existing in the same image is the Normalized Cross-Correlation (also a function in MATLAB, normxcorr2 as opposed to xcorr2, which is a classic matched filter). Theoretically, it is described here by Lewis...

http://s3.amazonaws.com/academia.edu.documents/3607630/10.1.1.21.6062.pdf?AWSAccessKeyId=AKIAIWOWYYGZ2Y53UL3A&Expires=1488301575&Signature=HpDWFoqWiZAQH7VPkr19hCbUrYQ%3D&response-content-disposition=inline%3B%20filename%3DFast_normalized_cross-correlation.pdf

But, although norm x-corr can deal with different types of noise, it is still useless under scaling and rotation, unless one repeats the function over and over for different combinations of rotations and scales. Even then, it can be difficult to find the target, unless you get very granular with the rotations and scales.

Hope this helps. I can describe more of these methods in detail if needed.

Note that although the theory of the matched filter might seem extreme, the implementation is very easy, especially if you just use xcorr2 or normxcorr2 in MATLAB.


How audio watermarks are captured I have read that we can use audio watermarks to "identify the original produces of the material or the consumers of the media."

How consumers of the media is measured by audio watermarks?

Lets say in TV we play an video which has audio watermark in it and 100's of people are watching it. How do we capture the audio watermark that got played in viewers TV?

In the U.S., a company called Nielsen has its watermarking technology encoded into broadcast content (Radio/TV/etc), while some 65k individuals called "panelists" wear a small pager-like meter called a portable people meter (PPM) which "listens" for and reports the watermarks that it "hears" (as the individuals go about their daily life). The "panelists" are essentially a statistically significant sample of the larger population of Radio/TV consumers.

  Nielsen’s watermarking involves codes, and each PPM transmits both
  complete and partial codes back to Nielsen. In situations where
  partial codes exist, Nielsen has edit rules that “fill in the blanks”
  by validating against nearby complete codes. Broadcasters are
  currently updating their station’s encoders with the Enhanced CBET (Critical Band Encoding Technology)
  technology. Prior testing and preliminary results from the roll-out
  show that there is an increase in the quality of the code detection
  and an increase in average quarter hour audience estimates.

Here are some links for you:

http://www.nielsen.com/us/en/press-room/2016/nielsen-to-increase-portable-people-meter-sample-size-by-10-percent-across-48-radio-metro-areas.html
http://sites.nielsen.com/newscenter/making-noise-nielsen-audio-gets-into-the-spirit-of-total-audience/


Remove watermark from video Is there a process to remove a watermark ever present across all frames of a video?
I was thinking that comparing frames and taking pairs poorly correlated, could give a clue about the watermark both share. Or perhaps merging or adding all the frames could reveal the strong watermark pattern. O maybe is a work for a neural network.
The presence of this fixed pattern between other variable information make me think is feasible. I will be glad to read your insights about this problem.

If the watermark is additive and static and, the content is smooth, and the vide long enough you could compute the watermark that minimizes $\sum ||S(F_i - W)||^2$, for the frames $F_i$ and a fixed watermark $W$. The solution $W$ is expected to converge to the watermark.
But it is always better to use the original video.


Problems computing the DFT of finite length sequence I am having trouble finding the same answer as the solution manual for this sequence.
The problem asks to compute the DFT of
$$
x[n] = \begin{cases}
            1 & \text{for even } n \in \{0\ldots N-1\}
         \\ 0 & \text{for odd } n \in \{0\ldots N-1\}
       \end{cases}
$$
Setting up the DFT equation I agree with the solution manual
$$
  X[k]=\sum_0^{N-1}x[n]W^{kn}_{N}=\sum_0^{\frac{N}{2}-1}W^{k2n}_{N}=\sum_0^{\frac{N}{2}-1}e^{\frac{-j2\pi k(2n)}{N}}
$$
However when apply the closed form geometric series formula $\sum_0^Na^k=\frac{1-a^{N+1}}{1-a}$ I get
$$
\frac{1-e^{-j2\pi k}}{1-e^{\frac{-j2\pi k}{N}}}
$$
while the solution manual says its
$$\frac{1-e^{-j2\pi k}}{1-e^{\frac{-j\pi k}{N}}}$$
Note their lower exponent does not have a 2. Could someone explain where and how the 2 goes away?  Or does someone agree that the provided answer is incorrect?
This is question 8.5c from Oppenheim's Discrete-Time Signal Processing 3e.


How do you measure the accuracy of a stereo system? I have a philosophical question of sorts for users familiar with a multi-view stereo setup. Imagine that I need to measure the dimensions of a large building and I take several photos of the building from different viewpoints. After calibrating the camera (intrinsic and extrinsic), I compute matches between points in these images and map them to 3D points. The output of this entire process is a cloud of 3D points. If I specify a scale, origin and co-ordinate axes, then I can perform measurements between points in this cloud which relate to the real-world scene (the building) that I captured. My question is this: how do I determine how accurate was my reconstruction.

Basically, this camera apparatus forms a measurement device. To compare a measurement, one typically needs to use another measurement device that is much more accurate that it. How do I measure the accuracy of a such a system?

From a more pragmatic than philosophical perspective, possibilities include:

Comparing with ground truth data--these data could come from surveying's instruments, or be assumed to be "reliable" (e.g. you could assess the system's accuracy on the Eiffel tower and other well-known buildings for which measurements have been curated);

Verifying physical properties are met (e.g., floors are horizontal, walls are perpendicular, distances between floors or between windows are constant);

Leaving one ore more views out of the 3D estimation process and assessing how far the projections of the estimated 3D points are from the actual 2D positions in these views (ie., jackknifing).


Reconstructing audio from a given spectrogram and phase data I have 2 seperate files called: spectrogram and phases. Spectrogram contains magnitude information of audio and phases contain the phase information of each windowed area. I need to reconstruct the audio using these information yet i am completely lost. Can you put me into the right direction. How can i do this?

A spectrogram is a graphical time-scroll representation of the frequency components in a signal. Sometimes it's called a "waterfall", because that's what it looks like as it is computed and plotted in real time. The core math behind this is the FFT.
A spectrogram generally shows the strength (or magnitude) of a range of frequency components in a signal (e.g. DC to 20kHz). The FFT has a complex output, one block of values for each block of input. The magnitude of these complex numbers is likely (somehow) stored in the "spectrogram" file, and the phase is stored in the "phases" file.
The trick to reconstructing the original signal from a spectrogram is to inverse the spectrum, and the math behind that is the IFFT.
You will need the magnitude and the phase data from the files to recreate the complex spectrum and then apply it to the IFFT to get the time domain signal back.
You will have to find out

what math was applied to calculate the spectrogram: FFT block size, any log-scale (dB?), any windowing, any sub-sampling etc... and how the data is formatted and organized in the file.
how the phase is calculated and stored (block size, degrees/radians, number format etc...)

As a commenter pointed out, you will have to find out all these details from whoever created the file. But beware that even then it might not be possible to recreate the original signal.
In that case, if you can specify how to generate the spectrogram with different parameters (no windowing, no sub-sampling etc...), and then (re)generate these files, you may stand a chance.
Good luck. Since this is a bit like forensics, you will need it.


What is the difference between Linear Interpolation factor and Sampling rate conversion factor? I have come across sampling rate conversion factor, which is given by:

                             S_factor = F_new/F_old

If F_new > F_old, then s_factor > 1

and

If F_new < F_old, then s_factor < 1

But, what is this "linear interpolation factor"?. Is it the same as sampling rate conversion factor?

Does the linear interpolation factors value should be in the range 0 to 1?, If so, why is that?

I'll move my comments to an answer:

With sample rate conversion factor $S$, you want to compute samples at "time" $\frac{i}{S}$ for $i=0,1,2,\dots$. "Time" relating to the index of the previous series. The linear interpolation factor is simply then the fractional part of $\frac{i}{S}$ for each $i$ with the integer part telling which sample index to use (together the next) for the linear interpolation.

Simple example to clarify: upsample a sequence a factor $\frac{10}{9}$. You would like to determine the samples at "time" $0, 0.9, 1.8, 2.7, 3.6, \dots$.  First, start with sample $0$, then, compute the linear interpolation between sample $0$ and $1$ with $\eta = 0.9$. Next sample should be at "time" $1.8$, so $\eta = 0.8$ for samples $1$ and $2$. "Time" $2.7$ gives $\eta = 0.7$ for samples $2$ and $3$ and so on.

Note that you conceptually can compute the "time" $1.8$ sample from sample $0$ and $1$ with $\eta = 1.8$, but using sample $1$ and $2$ with $\eta = 0.8$ is better. So while $\eta$ does not really have to be $0$ and $1$ the results will be better and there is no advantage to not have it in that range.


Fourier Transform of Impulse Train Why is the fourier transform of impulse train a impulse train? Is there a intuitive reason behind it?

Intuition can sometimes be misleading.  But here are some ideas that might help one move towards creating a mental picture.

An infinitely long pure sinewave in the time domain (consisting of just one frequency FT or DFT basis function) will be a single impulse in the frequency domain.

Distort the sinewave a little, but leave the waveform perfectly periodic, and the impulse will be followed by an evenly spaced harmonic series.  Usually the narrower and sharper the distortion (but keeping the waveform still perfectly periodic), the longer the harmonic series.  What might be considered a limiting case of maximum distortion, the narrowest waveform with the sharpest edge will have the longest harmonic series.  Or an infinitely long sine wave in the time domain maximally distorted into just an infinitely periodic impulse train, will produce a impulse followed by an infinitely long harmonic series, which looks a lot like another periodic impulse train.

Make it an even (cosine) function, and all the impulses in the FT will be real and thus symmetric around 0.  Add a DC offset to the distorted sine wave to complete the pulse train at 0.


What exactly happens with Spread Spectrum Modulation With regard to DSSS and CDMA(DS), what does it mean to spread a signal over a spectrum of frequencies in terms of individual frequencies? Is it about sending the same signal redundantly on multiple frequencies, or is the data stream actually divided and sent as multiple discreet signals over various frequencies? I've read explanations claiming either.

Direct-sequence spread spectrum (DSSS) is a technique that is used to generate a modulated signal that occupies more bandwidth than would be implied by its information content alone. A DSSS signal is generated via multiplication of a (typically digitally-modulated) baseband signal by another spreading code waveform. The spreading code waveform is constructed such that it is modulated at a higher rate than the baseband signal; this yields the "spreading" of the resulting signal's spectrum.

To see this, recall the multiplication property of the Fourier transform. Let the baseband signal $x(t)$ be multiplied by a spreading waveform $s(t)$:

$$
x(t)s(t)\Longleftrightarrow X(f) * S(f)
$$

where the $*$ denotes convolution. If $S(f)$ is broadband when compared to $X(f)$, then the result of convolving $X(f)$ with $S(f)$ is that $X(f)$'s spectral content is "smeared" or "spread" across a wider region of the spectrum. The amount of spreading is dependent upon $S(f)$ and therefore the choice of the spreading sequence $s(t)$. So, for example, you could take a 1 kbps BPSK signal, which would normally have a null-to-null bandwidth of 2 kHz, and spread it to instead have a null-to-null bandwidth of your choice: say, 1 MHz.

To answer your most pointed question, you're not really sending a single signal redundantly on multiple carrier frequencies, nor are you dividing the input signal into pieces that get transmitted at multiple carriers. In lieu of a discrete multi-carrier approach, DSSS instead is more of a continuous smearing of the input signal over a large region of spectrum.

If the receiver is also cognizant of the spreading sequence $s(t)$, then it is capable of despreading its observation of the transmission and therefore recovering the original baseband signal. As long as the spreading sequence $s(t)$ does equal zero for any $t$, then in principle, the spreading process is perfectly reversible. In practice, this isn't exactly the case, as there will be some distortion between what is transmitted and what the receiver sees, but it leads to an important point that is often misunderstood about DSSS:

  No benefit in bit error performance, measured in $\frac{E_b}{N_0}$ (the typical metric used for digital modulation performance evaluation) is offered by using direct-sequence spread spectrum. Since the spreading process is reversible, then the principle of reversibility states that there can be no improvement in performance of the optimal receiver for the spread signal when compared to the performance of the optimal receiver for the original, unspread signal.

So if there's no direct bit-error performance benefit, why would you want to do this? There are a few common reasons:

Multiple access: Code division multiple access schemes work using DSSS techniques. CDMA allows multiple transmitted signals to occupy the same bandwidth over the same time period. Transmitters are associated with separate carefully-selected spreading codes; these are also known by the receiver, which uses them to separate the previously-intermingled-in-time-and-frequency signals. CDMA can be made to seem like the greatest thing since sliced bread, but there are some practical issues (most notably the near-far problem) that limit its usefulness. It is, however, used in some cellular networks around the world, for one example.
Interference mitigation: Signals with high bandwidth are resistant to narrowband interferers on the AWGN channel. A hand-waving time-domain explanation: a wideband signal in the frequency domain will have a narrow autocorrelation function in the time domain, and vice versa. The optimal receiver will use a correlator (often implemented with a matched filter) matched to the transmitted spread-spectrum waveform, so its output (which ideally will contain copies of the signal's autocorrelation function at symbol-spaced intervals) will contain narrow, sharp peaks that will be easy to identify against the oppositely-characterized background interference.

A hand-waving frequency-domain explanation: the spreading process takes the original baseband signal's power and spreads it across a wide region of the spectrum. In the communication channel, narrowband interference finds its way into the signal that eventually is observed. During the despreading process, the DSSS signal is compressed back to the original baseband signal's bandwidth. However, at the same time, the despreading process actually spreads the narrowband interference across the same bandwidth as the original DSSS signal.

So effectively, the narrowband interferer has been transformed into a wideband interferer across the entire DSSS signal bandwidth, which as you recall, is much larger than the baseband (information-carrying) signal's bandwidth. Therefore, much of the interference power is pushed out of band with respect to the signal of interest, which yields a higher signal-to-interference power ratio (and better performance).
"Discreetness": No, not discrete, discreet. DSSS is often spoken of as being able to "pull signal out of the noise floor." This refers to the ability of a spread-spectrum receiver to transform a wideband DSSS signal back to the original, much narrower waveform. For a given transmitter power, as the transmitted spread signal bandwidth increases, the average power spectral density over the DSSS bandwidth will decrease, perhaps to the point where it is difficult to observe using spectral analysis at the receiver. This property can be useful for a number of reasons.

One aim might be to decrease the probability of the emitter being detected (and possibly intercepted) by an adversary. Another might be to mitigate interference with other, narrowerband systems that operate in the same region of spectrum. An example of this approach is in ultra wideband communications systems, which have increased in popularity in recent years.


Unilateral Z transform I tried to calculate the unilateral Z transform of x[n-2], is it right?

There is nothing wrong with your calculation. However, by definition, $x[n]=0$ for $n<0$, because otherwise we couldn't have used the unilateral $\mathcal{Z}$-transform in the first place. So you get

$$\mathcal{Z}\big\{x[n-k]\big\}=x^{-k}X(z),\quad k> 0\tag{1}$$

Eq. $(1)$ is valid for the unilateral as well as for the bilateral $\mathcal{Z}$-transform.

Things are different for a time advance. The unilateral $\mathcal{Z}$-transform of $x[n+k]$, $k>0$, is

$$\begin{align}\mathcal{Z}\big\{x[n+k]\big\}&=\sum_{n=0}^{\infty}x[n+k]z^{-n}\\&=\sum_{n=k}^{\infty}x[n]z^{-(n-k)}\\&=z^k\left(X(z)-\sum_{n=0}^{k-1}x[n]z^{-n}\right),\quad k>0\tag{2}\end{align}$$

For the bilateral $\mathcal{Z}$-transform we simply have

$$\mathcal{Z}\big\{x[n+k]\big\}=z^kX(z)\tag{3}$$


Reducing the number of data points of an FFT I'm in an odd situation where an FFT is giving me too much data to process. I'm performing an FFT on audio sampled at 500 Hz, over 30 seconds. The result is 7500 bins of useful information, however I don't need this sort of resolution. I still need to see frequencies up to 250 Hz so I can't reduce the sampling rate. I also need information from the entire sample, so I can't use less than 30 seconds of data either.

What can I do to reduce the number of points, but keep the magnitude and phase information relatively intact? I would prefer not to throw away points, is it possible to take an average of x points?

Do you want to know specific frequencies? Go for the Goertzel algorithm.

Do you need to know every 2nd frequency? This corresponds to downsampling in the frequency domain, which is aliasing in time domain. You can do it like this (where I have reduced the sampling rate and signal duration for better plotting, obviously it also works with 30s and 500Hz):

Fs = 50 # 50 Hz for less data
x = np.random.randn(2*Fs)  # The audio signal: 2s for less data for better plotting here
L = len(x)

# Calculate full FFT for reference
X = np.fft.fft(x)
f1 = np.linspace(0, Fs, L, endpoint=False)

# Calculate every 2nd sample of FFT
x2 = x[:L/2] + x[L/2:] # Perform the aliasing operation in time domain
X2 = np.fft.fft(x2)
f2 = np.linspace(0, Fs, L/2, endpoint=False)

plt.plot(f2, abs(X2), 'go-')
plt.plot(f1, abs(X), 'rx-')

plt.xlim((0, Fs/2))

Obviously, the approach of aliasing also works for every Nth sample in general, if $L$ can be divided by $N$ without remainder.


What is the meaning of $l_p$ norm in this model for sparse channel estimation? The signal $y(n)$ can be modeled in the time-domain as

$$y(n) = \sum_{i=0}^{L-1}h(i) x(n-i) + v(t) \tag{1}$$

where $\mathbf{h} = [h_1,h_2,\ldots,h_L]^T$ is the coefficient of length $L-1$, and $v(t)$ is the additive White Gaussian noise with variance $\sigma^2_v$.

In the Eq(20) found in the paper, http://www.eurasip.org/Proceedings/Eusipco/Eusipco2008/papers/1569101936.pdf

the eq has the term ${||h||}^p_p$ and the Authors mention that here $0<p \le 1$.

Consider an array $h = [1,0.2,0,0.5]$. What is meant by ${||h||}^p_p$ Is it the difference between the elements of $h$? What would be the meaning when p =0, p=1, and p=2?

The Authors in Eq(24)present the derivative of the cost function computed in Eq(20). The derivative is taken with respect to $h$ parameter. In that, they use the term $p\lambda \tilde{h}$ where $\tilde{h}$ is the sparse coefficient. I think here $p$ may not be zero.

The symbol ${||h||}_p$ refers to the p-norm of the vector. It is calculated as

$${||h||}_p=\sqrt[p]{\left|h_1\right|^p+\left|h_2\right|^p+...+\left|h_n\right|^p}$$

When they write ${||h||}^p_p$, it is just that multiplied by itself $p$ times, so that the root disappears. Mathematically:

$${||h||}^p_p=\sum\limits_{i=1}^n \left|h_i\right|^p$$


GSP as an extenstion of DSP I am a PhD. in pure mathematics.

Could you please illustrate the following statement: the eigenvectors of a graph Laplacian behave similarly to a Fourier basis, motivating the development of graph-based Fourier analysis theory.
I am reading the interesting  paper, but could not get how the  Fourier transform is extended to the graph Fourier transform as illustrated on page 23.  Indeed,  why should the matrix V be considered as the extension of (discrete) Fourier transform?

This is a valid question within the context of GSP although, it might not sound like much of a digital signal processing question.

To understand that point, you need to understand what the eigenvalues/eigenvectors tell us (about a matrix), then how is the Laplacian constructed (and what it means) and then relate it to DSP.

The Eigendecomposition of a matrix is a way of expressing the matrix as a product of eigenvalues and eigenvectors. To understand their geometrical interpretation, imagine that each row of the matrix describes a vector in some n-dimensional space. A vector has a beginning and an end and it points towards a direction. The eigendecomposition of a matrix decomposes the "most common" direction that these vectors are pointing towards as a sum of some elementary vectors (eigenvectors) multiplied by their strengths (eigenvalues).

If you were to fill a matrix with vectors that are pointing towards the same direction and asked for its eigendecomposition, you would get back 1 eigenvector at maximum strength. That eigenvector points towards the direction that all vectors in the matrix point towards. Try it here

If you were to fill a matrix with two subsets of vectors where one subset points in some direction and the other points to some other direction, you will get a different mixture.Try it here. You can try different values and explore this notion, I am just providing some simple illustrative examples here. For more information on how this is applied to image/signal processing more generally, you might want to see here.

Moving on to the Laplacian: The Laplacian, in Graph Theory is a very important matrix because it is the bridge to Spectral Graph Theory. The Laplacian's definition is of course one but there are different shortcuts depending on the type of graph you are dealing with. For example, for simple undirected graphs, you don't really have to evaluate the difference between the Degree and Adjacency matrices since that would result in a -1 when two nodes are adjacent and a 0 if they were not and the main diagonal filled with the degree of each node.

This is the construction. And the point is that the Laplacian can be used to broadly characterise graphs. It becomes very useful when you are dealing with graphs that are intractably (for a human brain) big. But, informally, the Laplacian is giving us a "proper" way (from the point of view of functional mathematics) to encode the connectivity of the graph and work with it algebraicaly.

If you couple this with the eigendecomposition, you can see how the pattern of node connectivity (how nodes are connected with each other) can be decomposed into sub-structures via the eigendecomposition we described earlier.

This is where the connection is made between the Fourier Transform, that is used to decompose a waveform into a sum of sinusoidal components and a graph's spectral decomposition, that is used to decompose the pattern of connectivity into a sum of sub-patterns. (All this, under specific constraints of course, so that these things do not break).

The last thing that remains now is to link all this with digital signal processing. DSP has its roots in Time Series analysis. Graphs are mathematical objects that capture relationships, how are these two things combined?

The answer to this is in the way you can map a signal on to a graph or, knowing that the signal originates from a process of which the graph is known, make some inferences about the signal. This latter case is very topical in biomedical signal processing and specifically brain connectivity, where you can infer which areas of the brain exchange information to achieve a certain cognitive task. (Those colourful matrices are basically weight matrices of graphs but you infer them through a series of correlations between various waveforms).

In the presentation (paper (?)) that you cite, this is outlined around slides 6 to 8, where some basic examples are given on trying to map certain processes on networks, then process the networks and conclude to something about the signal itself. And then slides 52 onwards provide some more concrete practical examples.

Hope this helps


How do I pick an optimal window size for overlap add/save method algorithmically Currently I've implemented an overlap-add method using resources on the internet, but I couldn't find a well documented way to minimize the cost of the method. In other words, how do I pick the window size so that the algorithm would work with the best performance possible? For now I use impulse response size and a constant multiplier to determine window size, but since the input signal size might vary, how to determine multiplier in relation to the input signal size?

My concern is performance, but the bottleneck is not addition, but the amount of performed FFT convolutions. Experimentally I found out that for each impulse response size there is an optimal window size, yet it varies. For example, if I pick 40000 samples of input signal and 128 samples of impulse response signal the optimal ratio between window and impulse response lengths is approximately 1:24.

It seems that Wikipedia suggests to use 1:8 ratio for both methods, as seen in the pseudocode of both overlap-add and overlap-save methods, then in the next section they provide measures of complexity for the methods, but I'm unable to wrap my head around how to utilize this information.

The optimal size is likely not just algorithmic, but depends on the cache sizes and levels and the compute-to-memory-bandwidth ratio for the particular processors and systems(s) on which you plan to run this code.  Also, your applications latency requirements (which filling big FFTs add).


Linear Prediction task $$x[n] = \big(A_1 r_1^n\cos(\omega_1 n) + A_2 r_2^n \cos(\omega_2 n) + A_3 \sin(\omega_3*n) \big)u[n],$$

where $\qquad \omega_1 \ne \omega_2 \ne \omega_3, \qquad   |r_1|,|r_2|<1 $

I need an linear prediction model that perfectly predicts the signal with zero error $\forall n \ge n_0 > 0$ and how to find the $n_0$

I am not sure how to proceed with this as my initial thought was to minimize the energy of the first two terms that include $\cos(\omega_1 n)$ and $\cos(\omega_2 n)$ but something seems off.

Is this a homework problem?

I figured out a way to solve this, but it isn't that simple.

First, you can split your function into three parts:

$$\begin{align}
x_1[n] &= A_1 r_1^n \cos(\omega_1 n) \\
\\
x_2[n] &= A_2 r_2^n \cos(\omega_2 n) \\
\\
x_3[n] &= A_3 \sin(\omega_3 n) \\
\end{align}$$

The $u[n]$ doesn't matter because we are only concerned with $n \ge 0$ and $x_3[0]=0$.

The cosines can be solved like this:

$$\begin{align}
x[n] &= r^n \cos(\omega n) \\
\\
x[n-1] &= r^{n-1} \cos(\omega n-\omega) \\
&= r^{n-1} \big(\cos(\omega n)\cos(\omega) + \sin(\omega n)\sin(\omega) \big) \\
\\
x[n+1] &= r^{n+1} \cos(\omega n+\omega ) \\
&= r^{n+1} \big( \cos(\omega n)\cos(\omega)  - \sin(\omega n)\sin(\omega) \big)\\
\\
r x[n-1] + \tfrac{1}{r} x[n+1] &= 2 r^n \cos(\omega n)\cos(\omega) \\
&= 2 \cos(\omega ) x[n]\\
\\
x[n+1] &= 2r \cos(\omega ) x[n] - r^2 x[n-1]\\
\\
x[i] &= 2r \cos(\omega ) x[i-1] - r^2 x[i-2] \\
\end{align}$$

The sine can be solved like this:

$$\begin{align}
x[n] &= r^n \sin(\omega n) \\
\\
x[n-1] &= r^{n-1} \sin(\omega n-\omega)  \\
&= r^{n-1} \big( \sin(\omega n)\cos(\omega)  - \cos(\omega n)\sin(\omega ) \big) \\
\\
x[n+1] &= r^{n+1} \sin(\omega n+\omega)  \\
&= r^{n+1} \big( \sin(\omega n)\cos(\omega)  + \cos(\omega n)\sin(\omega) \big) \\
\\
r x[n-1] + \tfrac{1}{r} x[n+1] &= 2 r^n \sin(\omega n)\cos(\omega) = 2 \cos(\omega ) x[n] \\
\\
x[n+1] &= 2r \cos(\omega) x[n] - r^2 x[n-1] \\
\\
x[i] &= 2r \cos(\omega) x[i-1] - r^2 x[i-2] \\
\end{align}$$

Notice, whether sine or cosine, the LPC equations are the same, what matters are the initial conditions.  Your parts are multiples of these, but that is also taken care of by proper selection of initial conditions.

Here is some Python code that demonstrates that these equations work:

(Edit:  The code is now in the followup below.)

What remains now is to combine these three parts into one equation.  This is where it gets messy.

What you have is your solution in the form of:

$$\begin{align}
x_1[n] &= a_{11} x_1[n-1] + a_{21} x_1[n-2] \\
x_2[n] &= a_{12} x_2[n-1] + a_{22} x_2[n-2] \\
x_3[n] &= a_{13} x_3[n-1] + a_{23} x_3[n-2] \\
\end{align}$$

You need to convert them to a common set of coefficients ($v_k$'s) in this form:

$$\begin{align}
x_1[n] &= v_1 x_1[n-1] + v_2 x_1[n-2] + v_3 x_1[n-3] \\
&+ v_4 x_1[n-4] + v_5 x1[n-5] + v_6 x_1[n-6] \\
\\
x_2[n] &= v_1 x_2[n-1] + v_2 x_2[n-2] + v_3 x_2[n-3] \\
&+ v_4 x_2[n-4] + v_5 x2[n-5] + v_6 x_2[n-6] \\
\\
x_3[n] &= v_1 x_3[n-1] + v_2 x_3[n-2] + v_3 x_3[n-3] \\
&+ v_4 x_3[n-4] + v_5 x3[n-5] + v_6 x_3[n-6]
\end{align}$$

That way they can be added together:

$$\begin{align}
x[n] &= x_1[n] + x_2[n] + x_3[n] \\
\\
x[n] &= v_1 x[n-1] + v_2 x[n-2] + v_3 x[n-3] \\
&+ v_4 x[n-4] + v_5 x[n-5] + v_6 x[n-6] \\
\end{align}$$

The way I can think of right now is a lot of work, but there may be another slick way to simplify this.  With a bunch of algebra, the following linear algebra problem can be set up:

$$
\left[
\begin{array}{c}
a_{11} \\
a_{21} \\
a_{12} \\
a_{22} \\
a_{13} \\
a_{23} \\
\end{array}
\right]
=
W
\cdot
\left[
\begin{array}{c}
v_1 \\
v_2 \\
v_3 \\
v_4 \\
v_5 \\
v_6
\end{array}
\right]
$$

Where $W$ is a 6×6 array and its elements are functions of the $a$'s.  The $v$'s can then be found my multiplying the left array by $W^{-1}$.

I haven't done the algebra, but there is a pattern there.

Hope this helps, phew.

Ced

Followup:

So there is a much easier way to figure out the $v$ values.  I didn't figure it out myself, I asked on the Math Exchange.  (see https://math.stackexchange.com/questions/2755034).

You simply have to multiply out:

$$ (S^2 - a_{11} S - a_{21})(S^2 - a_{12} S - a_{22})(S^2 - a_{13} S - a_{23}) \\  \qquad \qquad \qquad \qquad = S^6 - v_1 S^5 - v_2 S^4 - v_3 S^3 - v_4 S^2 - v_5 S^1 - v_6 $$

Messy, but it works.  Here it is in Python:

import numpy as np

#==========================================================
def MultiplyTuple( p, q ):

    Np = len( p )
    Nq = len( q )

    r = np.zeros( Np + Nq - 1 )

    for ip in range( Np ):
        for iq in range( Nq ):
            r[ip+iq] += p[ip] * q[iq]

    return r

#==========================================================

x1 = np.zeros( 100 )
x2 = np.zeros( 100 )
x3 = np.zeros( 100 )
x  = np.zeros( 100 )
y  = np.zeros( 100 )

w1 = .1
w2 = .18
w3 = .05

r1 = .9
r2 = .5

A1 = 1.0
A2 = 0.5
A3 = 0.4

a11 = 2.0 * r1 * np.cos( w1 )
a21 = - r1 * r1
a12 = 2.0 * r2 * np.cos( w2 )
a22 = - r2 * r2
a13 = 2.0 * np.cos( w3 )
a23 = - 1

x1[0] = A1
x1[1] = A1 * a11 / 2.0
x2[0] = A2
x2[1] = A2 * a12 / 2.0
x3[0] = 0
x3[1] = A3 * np.sin( w3 )

print 0, x1[0], x2[0], x3[0]
print 1, x1[1], x2[1], x3[1]

for i in range( 2, 20 ):
    x1[i] = a11 * x1[i-1] + a21 * x1[i-2]
    x2[i] = a12 * x2[i-1] + a22 * x2[i-2]
    x3[i] = a13 * x3[i-1] + a23 * x3[i-2]
    print i, x1[i], x2[i], x3[i]

print

for i in range( 0, 20 ):
    y1 = A1 * r1**i * np.cos( w1 * i )
    y2 = A2 * r2**i * np.cos( w2 * i )
    y3 = A3 * np.sin( w3 * i )
    y[i] = y1 + y2 + y3
    print i, y1, y2, y3

p1 = np.array( [ 1, -a11, -a21 ] )
p2 = np.array( [ 1, -a12, -a22 ] )
p3 = np.array( [ 1, -a13, -a23 ] )

p12  = MultiplyTuple(  p1, p2 )
p123 = MultiplyTuple( p12, p3 )

print p123

v1 = -p123[1]
v2 = -p123[2]
v3 = -p123[3]
v4 = -p123[4]
v5 = -p123[5]
v6 = -p123[6]

print

for i in range( 0, 6 ):
    x[i] = x1[i] + x2[i] + x3[i]
    print i, x[i], y[i]

for i in range( 6, 20 ):
    x[i] = v1*x[i-1] + v2*x[i-2] + v3*x[i-3] \
         + v4*x[i-4] + v5*x[i-5] + v6*x[i-6]
    print i, x[i], y[i]


Bandwidth of a composite signal How to find bandwidth of a composite signal that can be decomposed into five sine waves with given frequencies and having all peak amplitudes same ?

Unless I'm missing something here it would just be the difference in frequency between the highest frequency sine wave and the lowest:

BW = Fmax - Fmin


Get spectral picture from a wavelet transform According to my previous question, I have changed the generate command to:

y=generate1(100,1000,1);

and got the following picture:

Now I want to test a wavelet transform on the same signal, which of course exists in Matlab, but I have  one thing which I should  clarify: generally, as I know large scale values correspond to small frequencies and small scale values to large frequencies.  So, how should I determine if the given frequencies are small or large? Also, which wavelet transform should i use?

For example:

m=cwt(y,1:100,'sym2');

and plot it

plot(m)

How should I read the second picture? What kind of information can I get?

Try this one:

m=cwt(y,1:100,'sym2','plot');
colormap(pink)

Result (btw, I cannot reproduce your curve in time domain with the function generate1):

Low scale values compress the wavelet and correlate better with high frequencies. While high scale values stretch the wavelet and correlate better with the low frequency content of the signal.


Zigbee Transceiver Using OQPSK Modulation under Fading Channels I am working on transmitting ECG signal over wireless body area networks. According to IEEE802.15.4 standards , I am using ZigBee transceiver at 2.4 GHZ . Complex baseband equivalent channel model is used.

The steps are as follows :

1.      Signal Compression

2.      Quantization

3.      Coding

4.      OQPSK modulation using the matlab function (oqpskmod)

5.      Fading Channel plus AWGN

6.      Equalization

7.      OQPSK Demodulation using the matlab function (oqpskdemod)

8.      Decoding and Dequantization

*  According to the IEEE 802.15.4 standard, a pulse shaping step is performed in the transmitter after the OQPSK modulation step .

I didn`t perform this pulse shaping process, and I obtained reasonable results. Is it necessary to perform this step ?
If  yes , How should the receiver be modified ?

The pulse shaping is only required if you are concerned with spectral efficiency, as the pulse shaping will convert what would otherwise be a Sinc function in frequency to a spectrum that is much closer in overall bandwidth to your symbol rate. Also if the receiver assumes a specific pulse shape, the sensitivity will be optimized by transmitting with the matched pulse shape (as the matched filter in the receiver will then present the optimum SNR. Other than that, it is not necessary.

See my picture below demonstrating the time domain waveforms and frequency domain before and after pulse shaping. Note that the primary motivation to do additional pulse shaping (beyond the rectangular pulse you are otherwise using, which is indeed a pulse shaping in itself) is for spectral efficiency. The pulse shaping is required in IEEE 802.15.4 ultimately due to meeting regulations (such as FCC in the US) to ensure we all use spectrum efficiently and within regulated allocations and not interfere with users in other spectrum allocations. Observe the dramatic change in transmitted spectrum with and without the pulse shaping filter in the example below!


In the context of DFT, Where Does the Nyquist Frequency Sample Belong In a Double Sided Frequency Spectrum (Positive / Negative Side)? If we have an even number of data points $N$, after DFT in MATLAB, the output has the order:
$$(\text{DC}, f_1, f_2, \ldots, f_{N/2-1}, f_\text{Nyq}, -f_{N/2-1}, -f_{N/2-2}, \ldots, -f_1)$$
For real signals, the first output corresponding to $k$=0, is real and so is the Nyquist frequency. After that numbers are complex conjugates.
If we are interested in a single sided spectrum, the Nyquist frequency is shown on the positive side.
However, when a double-sided frequency spectrum is plotted, many authors put the Nyquist frequency on the negative side.
Some software like OriginPro, follow the opposite. Is there a fundamentally correct way or is it just a convention i.e.,
$$
\text { If } N \text { is even, } \quad k\quad\text { takes: }-\frac{N}{2}, \ldots,-1,0,1, \ldots, \frac{N}{2}-1
$$
Alternatively,
$$
\text { If } N \text { is even, } \quad k \text { takes: } -\frac{N}{2}-1, \ldots,-1,0,1, \ldots, \frac{N}{2}$$
where $k$ is the DFT index vector, which is used to construct the frequency axis as
$$\text {Frequency axis}=k/ N\Delta t$$
where $\Delta t$ is the sampling interval.
Many people say it is just a convention and both are correct. Thanks.

It's convention, they're equivalent:
$$
\exp{\left(j2 \pi \frac{N}{2}n/N \right)}
 = \exp{\left(j2\pi \frac{-N}{2}n/N\right)} \\
\Rightarrow e^{j\pi n} = e^{-j \pi n} \Rightarrow \cos(\pi n) = \cos(-\pi n)=(-1)^n,\ j\sin(\pi n) = j\sin(-\pi n) = 0
$$
MATLAB and Numpy go $[-N/2, ..., N/2-1]$, which is unfortunate for analytic representations (+ freqs only). Note also its value is doubled relative to other bins (but not manually; they correlate this way), so in a sense it's both a negative and positive frequency, so energy's preserved:

You can tell a library's preference by fftshift docs:


Total bandwidth in a network calculation I have data on the number of devices in a wired network and the data rates on links among these devices. I also know the positions of these devices from which I can compute the distances between the devices.

The idea is that if I were replace the wired network with wireless infrastructure, how much bandwidth am I looking for. Would let's say 200MHz be sufficient for such a case.

Assuming that the BER is very good, how should I go about in calculating the total bandwidth in such a network.

There are many variables involved, but you may proceed as follows. Assume the bit rate in a link is $R_b$. First of all, you choose a constellation; let's say you go with 16-QAM, which transmits four bits per symbol. Now you have a symbol rate $$R_1=R_b/4.$$ Next, you choose an encoder that gives you the BER you need. This encoder will have a certain code rate $R_c$. After encoding, you now have a symbol rate $$R_2=\frac{R_1}{R_c}.$$ Next, you need to take into account the overhead: you need framing and training data. Let's say the overhead is $x\%$. Then your final symbol rate is $$R_3=(1+x)R_2.$$ Now you can calculate the bandwidth you need: $$B=\gamma R_3,$$ where $\gamma$ is an "excess bandwidth" over the theoretical $B=R_3$. The excess bandwidth is typically between 1 and 2.

In summary: the symbol rate is given by the data rate, the constellation, the encoder, and the overhead. The bandwidth is determined by the symbol rate and the excess bandwidht factor.

Note that each of these variables affects the others. For example, if you choose QPSK instead of 16-QAM, you'll have better BER, but $R_1$ will be larger. However, now you need an encoder with smaller rate, so $R_2$ will be smaller. And so on. So, you need to find a combination of variables that you can implement within your constraints and that gives you the BER and data rate you need.


Train of impulses in frequency domain Why compressing a train of impulses in time domain produce a wide-spaced train of impulses in frequency domain and vice versa? i want to know the intuition behind that.

The intuition is simply the update rate in the time domain. If I repeat an impulse in time once per second, it should be quite intuitive that a 1 Hz tone will exist in frequency. When a waveform isn't exactly sinusoidal (such as a distorted sinewave), it will have higher frequencies that will only exist at harmonics of that waveform. Since the waveform repeats over a given time $T$, the only other frequency components that can exist must also repeat over that same time interval $T$, which are the harmonics. In general, whenever you repeat any waveform exactly the same in the time domain, as we see with the Fourier Series Expansion, the only non-zero frequencies that can possibly exist will be at DC, the repetition rate (fundamental frequency) and all higher harmonics. (Some of these may themselves be zero in the general case, but any non-zero frequencies can only exist at this location). If I said your name exactly the same way once per second, the only frequency content of my voice track will be at 1 Hz spacings.
This is the same with a train of impulses in time, just in this case we are repeating an impulse.
Further, you can take the Fourier Transform of any base signal (such as the impulse, which is a constant for all frequencies), and then take the property above showing where the only non-zero values can be, and we see that we end up "sampling" in frequency the Fourier Transform. This also works with any base signal where a Fourier Transform exists, and explains why a square wave has odd harmonics (which are simply samples of the underlying Sinc function as the Fourier Transform of a pulse).


How to mix two audio files recorded with different settings to obtain a realistic result I am a software engineer, currently working on a signal-processing related project, which I am new to. For evaluating the algorithm(which is written in R), I need to mix (merge) two different audio signals synthetically. The question is how to mix these two signals and obtain a “realistic” result?

To be more specific, consider the following scenario:

I have recorded Signal_A, and downloaded signal_B from the internet (therefore two different microphones/settings are used for recording).
The sound pressure level of both signals are known( I measured the sound pressure level of Signal_A in 1m distance, and also know the sound pressure level of Signal_B in the same distance).
As far as I know, the value recorded/represented in each waveform of any signal is relative to that corresponding signal, and since two signals are recorded using different settings we can expect that one might be more amplified than the other. Therefore a simple mixing of the two won’t give us a realistic result.

Now the question is: is there a way to manipulate recorded signal_A and signal_B programatically, assuming the fact that the SPL values of the original signals at 1m distance is known, and get the merged signal as they were playing simultaneously and recorded at 1m distance? i.e. to calculate factor ‘a’ and ‘b’ in a way that the realistic merged signal is formulated by:

Mereged_Signal = a*Signal_A  + b*Signal_B

Thanks in advance!

I have recorded Signal_A, and downloaded signal_B from the internet (therefore two different microphones/settings are used for recording)

You might be able to achieve a "seamless" mix (or transition) but the fact that you have two different microphones and settings might add artifacts in the final product. You can tell how this sounds like if you have ever noticed the differences between microphones. Accounting for such differences might be a bit more challenging.

  The sound pressure level of both signals are known( I measured the sound pressure level of Signal_A in 1m distance, and also know the sound pressure level of Signal_B in the same distance).

It is good that you have measured the SPL but most commonly you are unlikely to know the exact relationship that takes you from SPL to converted values.

  As far as I know, the value recorded/represented in each waveform of any signal is relative to that corresponding signal, and since two signals are recorded using different settings we can expect that one might be more amplified than the other. Therefore a simple mixing of the two won’t give us a realistic result.

I am not entirely clear on what you mean by "a simple mixing" but the first approach would be to simply adjust the amplitudes of the two recordings and then "add" them together as per your MergedSignal equation.

  Now the question is: is there a way to manipulate recorded signal_A and signal_B programatically, assuming the fact that the SPL values of the original signals at 1m distance is known, and get the merged signal as they were playing simultaneously and recorded at 1m distance?

Yes. BUT!, you have to make sure that the rest of the signal pathway was the same as well. Now:

Sound waves hit the transducer. The transducer converts them to voltage with a (hopefully) linear relationship of Volts per Pascal of pressure. In addition, this relationship is frequency dependent. But let's say that it is constant across the spectrum and call it $\alpha$.
Behind the transducer there is an amplifier whose job is to take the electrical signal and bring it up to some level. Again, this is achieved with a (hopefully) linear relationship. Let's keep things simple and assume that the job of this amplifier is to take the signal of the transducer and adapt it to a standard line level of approximately 2 Volts (peak-to-peak). Let's call this $\beta$.
The final stage is the Analog to Digital Converter (ADC). It converts the 2volts peak to peak to some range with a given word length (for example 0..255, -127..128 and so on).

So, assuming that $\alpha$ (mic), $\beta$ (amp/signal conditioning) are the same, you can indeed scale the recordings by the ratio of the SPLs. That is:

$$MergedSignal = (1.0 - \frac{SPL_A}{SPL_B}) \times Signal_A + \frac{SPL_A}{SPL_B} \times Signal_B$$

In the more general case and when you don't know the SPL, you can simply "match" the average level of the converted values. So, same relationship, instead of $SPL_A, SPL_B$ you have the averages of the signals.

Hope this helps.


multiplication in frequency domain equals convolution in time domain mismatch We know that a convolution in the time domain equals a multiplication in the frequency domain.
  As per the article below,

https://www.dspguide.com/ch9/3.htm

In order to multiply one frequency signal by another, (in polar form) the magnitude components are multiplied by one another and the phase components are added.
 To prove this, i created two sinusoidal signals in MATLAB,

Freq1 = 1000;
Freq3 = 3000;
Fs = 16000;
T = 1/Fs;
Nos = (0:128-1)*T;
Amp = 1.0;
Signal1 = Amp*sin(2*pi*Freq1*Nos);
Signal3 = Amp*sin(2*pi*Freq3*Nos);

In the time domain, i convolved these two signals,

Signal4 = conv(Signal1,Signal3);

Taking a 32 point FFT of the two signals & and a 32 point FFT of the convolved signal,

NFFT = 32;
freqdata1 = fft(Signal1,NFFT);
freqdata2 = fft(Signal3,NFFT);
freqdata4 = fft(Signal4,NFFT);

in the frequency domain, i multipled the magnitude components of the two individual signals and added the phase component of the two signals. I compared this Magnitude and phase value with the Convolved signal's phase and magnitude value. I expected the values, [Newmag' NewPhase'] & [Mag3' Phase3']  to be similar since the a convolution in time domain equals a multiplication in the frequency domain. But they are not. What am i missing here ? What have i done wrong ?

for ii = 2:((length(freqdata1)/2)+1)

   sig1_cc = real(freqdata1(1,ii));
   sig1_dd = imag(freqdata1(1,ii));
   Mag1(ii-1) = sqrt((sig1_cc^2)+(sig1_dd^2));
   Phase1(ii-1) = atan(sig1_dd/sig1_cc);

   sig2_cc = real(freqdata2(1,ii));
   sig2_dd = imag(freqdata2(1,ii));
   Mag2(ii-1) = sqrt((sig2_cc^2)+(sig2_dd^2));
   Phase2(ii-1) = atan(sig2_dd/sig2_cc);

   sig3_cc = real(freqdata4(1,ii));
   sig3_dd = imag(freqdata4(1,ii));
   Mag3(ii-1) = sqrt((sig3_cc^2)+(sig3_dd^2));
   Phase3(ii-1) = atan(sig3_dd/sig3_cc);

   Newmag(ii-1) = Mag1(ii-1)*Mag2(ii-1);
   NewPhase(ii-1) = Phase1(ii-1) + Phase2(ii-1);

end

[Newmag' NewPhase']
[Mag3' Phase3']

You missed the part where under the DFT, multiplication in one domain is equivalent to circular convolution in the other. Standard problem!

Just pad the signals you want to convolve with enough zeros on both sides, and things will look better.

Also, I hope you're just doing this as to exercise your matlab coding, because getting the phase out of two complex numbers just to multiply them is way, way, way more computationally intense than just directly multiplying the two numbers.


How much can we downsample by? and how much can we upsample by? e.g. let's say we have a signal with B = 100 Hz that is sampled at 400 Hz. The anti-aliasing filter will remove replicas, and there is nothing to alias, and we downsample. If we downsample by 2, we will get signal of 200 Hz. But what's stopping us? Can we downsample by 4, or 8?

(why can't we downsample by 4, get B = 400Hz, then have a new sampling rate at 800 Hz?)

Similarly, is there a limit to how much we want to upsample?

Nyquist. If your signal is real and has a BW of B = 100 Hz (meaning DC to 100 Hz), then the most you can downsample by is 2, in order to obey Nyquist and have a sampling rate that is at least 2B = 200 Hz.  Note you will typically want to sample more beyond this limit to allow room for reasonable filtering.


Architecture for high sample-rate FIR/IIR filter on FPGA or ASIC? Assume we have a hardware (FPGA or ASIC) which can operate at a fixed processing frequency of $100\textrm{ MHz}$. Assume we have an input data stream of $1\textrm{GSPS}$ (which has been sampled using any SERDES interface, but here the sampling technique does not matter). We would like to implement a filter which can NOT be implemented by a poly-phase decimation filter; for example assume we want to implement a band-pass filter which passes the frequency range $200-300\textrm{ MHz}$ in our example.

What architectures can be used for this purpose?
Can anyone provide a reference for designing such filters?

Short answer: If your sampling rate is higher than your throughput, you must use some kind of polyphase decomposition of your filters to work. So if you can't do that, you can't use a 100 MHz-clockable design to filter something with a sample rate of 1000 MS/s. It's as simple as that.

However, your example is exactly the standard example of a filter that can be implemented as a polyphase decimator: you get in 1000 MS/s, you filter to 100 MHz of that, and decimate by ten (note: assuming complex sampling here). You can do that, because your bandpass restricts the bandwidth that passes through to $\frac1{10}$ of the input bandwidth.

Anyway, even so, this will very probably not work out – there's simply zero margin for "real world" effects in your 100MHz clock to produce 100 MS/s, so you should really look into a design that can work faster.

Also, it's usually not that you buy an FPGA and can say "this works at 100MHz, no matter what I do on it", but you know its maximum clocking (which typically is much higher than 100MHz for devices that can take in 1GS/s via SERDES these days), and the maximum clock rate achievable depends on the DSP you do – so you can typically just ask your FPGA design tool to generate a FIR that is optimized for high rate (rather than e.g. for low power or small ressource consumption), and you'll get a different possible rate depending on what FIR architecture you chose.