Find a sinusoidal fit from data with fft I have a vector with raw data (120 samples). From the frequency spectrum I can tell there is one dominant frequency component. Now I am supposed to use the fft data (magnitude and phase) to find a good sinusoidal fit and offset to fit the original data.

I am completely stuck on this one. Can anyone help?

Depending on kind of signal you are dealing with I can think of the easiest method to do so.

Take the DFT of your signal (remember about proper scaling of values).
Find the dominant frequency $f_0$ and its amplitude $A$ at magnitude spectrum.
From phase spectrum find the corresponding phase shift $\theta$ for this given frequency - usually it is given in radians.
Reconstruct your signal as: $x(t)=A\sin(2\pi f_0t+\theta) $

Some quick and dirty code is included below. No warranty of working whatsoever.

#!/usr/bin/env python

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

if __name__ == "__main__":
    # Signal and analysis parameters
    fs = 4000.0    # Sampling frequency
    T  = 2         # Duration of a signal
    f0 = 100       # Fundamental frequency
    A  = 1         # Sine amplitude
    #th = np.pi/4   # Phase shift in radians (45 degrees)
    th = 0

    # Time vector
    t = np.arange(0, T*fs)/fs

    # Generate the signal with higher harmonics
    x = A*np.sin(2*np.pi*f0*t + th) + \
        A*np.sin(2*np.pi*2*f0*t + 2*th)/8 + \
        A*np.sin(2*np.pi*3*f0*t + 3*th)/16

    # Add some noise
    x += 0.5*np.random.randn(x.size)

    # Number of all samples
    N = x.size

    # Apply the window
    win = sig.gaussian(N, N)
    x = x*win

    # Perform the DFT
    X = np.fft.fft(x)
    # Get the frequency vector
    freq = np.linspace(0, fs-1/fs, N)

    # Get the phase and magnitude
    Xmag = np.abs(X)/np.sum(win)*2
    Xph  = np.angle(X)

    # Magnitude in logarithmic scale
    Xmagl = 10*np.log10(Xmag/Xmag.max())

    # Find the amplitude of the dominant frequency
    ind_max = np.argmax(Xmag[0:np.ceil(N/2)])

    # Extract the amplitude ...
    A_ = Xmag[ind_max]
    # ... frequency ..
    f0_ = freq[ind_max]
    # ... phase
    th_ = Xph[ind_max] + np.pi/2

    # Reconstruct the signal
    x_ = A_ * np.sin(2*np.pi*f0_*t + th_)

    # Do some plotting
    plt.subplot(311)
    plt.plot(t, x)
    plt.grid(True)
    plt.hold(True)
    plt.plot(t, x_, 'r')
    plt.title('Time domain signal')

    plt.subplot(312)
    plt.plot(freq[0:N/2], Xmagl[0:N/2])
    plt.grid(True)
    plt.title('Magnitude')

    plt.subplot(313)
    plt.plot(freq[0:N/2], Xph[0:N/2])
    plt.grid(True)
    plt.title('Phase')
    plt.show()


Windowing a time domain OFDM signal Windowing is a popular method of reducing the spectral sidelobes of OFDM. A popular window usually used for this purpose in the raised cosine (RC) window, because of its tapered and smooth edges. I understand the concept of windowing and why it reduces the sidelobes, however what I don't understand is this: "windowing a time domain OFDM signal is equivalent to every subcarrier being individually windowed by the window function (RC for example), basically the sinc sidelobes of every subcarrier will be reduced". This is true and you can actually see it when the window function is applied to a non-contiguous OFDM signal (some subcarriers are deactivated somewhere in the middle of the spectrum), the spectrum is not only reduced at the edges of the band, but also in the in-band. You can't get this if you apply a filter to the signal, with a filter you only see power spectrum reduction in the edges but not in the band. I'd appreciate if someone can explain why this is true, preferably with some equations.

If I understand the question properly, then it sounds like the difference between time-domain multiplication and frequency domain multiplication.

Windowing is a time-domain multiplication operation.  As such, it is equivalent to convolving the spectrum with the window kernel, which will operate all frequencies the same way.  On the other hand, filtering is (can be) a frequency domain multiplication, so the filter can be designed to operate specifically on chosen frequencies (or bands of frequencies if you like).

To look at it slightly differently: if you consider your signal as the sum of all your subcarriers (which you can do with linearity), then windowing your signal with an RC is the same as windowing each of your subcarriers with an RC before adding them together.


Meaning of negative frequencies in Baseband non sinusoidal, non periodic signal I can understand the meaning of negative frequencies in a sine or cosine signal, since by using Euler's identities, you have two complex phasors moving to different directions, which when added, give a real (cosine or sine) signal.
However if I have a random baseband (real, not periodic) signal, what does it mean for it to have negative frequencies?

It means the same thing.
The negative frequencies are there to cancel out the complex portion of the signal in the time domain, so that the time domain signal is constrained to the real axis.
The Fourier Transform decomposes a signal, $x(t)$ onto a set of complex basis functions: $e^{-j\omega t}$ for $\omega \in (-\infty, \infty)$.
If $x(t)$ is real, then the decomposition, i.e. the Fourier Transform $X(\omega)$, must be such that the complex portions cancel out, when the inverse transform is performed to recompose the signal from it's constituent parts.
This is a conceptual explanation of why the Fourier Transform of a real signal is symmetric about $0$ Hz.


THD: estimation using periodogram or fft? I want to estimate the THD of a voltage signal (output from a frequency inverter). At the begining, I need to estimate the power of the signal at frequency domain. To doing so, I should use a periodogram or a FFT.

Matlab uses periodogram in his THD function. I know that are differences between FFT and periodogram, as shown in this question. But, wich method should I use? and why?

Thanks.

Just for answer the question, basically there is no difference at all. My problem was that I've compared a modified periodogram against a FFT. Modified periodogram does have a difference; but simple periodogram doesn't.


Confusion implementing chebyshev type 2 band pass filter in matlab I want to design a band pass Chebyshev Type II filter using analog prototyping. The order of filter is 20 with a value of 60 dB stop band attenuation and 0.75 dB pass band ripple where,
Pass band edge = 800 Hz
Stop band edge = 2000 Hz
Sampling frequency = 6000
My Matlab code is below
clc;  close all;   clear all;
Rp = 0.75;    Rs = 60;
fp = 800;
fs = 2000;
Fs = 6000;

fn = Fs/2;     % nyquist frequency
wp = fp/fn;     % normalized Pass band corner frequency
ws = fs/fn;     % normalized Stop band corner frequency

[n,Wn] = cheb2ord(wp,ws,Rp,Rs);
[num,den] = cheby2(20,Rs,Wn);
freqz(num,den,512,Fs);

I have also attached snapshot of output which shows that output  frequency plot looks like low pass,but i think ,why output plot is not band pass?How can it look band pass?

You're using the function cheb2ord in the wrong way. Assuming that your lower stopband edge is $800$ Hz and your upper stopband edge is $2000$ Hz (note that a bandpass filter has two stopband edges, as well as two passband edges), and if the sampling frequency is $6000$ Hz, then the following command should do what you intended to do:
[b, a] = cheby2(10, 60, [4, 10]/15);
Note that for bandpass and bandstop filters, the resulting filter order is twice the number you use in the call to cheby2.


Do JPEG formats add significant overhead? For an extreme example, I have the following image that was a JPEG before stackexchange converted it to a PNG:-

Before you perform all sorts of funky image processing tricks on it, it's just a black square.  This is what happens with my web cam when there is insufficient light.  It just gets low pass filtered out and a black square results. The interesting thing is that when I pull up 100 consecutive images from the camera, they're all exactly 6616 bytes long. There must be an internal  light level threshold.

My question is, why was this square's original file size fixed at 6616 bytes? That seems like a lot of overhead to store just a single repeated value of "0", and given that extremely low frequency images should compress greatly. Also, given the random (quantum) nature of the various noise sources impacting the sensor, I would expect some natural variation in size.

For reference, the original image is available here, and it came from this:-

The picture you've posted most certainly didn't come directly from the image sensor – it was processed by the camera hardware, and software; it's pure dark grey.

Anyway; JPEG (to be more specific: JFIF as the JPEG-carrying file format) is not zero overhead, and you can expect a couple kB for this kind of image to work.

There's a couple bytes of headers (it's really negligible), but then things get a bit more complicated:

For a monochrome image, the actual JPEG quality setting won't do anything to the file size; it just sets a quantization mask, and you're always hitting the DC bin with a value of nonzero for anything. No matter how that turns out, Run-Length Encoding for the rest of the bins in each 8x8 pixel block will reduce the rest to the same number of bytes, regardless of quantization.
What does matter is the subsampling mask for the channels. Stanley's answer is a tiny bit wrong w.r.t. to the three channels: JPEG does support 4:2:0 subsampling, meaning that you actually omit one channel every row. But that's a story for a different day. Anyways, since your image doesn't change spatially at all, that saves space
Entropy encoding (I think it's a Huffmann method) is applied to the rest. Since the picture is built from identical 8x8 blocks, these should be reduced  to very little data at all.
There's the special JPEG/JFIF feature of interleaving the info in a way that allows the lower frequency parts of the image to be transferred first (e.g. for slow connections; if you grew up in the 90s, you'll remember pictures that got sharper gradually while you hoped for your dialup modem to take less costly time). While these practically make little difference for "typical" photos, for your monochrome picture, they'll make matters worse and increase the file size by some 20 to 25%, maybe.

So, 6 kB is actually pretty OK, and you should probably look into whatever gave you that picture – it's definitely heavily smoothing the image, to a degree that the data you're presented has nothing to do with what the sensor sees. The sensor sees noise.


can't find out the period of my signal I've a noisy signal, which is the sound of motor with a constant speed, so  the sound "should" be  periodic, I know that there is a way to use the autocorrelation function to get the period,I did it, but I can't figure out the period. Any idea how to do that bellow the signal and the result of the autocorrelation :
the signal

the autocorrelation result :

Here's an attempt to do what you need in scilab.

The top plot shows some data that I synthesized.  The second plot shows the auto-correlation of the raw data (zoomed in around the central peak of the auto-correlation).  The red circles show the peaks found using this find_peaks function.

The final plot shows the difference between all the peak locations. This will be an estimate of the period.  Because you're not really guaranteed that the underlying period will be an integer number of samples, you should probably find the mean of these values.

In this case, the "true" period is 1/f0 = 11.191996, and taking diffs = diff(peaks); and then mean(diffs(10:173)) yields 11.195122.

Code below.

N = 1000;
f0 = 0.0893495634;
phi = rand(1,1,'uniform')*2*%pi;
sigma = 0.5;

x = sin(2*%pi*[0:N-1]*f0 + phi) + sigma*rand(1,N,'normal');
XC = xcorr(x);

clf
subplot(311)
plot(x);

subplot(312)
plot(XC);

peaks=peak_detect(XC,0);

plot(peaks,XC(peaks),'ro')

a = get('current_axes');
a.data_bounds=[950 1050 -500 800];

subplot(313)
plot(diff(peaks));


What is the proof of the equation for the aliased frequency fa = |R*n - fs|? R (sampling rate)
fs (signal being sampled
fN (the Nyquist frequency)
fa (aliased frequency)

I was hoping someone could explain to me why this is true with a more mathematical proof. it seems reasonable for the aliased frequency to be related to these parameters, but I would like a more solid rationale.

This is incomplete but should point you in the right direction. Assume $k$ is an integer and $f_s=1/T_s$ is the sampling frequency. The signal being sampled is $$\cos(2\pi f_0t).$$ The sampled signal is $\cos(2\pi f_0kT_s)$. This signal has an alias at frequency $f_1=(f_s/2-f_0)+f_s/2=f_s-f_0$. To prove it:
\begin{align*}
\cos(2\pi(f_s-f_0)kT_s) &= \cos(2\pi f_skT_s-2\pi f_0kT_s) \\
&= \cos(2\pi k - 2\pi f_0kTs)\\
&= \cos(2\pi f_0kT_s).
\end{align*}
A similar argument can be made for all other aliases.


How to get coordinates of complex RTL-SDR signal? As I have already explained here (https://stackoverflow.com/questions/49358449/python-how-to-get-coordinates-of-complex-rtl-sdr-signal-peaks?noredirect=1#comment85722604_49358449), I need to find the coordinates of the peaks in my signal. While my original question is related more to the programming part, this one is about understanding what my signal has to look like for the peak finding to work at all.

It is my understanding that I'm getting a bunch of complex data from my RTL-SDR. For visualisation purposes I plot this data as a PSD, showing the frequencies on the x-axis and the power level in dB on the y-axis.

Now, I would like to use the Python peakutils.peak.indexes() method to get the y-coordinates of all peaks above a certain power level, e.g., -20.

This power level is specified as threshold in the method, which apparently means that the input data, i.e., the signal, somehow needs to represent all the power levels of the signal (so the method can search for those standing out)?! Is this assumption correct?

And if so, how do I get a representation of the signal that I could use to detect peaks by power level? Right now, the method seems to scan through frequencies (represented by my complex signal?!), looking for a power level (!) of -20, which probably is why it returns lots of useless results.

I already tried converting the complex signal array to a real array, which didn't work. I also tried converting the signal via fft and fftfreq, which didn't get me any further either.

As I said at the beginning: I can figure out the programming stuff myself. But I need someone to explain to me first, what I have to do, please. How do I need to convert the complex signal?

You are missing several steps.

You have (complex) time-domain samples.  You need to convert that to scalar frequency domain data to get something like your example plot.  Many steps are required to do that conversion before there are any frequency peaks to be found.

You first need to:  Pick a segment of complex IQ data from the RTL-SDR of length N (you may need to convert the raw RTL-SDR IQ samples from unsigned 8-bit to signed floating point), window it (von Hann or Hamming, etc.), convert it to the frequency domain via an FFT of length N, convert the FFT result to log magnitude, and label the FFT log magnitude result bins by frequency, which will be roughly

frequency(i) = sdr.center_freq + i * sdr.sample_rate / N

for bins 0 to N/2

frequency(i) = sdr.center_freq - (N - i) * sdr.sample_rate / N

for bins N/2 to N-1

Then you can search along that log magnitude array for peaks, and apply the frequency label to the peaks found.


Filtering signal with added noise (frequencies known) I have been assigned a project to take a given noisy sound file, with added broadband and narrowband noise. The goal is to reach a certain mean-squared error between the filtered signal and the original, non-noisy signal.

Here is my code for calculating the mean squared error (in MATLAB)

%  let x be the original signal, xn be the noisy signal, and xnf be the filtered
%   noisy signal. Let the prefix 'G' indicate a power spectral density,
%   a lowercase 'x' indicate time domain, and uppercase 'X' indicate FFT.
%   Compute the cross correlation between the original signal and the
%   filtered noisy signal
Rxxnf = N*real(ifft(ifftshift(Gxxnf))) ;
%   Find the number of samples of delay, I, for peak cross correlation
    [Rmax,I] = max(Rxxnf) ;
%   Circularly shift the original signal by I samples to make it have a
%   maximum correlation to the filtered noisy signal.  Call the shifted
%   signal "xs"
 xs = [x(end-I:end);x(1:end-I-1)] ;
%   Compute the mean squared error between the original signal (delayed)
%   and the filtered noisy signal.  Compute only for the samples that occur
%   after the delay.
   mse = mean((xnf(I:end)-xs(I:end)).^2) ;

The sound clip is about 10 seconds long, and is of a man speaking.

The locations of the narrowband noise are specified as ranges from

967.884 Hz to 1069.95 Hz
715.365 Hz to 810.994 Hz

Looking at the FFT of the original signal (left) and the noisy signal (right), we can see that there is indeed a significant amount of broadband noise.

To inspect this further, I look at the power spectral density.

I can confirm that the  broadband noise is the outermost section of the noisy PSD and that narrowband noise has indeed been added at the indicated frequencies.

I have applied a lowpass filter (6th order butterworth, $f_c=700\textrm{ Hz}$ (the broadband noise is actually introduced closer to $1000\textrm{ Hz}$, but $700\textrm{ Hz}$ gives the best mean-squared error), and 2 bandstops at the given frequency ranges (chebyshev type II, both second order, one with a passband ripple of $3\textrm{ dB}$ (low ripple gives best MSE), and one with a passband ripple of 18 (best MSE).

As of right now, my MSE is not as low as my professor's claimed MSE for this particular signal (he indicates a mean-squared error of .00143, I have only achieved .0023)

My questions are:

are there any other filter types I should be using rather than butterworth/chebyshev?
Should I place any other filters anywhere to perhaps mitigate ripple effects?

I very much appreciate the help.

You want a Wiener filter here.  Basically your frequency response at frequency $f$ should be $\frac{E |S|^2}{E |S|^2 + E |N|^2}$ where $S$ is the Fourier transform of the signal and $N$ of the noise (all dependent on $f$ of course).  So your filters should not aim for a stop band of zero but rather a stop band that will admit an amount of signal based on the expected SNR at that frequency.

There are comparatively few staircases in the total resulting frequency response here, so you might consider going for forward and backward in time IIR filtering (netting a zero phase delay, with normal delay the Wiener filter design becomes pointless for more than a 90degree phase shift) or some FFT with windowed overlap-add filtering.

Wiener filters like that minimize the expected mean square error, so they seem to be the right tool here.


LTI system, Fourier transform, unit step function, sine wave I just finished learning about Fourier Transforms and don't understand this signal:

$$x(t) = \cos(\omega t)u(t) $$

This is a cosine wave but only where $\omega$ is positive. My question is what can I do to figure out the output when the impulse response $h(t) = e^{-2\pi t}u(t)$ ?

I first try to break $x(t)$ up using Euler to get:

$$x(t) = \frac{1}{2} \left[e^{j2\pi t} + e^{-j2\pi t}\right]$$

then get FT of each component and use convolution property to get:

$$Y(j\omega) = H(j\omega)X(j\omega)$$

but kept getting an incorrect answer.

I can see three problems here:

1) $x(t)$ is a cosine wave but only where $t$ is positive, not $\omega$

2) You must keep the step function when applying Euler's identity, and I can't see the cosine frequency in your second expression for $x(t)$ either. Always make sure your identity is right, with WolframAlpha for example

3) Do not mix up your $\omega$ that is a parameter to $\cos(\omega t)$ and the $\omega$ variable of the frequency domain in $X(j\omega)$. Use different notations, like $\cos(\omega_0t)$


Detection periodic elements in image I am working on a university project in which I need to find a periodic net element in an image.

The net is a set of diagonal lines in a specific angle ( which I don't know in advance ) on a noisy image. They are hard to detect so I would like to use all the information I can get my hands on.

My current algorithm is based on using edge detection + hough transform to find lines at specific angle range.

I was wondering if there is any way to detect periodic signals in image? Something based on ff2 or something like that...

Thanks,

You could apply a 2d FFT. You will get a frequency spectrum as z values on a 2d plane:

 ← That's log(abs(fftshift(fft2(image))))) where image is from the question, but with white borders trimmed

From the z values you can determine the highest peak and calculate its position on the 2d plane. With some geometric calculation, you can derive wavelength, direction and phaseshift from this.


Determine surface homogeneity (luminance, colour, shape) by fourier transform Let's say I have a computer screen, or a TV, or a flat panel ceiling lamp, i.e. a certain luminance spread over an area.

If I wanted to determine, how homogeneous this source of light is, I'd take a picture of it and put it through a fourier transformation, e.g. a DFT.

If the panel were perfectly homogeneous, there'd only be one single peak at the zeroth frequency and nothing else. Of course, they never are, and the calculations are done in a discrete system, so it will resemble a 1/x function.
Depending on the amplitudes at frequencies > 0, the panel will be homogeneous to a certain degree. High frequencies might indicate sharp edges, circular blobs in the actual image might be represented by circles in the amplitudes spectrum as well.

From such an amplitude spectrum alone, it would be very difficult to get a sharp definition for "good" or "bad".

As far as I know, the phase spectrum of an image contains most of the information for reconstructing the actual image. But looking at the phase spectrum, I see just "noise". Would it be possible to use the phase information as well?

How would you determine, from a fourier spectrum alone, if the source of light is homogeneous or not?

I have been reading about fourier transforms for the last few months, but I feel like I am still missing a huge point when it comes to connecting fourier transformations and image processing. It seems everyone is just using the amplitude information of an image while being in the fourier space, and ignoring phase.

If I wanted to determine, how homogeneous this source of light is, I'd take a picture of it and put it through a fourier transformation, e.g. a DFT.

Why?

The ideal homogeneous area is a rectangle of some constant colour. A typical sensor has a response with well defined low and high limits. Therefore, excluding ceiling effects at those limits, if you were to shoot a gray card you could then obtain the standard deviation of its values.

A narrow standard deviation would point to a "good" level of homogeneity. A wide standard deviation would point to a "bad" level of homogeneity.

If the standard deviation seems too "coarse", too general a quantification, then the next level up (in complexity) would be to evaluate the co-occurence matrix of the image and through that, the inverse difference moment or homogeneity. The homogeneity is high in low contrast images and vice versa.

For more information please see this link

In any case, you would have to calibrate your metric to give meaning to "good" and "bad". Good homogeneity with a web cam is entirely different than good homogeneity with a cooled ccd at the end of a very fast lens (i.e. a lens with a large diameter) for shooting a grey card under the same lighting conditions.

  How would you determine, from a fourier spectrum alone, if the source of light is homogeneous or not?

If you absolutely had to use the Discrete Fourier Transform, then you would again be looking at how is a flat or "rough" surface "expressed" in the Fourier Transform. In general, the Fourier Transform is used for the "detection" of texture, especially in the case of repeating elements. Therefore, absence of texture leaves just an average value, full presence of texture is something like the spectrum of white noise and intermediate cases would have to be judged with something like spectral flatness but adapted to the two dimensional case.

Hope this helps.


Z-Transform of a^|n| I am wanting to compute the Z-transform of $f(n) = a^{|n|}$ . 'a' is a positive constant.

Looking at the transform table, I found that Z-transform for $a^n u(n)$
is available from the tables and is $\frac{Z}{Z-a}$. Where $u(n)$ is the unit step function.

I am trying to decompose $f(n)$ as --
$f(n) = a^n u(n)  +  a^{-n} u(-n) + \delta(n)$

Then using the table to find $Z( f(n) )$ as --
$\frac{Z}{Z-a} + \frac{Z^{-1}}{Z^{-1} - a} + 1$

Can anyone tell me if this approach is correct and/or suggest an alternate way.
Thanks!

I believe the right way to decompose this function is:

$f(n) = a^n u(n)  +  a^{-n} u(-n)$

Then the answer for your question is:

$Z( f(n) ) = \frac{Z}{Z-a} + \frac{Z^{-1}}{Z^{-1} - a}$
EDIT: Looks like it is wrong answer. :(

Ok, let's start from the beggining and try to find right answer step by step

Given sequence is $f(n) = a^{|n|}$ . 'a' is a positive constant.$

1) Let's  decompose it into sum of 2 sequences:

$f(n) = a^n u(n)  +  a^{-n} u(-n-1)$

2) If we define $f_1(n) = a^n u(n)$  and $f_2(n) = a^{-n} u(-n-1)$
Then $Z( f(n) ) = Z( f_1(n) )+Z( f_2(n) )$

3) $Z( f_1(n) ) = \sum_{n=-\infty}^{\infty} a^n u(n) z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n = \frac{1}{1-az^{-1}} = \frac{z}{z-a}  $

4) $Z( f_2(n) ) = \sum_{n=-\infty}^{\infty}  a^{-n} u(-n-1) z^{-n} = \sum_{n=-\infty}^{-1} (az)^{-n} = \sum_{n=1}^{\infty} (az)^{n} = -1+1+\sum_{n=1}^{\infty} (az)^{n} = -1+\sum_{n=0}^{\infty} (az)^{n} = -1 + \frac{1}{1-az} = \frac{-1+az+1}{1-az} = \frac{az}{1-az} $

5) $Z( f(n) ) = Z( f_1(n) )+Z( f_2(n) ) = \frac{z}{z-a} + \frac{az}{1-az} = \frac{z}{z-a} - \frac{z}{z-\frac{1}{a}} $

If there is anything unclear or wrong in this derivation, please comment.

As for wolfram answer - I don't know why it is so. I never used wolfram before.
But it perfoms symbolic computations and you must be sure that your assumptions about range of n (does $n\in[-\infty,\infty]$  or $ n\in[0,\infty]$ ) for example should coinside with wolfram assumptions.


Difference Between MFCC of speech and speaker recognition I am doing my project on speaker recognition. One of my friends is doing his project on speech recognition. Both of us need to calculate the MFCC for feature extraction.

Is the MFCC going to be the same one for both of us? Is there any difference between the algorithm of MFCC for speech and that for speaker recognition?

Both speech recognition and speaker recognition require some set of features to distinguish one speaker (or speech section) from another.

Suppose you have two speakers $S_1$ and $S_2$  and two words $w_1$ and $w_2$. Let's call $u_{ij}$ the utterance of speaker $i$ of word $j$.  Then the MFCC's of the speech are $M(u_{ij})$.

One way to pose the speaker recognition problem is to determine which speaker said a new utterance, $u_{\rm new}$. Assume that the utterance is word $J$. To do this, we just select the speaker such that
$$
I = \arg\min_{\forall i} \lVert u_{iJ} - u_{\rm new}\rVert
$$

How this works will be determined by the inter-speaker variation (i.e. how much the MFCCs of different speakers saying the same word differ).

One way to pose the speech recognition problem is to select:
$$
J = \arg\min_{\forall j} \lVert \bar{u}_{j} - u_{\rm new} \rVert
$$
where $\bar{u}_{j}$ is the average MFCC of all speakers saying word $j$.

How this works will be determined by the inter-word variation (i.e. how much the MFCCs of different words being said by different speakers differ).

So, yes, the MFCCs for both applications are (or can be) the same.


FIR Filter length Upper Bound Say I want to design a low pass filter with stop and pass band ripple constrained to $\epsilon$ and transition region $\delta$, what is the smallest FIR Filter length that suffices? Can you please provide a derivation?

Many people have asked themselves this question. Unfortunately, there is no exact solution, just approximations. One relevant paper is Practical Design Rules for Optimum Finite Impulse Response Low-Pass Digital Filters by Herrmann, Rabiner and Chan, where they basically fit a formula to their data obtained from lots of designs. A simplified formula is [1]

$$M=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\;\Delta\omega}$$

where $M$ is the filter order, and $\Delta\omega=\omega_s-\omega_p$ is the width of the transition band.

Also the Matlab function firpmord.m uses some formula to estimate the necessary filter order, but they also say that the required order is often underestimated.

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


Brownian noise is integration of gaussian or uniform white noise? On this wikipedia page it is written that

  Brown noise can be produced by integrating white noise

What kind of white noise is meant? The Gaussian or uniform white noise?

White noise is noise that has equal (uniform) amplitude across all frequencies. When we say "white" we're talking about the power spectral density (PSD) of the noise.

Saying something like "Gaussian noise" means the statistical properties of any one sample of the noise is distributed Gaussian. You can actually have Cauchy, Poisson, Gaussian etc. distributions that define any one sample of the noise. For any distribution, however, if the power spectral density of the noise is uniform, that noise is white.

So in the integral defining the Weiner process:
$$W(t) = \int_0^t\frac{dW(\tau)}{d\tau}{d\tau}$$

$W(t)$ has samples with Cauchy, Poisson, Gaussian, etc. distributions as long as the power spectral density is uniform, making it white noise.


Representation of frame image pixels In paper Human identification using gaits Section 3.2, the authors explain how to represent a 3D image (x,y,t) into 2D image with (x,t) ; (y,t) denoted 1D signals. I have a conceptual doubt which is in my case I have a matrix (x,y) represent the centroid of the motion- a pixel point from one frame to the next. So, if we want to do some kind of processing like finding the fourier transform, then is it theoretically correct to treat x- coordinate and y coordinate separately?Do we need to perform some trasnformations?I did not understand how the authors are representing an image pixel as 1D signal.Questions are :

It is not clear how the authors have converted to four 1D signals. How to use/represent the image pixel coordinates as 1D signals as done in the paper for further analysis?Should I perform autocorrelation between x and y to get 1D vector?
Phase space plot = it is plot of 2 variables. So, in this case, will it be a plot of x(t) vs y(t) or x(t-1) vs x(t) or the future_predicted_values for x cordinate vs the current measured values of x coordinate?Or the plot of features?
Do I need to smooth the data?If so then how to use bsplines as I found that this is popularly used for smoothening.

Example data points representing pixels (x,y) are

x      y
160    210
160    230
160    250
170    83

and so on for 1000 frames.

Can somebody please help in explaining the concept and the way to proceed with image signals so that a phase space plot can be created.

The paper that you are quoting assumes that the motions along the different coordinates are independant functions, i.e., $$P(t) = \left(x(t),y(t)\right) = (f_1(t),f_2(t)) \neq f(x,y)(t),$$ where $P$ stand for position.
The various 1D functions are obtained in this way.

A bit of smoothing is generally helpful.
However, when to comes to gait recognition, you will probably have to test different smoothing procedures: my intuition is that gait recognition relies on high frequency / low amplitude clues that can easily be cancelled out (but I may be wrong for this).


Background Removal and Object Tracking - Track Fish Jumping Out of the Water 0 down vote favorite

I am totally new to Image processing. I just started my Master thesis in the area of Computer vision and Machine learning. My background is Informatics. Now, my first task is to register the images of some fishes(Image registration) as they come out of water. I have got stream of images and I want to come up with a model of the fish by aligning the temporally different images of a fish. As, I can understand, I will have to fist of all remove the background and water from the images to work on just fishes, am I right?

Can anyone give me a brief idea of how I should go forward or the things I should read first before I can understand the things. For example, should I read the basics of Image Processing, Feature detection, Image segmentation....? and regarding the programming language I can use - one where I can find good libraries, forums, other help.

I would be really grateful, if anyone can help. Thanks

I`d recommend you to go through any of this courses (they include lectures, references and task for labs)

https://cw.felk.cvut.cz/wiki/courses/ae4m33mpv/start - course from Czech Technical University
http://cs.brown.edu/courses/cs143/ from Brown University


Equivalence of Convolution Through Polynomial and Matrix Products I am trying to understand the equivalence of convolution through polynomial (coefficient) multiplication, Fast Fourier Transforms, and matrix multiplication. I am using octave with the octave-signal package, and part of this question may possibly related to their implementation (which I believe would be the same for MATLAB).

So I create my original signal with length $N$, and my boxcar window length is $n$.

N = 48;
n = 6;

I use $\sin$ as the input function for this example.

dt = 2*pi/N;
t = [-N/2:N/2-1]'*dt;
f = sin(t);

Convolution by polynomial multiplication ($g = h \star f$):

h = ones(n,1);
g_conv = conv(h,f);
g_conv_same = conv(f,h,'same');

Convolution by FFT and FFT with zero-padding of the input (also $g = h \star f$):

g_ft = real(ifft(fft(h,N).*fft(f)));

f_padded = [f; zeros(n-1,1)]; % append zeros
g_ft_padded = real(ifft(fft(h,N+n-1).*fft(f_padded)));

Convolution by matrix multiplication ($\mathbf{g} = \mathbf{A} \mathbf{f}$):

A = convmtx(h,N);
g_mtx = A*f;

$\mathbf{A}$ looks like this:

Find indices to trim arrays:

idx_A = find(sum(A,2) == n);
idx_t = 1:length(idx_A);

Plot results (the signals convolved by different methods):

plot(t(idx_t),g_mtx(idx_A),'k','linewidth',3);
line(t-5*dt,g_ft,'color','g','linewidth',2);       %-> offset by 5*dt
line(t-7*dt,g_ft_padded(1:length(t)),'color','r','linestyle','--'); %-> offset by 7*dt
line(t(idx_t),g_conv(idx_A),'color','b','marker','o');
line(t-2*dt,g_conv_same,'color','b','marker','+'); %-> offset by 2*dt
legend({'matrix','FFT','FFT (padded)','conv','conv (same shape)'},...
       'location','southeast')
xlabel('x value')
ylabel('y value')

which looks like this:

I have several questions at this point:

Why are the results for the FFT methods offset by several units of dt?
If I pad the FFT results such that the input vector is no longer periodic (red dotted line), why do I miss part of the convolved function (on the right side) in the result?
How do I treat the ends for the FFT convolution so that the results are aligned with the polynomial and matrix multiplication results?

Thanks in advance. P.S. if there is a good textbook reference for these details I would appreciate your guidance.

(I have edited for clarification.)

Convolution and polynomial multiplication are equivalent by definition. Offsets are usually introduced through indexing. The treatment of ends is important.


Matlab Neural Networks: input size vs. number of inputs I'm trying to do time-series prediction with dynamic recurrent networks in Matlab (will try at least NARX and LRN architectures, if not one or two others). The question I have is, what is the functional difference between giving a network multiple inputs, versus giving it an input vector with multiple elements?

My inputs are values from multiple sensors. Can I create a single matrix (say, a 3xN matrix for 3 sensors that recorded N time steps), run 'con2seq' and use that as a single input to the network? What is the functional difference between doing that versus formatting it as 3 different sequential inputs?

As I understand it, each value (each row within the 3xN matrix, or alternatively each of the 3 single-row inputs) will get its own weight to every layer 1 neuron. Is that correct? What's the difference between the two options, then? Is the only purpose of creating multiple inputs for the case where you want to send some inputs to one subset of layers and other inputs to a different subset of layers?

After reading the first third of the Neural Network Toolbox User Guide, I looked around on here, stackexchange, and google with luck. Everyone says, "net.inputs{K}.size is the dimension of input K, while net.numInputs is the number of inputs," which is tautological more than clarifying.

Thanks much for your help!

This question is better fit in stackoverflow or cross validated site, yet I am trying to illustrate a simple example.

Suppose you have three measurements $a$, $b$, and $c$, with each measurement only $1$ time reading, the hidden layer node and output are scalar. Compose an input vector is shown below:

        w1
(a,b,c) => (H1)-------(O1) => Result

The value at output layer is $ \alpha ((a+b+c)w_1)$, where $\alpha$ is the activation function.

When you treat the measurements separately at the input layer, the neural network is:

  w1
a =>
     \
  w2
b => (H1)----------(O1) => Result
     /
  w3
c =>

and the output is $ \alpha (aw_1 + bw_2 + cw_3)$. Now you are dealing with three independent neurons, which indicates that updating one of the weights won't change the 'learned' weights on other input nodes.

As far as I know, a single neuron usually take in scalars rather than vector values (for the case of one time reading), or take in vectors rather than matrix arrays (for the case of a series of time readings).


Difference between filtering using Matlab and filtering using C code generated by Matlab hi，has anyone used Matlab coder to generate C code for a filter? Here is the thing: I tried to filter the same audio signal using Hamming window-based FIR bandpass filter, the filter is designed by Matlab DSP tool box. Then I used Matlab coder to generate the corresponding C code for deploying the filter on the MCU. But I found that the filtering in Matlab and that in MCU has differences, I think the result from Matlab is correct.
Is the filter C code generated by Matlab not reliable? Maybe I should try ARM DSP library

You have differences in the transient, but the steady-state looks the same therefore it is probably a memory initialization problem in the C-code

Here are a few things to check :

Are you sure your input is exactly the same?  Make sure your input buffer is properly initialized in the C-code.

You probably have some kind of data buffer to store the previous samples. Again make sure the data is properly initialized. If the data buffer is not properly initialized you will likely have garbage in the data buffer instead of having '0' in the data buffer which would explain why the transient is different.

The Matlab coder generates 3 user functions, Initialize, Step and Terminate. Make sure you call Initialize once, before make any call to Step.


Power Spectrum analysis arbitrary length signals I have signals belonging to diseased and healthy samples recorded from multiple electrodes. These signals are taken monthly and are of arbitrary length (say 4-5 mins). I am planning to do a power spectrum analysis (comparison between healthy and diseased). Now to do this two methods come to mind.

Take the normalized power spectrum of the 4/5 min signal as it is and average it over the electrodes.
Divide the 4-5 min recordings to smaller (say 5 seconds) segments, obtain the normalized power spectrum of each segment, average it over all the segments, and finally average it over the electrodes.

which of the following methods would give a better analysis?

The answer depends on the statistical characterization of the signal. If the random EEG signal is stationary over when seen over 5 second cycles, then you are better off averaging over 5 seconds cycles. If not, then averaging over these 5 seconds cycles might not yield the best estimate. In that case better off averaging over the entire duration.


How do you find the phase of the DFT of a discrete signal? My question is similar to this one and this one, but neither answers address my concern.

Suppose you have a signal x(n) = {-1,2,-3,2,-1} where we can assume x(0) = -3.  So this is an even signal.

How do you find the phase of the DFT of x(n)?

Here is a hint that will help you:

The DFT is cyclical in time and in frequency. For the sequence given by

$$x(n) = [-1,2,-3,2,-1]$$

With x(0) = -3 would be solved using the standard DFT equation that starts at n=0 using

$$x(n) = [-3, 2, -1,-1, 2]$$

From that you can solve for the DFT and then determine easily for each result what it's phase is.


How to detect overlapping peaks in FFT? I am working with a modular spectral processing STFT toolkit I have written in C.  I often need to carry spectral peak detection and analisys  (frequency estimation, peak displacement, etc) but I am into troubles when there are two or more peaks merging in a same apparent peak, i.e many overlapping peaks resulting in a single peak detectable as a single local maximum. I am aware that this problem has been addressed copiously in literature, in particular in chromatography or mass spectrometry. A commonly suggested approach is analysing the 2nd derivative looking for valleys, but this trick is not viable when the spectral resolution is insufficient, as when you have perhaps two or three peaks sharing few bins. Also I found many cases where a very low frequency sawtooth spectrum results in a smooth and large single peak, for which even the 2nd derivative of magnitude looks monotonous. So what to do ? I suspect that proper analisys of phases other than of magnitudes alone could give more clues of the actual underlieing peak structure but I am not a maths guru despite skilled enough. I would like to be able, once detected a whole peak "body", to decompose in its constituting overlapping peaks (if more than one) and possibly isolate them, or at least having a good estimate. Is there a way to do that, which is not too computationally intensive and which can be done in real time ? Thanks for any relevant answers!

A longer FFT on a longer window of actual data will better separate spectral peaks that share bins in a too short FFT.


obtaining "4 Hz modulation energy" Could anybody tell me a method to calculate the "4 Hz modulation energy" for a speech signal?

Thanks.

EDIT

I get some details : 4 Hz modulation energy: Speech has a characteristic energy modulation peak around the 4 Hz syllabic rate [3]. We use a portion of the MFCC algorithm [4] to convert the audio signal into 40 perceptual channels. We extract the energy in each band, bandpass filter each channel with a second order filter with a center frequency of 4 Hz, then calculate the short-term energy by squaring and smoothing the result. We normalize each channel's 4 Hz energy by the overall channel energy in the frame, and sum the result from all channels. Speech tends to have more modulation energy at 4Hz than music does.

The link in my comment suggests the following:

  I would recommend first extracting the envelope by using either a halfwave
  rectification (i.e., replace all the negative values in the time waveform
  with zeros) or a Hilbert Transform and then lowpass filtering the waveform
  at around 50 Hz (the lowpass is optional if you only care about the 4 Hz
  component).   Once you have the envelope, simply do an fft (making sure
  that your frequency resolution is at least 1 Hz) and look for the energy
  around 4 Hz.  If any of these steps don't seem obvious too you, let me know
  and I can send you some Matlab code.

So, what that means is you need to take your speech $s[n]$ and find the envelope:

$$
e[n] = \left | s[n] + j{\bf H}\left[s[n]\right] \right |
$$

where ${\bf H}$ is the Hilbert transform of your speech signal.

Then, take the FFT and look at the bin at 4Hz.


Looking for an algorithm that computes "Similarity" between two signals I'm new in this forum! ... and I'm not an expert in signal processing! ... So I hope to find help here:-)
I was wondering if someone could point me to an algorithm/technique that is used to compare time dependent signals. Ideally, this hypothetical algorithm would take in 2 signals as inputs and return a number that would be the percentage similarity between the signals.
In other terms I’m looking for a technique that allows to establish a degree of “similarity” between the two signals even if there is a delay or a shift between them. Obviously, if this technique could be implemented in a recursive way or with a small number of point ...it would be great! :-)

Thanks

Coherence? As in: http://en.wikipedia.org/wiki/Coherence_(signal_processing). By comparing the signal spectrums, you get real number numbers between 0 and 1, which each number is a measure of the similarity of the two signals at each frequency.  But it might be difficult for a beginner to understand and program.  You might look for simpler explanations and/or examples by doing a Google search for 'signal coherence', which results in about 17 million hits.


Flat vs frequency selective fading in MC-CDMA? I am trying to simulate an MC-CDMA system. I have NXK  matrix $S$ which consists of subcarriers ($s_n,_k$ represents subcarrier $n$ of user $k$) and KX1 vector $X$ consists of user symbols. So my received signal is $Y=SX+N$. To include the channel response I am confused in flat fading and frequency selective fading.

$\:$ For  flat fading  i.e., coherence bandwidth $>$ signal bandwidth,  and I have a 3 tap channel impulse response in time domain. So in frequency domain I have NX1 vector $h_k$ for user $k$.

Now to include the channel fading for user $k$,

Should I do  S(:,k) .* H(:,k) %% H is channel matrix i.e., each subcarrier is multiplied by a different channel coefficient while coeffecients have nearly same magnitude, or
S(:,k)*h(k) i.e., multiplying whole sequence by a single coefficient. and what that h(k) value should be ?

For frequency selective fading does only the values of channel coefficient will change or I will have multiple taps for each subcarrier and will need to convolve in frequency domain?

It would help if you mentioned where your references for MC-CDMA and your channel model come from.

The short answer is that it depends on the level of abstraction you want. For a flat channel, you would multiply by a complex-valued gain. For a frequency-selective channel, I would typically form the time domain signal and convolve with an impulse response (generated from some channel model that reflects your target environment). Perhaps if you provide us with a reference, I can provide additional assistance.


identify a specific table in document image I have some images of answered tests. I am working on recognizing the answers for the multiple choice questions on each test. A example:

I think the problem can be separate into two sub-problems:
    1.) Identify and segment the multiple choice answers from the rest of the test.  2.) Classify and recognize the answer as one of the letters in A, B, C, D, E.

For the second sub-problem, there exists libraries like Tesseract, also I can train a neural network to identify among the five letters which is not hard.

However, I am not sure how to approach the first sub-problem. The most important part in identifying the multiple choice answers would be identifying the region that contain the multiple choice answers (the red box region in the example). Since multiple choices answers are typically written in a table that has two rows and number_of_questions + 1 columns. My first step is to identify this table in the image.

What's the best way to identify such table? I am currently thinking about training a neural network to identify the table. Should I use template matching or other methods instead? Which one is better?

The simplest approach to this would be template matching via two-dimensional normalised cross correlation.

What this basically does is looking for a template: a small image representing the response letters, inside a (usually bigger) image. The result that is returned from this operation looks like an image but it essentially has peak values where the template matches the image.

This approach transferred to this problem means that one cross correlation would have to be carried out for each possible answer and each resulting image examined for "peaks".

Of course, this assumes that the letters appear at a single scale throughout the document. Otherwise, cross correlations would have to be applied (via an appropriate method) at different scales too.

There are two components to be carefully managed here:

Quality and scale of the image:

This can be roughly controlled by instructing users to align the top left corner of their camera's field of view with the paper. This will constrain the distance from the image but it still depends on the lens of the mobile phone.
Quality and scale of the template:

This depends on #1 and the hand-writing of the users. But on average, roughly equal size fonts can be used (as further below).

Here is an example based on Python, matplotlib and scikit-image:

from matplotlib import imshow, show
from skimage.io import imread
from skimage.feature import match_template

#Load the image
Q = imread("tst4.jpg")

#Load the template, here it is the letter A from the Sans font, approximately at the same size as the hand-written provided
A = imread("letter2.png")

U = match_template(Q, A, pad_input = True)

#Show the result
imshow(U);show()

(Please note, a more elaborate version of that is available via this link)

The result U looks like this:

For a Q that looks like this:

And an A that looks like this:

In the U image, please note that the red colour corresponds to high peaks (very close to 1.0). Just putting a threshold on the pixel values of U will return the centers of the 'A' symbols.

Of course, other formations seem to have high correlation as well in that image. These are "false-positives", that is, places where the algorithm THINKS there is an "A" but in reality there is not.

There are two ways by which these can be reduced (and increase the system's reliability):

Apply the template matching over a smaller area (this is now possible if the top left corner is roughly aligned with the edge of the paper)
Add more features to the classification operation than just the peak of the cross-correlation between the template and the image to make it more robust. This is probably the next level in complexity and more information can be provided if needed.

Hope this helps.


How to get the area(Pixels count) inside ring-like shape? I have the following image, and I do like to count the pixels inside the ring to get the area.

I did some morphological operations as a kind of post-processing to make the image as much as possible with clear smooth edges.
I tried to do that in different ways as you can see in the code down, but none of them was optimal.
Can you please advice me how to count the pixels inner area of the circle?
Note: some pixels inside are not totally black, they are with low intensity, that's why I was trying to do Otsu thresholding.
Thanks in advance

import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread, imsave
# import scipy.ndimage as ndi
from skimage import morphology, filters, feature

seg = np.squeeze(imread('prediction.png')[...,:1])
# meijering alpha=None,
# rem2 = morphology.remove_small_objects(seg, 4)
resf = filters.meijering(seg, sigmas=range(1, 3, 1),  black_ridges=False)

sobel = filters.sobel(resf)
# diam = morphology.diameter_closing(sobel, 64, connectivity=2)
gaussian = filters.gaussian(sobel, sigma= 1)
val = filters.threshold_otsu(gaussian)
resth = gaussian < val

# Morphology
SE = morphology.diamond(2)
# SE = np.ones((3,3))
# SE = morphology.disk(2)
# SE = square(7)
# SE = rectangle(3,3)
# SE = octagon(3, 3)

erosion  = morphology.binary_erosion( resth, SE).astype(np.uint8)
dilation = morphology.binary_dilation(resth, SE).astype(np.uint8)
opening  = morphology.binary_opening( resth, SE).astype(np.uint8)
closing  = morphology.binary_closing( resth, SE).astype(np.uint8)
#thinner = morphology.thin(erosion, max_iter=4)

rem  = morphology.remove_small_holes(resth, 2)

# entropy  = filters.rank.entropy(resth, SE)
# print(seg.shape)

plt.figure(num='PProc')
# 1
plt.subplot('335')
plt.imshow(rem,cmap='gray')
plt.title('rem')
plt.axis('off')
# 2
plt.subplot('336')
plt.imshow(dilation,cmap='gray')
plt.title('dilation')
plt.axis('off')
# 3
plt.subplot('337')
plt.imshow(opening,cmap='gray')
plt.title('opening')
plt.axis('off')
# 4
plt.subplot('338')
plt.imshow(closing,cmap='gray')
plt.title('closing')
plt.axis('off')
# 5
plt.subplot('332')
plt.imshow(seg,cmap='gray')
plt.title('segmented')
plt.axis('off')
# 6
plt.subplot('333')
plt.imshow(resf,cmap='gray')
plt.title('meijering')
plt.axis('off')
# 7
# 8
plt.subplot('334')
plt.imshow(resth,cmap='gray')
plt.title('threshold_otsu')
plt.axis('off')
# 9
plt.subplot('339')
plt.imshow(erosion,cmap='gray')
plt.title('erosion')
plt.axis('off')
#
plt.show()

The answer is down from the question on StackOverflow website, thanks for Paul Brodersen the credit of the answer is for him.

#!/usr/bin/env python
"""
Determine areas in image of ring.

SO: https://stackoverflow.com/q/61681565/2912349
"""
import numpy as np
import matplotlib.pyplot as plt

from skimage.io import imread
from skimage.filters import threshold_otsu
from skimage.measure import label, regionprops
from skimage.color import label2rgb

if __name__ == '__main__':

    raw = imread('prediction.png', as_gray=True)
    threshold = threshold_otsu(raw)
    thresholded = raw > threshold
    # Label by default assumes that zeros correspond to "background".
    # However, we actually want the background pixels in the center of the ring,
    # so we have to "disable" that feature.
    labeled = label(thresholded, background=2)
    overlay = label2rgb(labeled)

    fig, axes = plt.subplots(1, 3)
    axes[0].imshow(raw, cmap='gray')
    axes[1].imshow(thresholded, cmap='gray')
    axes[2].imshow(overlay)

    convex_areas = []
    areas = []
    for properties in regionprops(labeled):
        areas.append(properties.area)
        convex_areas.append(properties.convex_area)

    # take the area with the smallest convex_area
    idx = np.argmin(convex_areas)
    area_of_interest = areas[idx]
    print(f"My area of interest has {area_of_interest} pixels.")
    # My area of interest has 714 pixels.
    plt.show()
```


Proof of the convolution property of Fourier Series in continuous time I am facing problem in understanding the proof of Convolution property of Fourier Series (FS) in continuous time CT;
that is:  $$\mathrm{FS} \big\{x_1(t)\star x_2(t)\big\}=T\sum_{n=-\infty}^{\infty}c_n d_ne^{j \Omega_0 n t}$$ where $x_1(t)$ and $x_2(t)$ are periodic with $T$, and $c_n$ and $d_n$ are the Fourier series coefficient of $x_1(t)$ and $x_2(t)$.

Below I am showing the proof and my confusion,

PROOF:
\begin{align}
x_1(t)\star x_2(t) &= \int_T x_1(\tau)x_2(t-\tau)d\tau\\
\therefore \mathrm{FS}\big\{x_1(t)\star x_2(t)\big\}&=\int_T FS[x_1(\tau)] FS[x_2(t-\tau)]d\tau\\ &=\int_T \sum_{n=-\infty}^{\infty}c_n e^{j \Omega_0 n \tau} \sum_{p=-\infty}^{\infty}d_p e^{j \Omega_0 p (t- \tau)} d\tau\\
&=\sum_{n=-\infty}^{\infty} \sum_{p=-\infty}^{\infty} c_n d_p \int_T e^{j \Omega_0 (n-p) \tau} e^{j \Omega_0 p t} d\tau\\
&=\sum_{n=-\infty}^{\infty} c_n d_n T e^{j \Omega_0 n t}
\end{align}
CONFUSION:
Now I am confused  how we are getting
$$ \sum_{n=-\infty}^{\infty} c_n d_n T e^{j \Omega_0 n t}  \quad\text{from}
\quad\sum_{n=-\infty}^{\infty} \sum_{p=-\infty}^{\infty} c_n d_p \int_T e^{j \Omega_0 (n-p) \tau} e^{j \Omega_0 p t} d\tau \quad ?$$

Please explain.

HINT:

Show that with $\Omega_0=2\pi/T$

$$\int_T e^{j\Omega_0(n-p)\tau} d\tau=\begin{cases}T,&\quad n=p\\0,&\quad\text{otherwise}\end{cases}$$

So the sum over $p$ reduces to a single element $Td_ne^{j\Omega_0 nt}$, and the result follows.


How to find first 25 coefficients of DCT in Matlab? I need first $25 \times 25$ coefficients of DCT to encrypt an image in Matlab.

How do I obtain them ?

I used this code to compute DCT of an image :

I = imread('cameraman.tif');
I = im2double(I);
T = dctmtx(8);
dct = @(block_struct) T * block_struct.data * T';
B = blockproc(I,[8 8],dct);
mask = [1   1   1   1   0   0   0   0
    1   1   1   0   0   0   0   0
    1   1   0   0   0   0   0   0
    1   0   0   0   0   0   0   0
    0   0   0   0   0   0   0   0
    0   0   0   0   0   0   0   0
    0   0   0   0   0   0   0   0
    0   0   0   0   0   0   0   0];
B2 = blockproc(B,[8 8],@(block_struct) mask .* block_struct.data);
invdct = @(block_struct) T' * block_struct.data * T;
I2 = blockproc(B2,[8 8],invdct);
imshow(I), figure, imshow(I2)

The statement in the paper I read is :

  Basic constructions: Let $V$ be a feature space of dimension $d$ (for
  instance, the space $R^d$ corresponding to the entries in the top left
  corner of the DCT of a digital image excluding the DC component) and
  consider an original feature $φ = (φ_1, . . . , φ_d)$.

The paper gives also some practical details like this :

  We choose as a feature space the upper left $25 \times 25$ coefficients of
  the DCT (excluding the DC component).

In JPEG compression, DCT coefficients are generally parsed not as sub-blocks, but in an order:

  corresponding to the entries in the top left corner

This corresponds to a zigzag indexing:

This can be obtained by the following code. Here I use the $8\times8$ JPEG quantization table, and extracted the $5\times5$ top-left coefficients. On a light note, the bijection between the $2$D coordinates and the zigzag  $1$D coordinate can be obtained by Cantor pairing polynomials (Matlab implementation).

% ExtractComponent
nCoomponent = 25;
quantTable = [    16,  11,  10,  16,  24,  40,  51,  61;
    12,  12,  14,  19,  26,  58,  60,  55;
    14,  13,  16,  24,  40,  57,  69,  56;
    14,  17,  22,  29,  51,  87,  80,  62;
    18,  22,  37,  56,  68, 109, 103,  77;
    24,  35,  55,  64,  81, 104, 113,  92;
    49,  64,  78,  87, 103, 121, 120, 101;
    72,  92,  95,  98, 112, 100, 103,  99];

index = zigzag(length(quantTable));
indexCrop = index(2:26);
quantTableCrop = zeros(size(quantTable));
quantTableCrop(indexCrop) = quantTable(indexCrop);

figure(1)
image(quantTable);
figure(2)
image(quantTableCrop);

and the function

function index = zigzag(size)

% ZIGZAG(size) returns the zigzag order index
% for a size x size square matrix

    index = 1:size*size;
    [j,i] =  meshgrid(1:size);
    i = i(:);
    j = j(:);
    iplusj = i + j;
    even_ind = (rem(iplusj,2) == 0);
    odd_ind = ~even_ind;
    j(even_ind) =   j(even_ind);
    j(odd_ind) =    i(odd_ind);
    [nouse,tags] = sortrows([iplusj,j]);
    index = index(tags);
end


Jamming resistance in spread spectrum systems In direct sequence spread spectrum, which is the basis of CDMA systems, the baseband signal is multiplied with a pseudo noise (PN) sequence which has a high chip rate, thus the resulting modulated signal will have high rate and thus larger bandwidth. But the power of the original baseband signal and the power of spreaded signal resulting from the multiplication process is the same.
Then, why this spreading process provides jamming resistance property?

Viewed in time domain, they have the same power. Does the jammer depend on the power value to detect the signal or not? I read in a reference that the spreaded signal appears as noise for the unintended receiver. But this only when viewed in frequency domain. Does the jammer depends on frequency domain measurements for detecting the signal, not on time domain power measurements?

In order to work well, a jamming signal must be well-correlated with the signal it is trying to jam.  In the case of sine-wave modulated signals, it's easy to detect the signal one wishes to jam, to infer that it's narrowband (and, hence, riding on a real or suppressed sine-wave carrier), to infer its bandwidth, and from there to choose a suitable jamming signal.

To jam a CDMA spread spectrum signal you'd either need to overwhelm it with an uncorrelated signal (random, single tone, pseudo-random noise) or you'd have to find it, find its PN sequence, synchronize to its PN sequence, and then in most cases adjust your synchronization to match the time offset between your reception of the signal and the receiver's, then you'd have something to send -- but if I were trying to receive the signal, I could defeat you by having multiple receivers, located more than a wavelength of the chipping rate away from each other, each trying to receive the signal.  One receiver may get the jamming signal and the intended signal on top of each other, but unless your jamming signal is co-located with the intended transmission, other receivers would be able to pick out the good from the bad.

But keep in mind that civilian CDMA is intended to minimize unintentional jamming, and does pretty well.  It's not intended to defeat malicious jamming.

I suspect that there are smart military types out there that know a lot more than I do about jamming CDMA -- I'm just making this up as I go from first principles.  I also suspect that military systems will use something more complicated than civilian CDMA, for just that reason.


What does standard deviation of noise 1 mean? In theory of noise,Typically, the standard deviation of noise has 1 and mean has 0.
I think that the reason of noise mean has 0 that we can assume that all noise signal go to zero when we sum it all.
But I can't understand standard deviation of noise has 1 in image noise.
They are always 0~255. There is no negative value.

What does zero-mean noise mean?

Especially, but I have got an idea as following. I can't catch quite a bit. So would you please let me know how to illustrate following answer?

"for an image intensity which is always non-negative, you can simply modify your added N(0,1) standard gaussian noise in 3 steps: 1-scale it so that your noise has whatever std you wish. 2-add whatever proper mean you need such as 128. 3-clip the noise for below zero and above 255 samples. This is simple bt It will help you get what you need at least initially."

Your question is too basic, and it gives the impression that you haven't done enough study of the subject on your own. I'll provide a few hints:

To understand noise, you first need to specify its probability distribution. I assume you're refering to Gaussian noise in your question, but not all noise is Gaussian. In particular, in digital capture of an image, there are several kinds of noise present.
Gaussian noise is specified by its mean and its variance ($\sigma_n^2$) or its standard deviation ($\sigma_n$). The physical interpretation of the variance is that it is related to the noise power. There is nothing special about $\sigma_n=1$.
Most computer random generators (for instance, Matlab's randn) generate samples of a Gaussian random variable with zero mean and variance 1. It is easy to convert this random variable to another Gaussian random variable with any specified mean and variance.
You need to understand the effects of clipping Gaussian noise. Think about a sensor photosite being pulled or pushed beyond its dynamic range by the noise. However, as I said in another answer, the resulting noise is no longer Gaussian.


Show Equivalence Between Multiplication in Time Domain to Convolution in Frequency Domain My goal is to compute the Fourier of the product between two discrete time signals, y1 and y2. This can be done by computing the convolution between the fourier transform of y1, f1 and the fourier transform of y2, f2 (This is what I understood from the Wikipedia Page).

I tried the following code in MATLAB, it shows that the convolution between y1 and y2 is equal to ifft(f1.*f2) as shown on the Wikipedia page, the code worked perfectly:

y1=sin(x);
y2=sin(3*x);
n=length(y1);
convol=conv(y1,y2);
f1=fft(y1,n*2-1);
f2=fft(y2,n*2-1);
plot(ifft(f1.*f2))
hold on
plot(convol,'x')

Now I want to be able to compute the fourier transform of the product of y1 by y2 using the convolution theorem, this can be done by computing the convolution between f1 and f2, so I wrote the following code:

x=0:0.01:10;
y1=sin(x);
y2=sin(3*x);
n=length(y1);
f1=fft(y1,n*2-1);
f2=fft(y2,n*2-1);
convol=conv(f1,f2);
ffty1y2=fft(y1.*y2,length(convol));
plot(abs(ffty1y2))
hold on
plot(abs(convol)/length(convol),'x')

As you can see the results do not match (neither in location of the peaks nor in their values), what could be wrong here? How do I fix this?

Let's assume you have 2 signals: vX and vY.
So:

clear();

numSamplesX = length(vX);
numSamplesY = length(vY);

numSamplesConv = numSamplesX + numSamplesY - 1;

vTimeDomainConv = conv(vX, vY);

vFrequencyDomainConv = ifft(fft(vX, numSamplesConv) .* fft(vY, numSamplesConv), 'symmetric');

max(abs(vTimeDomainConv - vFrequencyDomainConv)) %<! Should be < 1e-12

This is the Convolution Theorem for Discrete Signals to show convolution in time domain is equivalent to element wise multiplication in frequency domain.

If you want to show element wise multiplication in time domain can be done using the convolution in frequency domain you need to either interpolate the time domain signal to length of linear convolution or just use circular convolution in frequency domain:

clear();

numSamples = 6;

vX = randn(numSamples, 1);
vY = randn(numSamples, 1);

vXY = vX .* vY; %<! Time Domain Mulaiplication
vDxy = fft(vXY); %<! The DFT of the Time Domain Multiplication

vDx = fft(vX);
vDy = fft(vY);

% Applying circular convolution
vXYC = cconv(vDx, vDy, numSamples) / numSamples; %<! Normalization

% Comparing the result in Frequency Domain
max(abs(vXYC - vDxy)) %<! Should be < 1e-12

% Time Domain
vXYFromFrequency = ifft(vXYC, 'symmetric');

% Comparing the result in Time Domain
max(abs(vXY - vXYFromFrequency)) %<! Should be < 1e-12


Fourier transform of an integrator filter I have to find the Fourier transform , and $y(t)$ of an $ x(t) = e^{- \frac {t}{T} } u(t) $ that passes into a integrator filter. I know that $ Y(f) = X(f) H(f) $ so I first calculate the Fourier transform of $x(t)$ and I found that $ X(f)= \frac{T}{1+ i2\pi fT} $. Now I don’t know how to calculate the step response of the integrator filter. After, using inverse transform I should obtain $y(t)$.

As you know, we don't solve homework problems here, but you can get some hints. I don't know if you're really required to solve the problem in the frequency domain, but I would suggest to compute the output of the integrator in the time domain:

$$y(t)=\int_{-\infty}^{\infty}x(\tau)u(t-\tau)d\tau=\int_{-\infty}^{\infty}e^{-\tau/T}u(\tau)u(t-\tau)d\tau\tag{1}$$

This integral turns into something very simple after noting how the infinite integration limits change to finite values when taking into account the two step functions. I'm sure you can continue from here.

If you really want to know the Fourier transform of the unit step $u(t)$, it can be found in most text books on signal processing, and in most Fourier transform tables:

$$\mathcal{F}\{u(t)\}=U(f)=\frac12\delta(f)+\frac{1}{j2\pi f}\tag{2}$$

It's not straightforward to derive this result yourself, so I don't think that this was expected of you.

EDIT:

If you solve the problem in the frequency domain, then with $(2)$ and with $X(f)=\frac{T}{1+j2\pi fT}$ you obtain

$$\begin{align}Y(f)&=X(f)U(f)\\&=\frac{T}{1+j2\pi fT}\left[\frac12\delta(f)+\frac{1}{j2\pi f}\right]\\&=\frac{T}{2}\delta(f)+\frac{T}{1+j2\pi fT}\cdot\frac{1}{j2\pi f}\tag{3}\end{align}$$

where the last equality follows from the fact that $X(f)\delta(f)=X(0)\delta(f)$ is $X(f)$ is continuous at $f=0$.


Showing a system is always controllable? I need to show that the following system is always controllable:

\begin{align}A &= \begin{bmatrix}
-\alpha_1I_{k\times k}& -\alpha_2I_{k\times k}& \cdots &-\alpha_{n-1}I_{k\times k}&-\alpha_nI_{k\times k}\\
I_{k\times k}&0_{k\times k}&\cdots&0_{k\times k}&0_{k\times k}\\
0_{k\times k}&I_{k\times k}&\cdots&0_{k\times k}&0_{k\times k}\\
\vdots & \vdots & \ddots&
\vdots&\vdots\\
0_{k\times k}&0_{k\times k}&\cdots&I_{k\times k}&0_{k\times k}\end{bmatrix}_{nk\times nk}\\
B&=\begin{bmatrix}
I_{k\times k}\\
0_{k\times k}\\
\vdots\\
0_{k\times k}\end{bmatrix}_{nk\times n} \quad C=\begin{bmatrix}N_1&N_2&\cdots&N_n\end{bmatrix}_{m\times nk}.
\end{align}

Now, I think this is actually the controllable canonical form, but I'm really confused about how to show this is always controllable. Can anyone help? Sorry if this is a dumb question.

There's a very simple way to check controllability, indeed if you define the reachability matrix
$$
R = \begin{pmatrix}B & AB & \dots & A^{n-1}B\end{pmatrix}
$$
then the reachable subspace is the image of R. Hence to check complete controllability you just have to check that $R$ is full rank.
First, I think there's an error in the question, $B$ should be $nk\times k$ as you defined it.
Anyway, in your case $R$ will be a $nk\times nk$ matrix, the first $k$ columns are given by the columns of $B$, which are linearly independent by definition of $B$ (in fact you have the identity as the first block), thus $rankR\ge k$. The second $k$ columns are given by $AB$ which has the form
$$
AB = \begin{pmatrix}\star\\I_k\\0_{k\times k}\\\vdots\\0_{k\times k}\end{pmatrix}
$$
where $\star$ means a block you dont care about. Note that the columns of the matrix
$$
\begin{pmatrix}B & AB\end{pmatrix} = \begin{pmatrix}I_k & \star \\ 0_{k\times k} & I_k \\ 0_{k\times k} & 0_{k\times k} \\ \vdots & \vdots \\ 0_{k\times k} & 0_{k\times k}\end{pmatrix}
$$
are linearly independent reguardless of $\star$. Hence $rankR\ge 2k$. In the same way, for the term $A^k B$ you get the form
$$
A^k B = \begin{pmatrix}\star_1 \\ \star_2 \\ \vdots\\ \star_{k}\\ I_k\\0_{k\times k}\\ \vdots \\ 0_{k\times k}\end{pmatrix}
$$
where the Identity block is at position $k+1$. hence for any $0\ge k<j\ge n-1$ $A^{k}B$ and $A^jB$ have linearly independent columns, hence you obtain thart $R$ is full rank an the system is completely controllable


Biquad cookbook formula for broadcast FM de-emphasis? Can a simple IIR biquad filter be used as an audio de-emphasis filter for broadcast FM radio?  If so, how does one translate the filter requirements (50 uS time constant in the U.S.?) into parameters for one of the common (e.g. RBJ's) biquad cookbook algorithms, given typical audio sample rates?  Is the formula for a low shelving filter appropriate?

You could also implement a 2nd order IIR filter which has response of 1st order filter by using some other TF but one used in RBJ's cookbook.
Here's one way to do it using MZT (C++).

double Fs = 44100; // higher sample rate results better filter accuracy
double T1 = 0.000075; // time constant
double p1 = -exp(-1.0/(Fs*T1));
double z1 = 1 + z1;

a0 = 1.0;
a1 = p1;
a2 = 0;

b0 = 1.0;
b1 = z1;
b2 = 0;

% swap between a1, b1 to select pre- or de-emphasis

Coefficients from above calculations adds +13.63dB gain so you need to normalize before use.

Same as Octave code:

% PACKAGES
pkg load control
pkg load signal

clear all;
clc;

fs = 44100;
samplingtime = 1/fs;

% analog prototype
A2 = [1];
B2 = [0.000075 1];
Ds = tf(B2, A2); % swap between A2 and B2 to select pre- or de-emphasis
Ds = Ds/dcgain(Ds);

% MZT
T1 = 0.000075; % 75us
z1 = -exp(-1.0/(fs*T1));
p1 = 1+z1;

a0 = 1.0;
a1 = p1;
a2 = 0;

b0 = 1.0;
b1 = z1;
b2 = 0;

% swap between a1, b1 to select pre- or de-emphasis

Bmzt = [b0 a1 b2]
Amzt = [a0 b1 a2]

DzMZT = tf(Amzt, Bmzt, samplingtime);
DzMZT = DzMZT/dcgain(DzMZT);

%% Plot
wmin = 2 * pi * 20.0; % 20Hz
wmax = 2 * pi * ((fs/2.0) - (fs/2 - 20000)); 20kHz

figure(1);
bode(Ds, 'b', DzMZT, 'c', {wmin, wmax});
legend('Analog prototype', 'MZT 2nd order','location', 'northwest');
grid on;

There are few other ways to calculate the coefficients as like what mentioned here.

If you don't want to code then here are (75us) IIR coefficients for most common Fs:

44100
den = [ 1 -0.2430506562509079 -0.3653116688328392 ]
num = [ 0.226236781592716 0.1524695626369773 0.01293133068655965 ]

48000
den = [ 1 -0.2625656192085276 -0.3736509425970024 ]
num = [ 0.2102915623850396 0.1415200473180778 0.01197182849135262 ]

88200
den = [ 1 -0.3702185951814965 -0.419953874650389 ]
num = [ 0.1217100727211314 0.08132561501764404 0.006791842429339071 ]

96000
den = [ 1 -0.3813262768034237 -0.424764645791489 ]
num = [ 0.1125121591800709 0.07513002079134688 0.006266897433669414 ]

176400
den = [ 1 -0.4405600251706213 -0.4505357896180226 ]
num = [ 0.06328932872397297 0.04212236025403318 0.003492496233349943 ]

192000
den = [ 1 -0.4464885169088165 -0.4531262255096845 ]
num = [ 0.05834721577148291 0.03882109525679095 0.003216946553225094 ]

352800
den = [ 1 -0.4775710488916642 -0.4667412351991018 ]
num = [ 0.03239176683609103 0.02151798489625147 0.001777964176891494 ]

384000
den = [ 1 -0.4806342452619782 -0.4680860610215655 ]
num = [ 0.02982988806358219 0.01981315821469544 0.00163664743817871 ]

Magnitude error less than 0.1 dB (20Hz-20kHz)


Ultrasonic Interference with Electromagnetic Waves? I am currently working on ultrasonic frequencies (20 kHz-60 kHz). According to a news snippet I found:

  The state government has also tried employing ultrasonic guns that can scare monkeys away and create van vatika or small sanctuaries for different troops. The plan to use ultrasonic guns, which cost Rs 20,000 each, was shelved after municipal officials said the ultrasonic frequencies jammed phone and internet services.

stated that the Ultrasonic Guns jammed up phone and internet services.

Aren't Ultrasonic (Mechanical energy) and Phone / Internet transmissions (Electromagnetic energy) completely different in nature and operate on different spectrum? Ultrasonic waves they might have used might go up to 100 kHz, while the Internet (2.4 GHz) and mobile systems (GSM: 900 MHz and above) are in completely different spectrum. Is this possible?

Not sure about up into the GHz ISM bands, but unfiltered, unshielded switching power amplifiers and power supplies running at high audio frequencies can generate and radiate harmonics well into the RF HF bands.


How can extract the cosine transform formula used for 2D by scipy.fft.dct For the 1D cosine transform the documentation is clear in here, and I can reproduce it easily:
The formula is:
$$y_k= 2 \sum_{n=0}^{N-1}x_n \cos \left( \frac{\pi k(2n + 1)}{2N}\right)$$
and here it is reproduced manually, for example for the third harmonic:
import numpy as np
from scipy.fft import fft, dct

y = np.array([10,3,5])

# dct() call:
print(dct(y, 2))

# manual reproduction for k = 2:
N = 3
k = 2
n = np.array([0,1,2])
2 * (
    y[0] * np.cos((np.pi*k*(2*n[0]+1))/(2*N)) +
    y[1] * np.cos((np.pi*k*(2*n[1]+1))/(2*N)) +
    y[2] * np.cos((np.pi*k*(2*n[2]+1))/(2*N)))

in both cases I get a 9.
But there is no documentation for the 2D DCT, and my attempts at hacking the formula with a toy matrix are not working out:
Compare:
z = np.array([[ 2,  3       ],
             [ 10, 15]])
dct(z, axis=0)  # dct() call

to for instance:
N = 2
M = 2
k = 0
l = 0
n = np.array([0,1])
m = np.array([0,1])
M*N * (
    z[0,0] * np.cos((np.pi*k*(2*n[0]+1))/(2*N)) * np.cos((np.pi*l*(2*m[0]+1))/(2*M)) +
    z[0,1] * np.cos((np.pi*k*(2*n[0]+1))/(2*N)) * np.cos((np.pi*l*(2*m[1]+1))/(2*M)) +
    z[1,0] * np.cos((np.pi*k*(2*n[1]+1))/(2*N)) * np.cos((np.pi*l*(2*m[0]+1))/(2*M)) +
    z[1,1] * np.cos((np.pi*k*(2*n[1]+1))/(2*N)) * np.cos((np.pi*l*(2*m[1]+1))/(2*M))
)

for the first coefficient.
Can anyone help me reconcile the output of dct() with my attempted manual calculation?
I guess the formula is not...
$$y_{k,l}= 4 \sum_{n=0}^{N-1}\sum_{m=0}^{M-1}x_n y_m \cos \left( \frac{\pi k(2n + 1)}{2N}\right)\cos \left( \frac{\pi l(2m + 1)}{2M}\right)$$
but it would be really easy to correct if I could get the same output manually for one of the coefficients on the example matrix above.
The methods applied here make explicit reference to this author, but they don't help me understand what the function is doing.

PS The formula turns out to be spot on if it is called twice as in:
dct(dct(z.T, axis=0).T, axis=0)

I'd take bets that the 2D-DCT is just separable into identical row-wise and column-wise DCTs. Anything else would be surprising.
Anyways, you know the 2D-DCT is a reversible linear operation: It maintains dimensionality, and has only a trivial null space. Therefore,  you can analyze what it does directly from applying it to an arbitrarily chosen basis of your matrix space.


Filtering of a sampled signal Let's say we have 2 signals of frequencies 4 kHz and 8 kHz.

If we sample both signals with sampling frequency 10 kHz and pass the signal through a low pass filter with cut off frequency of 10 kHz, what frequencies will be present in the output signal?

The 8kHZ signal causes aliasing as we should sample at 2fmax = 2*8kHz = 16kHz (normally we should sample faster than 2fmax). It is also strange that you have 2 signal @4kHz and 8kHz but you want to use a LP at 10kHz.


Improving spectrogram resolution in Python? I'm using the specgram() function in matplotlib to generate spectrograms of speech wave files in Python, but the output is always of vastly inferior quality to what my normal transcription software, Praat, can generate. For example, the following call:

specgram(
    fromstring(spf.readframes(-1), 'Int16'),
    Fs=framerate,
    cmap=cm.gray_r,
)

Generates this:

While Praat, working on the same audio sample with the following settings:

View range: 0-8000Hz
Window length: 0.005s
Dynamic range: 70dB
Time steps: 1000
Frequency steps: 250
Window shape: Gaussian

Generates this:

What am I doing wrong? I've tried fiddling with all the specgram() parameters, but nothing seems to improve the resolution. I have virtually no experience with FFTs.

Here are the matplotlib.specgram parameters

matplotlib.mlab.specgram(x,
                         NFFT=256,
                         Fs=2,
                         detrend=<function detrend_none at 0x1dd6410>,
                         window=<function window_hanning at 0x1e0b1b8>,
                         noverlap=128,
                         pad_to=None,
                         sides='default',
                         scale_by_freq=None)

The parameters provided in the question description need to be converted to comparable mpl.specgram
parameters.  The following is an example of the mapping:

View range: 0-8000Hz            Fs=16000
Window length: 0.005s           NFFT = int(Fs*0.005) = 80
                                noverlap = int(Fs*0.0025) = 40
Dynamic range: 70dB             n/a
Time steps: 1000                n/a
Frequency steps: 250
Window shape: Gaussian          default window is hanning change to gaussian

If you use 8ms you will get a power of 2 FFT (128).  The following is the description of the Praat settings from their website

  View range (Hz) :
  the range of frequencies to display. The standard is 0 Hz at the bottom and 5000 Hz at
  the top. If this maximum frequency is higher than the Nyquist frequency of the Sound
  (which is half its sampling frequency), some values in the spectrogram will be zero,
  and the higher frequencies will be drawn in white. You can see this if you record a
  Sound at 44100 Hz and set the view range from 0 Hz to 25000 Hz.

  Window length :
  the duration of the analysis window. If this is 0.005 seconds (the standard),
  Praat uses for each frame the part of the sound that lies between 0.0025 seconds
  before and 0.0025 seconds after the centre of that frame (for Gaussian windows,
  Praat actually uses a bit more than that). The window length determines the
  bandwidth of the spectral analysis, i.e. the width of the horizontal line in
  the spectrogram of a pure sine wave (see below). For a Gaussian window, the -3 dB
  bandwidth is 2*sqrt(6*ln(2))/(π*Window length), or 1.2982804 / Window length.
  To get a broad-band' spectrogram (bandwidth 260 Hz), keep the standard window
  length of 5 ms; to get anarrow-band' spectrogram (bandwidth 43 Hz), set it to
  30 ms (0.03 seconds). The other window shapes give slightly different values.

  Dynamic range (dB) :
  All values that are more than Dynamic range dB below the maximum (perhaps after
  dynamic compression, see Advanced spectrogram settings...) will be drawn in white.
  Values in-between have appropriate shades of grey. Thus, if the highest peak in
  the spectrogram has a height of 30 dB/Hz, and the dynamic range is 50 dB (which
  is the standard value), then values below -20 dB/Hz will be drawn in white, and
  values between -20 dB/Hz and 30 dB/Hz will be drawn in various shades of grey.

Link to Praat settings

The OP's question might be concerning the contrast difference between the Praat
specgram and the mpl (matplotlib) specgram.  Praat has a Dynamic Range setting
which affects the contrast.  The mpl function does not have a similar setting/parameter.  The mpl.specgram does return the 2D array of power levels (the spectrogram) the
dynamic range could be applied to the return array and re-plotted.

The following is a code snippet to create the plots below.  The example is ~1m15s speech with a chirp from 20Hz-8000Hz.

import numpy
import pylab
import wave
import array
pylab.close('all')
w1 = wave.open('example_no_noise.wav')
w2 = wave.open('example_noise.wav')
# hmmm, probably a better way to do this, scipy.io function?
x1 = numpy.array(array.array('h', w1.readframes(w1.getnframes())))
x2 = numpy.array(array.array('h', w2.readframes(w2.getnframes())))
x1 = x1 / (2.**(16-1))  # normalize
x2 = x2 / (2.**(16-1))  # normalize
Fs = 16000.
NFFT = int(Fs*0.005)  # 5ms window
noverlap = int(Fs*0.0025)
pylab.figure(1)
pylab.specgram(x1, NFFT=NFFT, Fs=Fs, noverlap=noverlap,
               cmap=pylab.get_cmap('Greys'))
pylab.title('Full 1m15s example min noise')
pylab.figure(2)
pylab.specgram(x2, NFFT=NFFT, Fs=Fs, noverlap=noverlap,
               cmap=pylab.get_cmap('Greys'))
pylab.title('Full 1m15s example more noise')
pylab.figure(3); n=2100*176;
pylab.specgram(x2[n:n+256*256], NFFT=NFFT, Fs=Fs, noverlap=noverlap,
               cmap=pylab.get_cmap('Greys'))
pylab.title('Full ~4s example min noise')
pylab.figure(4); pylab.plot(x1[n:n+256*256])


Error Propagation through an IFFT I'm not sure how to approach this problem.  I will describe what I am trying to do, and what some of my matlab outputs look like.  I'm not trained in signal processing or anything, so please be patient reading.

I have a fourier transform of some signal, and its frequency domain covariance matrix.  For example, say my signal $x[f]$ is $100$ points long.  Its covariance matrix $\Sigma$ is $100 \times 100$.  I noticed that its diagonal is real, and when I plot it, it looks like a PSD.  All the off diagonal elements are complex.

What I wish to do, is take the IFFT of the signal $x[f] \rightarrow x[t]$, and then plot $x[t]$ with error bars from the covariance $\Sigma$.  Is this possible?  After I take the IFFT of $x[f]$, what do I need to do to $\Sigma?$

Thank you for reading!

Maybe it helps you to note that the FFT or IFFT with length $N$ can be realized by multiplication with the Fourier matrix $F$, or its conjugate transpose divided by $N$ respectively.

If you have the covariance matrix $\Sigma$ of your signal in frequency domain, the covariance matrix of the signal after IFFT (or multiplication with $F^H/N$) is then given by the usual way to compute the covariance matrix of a linear transform, i.e.

$$\Sigma' = 1/N^2 \; F^H \Sigma F$$

Note the prefactor $1/N^2$ is just because the IFFT is given by $F^H / N$.


Question About Sampling White Noise Assuming that we have a continuous signal $v(t)$

$$v(t) = \int_0^t g(u)\, \mathrm du\tag1$$

where $g(t)$ is white noise.Then take the derivative of it.

$$\frac{\mathrm d}{\,\mathrm dt} v(t) = g(t)\tag2$$

Then sample $v(t)$ so that we get discrete-time model, and calculate the finite-difference instead of derivative

$$ v\left( (n+1) T_s\right) - v\left(nT_s\right) = \tilde g \left(nT_s\right) \cdot T_s \tag3$$

where $T_s$ is sampling interval.

I'm told that $\tilde g(nT_s)$ is the low-passed white noise, so how to prove it?

A continuous-time white noise process $\{X(t)\colon -\infty < t < \infty\}$ is a hypothetical construct that we can treat (in the simplified versions that we use on dsp.SE) as a zero-mean wide-sense stationary process with autocorrelation function $K\delta(\tau)$ where $\delta(\cdot)$ is the Dirac delta. More strongly, all the random variables $X(t)$ are mutually independent, not just uncorrelated. Note that the power spectral density is $S_X(f) = K, \infty < f < \infty.$ Now, the stochastic integral
$$Y(t) = \int_0^t X(u) \,\mathrm du$$ gives us a random variable whose mean we can find as
\begin{align}
E[Y(t)] &= E\left[\int_0^t X(u) \,\mathrm du\right]\\
&= \int_0^t E[X(u)] \,\mathrm du\\
&= 0 & \scriptstyle{E[X(u)] = 0~\text{as per definition of white noise}}
\end{align}
and whose variance is
\begin{align}
\operatorname{var}(Y(t)) &= E[(Y(t))^2]-\left(E[Y(t)]\right)^2\\
&= E[(Y(t))^2]\\
&= E\left[\int_0^t X(u) \,\mathrm du \int_0^t X(v) \,\mathrm dv\right]\\
&=  \int_0^t \int_0^t E[X(u)X(v)] \,\mathrm du   \,\mathrm dv\\
&= \int_0^t \left[\int_0^t K\delta(u-v) \,\mathrm du \right]  \,\mathrm dv\\
&= \int_0^t K  \,\mathrm dv\\
&= Kt \tag{1}
\end{align}
which shows us immediately that if we collect the $Y(t)$'s to form a random process $\{Y(t)\colon 0 \leq t < \infty\}$, then this process is not wide-sense stationary; its variance increases with time instead of being constant as is needed for wide-sense stationarity.

Turning to sampling from continuous-time processes to form a discrete-time process, note that it is not possible to sample a white noise process $\{X(t)\colon -\infty < t < \infty\}$ to extract a single random variable, say $X(5)$, from it in part because white noise is an abstraction and in part because we (as engineers) must model reality and the reality is that an actual honest-to-goodness physical sampler doesn't sample instantaneously at $t$ but grabs a small chunk (say from $t-\varepsilon/2$ to $t+\varepsilon/2$) of whatever it is sampling, and it might be set up to report the average value $$ \frac{1}{\varepsilon} \int_{t-\varepsilon/2}^{t+\varepsilon/2} x(u) \,\mathrm du$$ of the waveform as the sample value.  For deterministic signals, this is of little importance as long as $\varepsilon$ is small in comparison to how fast $x$ can change (for small $\varepsilon$, $\int_{t-\varepsilon/2}^{t+\varepsilon/2} x(u) \,\mathrm du \approx x(t)\cdot \varepsilon$ for continuous signals $x$) but for white noise,
$$ \frac{1}{\varepsilon} \int_{t-\varepsilon/2}^{t+\varepsilon/2} X(u) \,\mathrm du$$ is a zero-mean random variable with variance $\displaystyle\frac K\varepsilon \gg K$.

But the OP is not interested in sampling the white noise process but rather in looking at the random variables
$$Z[n] = Y((n+1)T_s) - Y(nT_s) = \int_{nT_s}^{(n+1)T_s} X(u) \,\mathrm du$$ which can be easily shown to be independent zero-mean random variables of variance $KT_s$. That is,
$\{Z[n]\colon 0 \leq n < \infty\}$ is a discrete-time white noise process: a collection of independent identically distributed random variables with mean $0$ and variance $KT_s$.

Now consider the random process $\{\hat{X}(t)\colon -\infty < t < \infty\}$ that is obtained by passing our original white noise process $\{X(t)\colon -\infty < t < \infty\}$ through an ideal lowpass filter (unit gain in the passband) of bandwidth $W = \frac{1}{2T_s}$ so that the power spectral density $S_{\hat{X}}(f)$ is given by
$$S_{\hat{X}}(f) = S_X(f)\cdot \operatorname{rect}\left(\frac{f}{2W}\right) =\begin{cases}K, & -W < f < W,\\ 0, &\text{otherwise},\end{cases}$$
and autocorrelation function
$$R_{\hat{X}}(\tau) = K(2W)\operatorname{sinc}(2W\tau) = KT_s\operatorname{sinc}\left(\frac{\tau}{T_s}\right).$$
This process $\{\hat{X}(t)\colon -\infty < t < \infty\}$ is called a low-pass or band-limited white noise process (meaning that the power spectral density is constant within the frequency band of interest), and is not an abstraction as white noise (of infinite bandwidth and infinite power) is.  So, we can talk of sampling this process instantaneously if we like  --effectively assuming that the real-life sampler is returning exactly $\hat{X}(t)$ as the sample value, and the samples $\hat{X}(nT_s)$ all are zero-mean independent random variables with variance $KT_s$. Note that this is precisely the description of the $\{Z[n]\}$ process, that is, we can consider $\{Z[n]\}$ to be samples from a white noise process that has been low-pass filtered to bandwidth $\frac{1}{2T_s}$ and then sampled at $T_s$ second intervals (i.e. the Nyquist rate).


Gain of multipath fading components For the Rayleigh fading model, I want use four multipath components,  each component has its delay and gain.

Can the gain in $\textrm{dB}$ of any of the individual four components be positive value?

I can understand that it can be zero or maybe negative indicating reduction in its amplitude due to delay and effects of channel, but I can not understand how it can be positive.

The overall signal resulting from the four components at the receiver may be with positive gain. This is reasonable, due to constructive addition of maybe inphase components, but how can a single component be with positive gain?

It can have a positive gain, it all depends on your definition of 0dB. Sometimes, the first path is assumed to be 0dB, then it can happen that later paths have a higher gain than the first one, leading to positive path gain.

However, in reality, there is always the path loss between the transmitter and the receiver, which is normally added to all paths similarly.

So, for modeling and simulation, it is possible to have positive path gain. However, in reality the signal at the receiver will always have less power than the transmitted signal.


Why Low Frequency called Approximation and High Frequency Detail? While we are decomposing a signal using Wavelets into levels why do we call our high frequency components as Detail and Low Frequency Components as Approximation ?

Intuitively speaking, anything that is 'high frequency' is something that is 'rapidly changing in time'. Anything that is 'low frequency' is something that is 'slowly changing in time'.

If you think about it, any time you have 'detail' in a signal or image, it means that you have, very quick, rapid variations in time or space. This then becomes the 'detail' of your signal/image.

In contrast, if you just care about how something changes generally, you care about long-term trends - you dont care about how things change with every sample, you just care about how things change over a long time, or over long space. In this case, you just want to 'approximate' your signal.


How to define SNR for multiple signals in noise? I have to estimate the number of signals present in a measurement contaminated by additive noise given $n$-dimensional snapshot vectors $\bf x$, modeled as $ \bf x = \bf A \bf s + \bf z $ where $\bf s$ is a $k \times 1$ vector representing $k$ different signals. $\bf A$ is an $n \times k$ non-random matrix and $\bf z$ is an $n \times 1$ noise vector. (This model is common in array processing  problems.)

How should I define the SNR—per signal, averaged or otherwise?

The measurement that is important to you will depend on the application. If you are looking for an over all measurement of all signal power to noise power then you define signal to be the power in all signal bands and the noise all of the powers in the noise bands.

However, if you are trying to use an sub-band adaptive filter to correct for some sort of distortion then the SNR in the particular band that you care about would be important not an overall SNR.

In the same application both measurements could be important at different stages of the system. If you are working in wireless communications and the frequency range that is being digitized contains multiple signals then the over all SNR into the receiver needs to be high enough to limit the noise introduced by amplifiers before digitization. But if after I digitize the band it is then split into multiple sub sections for demodulating the signals I would only care about the SNR of the signal currently being demodulated.

For your application it sounds like an over all measure is important because you are treating the range of frequencies as a whole.

In the application of mismatch correction algorithms for time-interleaved ADCs (my topic) we sometimes use multiple sinusoids input to the interleaved converters to measure the performance increase the correction algorithms. This makes it easier to calculate and visualize things like SNR and SFDR when the desired signal locations are known.

I hope this helps,

Charna


remove spatial sinusoid noise from 2d image Given a 2D image where each row of pixels is shifted by some amount horizontally in a sinusoid fashion, is it possible to recover the original image?

Imagine a camera that can capture one row at a time, and is moved vertically to capture the whole image. While it is moved, the camera vibrates, resulting in noise similar to above picture.

If there is a known reference object that spans the entire height of the picture, such as a vertical bar, this is simple enough to solve.
But can you detect and correct this sort of noise if there is no fixed object in the picture?

If that is not possible, what if you assume some limits for the noise: the approximate amplitude and frequency can be known, but not the phase. Is there a solution under these constraints?

Here is one approach to try:

This problem is about finding the correct re-indexing algorithm for pixel values at every raw.

For each raw we can search for possible shifts optimizing some cost function that represent the correct alignment. The cost function will take the raw intensities and two neighboring raw's (upper and lower raw's) as input parameters.  The output of the cost function should represent a smoothness/continuity measure.

Possible cost function to experiment with:

Minimizing total variation
$l_{1}$ norm minimization in some transform basis.

Constraints on the sine function can be used to improve the optimization algorithm. If you assume a narrow range for pixel shifts you should be able to use a brute-force approach ( searching for all possible shifts).


Effective field in MRI when solving for Bloch equations I have a confusion regarding MRI signal formation. So, as I understand it, we need to solve the Bloch equations for the excitation and the relaxation stages:

So, in the excitation stage, let us assume that we have an RF field excitation which is applied along with the slice selection gradient. The bloch equations are given by:

$$
\frac{dM}{dt} = \gamma M \times B - R
$$

where $R$ is the vector accounting for the relaxation effects. Now, during excitation, the B1 field is on and as I understand it, the B vector is given as:

\begin{align}
B_x &= B_1 \cos(\omega t)\\
B_y &= B_1 \sin(\omega t)\\
B_z &= z G_z
\end{align}

The $B_z$ field is basically the gradient times the position offset from the origin and the $B_x$ and $B_y$ fields are the real and imaginary parts of the B1 waveform.

Now, when the B1 field is turned off, I think $B_x$ and $B_y$ components should be 0. However, let us assume that we have a readout gradient on when we are sampling the MR signal.

What will contribute to the $B_z$ field now? Assuming that the RF excitation is with a 90 degree pulse and after the excitation the magnetization vector lies completely along the $Y$ axes.

In the rotating reference frame, $B_z$ from the static magnetic field $B_0$ is transformed away, hence $B_z = 0$ if no gradient field is applied. All gradients change the $B_z$ value, as you described it with the slice selection gradient. This gradient causes $B_z$ to increase when you move along $z$. Note that it is irrelevant if you do this from the center of the magnet or somewhere close to the bore - you will have the same $z$-dependence at each position $(x,y)$.

For the readout-gradient (assuming it is along left-right direction that we'll call $x$), the same is true: $B_z = G_Rx$. If you go from left to right, $B_z$ changes. In this case, it is independent of the position along the bore, so moving along $z$ does not have an influence on the value, if you stay at the same $x$ coordinate.

A gradient field always acts along the $z$-component of the magnetic field. This is always a bit tricky, since gradients can be "applied along each direction".