Does chirp have constant magnitude frequency response? Pg. 223 claims so, yet my results via DFT differ:

Is the textbook wrong?

My attempted explanations: (code)

DFT vs DTFT: "frequency response" is computed via latter. Still, DFT should resemble a sampled DTFT.

DFT time-domain periodicity, whereas DTFT assumes aperiodic, or "repeats at infinity" with infinite zero-padding.

To address each, I try greater N, and zero-padding - below. Zero-padding appears to correct phase (quadratic if unrolled), and more samples tend to flatten the magnitude for an ever-growing portion of frequencies to the right.
I figure, in the limit N -> inf, the amplitude spike has zero width (like in Gibbs) - but this appears contradicted in the "large N long padding" case, where a nontrivial portion of the amplitude decays with oscillations. Further, the left peak appears to scale with N, behaving more like an impulse in the limit, which won't yield zero energy as in Gibbs phenomenon.

Update: turns out magnitude doesn't spike, but rather decays exponentially, which is far from the expected horizontal line - and, the phase is linear:

Note: see comments below accepted answer for further info.

The book is not wrong, but it does present the concepts on LFM in a clunky manner and can be misleading. The book presents the analytical expression for the LFM spectrum, which is an approximation. It also plays with the plot views and most likely unwraps the phase angles, which is usually required to see the phases you expect.
Usually when you're introducing LFM, you'll show the modulated pulse itself as well as the phase progression in the time-domain. The analytical expressions in the time domain is all you need to observe the linear frequency and therefore quadratic phase nature of LFM. Doing this in the frequency-domain just tends to introduce more confusion. An example of the time-domain LFM pulse and its phase is shown below.

When deriving the expression for the Fourier transform of an LFM pulse, you do indeed yield a magnitude of 1 over the bandwidth of the pulse. This is intuitively satisfying because you have the same contribution from each frequency over the bandwidth.
Confusion does occur however when one goes to plot this if they expect a constant frequency response. With any practical LFM spectrum plot, even with very long pulse widths, one should expect a ripple effect which you have already identified. The quadratic nature of the phase is still captured in the DFT. The spectrum of the LFM pulse from above is shown below.

I haven't tried your code, but it might be that you just need to zoom in the proper areas and unwrap the phases to see what you want. Provided is the MATLAB code to produce the plots to help you in converting it to python.
%% LFM - Time and Frequency Domain

% Sampling
Fs = 50e6;

% Pulse parameters
tau = 50e-6;
bandWidth = 10e6;
alpha = bandWidth/tau;

% Define waveform
t = 0:1/Fs:tau - 1/Fs;
fmcwPulse = exp(1i*pi*alpha.*t.^2); % Complex transmitted LFM waveform

% Plot
figure;
subplot(2, 1, 1);
plot(real(fmcwPulse));
xlabel("Samples");
title("LFM Pulse - Real Part");

subplot(2, 1, 2);
plot(unwrap(angle(fmcwPulse)));
xlabel("Samples");
title("LFM Phase");

figure;
subplot(2, 1, 1);
plot(abs(fftshift(fft(fmcwPulse))));
title("LFM Spectrum Magnitude");

subplot(2, 1, 2);
plot(unwrap(angle(fftshift(fft(fmcwPulse)))));
title("LFM Spectrum Phase");

Update
Modifying the code above so that $\tau = 1 \space s$, which is relatively long, yields a spectrum closer to the ideal flat spectrum that one would expect analytically. The spectrum is shown below.

The textbook is wrong, or very misleading depending on interpretation. The frequency response is square - or, constant over the swept frequency. Assuming an ideal chirp system which can scale up $f$ to any arbitrary duration $t$ (infinite bandwidth);, then the magnitude is a square of infinite width - which happens to be the same as plain "constant".
What follows are a number of related descriptions; thanks to Envidia for helpful discussion - refer to his answer for complex chirp plots.

Phase (of chirp): quadratic. Must be careful with what's meant by "phase"; the instantaneous, i.e. the complete input, $\phi(t)$: $\cos(\phi(t))$. The instantaneous frequency of $\cos(t^2)$ isn't $t$, but $t/2$; see here.

Phase (of freq response): parabolic, with half-parabolas symmetric about and adjacent to the spike at $f_{\text{max}}/2$. So is the textbook wrong? For any physical, finite chirp system, yes - but book again assumes an ideal system with an infinite frequency swept, so the "peak" sits at infinity, rendering the phase a true parabola.

Magnitude (of chirp): for a complex chirp, it's unity / peak amplitude, as for any complex signal with shared inst. phase ($\sin^2 + \cos^2$). For a real chirp, it's the abs. value of amplitude.

Magnitude (of freq response): why the ripples? Because $f$ sweeps fractional frequencies, which are represented as combinations of nearby and higher frequencies when a given frequency doesn't last long enough (for a chirp, each lasts for a mere one sample).

Symmetry: real -> even about DC; complex -> asymmetric, no negative frequencies.

I've used mainly this code for visualizing various scenarios. As a puzzle (to readers; I already understand) for understanding the topmost bullet point, I pose the following: how does a frequency swept from 1 to 2 complete more cycles than a signal with $f=10$ over the same period? What if instead $t$ is from 11.9 to 12.0?

Don't look at the Wiki, or especially this code, as both answer the question (for most part).


A few questions about Slepian and generalized gaussian windows I'm trying to add documentation for all the window functions in scipy.signal, and I'm stuck on the Slepian (same as DPSS?) and Generalized Gaussian windows, which I'd never heard of before.

There are two variables that are shape parameters of some type, p in the generalized Gaussian, and width in the Slepian.  (sig appears to be sigma, the standard deviation.)

2 questions:

Instead of me reverse-engineering and guessing, can anyone explain what these variables are called and what they do?
Can you explain what these windows are useful for or where they are used?

def general_gaussian(M, p, sig, sym=True):
    """Return a window with a generalized Gaussian shape.

    The Gaussian shape is defined as ``exp(-0.5*(x/sig)**(2*p))``, the
    half-power point is at ``(2*log(2)))**(1/(2*p)) * sig``.

    """
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    n = np.arange(0, M) - (M - 1.0) / 2.0
    w = np.exp(-0.5 * (n / sig) ** (2 * p))
    if not sym and not odd:
        w = w[:-1]
    return w

def slepian(M, width, sym=True):
    """Return the M-point slepian window.

    """
    if (M * width > 27.38):
        raise ValueError("Cannot reliably obtain slepian sequences for"
              " M*width > 27.38.")
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1

    twoF = width / 2.0
    alpha = (M - 1) / 2.0
    m = np.arange(0, M) - alpha
    n = m[:, np.newaxis]
    k = m[np.newaxis, :]
    AF = twoF * special.sinc(twoF * (n - k))
    [lam, vec] = linalg.eig(AF)
    ind = np.argmax(abs(lam), axis=-1)
    w = np.abs(vec[:, ind])
    w = w / max(w)

    if not sym and not odd:
        w = w[:-1]
    return w

Possible matches:

nipy's dpss_windows function uses NW, the "standardized half bandwidth corresponding to 2NW = BW*f0 = BW*N/dt but with dt taken as 1"

Matlab's dpss uses time_halfbandwidth Is this the same window?  Is time_halfbandwidth the same thing as width?

This DPSS definition has $ \omega_c$ "the desired main-lobe cut-off frequency in radians per second".

Generalized normal distribution has β (equal to twice p?) which is just called a shape parameter, with normal distribution for β = 1 and Laplace distribution for β = 2.

The Slepian sequences are a family of functions.  Most algorithms calculate 2*NW - 1  sequences at once for a given NW.  N is the number of points in the sequence and W decides half the width of the mainlobe in the frequency domain for the Fourier transform of a given Slepian sequence.  Typically you would use an NW of 3 or 4 for your signal processing.

In scipy, they are asking for the $N$ (m in the Python code) and the $W$ (width in the Python code) parameter separately while in matlab you input the time-bandwidth product $NW$ as a single parameter.  This makes sense because you're typically calculating the Slepian sequences for a fixed window size, $N$.

If you are estimating the power spectrum of a stationary time series the DPSS are the set of windows you should use.

The generalized Gaussian function returns a Gaussian-like function raised to successively higher powers depending on the p parameter.  As p is raised to successively higher powers the generalized Gaussian becomes narrower in the time domain.  The nice property of a Gaussian is that it is it's own Fourier transform and it is the function that achieves the limit concerning the uncertainty principle.  A Gaussian function may be useful if you would like to compute a short-time Fourier transform or spectrogram as an estimator for a time-varying power spectrum of a non-stationary time series.

Fourier transform of GG is also a Gaussian. Using the convolution theorem, the FT(Gaussian \times Gaussian) = FT(Gaussian) \conv FT(Gaussian) = Gaussian \conv Gaussian = Gaussian. Hope that helped!


Representing data in combination of QAM and BPSK I have a question regarding the data representation in different modulations. Assuming we are using MIMO system whith $4$ $T_x$ x $4$ $R_x$  which means four transmission antennas by four receiving antennas. My question, if three antennas are using QAM modulation, and the other antenna uses BPSK modulation. Is it possible the Receiving antennas can distinguish the data sent from each antennas ? I mean is it possible to detect the data sent from transmission antennas?

I'm asking that because in QAM modulation, data will be represented with real and imaginary parts, where in BPSK, data will be represented in real part only. So I see it logic that receiving antenna can detect the data perfectly. But I'm not sure.

Let's us take the simple case $2\times 2$ MIMO. Let $y_1$, and $y_2$ be the received signal at the first and second receive antennas, respectively. Then we have

$$y_1 = h_{11}x_1 + h_{21}x_2 + n_1$$

and

$$y_2 = h_{12}x_1 + h_{22}x_2 + n_2$$

where $h_{ij}$ is the channel coefficient from transmit antenna $i$ to receive antenna $j$, $x_i$ is the transmitted signal from transmit antenna $i$, and $n_i$ is the AWGN at receive antenna $i$.

We can write the above signals in matrix-vector form as

$$\mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{n}$$

Apparently, there is interference here. So, we need to do some sort of equalization. Let's use zero forcing (ZF) for simplicity. In this case you get

$$\mathbf{H}^{-1}\mathbf{y} = \mathbf{\tilde{x}}=\left[\begin{matrix} \tilde{x}_1\\\tilde{x}_2\end{matrix}\right]$$

where $$\tilde{x}_i = x_i + \tilde{n}_i$$

for $i = 1,\,2$.

Now, suppose that $x_1$ is 4-QAM, and $x_2$ is BPSK, both have the same symbol time. Will that be a problem? and why? They seem separable to me, and can be detected given that the receiver knows in advance that $x_1$ is drawn from a 4-QAM constellation, while $x_2$ from a BPSK constellation. Actually we do this in OFDM. It's called bit loading, where different subcarriers are modulated with different modulation sizes, and schemes when the CSI is know at the transmitter. Of course in OFDM the symbols are separable because of the orthogonality of the subcarriers, but we can somehow separate the signals in spatial multiplexing, although not perfectly.

No, they can't simply tell things apart. How should they?

On the air, the transmitted QAM-based signal adds up with the BPSK signal, and even in the exact symbol times, you get numbers with real and imaginary part.

Also, you forget that signals still are pulse shaped; a sensible pulse shaper can have a complex part when going from one BPSK symbol to the other.

What you describe has nothing to do with how MIMO works! MIMO works because for uncorrelated channels, you can find a matrix decomposition of your $\mathbb C^{N_{TX}\times N_{RX}}$ (flat channel case) matrix, not because the signals are different. With that matrix decomposition, you can find "separable" channels, even though the signals add up "on the air" and you just receive a mishmash of all TX signals at each RX antenna (though the amount of the TX signals in that mishmash has to be different for each RX antenna).

Because Zeyad_Zeyad introduces SM but doesn't agree on whether in a MIMO channel, signals from different antennas sum in the air:

The typical MIMO channel is modeled like this:

With each RX antenna seeing a linear combination of all TX signals; hence, the representation as channel matrix $\mathbf H$

$$\mathbf H=\begin{pmatrix}
h_{11} & h_{12} & h_{13}\\
h_{21} & h_{22} & h_{23}\\
h_{31} & h_{32} & h_{33}
\end{pmatrix}$$

so that you can get the first receive signal as the first entry, the second as the second entry, … of

$$\vec r =\mathbf H \vec t\text,$$

where $\vec t$ is the transmit signal vector (i.e. the first entry is the symbol at TX antenna 1, the second at TX antenna 2 …).

As you can see, if say $t_1$ is real but $t_2$ or $t_3$ have an imaginary part, $r_1$ can very well be complex.

So, as you can see, for MIMO to work, you don't need the signals to have different characteristics – you need the channel coefficients $h_{i,j}$ to be independent enough!

Please allow me a personal note from the perspective of someone who helps students survive exams in "Digital Communications Basics" courses at university, @Eng. Badr:

In your previous questions, you were wildly jumping from extremely involved topics like FBMC to more entry-level topics like CDMA then down to the very, very basic of QAM vs PSK (which you had a great misunderstanding about!). Now you're back to a topic that is relatively complex (MIMO) and requires solid understanding of the things like constellations.

You sound like a student who's either preparing for an exam for which they are missing the basic course, or like someone doing a thesis on a topic where they have to acquire basic knowledge by themselves. In either case, jumping into topics, asking a single question, and then jumping to the next topic will not enable you to understand much, and cost very much of your time, without enabling you to do your own deductions. Especially in the thesis case, that's a hazard, because at some point, someone will point at a blaring mistake in what you're doing and say "ok, every person who works on {topic} should know that this can't work, it's really basic".

So, do yourself the favor of getting a digital communications textbook that takes you from the basics of complex baseband over constellations, pulse shaping, a tiny bit of channel coding, over to multi-access schemes like FDMA and CDMA and finally MIMO. I promise that in the medium run, reading that book from start to finish will save you time and effort and, if my hypothesis on the situation you're in is right, will definitely improve the grade you'll be getting.


Problem getting DTW to work for speech recognition I used my ipad to record the words "This. Is. A. Test. This. Is. A. Test.", then downloaded it to my computer as a .MOV file and used ffmpeg to extract audio as WAV. I then used the tuneR package to load the WAV file into R.

library(tuneR)
sndObj <- readWave('myinput.wav')
plot(sndObj@left, type = "l")

It loads fine and the plot looks ok. The problem comes when I segment the words and calculate DTW distances between all pairs of words (using the DTW package in R).

I can't find any evidence that DTW distance is smaller between the correct pairs than between any other pairing (i.e "this->this" is about the same DTW distance as "this->a". Am I missing any important preprocessing?

Below: sample heat map of actual vs expected DTW distances.

Usually DTW includes mel-cepstrum feature extraction and cepstrum mean normalization. I'm not sure if you have it.


Peak to peak amplitude of sum of sinusoidals (harmonic frequencies) I would like to estimate the peak-to-peak amplitude of a periodic signal whose frequency components are known. This is, I have the frequency spectrum (a peak in the fundamental frequency and other peaks in its harmonics) and I would like to compute the peak-to-peak amplitude. Is that possible without reconstructing the signal in the time domain and then detecting peaks?

Suppose that we have a continuous-time periodic signal
$\displaystyle s(t) = a_0 + \sum_{n=1}^N a_n \cos(n\omega_0 t + \theta_n)$.

What does this mean? Do we have a  trace of $s(t)$ on some recorder chart
and the trace looks periodic?  If we did, the question to be solved would be
trivial since we could simply measure the maximum and minimum values of $s(t)$ on
the chart.  So the problem of interest is:

  Given the values of
  $\omega_0, a_0, a_1, a_2, \cdots, a_N, \theta_1, \theta_2, \cdots, \theta_N$,
  find the values of
  $$M_\max = \max s(t)\quad \operatorname{and} \quad M_\min =  \min s(t).$$

The OP wants to find the peak-to-peak amplitude of $s(t)$ and this is, of
course, just $M_\max-M_\min$.  It is worth noting that since $s(t)$ is periodic
with period $T = \frac{2\pi}{\omega_0}$, the maximum value and minimum value
occur (at least once) in each interval of length $T$ on the time axis.  The key
issue here is since we know the $a_i$ and the $\theta_i$ etc., we can write
down the formula
$$s(t) = a_0 + \sum_{n=1}^N a_n \cos(n\omega_0 t + \theta_n)\tag{1}$$
(whether this constitutes "reconstruction" of the time-domain signal
is a question that I will leave to others to answer) but the
standard method for finding the maxima and minima of $s(t)$
requires us to find solutions to the nonlinear equation $(2)$ below
(that is, values of $t$ for which $(2)$ holds):
$$\frac{ds(t)}{dt} = -\sum_{n=1}^N a_n\cdot n\omega_0 \sin(n\omega_0 t + \theta_n) = 0.
\tag{2}$$
Once we have found the values $t_1, t_2, \ldots, t_k, \ldots$ for which
$\displaystyle \sum_{n=1}^N a_n\cdot n\omega_0 \sin(n\omega_0 t_k + \theta_n) = 0$,
we can simply calculate $s(t_1), s(t_2), \ldots, s(t_k), \ldots$ and search
through this list of numbers (it is not necessary to sort the list
as recommended in hotpaw2's answer) to find $M_\max$ and $M_\min$.  But,
ignoring the difficulty of finding the $t_k$'s makes hotpaw's answer not
particularly useful in solving the OP's problem.  The answer by ethereal
is even worse in this regard since all it boils down to is the assertion
that (for $\alpha > 0$)
the maximum and minimum values of $\alpha\cdot s(t)$ are $\alpha\cdot M_\max$
and $\alpha\cdot M_\min$ respectively. Finding the maximum and minimum
of $s(t)$ from knowledge of its Fourier series is a nontrivial task, not at
all as easy as it is made out to be by ethereal or by hotpaw2.

A related problem has received a lot of research attention in the past
thirty years or so.  For a periodic signal such as $s(t)$, the average power
in the signal is readily computed as
$$\bar{P} = \frac{1}{T}\int_0^T |s(t)|^2\,\mathrm dt = a_0^2 + \frac{1}{2}\sum_{i=1}^N |a_i|^2.$$
On the other hand, the peak power is $P_\max = \max\{M_\max^2, M_\min^2\}$,
and is also of interest, especially to system and power amplifier designers,
and the ratio $\frac{P_\max}{\bar{P}}$, aptly named the _peak-to-average-power
ratio (PAPR) has received much attention. It can be calculated exactly for
a given $s(t)$ but only with great computational effort, and so a lot of effort
has gone into finding bounds on the PAPR.  Some of these could help in getting
bounds for the OP's problem.
only with a great deal of effort.

Consider the linearity property of the DFT:

This indicates that the DFT amplitude scales directly with the time domain amplitude of your signal.  Therefore, you should be able to deduce the time domain amplitude directly from the frequency domain amplitude empirically; find the DFT amplitude of a baseline or unity TD signal, and from that you can easily compute a scaling factor that can be applied to your spectral peaks to derive the TD amplitude.


Is it customary to correct for the gain of a window? Consider how the Hanning window is defined:

0.5 - 0.5 * cos(n*2*Pi/(N-1))

By this definition, it has a gain of 0.5, which is simply the average value of the coefficients. By contrast, Flattop windows, as defined, have unity gain, presumably by design.

It would seem appropriate to scale the Hanning window by a factor of 2, but I have never seen this discussed anywhere. It would seem that all windows should be scaled for unity gain.

In practice, are windows typically corrected for their gain?  If not, why not?

EDIT:

Since nobody has given an answer, I'll elaborate a bit.

It is quite easy to find papers that report the gain of the more common windows. But nowhere have I seen anyone refer to correcting the gain before using it for spectral analysis. Maybe I have always missed that statement, or everyone assumes gain correction to be an obvious requirement.

It seems like common sense to set the gain of a window to unity so that signal's energy level is preserved. Furthermore, how can one compare the various windows for amplitude accuracy if one has 0 dB gain, as a flattop does, and the other has nearly 10 dB loss, as the Gauss does.

Windows are also widely used for FIR filter design. In this application, it should be clear that the signal to be windowed, a sinc pulse, has most of its energy in the center of the window. Consequently, the window does little to reduce the sinc pulse's total energy. Thus, when used for filter design, we don't want unity gain, but rather unity peak amplitude, as most windows have, except the flattops. Something other than unity peak amplitude would affect the gain of the resulting FIR filter.

Yes, it is customary to correct for the gain of a window, except for some cases I refer to later. (If you are interested only in the relative amplitude, of course you do not need to correct for the gain.)

Because the window reduces the gain of the original signal (time domain), the amplitude obtained through FFT need to be corrected. For example, if you use the Hanning window, you need to multiply all the amplitudes by 2 (the reciprocal of 0.5). As I understand it, most of the software packages for FFT automatically correct for the window used.

However, such correction is good only when all the frequencies of interest distribute throughout the time domain window. For example, suppose you have 1024 data with all signal levels of zero except for #512 point which has a value of 1 (impulse signal). Obviously, any windows do nothing to the data. So, if you correct the amplitudes for the window gain (multiply by 2), then you will end up with overestimation of the amplitude. If your 1024 data are all zero except for the very 1st point with a value of 1, then every point has a value of zero after windowing, and you lose the signal.

So, if you are dealing with random signals, with all the frequency components expected to lie almost evenly over the length of the signal, you need (or should) correct for the gain of the window you use.

The half factor normalizes to unit amplitude.


Why is it a bad idea to filter by zeroing out FFT bins? It's very easy to filter a signal by performing an FFT on it, zeroing out some of the bins, and then performing an IFFT.  For instance:

t = linspace(0, 1, 256, endpoint=False)
x = sin(2 * pi * 3 * t) + cos(2 * pi * 100 * t)
X = fft(x)
X[64:192] = 0
y = ifft(X)

The high frequency component is completely removed by this "brickwall" FFT filter.

But I've heard this is not a good method to use.

Why is it generally a bad idea?
Are there circumstances in which it's an ok or good choice?

[as suggested by pichenettes]

Zeroing bins in the frequency domain is the same as multiplying by a rectangular window in the frequency domain.  Multiplying by a window in the frequency domain is the same as circular convolution by the transform of that window in the time domain.  The transform of a rectangular window is the Sinc function ($\sin(\omega t)/\omega t$).  Note that the Sinc function has lots of large ripples and ripples that extend the full width of time domain aperture.  If a time-domain filter that can output all those ripples (ringing) is a "bad idea", then so is zeroing bins.

These ripples will be largest for any spectral content that is "between bins" or non-integer-periodic in the FFT aperture width.  So if your original FFT input data is a window on any data that is somewhat non-periodic in that window (e.g. most non-synchronously sampled "real world" signals), then those particular artifacts will be produced by zero-ing bins.

Another way to look at it is that each FFT result bin represents a certain frequency of sine wave in the time domain.  Thus zeroing a bin will produce the same result as subtracting that sine wave, or, equivalently, adding a sine wave of an exact FFT bin center frequency but with the opposite phase.  Note that if the frequency of some content in the time domain is not purely integer periodic in the FFT width, then trying to cancel it by adding the inverse of an exactly integer periodic sine wave will produce, not silence, but something that looks more like a "beat" note (AM modulated sine wave of a different frequency).  Again, probably not what is wanted.

Conversely, if your original time domain signal is just a few pure unmodulated sinusoids that are all exactly integer periodic in the FFT aperture width, then zero-ing FFT bins will remove the designated ones without artifacts.

FFT gives poor time resolution i.e it doesn't give information at what time that particular frequency exist. It gives information on existing frequency components for given signal duration.

By zeroing bins in FFT gives poor resolution after IFFT in time domain.


How to find phase values on precise pixel locations of a given Image? I have computed magnitude and phase spectrum of cameraman Gray scale image  using fft function. Here, we get magnitude and phase spectrum of the whole image. .Please see the reconstruction of image from phase only part from the accepted answer of Image Reconstruction:Phase vs Magnitude  which says  that in an image neighboring pixels of the edge,lines and corner are in phase . So  I want to find phase values on precise pixels

Given an image of dimensions 256*256, how can I find phase values of pixels at locations(100,121),(100,122) and (100,123) ?

Whether it is possible using Fourier Transform or not? If not ,what is the other way to do it? Can anybody explain it with (or without )  code?

clc;
clear all;
close all;

i=imread('C:\Users\RK\Desktop\cameraman.gif');
%i=rgb2gray(i);
i=uint8(i);

figure,
subplot(1,3,1);imshow(i);
title('Cameraman Gray scale Image');

f1=fft2(i);
f2=log(1+fftshift(f1));

m=abs(f2);
subplot(1,3,2);
imshow(m,[]);
title('Magnitude Spectrum');

phase=angle(f2);
subplot(1,3,3);a
imshow(phase,[]);
title('Phase Spectrum');

You can't find the precise phase value of a given pixel because of the uncertainty principle. You can only find the phase value for the pixel's neighborhood using short-time Fourier tranform.

If you want to know the phase for pixel (100,100), ou have to take a neighborhood, for instance (50:150,50:150) and calculate the 2D Fourier transform on this window.


Is Ideal LPF BIBO unstable? In one of other discussions : How to find frequency response, stability, and causality of a linear system?

I found a comment which was quite strong and definitely caught my attention.

  An ideal low-pass filter is an example of a system that is not BIBO
  stable even though its frequency response is bounded for all $f$

I am following the definition of stability as per here in wiki http://en.wikipedia.org/wiki/BIBO_stability

Can anyone give me a proof that ideal LPF can indeed be BIBO unstable?

Of course, ideal LPF with infinite gain can produce unbounded output. The question is restricted to LPF when gain is finite.

This answer is a response to a comment by the OP on on yoda's answer.

Suppose that $h(t)$, the impulse response of a continuous-time
linear time-invariant system, has the property that
$$\int_{-\infty}^{\infty} |h(t)| \mathrm dt = M$$ for
some finite number $M$.  Then, for each and every
bounded input $x(t)$, the output $y(t)$ is bounded also.
If $|x(t)| \leq \hat{M}$ for all $t$ where $\hat{M}$
is some finite number, then $|y(t)| \leq \hat{M}M$ for all $t$
where $\hat{M}M$ is also a finite number.
The proof is straightforward.
$$\begin{align*}
|y(t)| &= \left |\int_{-\infty}^\infty h(\tau)x(t - \tau)\mathrm d\tau\right |\\
&\leq \int_{-\infty}^\infty |h(\tau)x(t - \tau)|\mathrm d\tau\\
&\leq \int_{-\infty}^\infty |h(\tau)|\cdot|x(t - \tau)|\mathrm d\tau\\
&\leq \hat{M}\int_{-\infty}^\infty |h(\tau)|\mathrm d\tau\\
&= \hat{M}M.
\end{align*}$$
In other words, $y(t)$ is bounded whenever $x(t)$ is bounded.

  Thus, the condition
  $\displaystyle\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$
  is sufficient for BIBO-stability.

  The condition  $\displaystyle\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$
  is also necessary for BIBO-stability.

Assume that every bounded input
produces a bounded output.  Now consider the input
$x(t) = \text{sgn}(h(-t)) ~\forall~ t$.  This is clearly bounded,
($|x(t)| \leq 1$ for all $t$), and at $t=0$, it produces output
$$\begin{align*}
y(0) &= \int_{-\infty}^\infty h(0-\tau)x(-\tau)\mathrm d\tau\\
&= \int_{-\infty}^\infty h(-\tau)\text{sgn}(h(-\tau))\mathrm d\tau
&= \int_{-\infty}^\infty |h(-\tau)|\mathrm d\tau\\
&= \int_{-\infty}^\infty |h(t)|\mathrm dt.
\end{align*}$$
Our assumption that the system is BIBO stable means that $y(0)$ is
necessarily finite, that is,
$$\int_{-\infty}^{\infty} |h(t)| \mathrm dt < \infty$$

The proof for discrete-time systems is similar with the obvious
change that all the integrals are replaced by sums.

Ideal LPFs are not BIBO-stable
systems because the impulse response is not absolutely integrable,
as stated in the answer by yoda.  But his answer does not really answer the question

  Can anyone give me a proof that ideal LPF can indeed be BIBO unstable?

A specific example of a bounded input signal that produces an unbounded output
from an ideal LPF (and thus proves that the system is not BIBO-stable)
can be constructed as outlined above (see also my comment on the main question).

The fourier transform of ideal lpf is sinc function in time domain which exist from -infinite to +infinite so it is noncausal and area within it is infinite so unbounded. ..


window for short time fourier transform When we implement the 1D-short time fourier transform, the formula is as such:

$$
Sf[m,l]=\sum_{n=0}^{N-1}f[n]g[n-m]\exp\left(\frac{-i2\pi l n}{N}\right).
$$

the window $g$ used is normally symmetric, why is this so? (Hint: Use the phase of Short time fourier transform.

You are talking about symmetry in two contexts :-

Gaussian
Window

I will talk about these two contexts one by one

Gaussian

A symmetric 2-D gaussian (I assume 1-D as in 1-D the question of symmetry does not arise) means that it does not have any directional preference (dictated by equal variances $\sigma_x$ and $\sigma_y$) . Now the job of Gabor Transform is to give us a joint time-frequency distribution (i.e telling us about frequency and phase content of a local image patch)
Now consider the following two cases :-

Symmetric Guassian ($\sigma_x$ = $\sigma_y$)
In this case the gaussian does not have any directional preference.

In the above image, suppose you are centered on the central pixel and consider a patch around it which is indicated in the above image by a black rectangle in the center. Following points are important for this image :

Doing a Gabor Transform over that patch would yield a time-frequency representation, which involves all the pixels in the patch accordingly weighted by the symmetric gaussian.
All the directions at a certain distance from the origin get the same weigthage.

However, if a gaussian oriented at an angle of $135^{\circ}$, then :-

the white principle diagonal pixels would get almost all the weightage while the other pixels will get little weightage.
The time-frequency representation will be mainly a time-frequency representation around the principle diagonal.
If the values of $\sigma_x$ and $\sigma_y$ are sufficiently small, the gabor transform will be very very close to the time-frequency representation for the straight line which marks the principle diagonal.

Thus, for the Gaussian function,

Symmetric gaussian is non-directional.
Non-Symmetric gaussian is directionally oriented.
Non-Symmetric gaussian will be useful in cases where the time-frequency analysis of a directionally oriented feature in a patch has to be carried out.

Window

The Window size basically represents the time-resolution (or space-resolution in the case of 2-D signals like images). A narrow window means that the focus on some local patch of the image is high and hence a more faithful time resolution of the Time-frequency analysis will be possible.

However, this is limited by the Heisenberg's Uncertainty principle. Although, Gabor Transform attains the theoretical lower bound for the uncertainty principle, it stays within the bounds of the principle. Thus, changing the window size is akin to changing your time resolution.

A Non-symmetrical window (like $3\times7$) is generally not chosen because generally we do not have a reason or a preference for changing our space-resolution in one direction (x or y) relative to other. In case there is a concern for the same, it can be done. However, I have never till now seen such a situation.


FastICA - Purpose of centering and whitening the Data I'm learning about ICA, but dont know whats the purpose of centering the data (making each component have zero mean)and whitening the data (not sure what is this for, is it the same as PCA?), I'm following FastICA algorithm for wikipedia
FastICA

The goal of FastICA is to rotate your data (unitary transform) so that each axis looks as non-Gaussian as possible.  Gaussian data still looks Gaussian when you rotate it.

If you don't "sphere" the data, all the algorithm can really do is rotate the whole block to one axis.  By bringing the mean to zero (centering), and normalizing the variance in all directions (whitening), you give the algorithm freedom to rotate in all directions.

This lets FastICA find the rotations that correspond to non-gaussian data.  All of the rotations it tries will leave the Gaussian components of your data as zero-mean and unit-variance, but non-gaussian components can be optimized.

ICA finds the independent components and PCA find the uncorrelated components. ICA cannot distinguish the strength of the components i.e. ambiguity so whitening does not loose any information but fast to converge.


identification of sounds within audio signals I am looking for methods to identify sounds within audio files.  I was hoping to find some ready-made user-friendly software to do this but haven't found it yet.  I am familiar with Matlab and some C++ and have mech. eng. background so I can do some coding work myself to implement a solution but I'm not a serious programmer.

From my online searches I seem to have found there are 2 potential methods; cross correlation and fingerprinting.  It seems to me that the audio fingerprinting method is potentially more accurate than raw cross correlation.  At this point this investigation is for my own hobby purposes to be implemented with a robot.

Can anyone with related experience comment which approach may be more successful?    It seems there are several audio fingerprinting methods and some open source codes even for matlab so I am intrigued, however, these methods are mostly developed to identify entire songs.

The sounds I want to identify are short (from < 1s to ~5s).  I was thinking if the analysis sample rate was increased to fill the short sound file that maybe the method could work with very short files such as these sounds, but I do not know the implications for the rest of the techniques of fingerprinting.  Does this seem like a reasonable / possible approach to anyone?  If not are there references for other potentially applicable DSP techniques for this I can research?

Thanks.

Audio Fingerprinting can be used to match near identical recordings of sounds with a bit of additive noise. It will not work for similar sounds, cover versions etc.

Without knowing your application, I'd recommend giving Dan Ellis' audio fingerprinting code a try. You'll find MATLAB and Python versions here:

http://labrosa.ee.columbia.edu/matlab/audfprint/


BER Analysis of Dual hop Decode and Forward Relay I want to perform fair BER analysis between Two-hop Decode and Forward relay communication (S-->R-->D) and Direct Transmission from S-->D. The links have Rayleigh fading.
Now for fair analysis 1) Total Transmit Power is kept same 2) Rayleigh Channels with  distance are considered for each link.

In the literature,BER vs (Average SNR) is plotted.

How do I calculate Average SNR for dual hop transmission?

and Can I plot BER vs Transmit power?

I am working with a similar kind of network. I don't know if it is correct, but I think the most proper way is to consider the SNR from the first transmitter, not considering what the relay do to the signal. This because you want to analyze the performance of the end-to-end communication. SO you send the signal with some energy, and you want to know the performance of your signal once you have sent it with that energy.

That's the way I am doing, and comparing it to the direct transmission, I am getting a very similar performance: a difference of less than 1dB for the same BER.

I am not amplifying the signal in the relay, so I am just doing the decode and forward (DF). I think it is reasonable to get this "worst" performance in the DF relay compared to the direct transmission because you have two channels, therefore you have two components of noise. You perform the DF in the relay, but it can decode wrong.

Let me know if this answer was useful, and if you have any more questions. I am not an expert in this area, but I am working on it.

You will need to think about what model you want for the dependence of errors on one hop on the second hop. For instance, the simplest model might be to assume that errors are independent on the two hold. In that case, think about how the BER for the entire link can be determined by the BER in each link.

Hopefully, this gives you enough information to move forward.


FFT the whole input vs FFT segment suppose I have a 4 seconds signal captured at $f_s$ = 1024Hz, which gives 4,096 data array. What is the different between these 2 FFT methods (MATLAB code)

[1] FFT the entire input with 4096-points FFT

N = length(input);
Y = fft(input)/N;

[2] Divide the input into 8 segments (4096/8=512 data point each), and feed each segment to 512-points FFT, and then average the output:

N = length(input);
Y = zeros(1,512);
for i = 1:8
  segment = input(512*(i-0)+1 : 512*i);
  Y = Y + abs(fft(segment,512))/512/8;
end

Question #1: Are these 2 methods equivalent?

Question #2: Is the frequency resolution for both methods equal to $\frac{f_s}{N}$ = 0.25Hz per bin?

Question #1: Are these 2 methods equivalent?

No, method #1 is taking an FFT of the entire sequence resulting in frequency bin widths of fs/N = 1024/4096 = 0.25 Hz. Method #2 is summing shorter FFTs at a frequency resolution (bin width) of 1024/512 = 2.0 Hz. For #2, if you divide by the number of FFTs after you sum them you are then averaging them. This is often times used to reduce the variance of the result.

  Question #2: Is the frequency resolution for both methods equal to 0.25Hz per bin?

No, they are 0.25 & 2.0 Hz. The resolution, or bin width, of the FFT is fs/N.

The frequency resolution of a DFT operation is solely dependent on the #of samples ($N_s$)only. Taking a larger than $N_s$ length DFT implicitly makes use of zero-padding in matlab. This does not increase the frequency resolution (though it may appear so when you see the plots) and rather makes use of spectrum obtained without zero-padding to interpolate to the increased DFT length sample values. This operation helps in smooting the appearance of the spectrum and resolving potential ambiguities between the peaks.

In your case you are dividing the signal into frames and DFT of each frame will have thus a lower fundamental resolution (dependent on the number of sample you use 512 now). You are seeing good resolution (though the second case has poor resolution than first) because the sampling fvequency is much higher than the signal content frequencies hence you have captured enough number of cycles of the constituent components in both cases.

The approach of dividing the signal into frames and doing then averaging the resulting spectrum^2 (energy in the frequency domain) is used in power spectral density computation of stochastic (deterministic (sine waves)+noise) signals where number of available samples is always a limitation, hence it is appropriate to deduce interpretation from an average PSD than a single PSD estimate.


What Are Linear and Circular Convolution? I have some basic understanding of signals and convolution. As far as I know it shows the similarities of two signals. Could I get some explanation in plain English of:

what are the linear and circular convolution
why they are important
practical situation where they are used

Linear convolution is the basic operation to calculate the output for any linear time invariant system given its input and its impulse response.

Circular convolution is the same thing but considering that the support of the signal is periodic (as in a circle, hence the name).

Most often it is considered because it is a mathematical consequence of the discrete Fourier transform (or discrete Fourier series to be precise):

One of the most efficient ways to implement convolution is by doing multiplication in the frequency.
Sampling in the frequency requires periodicity in the time domain.
However, due to the mathematical properties of the FFT this results in circular convolution.

The method needs to be properly modified so that linear convolution can be done (e.g. overlap-add method).

Linear convolution:
     For aperiodic and infinite sequence.
Circular convolution:
     For periodic and finite sequence.


Determine minimum sample rate for continuous sinusoid Consider a signal $$ x(t) = \cos(175\pi t) $$ which is sampled to produce discrete time signal $$ x[n] = x(nT_s) $$ The fundamental period of $x[n]$ is $$ N_0 = 7 $$

Given this, what is the smallest possible sample rate $T_s$? (Ans: 1.6327 ms).

I would assume that this is related to finding the Nyquist frequency. I was thinking:

Since, $$N_0 = 7 \implies f_0 = \frac17 \implies f_{\mathrm{Nyquist}} = 2 \frac17\implies T_s = \frac72 $$ However, this is obviously incorrect. I am not even using any information of the original signal. Any suggestions on what I could be doing wrong here?

Remember that the cosine function repeats itself every $2\pi$. The discrete signal can be written as $$x(n)=\cos(175\pi n T_s)$$

So when the argument reaches $2\pi$, we'll have a period. As we know that the fundamental period is $7$, then

$$175\pi \cdot7T_s=2\pi \implies T_s =\frac{2}{7\cdot 175}=1.6327\cdot 10^{-3}$$

For a pure sinusoid in zero noise, just 3 or 4 non-aliased samples (at almost any sample rate) will work.  Here are derivations of two closed form solutions by C.Turner:

http://www.claysturner.com/dsp/3pointfrequency.pdf

http://www.claysturner.com/dsp/4pointfrequency.pdf

For finite sampling in noise, the required sample rate will depend on the duration of sampling, the S/N ratio, and the parameter (frequency, etc.) estimation accuracy required.  It's not just 2X the frequency or slightly over (except for an infinite length of sampling).


Finding linear convolution of two time series PROBLEM

Two real, causal time series

$$f_k = a^k \quad \text{and} \quad g_k = b^k$$

where $a \neq b$, $|a| < 1$, and $|b| < 1$, are given for $k = 0,1,2,...$

Find the linear convolution $h_n$ of the time series $f_k$ and $g_k$.

ATTEMPTED SOLUTION

I think perhaps my attempted solution here can be improved upon. We have from the definition:

$$h_n = \sum_{k=0}^{n} f_k g_{n-k} = \sum_{k=0}^{n} a^k b^{n-k} = b^n \sum_{k=0}^{n} \left(\frac{a}{b} \right)^k$$

If we let $R = \left(\frac{a}{b} \right)$, the last term can be written as:

$$h_n = b^n \cdot \frac{1 - R^{n+1}}{1 - R}$$

or

$$h_n = \frac{b^n - \frac{a^{n+1}}{b}}{1 - \frac{a}{b}}$$

$$h_n = \frac{b^{n+1} - a^{n+1}}{b-a}$$

I don't see how I can improve upon this answer though, but since we are given certain conditions in the problem, I have a feeling that it is possible to get a more elegant solution.  I see that as $n \to \infty$, $h_n \to 0$, but is there any way I can "improve" on the answer for any $n$?

There is a small mistake in your solution. Here $f_k$ and $g_k$ are infinite duration signal, ie. $n \to \infty$. You made a mistake that you took  n in $h_n$ and number of coefficients in $f_k$ or $g_k$ as same, but it is different.

$$h_n = \sum_{k=0}^{\infty} f_k g_{n-k} = \sum_{k=0}^{\infty} a^k b^{n-k} = b^n \sum_{k=0}^{\infty} \left(\frac{a}{b} \right)^k$$

$$h_n = b^n \cdot \frac{1}{1 - \frac{a}{b}}$$

$$h_n = \frac{b^{n+1}}{b-a}$$

If you want to finite the number of coefficients in $f_k$ and $g_k$, use other variables other than n.

Then n in the RHS of your equation will replace with the new variable.


When does a Window function's sidelobe level become out of band attenuation? Another question I'm struggling to answer by myself.

For a Window function, say Blackman (response pictured), what is the out of band attenuation? If I'm designing a filter with out of band attenuation of 40dB, is it the level of the first sidelobe that I'm interested in? Or is the out of band attenuation dictated by something else (like 3rd sidelobe level for example?)

Thanks!

Windows with low sidelobe levels have large transition
bandwidth's. Transition bandwidth is inversely proportional to N for a
given window.  Indeed, the ratio of transition width over sidelobe level
in dB for fixed length N is virtually a constant.

We can use the following empirical formula to estimate the required filter length

where,

N is filter length,
A is the out of band attenuation in dB
fs is the sampling frequency


audio analog input + usb data cable = sound noises I connect my android smartphone from audio jack output to my pc audio input and get a good analog sound signal. But when I connect this smartphone to  the pc with usb data cable the analog sound becomes hissing e.g. music from smartphone overlaps with noises on pc.
Why does this happen?
How can I eliminate these noises? maybe I have to use some shielded cables? If yes - then audio or usb? or I need some sound direct box or something like that?

The problem was resolved by adding ground loop noise filter like this


Determining channel frequency response from measured data via IFFT I've been working on a small system that transmits a linear chirp over the 35Hz-20kHz frequency range in a room and simultaneously records the transmitted signal with echoes (from the floor, walls, etc.), using two microcontrollers.
Using the received data $y[n]$ and reconstructing the transmitted chirp $x[n]$ in the receiver, I obtain the channel frequency response by $H[k]=Y[k]/X[k]$, $X$ and $Y$ being the FFTs of $x$ and $y$ respectively. I then attempt to compute the channel impulse response by computing $h[n]={\tt IFFT}(H[k])$. I then try to use time-gating to get rid of the echoes.

The fact that I am not using external memory limits my computational ability to using float arrays of length N=2048 at most, so the signals I'm dealing with are in the single milliseconds regions (this of course limits my resolution, but that's a problem I have to come to terms with).

The problem that's bothering me is, when I look at the calculated impulse response I get something like this:

It seems that the main response is seemingly moved cyclicly a bit to the left, so that its beginning is at the end of the array and there is a large "spike" there (either that or the response has a big delay, which doesn't fit the physical system as the speaker and the microphone are just a few centimeters apart for the tests I'm performing).
What might be the reason for this?
I've read some posts around here which looked similar and there was talk about pre-ringing; could this be the case here?

I want to compute the FFT of the time-gated (and ideally, "fixed") impulse response to get some data out of it, and a part of it is obtaining the delay of the main response- can this be done in this case?

Insight would be greatly appreciated!

Avner

Edit: Attempting to bump this along with a follow-up question I had to Peter's response.
If I should attempt to ask this in a separate question instead- please let me know.

"Thank you very much! So far I've been transmitting a chirp that lasts 25ms, and the reception goes on for 32ms (so there's an underlying assumption that the channel is 7ms long due to short distances between speaker and mic). I've padded the chirp samples with zeroes to match the reception length (in samples) and using this scheme gave me the graph attached in the question.

I've read that using a sine-sweep and performing a deconvolution results in separating nonlinear harmonics to "behind" the channel response in time, so that might be what's moved cyclicly and runs over the channel. Do you think there is a solution to this problem, given the limitations? I'm pretty much filling up the 32kB memory on the microcontroller I'm using with two 2048-float arrays that constitute a general array I'm working on and take up 16kB (doing things in-place) and two more int-arrays that take up 12kB more, so there's not a lot of memory left to use."

Without seeing your code directly, my guess is that you're doing something like:

N = length(x);
X = fft(x,N);
Y = fft(y,N);
H = Y./X;
h = ifft(H);

The problem with doing that is this operation performs circular convolution over a length of N.  The trouble is, the "true" length of the convolution of y with 1/X's impulse response will be of length 2N - 1 (assuming x and y are both of length N).

So you are getting time domain aliasing.

The solution is probably to set

N = length(x) + length(y) - 1

or greater.

You should use cross-correlation between the received signal and chirp. Your approach won't work in real time.


Current harmonics via FFT -- do peaks matter? Typically, one would perform an FFT on a current waveform and measure it's magnitude at multiples of the line frequency (e.g. 50Hz, 60Hz) to get the harmonics, and subsequently the Total Harmonic Distortion.

If I know that my line frequency should be 50Hz but the FFT peak is not exactly at 50Hz (say, it is 55Hz), then what is my fundamental?  50Hz or 55Hz?
How about the harmonics?  Should I then measure at multiples of 50Hz or 55Hz?
If say my fundamental peak is exactly 50Hz, but my second FFT peak is at 110Hz, where should I measure the 2nd harmonic?  100Hz or 110Hz?  In other words, should I seek out the peaks?
Lastly, THD is usually defined to take the sum at multiples of the fundamental.  What if my FFT peaks are not exactly at multiples of the fundamental?

Update:
- I was trying to measure THD for my power supply with line frequency of 50Hz. If measure the peak of the FFT, it's 50.35Hz. If I record that as the fundamental, the fractional points add up in the subsequent harmonics (e.g. for the 29th harmonic, the value becomes 1460.26 vs 1450 if 50Hz was used as the fundamental).  From the FFT plot, clearly 1450 is the peak, while 1460 is not.
- Should I look for peaks at each harmonic?

You can interpolate the FFT bin results to get better estimates of the actual peaks in the spectrum.  Use a high quality interpolator, such as a windowed Sinc kernel.


Modulating ofdm in audio spectrum Trying to modulate audio with OFDM data.

My fs = 44.1khz, i'm using 64 sub-carriers, but i'm only using 5 sub-carriers for data, so my bandwidth is about 3.445khz. I'm not sure how to properly write the all real data to a wav file and then reconstruct the complex signal when I read it back. Can someone point me where I'm messing up?

Matlab code:

 close all
clear all

%fft length
N=64;
NN=1024;

%carrier frequency
fc = 15000;

%create the data to send
qpsk1 = (floor(2*rand(1,2))-.5)/.5 + 1j*(floor(2*rand(1,2))-.5)/.5;
qpsk2 = (floor(2*rand(1,3))-.5)/.5 + 1j*(floor(2*rand(1,3))-.5)/.5;

inputiFFT = [zeros(1,29) qpsk1, 0, qpsk2 zeros(1,29)];
outputiFFT = ifft(inputiFFT,N);

%spin data to carrier
y3 = outputiFFT.*exp(1j*2*pi*15000/44100*(0:length(outputiFFT)-1));
%despin to baseband for sanity check
y4 = y3.*exp(-1j*2*pi*15000/44100*(0:length(y3)-1));

figure
subplot(211)
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(outputiFFT,NN)))))
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))))
hold off
title('Data - blue, Data at carrier - r')
subplot(212)
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y4,NN)))))
title('sanity check down mix')

%upconvert to 15khz and write to wav (all 'real')
towav = real(outputiFFT).*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1)) + imag(outputiFFT).*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1));

%write to wav then read from wave
wavwrite(towav, 44100, 'ofdmTest.wav');
[y_wav,FS] = wavread('ofdmTest.wav');

figure
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(towav,NN)))),'color','b')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))),'color','r')
hold off
title('Theoretical data up converted to carrier -r , data written to wav - b')

%convert to I-Q
yRx = (y_wav.*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1)).'+ 1j*y_wav.*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1)).');

figure
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(yRx.',NN)))),'color','b')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y_wav,NN)))),'color','r')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))),'color','g')

hold off
title('yRx -b, data written to wav -r, theoretical - g')

Typically you use Orthogonal Frequency Division Multiplexing (OFDM) with communications. In this context having complex data in the time domain makes sense because the real and imaginary parts of the number represent the inphase and quadrature parts of the signal.

In an audio signal which is what a wav file is used for is only real. As sound is a naturally occurring signal.

I am not sure why you'd want to store the inverse of an FFT in a wav file (listening to it would sound awful). As you've discovered there's no way of guaranteeing that this data is real. You seem to get around this by taking the real part of the output and only writing that to your file.

However, if I assume you need this data stored in a wav file you could store the real and imaginary part in the wav file by interleaving them. You could write something like:

towav = zeros(length(outputiFFT)*2);
towav(1:2:end-1) = real(outputiFFT).*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1));
towav(2:2:end) = imag(outputiFFT).*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1));

in your code above.

Bear in mind this "halves” your sample rate when you write the 'towav' variable to your wav file. Although I am assuming you don't need to listen to the file, so this is unlikely to matter.

EDIT: as the comment below points out if you weren’t generating random QPSK signals and you ensured that the input sequence in the frequency domain was conjugate symmetric you would have a real sequence in the time domain and you would not need to store the imaginary part.


Performing a Deconvolution So I have the following data:

[[-344    13771   4600 ]
 [-275.2  12478   6410 ]
 [-206.4  19443    830 ]
 [-137.6  69392   3830 ]
 [ -68.8  143737  3780 ]
 [   0    189278 16870 ]
 [  68.8  184486  5090 ]
 [ 137.6  188466  9380 ]
 [ 206.4  185023 21680 ]
 [ 275.2  128133  1460 ]
 [ 344    51288   1950 ]
 [ 412.8  10854   4290 ]]

First column is the x value (position in microns).  Second value is the recorded data point (y), and the third column is the error (noise).

I am wondering how I can find the response of my device if the data I record is known to be a convolution of the response and a slit of width 83.6666 (plus noise).

My attempt was to use the convolution theorem of fourier transforms and use MatLab's fft() to solve for the desired function, but I could not figure out how to get everything the same length.  I also thought about being blunt and using deconv() but that gave me something very strange, (with my rect function defined as

>> x=-344:1:412;
>> h=zeros(size(x));
>> h(300:386)=1;

So I used deconv(h,y) and got a very strange looking plot that was not a function.

So I was hoping someone could give me some advice on how to perform this mathematically (code or no code).  Especially with the addition of noise.

Thanks!

en.wikipedia.org/wiki/Wiener_filter

what you are looking for is a Weiner  filter  that takes into account the standard deviations of your data. See http://www.astroml.org/book_figures/chapter10/fig_wiener_filter.html
 and the description in Wikipedia  for implementation.  Basically you will do what you were suggesting,  but with an additional term in the denominator of the ratio of fft's  that is the fft  of the noise.

The functions you want to process have different sample rates. The measured data is sampled every 68.8 microns and the slit data is sampled every 1 micron. So you have to downsample the slit data to the same sample rate, or even better upsample the measured data, to prevent loss of information. Then both signals will have equal lenght and the deconvolution should look something like:

ifft(fft(y)./fft(h))


What is bit rate? How is it related to bandwidth? I am very confused about when to use bandwidth as bit rate and where not. Are these terms one and the same?

As hinted by Mr Duval, the term "bandwidth" has been adopted by networking engineers to mean something similar to "transfer rate" (just as "wideband" in networking tends to mean "high speed"). You've tagged your question "digital communications", though, where the term is very precisely defined.

Bit is the unit of information in information theory. Intuitively, a bit is the information conveyed by the answer to a yes/no question, where each answer is equally likely.

In communications we want to transfer information from a transmitter (information source) to a receiver (information sink). The information transfer rate is measured in bits per second. Informally, you can think of this as the number of zeros and ones that a system can transmit per second.

However, it turns out that information is transmitted and received using physical, analog signals (voltages or electromagnetic (radio) radiation). These signals pass through a physical medium, called a "channel", that acts like a filter: only certain frequencies make it from the channel input to its output. So, the channel is said to have a bandwidth, just like any filter.

One of the main jobs of communications engineers is to transfer information as fast and reliably as possible over a channel of given bandwidth. Basically, the problem is how to design a physical signal that:

Fits over the channel bandwidth
Conveys bits at a high rate
Does so reliably (with few bit errors at the receiver)
Does so at feasible complexity at both transmitter and receiver.

The design of such signals is, in a broad sense, covered by the term "modulation".

I hope this is clear -- please comment if you need clarification.

Bit rate is the number of bits transferred/processed per second. Bandwidth in measured in bits per second. So, you confuse distance with meters. In fact, you should confuse the bandwidth with throughput.


matlab How to design lpf\bpf\hpf without builtin functions related to Filtering voice signal from randomly added frequencies

I'm trying to create a simple low pass filter but i'm not sure what is wrong.

my code:

[y,Fs]=wavread('sound.wav');
m = length(y); %30586
Fc = 500;               %% frequency cut
n = [0:m-1];            %% kernel length / samples
Fs = Fs;    %22000      %% sample frequency / (nyquist - %2*m)
Tc = 2*pi*(Fc/Fs);      %% theta c
Hm = sin(Tc*n) ./ (pi*n);
Hm(1) = Tc/pi;
figure,plot(Hm);

1.

LPF: How to get to the 0-500 range?

We used this code in the class but with Fc=45,m=100,Fs=200 and it looks ok.
but when I try this on this wav file it's not working as expected. I thought it will pass from 0-500 but instead it returns some crowled plot.

only when I set Fc to be a low number I can see the sin wave. so i've try to manipulate by setting Fc=real(ceil(m/500)). it better than nothing but it's not it.

2.

hpf = 1-lpf, right? well is seems not to be. my guess is that my lpf is wrong

3.

I will deal the ripple after i'll manage to do this

You seem to be confused with the time and frequency domain. Your code for an LP filter is the inverse Fourier transform of an Ideal low-pass filter. So by getting a frequency response for the function you wrote you'll get an ideal LP response you want.

However, there are two mistakes:

1) Your frequency for the ideal LP filter needs to be $[-\frac{m-1}{2} : \frac{m-1}{2}]$

2) You need to change the value for the 0-frequency only if M is an odd number.

Now, you can use the function freqz to get your freqency response, so your final code would be like this:

%[y,Fs]=wavread('sound.wav');
%m = length(y); %30586
m = 30586;
Fc = 500;               %% frequency cut
n = -(m-1)/2 : (m-1)/2;            %% kernel length / samples
%Fs = Fs;    %22000      %% sample frequency / (nyquist - %2*m)
Fs = 22000;
Tc = 2*pi*(Fc/Fs);      %% theta c
Hm = sin(Tc*n) ./ (pi*n);
if (rem(m, 2) == 1)
    index = (m+1)/2;
    Hm(index) = Tc/pi;
end

[H w] = freqz(Hm, 1, 1024); % 1024 - or how many points you want for your Fourier Transform

% Now your H is your frequency response, and w is digital frequency in [0,pi].
% You can now plot the amplitude response:
figure,plot(w, abs(H));

And you have your Low-Pass filter. However, I'd like to suggest using the many ways to create filters in Matlab, they were made for a reason, and allow for much easier filter design :)

One more remark, I have to suggest that you know the difference between analog and digital frequency, as they have different notation. Analog frequencies you write in uppercase, e.g. $F, \Omega$, and digital frequencies in lower-case, $f, \omega$. Here's a version of the code which is a little prettier:

M = 30586;
Fc = 500;
Fs = 22000;
wc = 2*pi*Fc/Fs;

n = -(M-1)/2:(M-1)/2;
h = (sin(n*wc))./(n*pi);
if ( M-2*fix(M/2) ) > 0
        index = (M+1)/2;
        h(index) = wc/pi;
        display(index);
end

[H w] = freqz(h, 1, 1024);
plot(w, abs(H));

Cheers :)

EDIT:

When you use freqz, the digital frequency w can be somewhat hard to control, because it's not created in the same way, or more clearly, the points are not selected like a normal frequency we create with, for example, F = 0 : Fs/(M-1) : Fs, even though it should be the "same" frequency.
The result is that it's not really possible to get the output filtered signal with Y = H .* X, and  y = ifft(Y).

However, I would suggest using the function filter(b, a, x) for filtering your signal, which works good with freqz. You can proceed to do just that:

% ------
% Create a random signal x
n = 0:M-1;
F = -Fs/2:Fs/(M-1):Fs/2;
Fx = [2000 7000 10000];
x = cos(2*pi*Fx(1)/Fs*n) + cos(2*pi*Fx(2)/Fs*n) + cos(2*pi*Fx(3)/Fs*n);
X = abs(fftshift(fft(x, Nfft))); % fftshift because we need domain [-Fs/2, Fs/2]
% ---------

% -------
% Now filter the signal
y = filter(H, 1, x);
Y = abs(fftshift(fft(y)));
% -------

That's it :)

fs=70MHz % Sampling freq = 70 MHz
fc=fs/(10);   % pass band corner frequency
fc=fs/(10);   % pass band corner frequency
fc1=fs/(8);   %stop band corner frequency
%change the scaling factor according to ur cutoff frequency
lpfLength=127; % Order/Number of Filter coefficients

% Low Pass Filterlpf=fir2(lpfLength,[0 fc fc1 fs/2]/(fs/2),[1 1 0 0]);

FFThn=fft(lpf,NFFT);    % FFT of filter h(n) i.e H(w). NFFT can be any value like 2048

magfhn=20*log10(abs(FFThn)).';    % Magnitude in dB
magfhn=circshift(magfhn,NFFT/2).';  %circshift is used if u want to plot it from -fs/2 to fs/2

figure;
plot((-NFFT/2:NFFT/2-1)*fs/NFFT,magfhn); %Freq. plot of LPF, If u want to plot from 0 - fs than use 0:NFFT
xlabel('Frequency (Hz)');
ylabel('Suppression (dB)');
title('Low Pass Prototype Filter')
axis([-fs/2 fs/2 -100 5]);


How do I add AWGN to an I and Q representation of a signal? I have a wireless communication system that I am simulating in Matlab. I am performing some watermarking through slightly adjusting the phase of the transmitted signal. My simulation takes the original I (inphase) and Q (quadrature) values and adds in the watermark. I then have to simulate the resulting bit error rate after being transmitted. For now I just need to add varying amounts of thermal noise to the signal.

Since I have the signal represented as its I and Q channel it would be easiest to add AWGN(additive white Gaussian noise) to the I and Q directly. One thought was to add noise to both channels independently, but my intuition tells me that this isn't the same as adding it to the signal as a whole.

So how can I add noise to it when it is in this form?

Yes, you can add AWGN of variance $\sigma^2$ separately to each of the two terms, because the sum of two Gaussians is also a Gaussian and their variances add up. This will have the same effect as adding an AWGN of variance $2\sigma^2$ to the original signal. Here's some more explanation if you're interested.

An analytic signal $x(t)=a(t)\sin\left(2\pi f t + \varphi(t)\right)$ can be written in its in-phase and quadrature components as

$$x(t)=I(t)\sin(2\pi ft) + Q(t)\cos(2\pi ft)$$

where $I(t)=a(t)\cos(\varphi(t))$ and $Q(t)=a(t)\sin(\varphi(t))$. If you wish to add AWGN to your original signal as $x(t)+u(t)$, where $u(t)\sim\mathcal{N}(\mu,\sigma^2)$, then you can add AWGN to each of the terms as

$$y_1(t)=\left[I(t)\sin(2\pi ft) + v(t)\right] + \left[Q(t)\cos(2\pi ft) + w(t)\right]$$

where $v(t), w(t)\sim\mathcal{N}(\mu/2,\sigma^2/2)$

Also note that because the in-phase and quadrature terms are additive, the AWGN can also be simply added to either of the two terms in the $IQ$ representation of $x(t)$ above. In otherwords,

$$y_2=I(t)\sin(2\pi ft) + \left[Q(t)\cos(2\pi ft) + u(t)\right]$$
$$y_3=\left[I(t)\sin(2\pi ft) +u(t)\right]+ Q(t)\cos(2\pi ft)$$

are statistically equivalent to $y_1$, although I prefer using $y_1$ because I don't have to keep track of which component has noise added to it.

Kellenjb,

The noise in both the I and Q are not in fact going to be gaussian. In fact they are going to originate from the same original noise vector. This is because there was only one noise vector to begin with at the receiver. So what is happening, is your signal comes into the receiver, where AWGN is added of course. Soon afterwards though, the receiver is going to project that (signal + noise) onto a sin basis, and onto a cosine basis, thereby giving you your I and Q components.

So now the noise in either branch is no longer gaussian, but are in fact, the product of a sin basis times orignal noise vector, and product of cosine basis times original noise vector.

The way I would recommend to simulate this, (are you doing all of this in baseband?), is to simply construct a sin and cosine basis, and simply multiply against (signal+noise), where 'signal' is your original signal of course, and then of course take it down to baseband after that. In fact once you filter for taking it down to baseband, your noise vectors are going to be non-white, and non-gaussian.

Hope this helps! :)


Why more poles than zeroes? I read that an "improper system" "has more zeros than poles; it is not causal, cannot be implemented, has a strictly proper inverse and has infinite high-frequency gain."

Does causality fail due to instant signal change? I try to make up a recurrence, where $x_n$ depends on both $x_{n-1}$ as well as $x_{n+1}$ but they all can be interpreted as $x_{n+1}$ dependence on $x_{n}$ and $x_{n-1}$ and, thus, are causal.

Is it related with solving recurrent relations saying that P(x) in generating function $f(x) = f_0 + f_1 x + f_2 x^2 + \cdots = P(x)/R(x) \ \ \ \ \ $ must have degree less than R(x) as Miguel A. Lerma says in Generating function of Linear Homogenous Recurrence. Here, $R(x) = 1 + r_1 x^1 + r_2 x^2 + \cdots$ is a reciprocal of the characteristic polynomial of our recurrence $x_n + c_1 x_{n+1} + \cdots = 0.$ This indeed seems to happen whenever which R(x) I try.

On the other hand, https://ccrma.stanford.edu/~jos/fp/Existence_Z_Transform.html seem to identify causality with pole values (if your pole > 1 then series does not converge, which means anticausality) rather than their number and others say that causality is $x_n = 0$ for all $n < 0.$ I do not see how this is related with poles/zeros ratio and my guess that poles > zeroes = causality may be wrong. I have asked to relate various kinds of causality here.

I think the system is always non-causal if it has more zeros than poles. @Hilmar, your implementation might be possible digitally but if you see a transfer function like this:
$H(z) = \frac{(z-1)(z-2)}{(z-3)}$, taking the inverse Z-Transform,
$y(n-1)-3y(n-2) = x(n)-3x(n-1)+2x(n-2)$ is clearly non-causal and cannot be implemented real-time(of course you can store the signals and do processing offline, but real time implementation is impossible). So I guess Non-Causality is the answer, apart from the infinite high-frequency gain.

Update by Val
Ok, I understand now. In wikipedia, we see that $$a_0 y_n + a_1 y_{n-1} + a_2 y_{n-2} + \cdots = b_0 + b_1 x_{n-1} + b_2 x_{n-2} + \cdots$$ is identical to $$H(z) = {b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots \over a_0 z^{0} + a_1 z^{-1} +\cdots}.$$ I have finally noted that there are no positive powers. That is, $a_0$ corresponds to the current output, $z^0 = y[n]$, and $\cdots = a_{-2} = a_{-1} = \mathbf{0} = b_{-1} = b_{-2} = \cdots$. If we had more zeroes then we have more power in the nominator and, after normalizing $H(z)$ so that $a_0$ has $z^0$, we still have a positive power of z in the nominator, which means that output depends on future input.

In the simplest case, $Y(z) = zX(z)-x_0$ corresponds to $y[i] = x[i+1]$, that is output depends on future input (we can also drop $x_0$ assuming our sequences are causal, i.e. $0 = y_{-1} = x_0$ aka $z^{-1}Y = X$ or $Y = zX$). It is curious that causal is an example of non-causal. Could anybody clarify it further?

It also worth noting that $(1-az^{-1})(1-bz^{-1})$ is not a zero of order two. It has equal number of zeroes and poles, since it is equal to $(z-a)(z-b)\over z^2.$

EDIT: based on the comments, the answer below only applies to discrete time systems and poles and zeroes in the Z-plane. NOT to time continuous system and poles and zeros in the s-plane.

A time discrete system with all zeroes and no poles is simply an FIR filter which is always stable and causal. A system with more zeros than poles can be easily split into a cascaded FIR and IIR filter where the excess zeros are end up in the FIR filter.

I don't see any problem with that. It's easily possible to create a system with more zeros then poles that' stable and causal. Here is a simple MATLAB example

z = [.1 .2 .3 .5 .6 .7]; p = [.15 .25 .85]; % 6 zeros & 3 poles more or less random
% plot the impulse response
x = zeros(128,1); x(1) = 1;
[b,a] = zp2tf(z',p',1);
y = filter(b,a,x);
plot(y);


Given an experimental FFT plot, how to find the best match from a set of generated FFT plots? What I'm trying to do is, given some experimental data and it's (magnitude normalized to 1) FFT, match the FFT with one of say 1000 other (magnitude normalized to 1) FFT plots. Please note that the experimental/generated plots all have the same bin width/count.

Initially what I'm doing is taking the normalized cross correlation between the experimental FFT plot and each of the generated FFT plots, and then saving only the zero lag value. Then the highest of these values is found at the end which tells which plot correlated best..

This actually works pretty well, but it doesn't seem like something that is done often.. so

Does this make sense? Is it something that's ever done in the DSP community? Are there better methods of doing this which I'm oblivious to? (Note I have tried an MSE metric with limited success).
If what I'm doing is okay, are there any improvements I can make to it. I have two small problems. The first is computational time, I can compute the zero lag correlation by hand in Matlab as sum(x.*y)/scalefactor (scalefactor = sqrt(sum(x.^2).*sum(y.^2)) ) or I can use xcorr(x,y,0,'coeff') to grab the zero lag value. I'm not sure if these are optimized as well as they could be...  The second issue is that if any of the set of generated FFT's are large, say a constant 1 over f=0 to 0.5, the cross correlation metric will be large despite the experimental FFT being much "different".  This actually makes me think cross correlation is not the right metric to use, so, any thoughts?

UPDATE: So thinking about this some more, I determined what I was doing was equivalent to just taking the two sequences and performing a circular convolution on them and saving the maximum value. This is the same as computing the FFT of the sequences and then cross correlating them and saving the zero lag value.

With this method I end up doing something like ifft(fft(x).*conj(fft(y))) to do the circular correlation quickly, although I'm still trying to determine if this is more efficient than what i was doing before...

Regarding the problem with high correlation values, this happens when both signals are not zero-mean.

Assume two signals $x[k]$ and $y[k]$ of length $N$ each. Each signal can be decomposed into the sum of a constant and a mean-free signal as follows:

$$x[k] = \overline{x} + x_0[k]$$

with $\overline{x}$ being the mean

$$\overline{x} = \frac{1}{N} \sum_{k=0}^{N-1}x[k]$$

and $x_0[k] = x[k] - \overline{x}$. Do the same to get

$$y[k] = \overline{y} + y_0[k]$$

Correlating both sequences at zero lag yields (ignoring scaling here, but this does not change the issue)

$$c = \sum_{k=0}^{N-1}x[k] \cdot y[k] = \sum_{k=0}^{N-1}\overline{x} \cdot \overline{y} + \overline{x} \cdot y_0[k] + \overline{y} \cdot x_0[k] + x_0[k] \cdot y_0[k]$$

Since $x_0[k]$ and $y_0[k]$ are zero-mean, the cross-terms cancel out so the final result is

$$c = N \cdot \overline{x} \cdot \overline{y} + \sum_{k=0}^{N}x_0[k] \cdot y_0[k]$$

So if both correlated signals are not mean-free, the correlation will always contain the product of both means times the amount of samples. If you don't want this and the mean is not relevant, you should subtract it prior to correlating.

Regarding the Matlab question. You can pack all reference signals next to each other in a matrix and then correlate with all at once at zero-lag. If you have $M$ reference signals of length $N$ each, you could do (assuming 'signal' as a column vector and the reference signals saved in columns of 'ref' as well)

c = sum(repmat(signal, 1, M) .* ref)

to obtain a vector row vector c containing the correlations of your experimental signal with all the reference signals ref. No need to calculate the whole correlation or convolutions which is indeed a huge overhead. If you do not want to use repmat, you can also have a look at bsxfun which can achieve the same result and is generally interesting when you want to do element-wise operations with a vector on each row or column of a matrix.

If the signals do not have any dominating harmonics in them and they all look similar in envelope and smooth enough, maybe you can try sampling the FFTs and performing the same(correlation metrics). It might save you some amount of computation.

And as for your question as to if you're doing something out of the ordinary, I don't think so. Cross correlation at zero lag is the best metric that can give you the closest answer.


Finding the Harmonics contained in an Image? I am beginer in image processing. Any signal whether it is 1D,2D or any multidimensional signal can be represented using combination of number of sine and cosine wavesforms (harmonics). Similerly any image can be termed as a function of sinusoidal signals.

I want to see individual pattern for the number of waveforms/harmonics present in an image.for example the following  link contains few of them http://www.cs.toronto.edu/~rfm/factored/filters_out.png . So is there any code in Matlab/ OpenCV  or results to understand easily these components(sine and cosine functions(harmonics)) present in any image?

this questions is asked fairly often. If you don't have a grasp of fft in 1d, higher dimensions may be difficult to grasph. But logically it makes a lot of sense once you see it What does frequency domain denote in case of images?

basically sharp changes (high contast noise) usually has high frequencies. and smoother gradients have low frequencies. just as a quick crash course

edit

sorry i thought the other post had code, its just code agnostic conceptual stuff. anyway ill tell you matlab (since your tag says matlab in it) and it is more familiar to me than openCV

function fft_im = createFFTIm(grayscale_im)
    fft_im = fft2(grayscale_im);

    figure()
    subplot(2,2,1:2);imshow(grayscale_im); title('original');
    %fftshift puts the 0 frequency in the middle of the plotting area, normally
    %it would be split, part of it in the 4 corners of the image
    %we add 1 because log(0) is -inf, which we dont want
    %we scale it jsut to show the colors better
    subplot(2,2,3);imagesc(log(1+abs(fftshift(fft_im)))); colormap(gray);
    title('magnitude spectrum');

    subplot(2,2,4);imagesc(angle(fft_im));  colormap(gray);
    title('phase spectrum');
end

edit #2 i slightly modified the code

and here are two examples

notice with the chess example, there are visible peaks in the magnitude spectrum. Just like the regular fft, these peaks correspond to a sing wave of a specific frequency. Since out image is periodic there is are strong sines at the corresponding frequencies. You also notice since the image is near symmetric, the phase is also periodic.
The phase is responsible for modeling the fact that an image is not symmetric

edit #3
this decomposes the figure into individual elements. In the code I limited it to a 4x4 subblock of pixels (the upper left portion of the figure) but you can change this all the way up to the size of the image (though for any large image this would be incredibly slow)

As in my comment, I processed the image one pixel at a time, and did an FFT2 on each of those 1 pixel images. The images themselves are still the same size as the original, it just has only one pixel value populated

function fft_im = decomposedcreateFFTIm(grayscale_im)
    %i only did a small 16 pixel sub block of the image, 4x4 pixels from
    %the upper left corner of the image
    max_rows_to_process = 4;
    max_cols_to_process = 4;

    [m,n] = size(grayscale_im);

    %this is used to isolate pixels, for individual analysis
    zer_array = zeros(m,n);

    figure(1);title('magnitude')
    figure(2);title('phase')

    %i only did a small 16 pixel sub block of the image, 4x4 pixels from
    %the upper right corner of the image
    m=4;n=4;

    %for all rows
    for ii=1:1:max_rows_to_process
        %for all columns
        for jj=1:1:max_cols_to_process
            curr_idx = (ii-1)*max_cols_to_process + jj;

            %creates an image with only one pixel
            indiv_pix_pic = zer_array;
            indiv_pix_pic(ii,jj) = grayscale_im(ii,jj);

            %does fft
            fft_im = fft2(indiv_pix_pic);

            %fftshift puts the 0 frequency in the middle of the plotting area, normally
            %it would be split, part of it in the 4 corners of the image
            %we add 1 because log(0) is -inf, which we dont want
            %we scale it jsut to show the colors better
            figure(2)
            subplot(max_rows_to_process,max_cols_to_process,curr_idx);
            imagesc(angle(fft_im));  colormap(gray);

            figure(1)
            subplot(max_rows_to_process,max_cols_to_process,curr_idx);
            imagesc(log(1+abs(fftshift(fft_im)))); colormap(gray);
        end
    end
end

magnitude

phase

If you want to see the results of a Fourier transform, then it's not that difficult.

The FFT of an MxN image is an MxN array of complex numbers.  Each complex value represents the magnitude and phase of a sine wave at one particular frequency - or more precisely one pair of frequencies (horizontal and vertical), since it's a 2D FFT. (Aside: a cosine wave is just a sine wave with a different phase, so we can treat the result as just being a sum of sine waves)

To see the magnitudes, take the absolute value of those complex numbers.  Then re-scale the results to either 0 .. 255 or 0.0 ... 1.0 (the latter is probably better under Matlab).  You now have an MxN array of numbers, which is just the same as an image.  So display it as if it were an image.

If the results are pretty much black with a few white spots, you may need to apply some form of contrast enhancement (such as histogram equalisation) to see anything else.

The resulting "image" will be symmetrical top-to-bottom, and left-to-right.  The value in the centre is the "DC" or zero frequency component.  As you move away from the centre, you're looking at higher frequency sine waves.  The ones at the edges are the highest possible frequencies - corresponding to alternating black and white pixels in the original image.

For most real-world images, the FFT result viewed as an image will be a largely featureless grey mush.  That just tells you that the image contained a mixture of many different frequencies, or in other words that the objects in the image were of a wide range of different sizes and orientations.


Is there two Types of Scale Space in SURF algorithm? I couldn't figure out below paragraph on SURF paper and hope that someone can help me to understand it.

  Bay H., Ess A., Tuytelaars T. Van Gool L. - Speed-Up Robust Features
  (SURF), page 4, column 2

  The large scale changes, especially between the first filter within these octaves (from 9 to 15 is a change of 1:7),
  renders the sampling of scales quite crude. Therefore, we
  have also implemented a scale space with a finer sampling
  of the scales. This first doubles the size of the image, using
  linear interpolation, and then starts the first octave by filtering with a filter of size 15. Additional filter sizes are 21,
  27, 33, and 39. Then a second octave starts, again using filter which now increase their sizes by 12 pixels, after which
  a third and fourth octave follow. Now the scale change be-
  tween the first two filter is only 1.4 (21/15). The lowest
  scale for the accurate version that can be detected through
  quadratic interpolation is s = (1.2(18/9)/2)=1.2.

is that mean they used two different scale space at same time or they used each one individually and choose the one with better result?

They implemented each one separately, then compared their results to each other.

FH-9  --->  it's a detector that it's initial filter size is 9x9

FH-15 -->  it's a detector that it's initial filter size is 15x15


Can the FFT representation of a signal be modified w/o creating glitches? I've read that modifying the bin values in the frequency domain after a FFT poses problems that lead to glitches when one takes the signal back to time domain using iFFT. This is caused iirc because the modified data does not align "continuously" anymore to the bins (and discontinuities $\implies$ artefacts). Because of this it's e.g. not recommended to perform ("ideal filter") equalisation in the frequency domain.

But since there are e.g. interpolation methods,

Are there ways for doing processing such as shuffling or removing bins in the frequency domain such that one can retain sound quality? E.g. by "synthesising" new data using interpolation methods.

In the context of audio something seems to be feasible, since there exists plug-ins such as this that seem to do filtering in the frequency domain:

The way you do it is you use a window that tapers to 0 at the ends, both before and after the FFT (analysis and re-synthesis windows), along with 50% overlap (needs to be a COLA window), to fade out the glitches and keep the signal continuous.

I think as long as bin modifications represent a continuum in time you'll be OK. In other words, bin modifications (scale, offset, etc.) should be close to that of the same bin for the previous window.


Basic Confusion About the DFT and Convolution I am learning DSP (with Digital Images) and I have some elementary confusion about the convolution between two discrete periodic signals. Specifically, I have learnt that when filtering an image, we apply convolution between the image (as discrete 2D signal) and a kernel filter. To visualize what happens in the frequency domain, we calculate the DFT of the image and the kernel and since we are convolving in the space domain, in the frequency domain this corresponds to a point per point multiplication of the two spectrums.

Now, the DFT assumes that the signal to transform is discrete AND periodic, resulting in a spectrum which is discrete and periodic. So here we are assuming that both the signal and the filter are discrete and periodic signals which repeat (usually) after the same M samples. My confusion is, how can convolution work for two discrete periodic signals? Since they repeat infinitely, applying the convolution usual definition as a flipped-shifting multiplication and sum, we would end up with an infinite sum in some single points if the signal to be filtered and the filter overlap. (Because they would overlap in the current repetition, in the next, and the next ad infinitum)

As an example in 1D, let $x[n]$ be our discrete function, which is periodic with a period of $N$ samples in the spatial domain. And let $h[n]$ be the filter impulse response in the spatial domain, also periodic of $N$ samples. We want to filter $x[n]$ with $h[n]$ and get the filtered signal $y[n]$. To do this we convolve these two periodic discrete signals:

$$y[n] = \sum_{k = -\infty}^\infty x[k]h[n - k]$$

However this convolution for a certain $n$, e.g. $n=0$, can lead to an $y[0]$ which is infinite, as both $x[n]$ and $h[n]$ are periodic.

When we convolve between the image and the kernel, we are only considering one repetition of $h[n]$ and one repetition of $x[n]$. Like this:

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

Why can we do this? Will this still give us a periodic discrete $y[n]$ whose representation in the frequency domain is the point-per-point multiplication of the original signal and filter?

As I understand your question, it boils down to: why is it ok, to take the convolution of a periodic signal only on an interval rather than over infinity? The DFT or discrete nature of the signal does not really play a role in your question. It merely provides the basis of assumung your two signals $f$ and $g$ as periodic, on the interval $[0,T]$ such that $$f(t+T)=f(t-T)=f(t)$$ and obviously $g(t+T)=g(t-T)=g(t)$.

Now assume that you define the convolution over the whole $(-\infty,\infty)$ like so:
$$y(t):=(f\ast g)(t)=\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau)$$

From the periodicity above follows that
$$\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau) = \int_{-\infty}^\infty d\tau\, f(t+T-\tau)g(\tau)$$
which you can intepret intuitevely to mean that the convolution should also be periodic. I write "intuitively" because the integral is divergent as we will see below.

Furthermore you can split the integral into a sum of integrals where $\tau$ goes over intervals of length $T$, such as $[n\cdot T,(n+1)\cdot T]$ with $n$ a signed integer number in $(-\infty,\infty)$. Doing that we have (after some simplifications)

$$\begin{eqnarray}
y(t)&=&(f\ast g)(t)=\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau) \\
&=&\sum_{n=-\infty}^\infty \int_0^Td\tau\,f(t-\tau)g(\tau)
\end{eqnarray}$$
which would be divergent just because of the sum. Thus, the convolution of periodic functions produces a divergent integral. Sticking to this definition of the convolution would make it impossible to define the convolution like that. At least one of the functions must go to zero sufficiently fast, to produce a finite integral. Your periodic functios, however, will violate this decay condition. Thus it makes sense to window the functions with a rect function that captures all the information in one period.

So you see why it would make sense to define the cyclic convolution as just the integral

$$y(t)=\int_0^Td\tau\,f(t-\tau)g(\tau)$$

and use the periodic extension, since we already provided some intuition as to why setting $y(t+T)=y(t)$ makes sense. Thus the cyclic convolution carries all the information you would have from the regular definition but avoids the divergence due to the infinite sum. I hope this helps.


Effects of averaged periodogram (Welch's Method) on resolution? I'm trying to make an app(running on an android phone) that can detect the frequency of a guitar string being played. I know it'll be hard to detect that but that's a problem I'll deal with later.

Right now I'm trying to reduce the effects of noise on my FFT of an audio signal. I'm sampling at $44.1\textrm{ kHz}$ and have a sample size of $4096$. This gives me a bin size of $10.7\textrm{ Hz}$.

I've read good things about averaged periodograms killing noise. I want to try to use an averaged periodogram, Welch's Method, to reduce the effects of noise.

What I'm wondering is, if I window it into $1024$ sample windows, then will the resolution of my final product become $\approx 40\textrm{ Hz}$?
If I have overlapping windows to keep my resolution high, will I still have effective noise reduction?

Sorry if this is a bad question. This is my first time doing DSP outside MATLAB so I'm struggling.

Yes, if you reduce the FFT window from 4096 to 1024 samples, then you will have 4 times less frequency resolution. Also, if you have significant overlap between the windows, your resolution will not change. It only depends on the length of the FFT window. Furthermore, with Welch's method, a window for each segment is usually used, which again reduces the spectral resolution due to convolution of the signal spectrum with the window spectrum.

So, yes Welch's method reduces noise, but also spectral resolution.

If you want to have both, why dont you increase the recorded signal length to lets say 4*4096? This would amount to 360ms, which, depending on your requirements, might still be fine for guitar pitch detection. THere's always the tradeoff between spectral resolution and required amount of time.

The answer by Maximillian is correct. I would only add that technically, I think noise is reduced must when the samples are independent, and as such, one should not overlap windows. Overlapping is a trade off between noise reduction and rapid response when a signal is changing.


How to add pitch-shifted version of signal to itself? I want to write an algorithm to add pitch-shifted version of audio signal to itself (like Electro-Harmonix POG does). And I need to do it in realtime. How to implement it?

My guess is: break signal into pieces, get FFT of each piece, copy each of the frequencies in the original signal (for example, for 1 octave up shift copy 20hz to 40hz, 40hz to 80hz and so on), get the iFFT, play it back. But I don't know if I got the right idea. Will it work?

okay, try googling "pitch shifting AMDF" or "pitch shifting TDHS" and see what you hit there.  also check out DSP Dimension and maybe music-dsp and the mailing list.

you need three concurrent processes:

a good pitch detector that gives you a good estimate of the period (to a precision of a fraction of a sample).  this period estimate should be updated as least as often as every 15 or 20 ms.
a real-time resampling and interpolation alg with a pointer that has precision at the fractional sample.  this pointer will have both integer and fractional parts, the integer part tells you where in your delay buffer you are interpolating, the fractional part tells you how you're interpolating (like the interpolation coefficients).  try linear interpolation at first, it's easier and you can improve it later.  now this pointer will advance faster than the write pointer for upshifting and slower for downshifting.  either way, that pointer will run into the edge of the delay buffer (and glitch or click or pop when it wraps around) unless you do something about it.
a splicing and crossfading alg that allows you to splice or seamlessly "jump" from one instantaneous pointer position in the delay buffer to another pointer position.  the difference between those two pointer positions should be one period as determined by your pitch detector.

if i hint more than that, i might be spilling beans (internal trade secrets) from one manufacturer or another.  but this is the main idea.  think of splicing as necessary to time-scale the input (say time-expand if you're upshifting), and then step 2 is playing that time-expanded audio back faster which both returns it to the original tempo, but bumps it up in pitch.

There are many ways to do pitch shift.
here is sample way. You can do it by re sampling (up and down).
As you know down sampling and up sampling , they make your data frequency change too.

You can use lib Sox.


How do I plot the spectrogram of a wav file? I am dabbing a little bit into signal processing on MATLAB. And I have the following question that I would like to know the answer to.

What has to be done : Plot the Spectrogram of an audio file.
What I know : Since its an audio file, I use wavread to read it and use the array returned to plot any of the graphs that I would like to. Only recently, I found that wav files store data on multiple channels. All the files that I have, have two channels. I know that the Spectrogram is supposed to give you an idea of the fundamental frequencies present in your audio file.
What I would like to know : What does the Spectrogram exactly do? What do the parameters mean? This is the structure of the Spectrogram function call : S=spectrogram(x,window,noverlap,nfft,fs). What do the window and noverlap parameters mean? How do I set these parameters? Are they input file dependent? After all this, how do I plot the Spectrogram of the audio file since it contains two channels? The function returns an error when I try to pass the two dimensional(since there are two channels) array returned by wavread. It says that I can pass only a row or a column at a time. How do I go about plotting the Spectrogram of my wav file?

Thanks a million, in advance! : )

This isn't a great question, but I'll try to answer:

  All the files that I have, have two channels.

A spectrogram can only show one channel at a time.  Either make one spectrogram for each channel, or average your channels together into one.  Here's an example of the syntax I found with a quick google search:

[y, fs]=wavread(fileName); % Read wave file
left=y(:,1); % Left channel
right=y(:,2); % Right channel

then you can do spectrogram(left)

  I know that the Spectrogram is supposed to give you an idea of the fundamental frequencies present in your audio file.

Not quite.  The term "fundamental frequency" means something different.

  What does the Spectrogram exactly do?

It shows the frequency spectrum of the signal and how it changes over time.  If the beginning of the signal is white noise (flat spectrum), and the end of the signal is a tone (spikey spectrum), the spectrogram will show how it changed from one spectrum to the other over time.  It does this by dividing the signal up into small chunks and calculating the spectrum of each chunk.

  What do the parameters mean? This is the structure of the Spectrogram function call : S=spectrogram(x,window,noverlap,nfft,fs). What do the window and noverlap parameters mean?

Did you read the documentation?

window is a window function that is applied to each chunk of the file before the frequency transform.  If you didn't use a window function, there would be abrupt changes at the beginning and end of each chunk, which produce wide frequency components that you don't want.  The window function is like a "fade-in, fade-out" envelope, so that you only see the spectrum of the middle of the chunk, and not the artifacts caused by chopping the chunk abruptly at the endpoints.

noverlap is the number of samples that each chunk overlaps the next.  The windowing function attenuates the signal to near 0 at the boundaries of each chunk, so if you didn't overlap the chunks at all, you would be emphasizing parts of the signal and ignoring other parts.  The overlap lets you measure the same part of the signal twice, once in the middle of the chunk where it's important, and once at the ends where it's mostly ignored.

  How do I set these parameters?

Basically you try different ones and see which is best at showing you what you want to see.  Different window functions are designed to do different things, so read about them and what they are meant to do.

  Are they input file dependent?

Only if you're looking for different types of things in different files.

I guess it should be this (in frequency domain, not time domain):

[y, fs]=wavread(input);             % read wave file

left=y(:,1);                        % Left channel
right=y(:,2); % Right channel

fy = fft(left);                     % transform waveform to frequency domain

figure;
spectrogram(fy);